Widespread exposure to SARS-CoV-2 in wildlife communities

Domain: earth_science
Published: 2024-07-29
Editorial outcome: major_revision

Abstract
Pervasive SARS-CoV-2 infections in humans have led to multiple transmission events to animals. While SARS-CoV-2 has a potential broad wildlife host range, most documented infections have been in captive animals and a single wildlife species, the white-tailed deer. The full extent of SARS-CoV-2 exposure among wildlife communities and the factors that influence wildlife transmission risk remain unknown. We sampled 23 species of wildlife for SARS-CoV-2 and examined the effects of urbanization and human use on seropositivity. Here, we document positive detections of SARS-CoV-2 RNA in six species, including the deer mouse, Virginia opossum, raccoon, groundhog, Eastern cottontail, and Eastern red bat between May 2022–September 2023 across Virginia and Washington, D.C., USA. In addition, we found that sites with high human activity had three times higher seroprevalence than low human-use areas. We obtained SARS-CoV-2 genomic sequences from nine individuals of six species which were assigned to seven Pango lineages of the Omicron variant. The close match to variants circulating in humans at the time suggests at least seven recent human-to-animal transmission events. Our data support that exposure to SARS-CoV-2 has been widespread in wildlife communities and suggests that areas with high human activity may serve as points of contact for cross-species transmission.

Introduction
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the causative agent of coronavirus disease 2019 (COVID-19), has resulted in over 771 million human cases and over six million deaths worldwide1. As SARS-CoV-2 becomes endemic in humans, more virulent and transmissible variants. The considerable pathogen pressure imposed by the pandemic has caused concern as to whether SARS-CoV-2 will spill into wildlife populations, establish a sylvatic cycle, and potentially serve as a source for new variants.

Transmission of SARS-CoV-2 to captive animals has been well documented2–4, but detections in free-ranging wildlife are currently limited to only a few species including white-tailed deer (Odocoileus virginianus5–7), feral mink (Neovison vison8), and Eurasian river otters (Lutra lutra9). Experimental infections and modeling of the functional receptor for SARS-CoV-2 (angiotensin-converting enzyme 2: ACE2) have shown that numerous wildlife species may be competent hosts10–15. However, it remains unexplored whether a diversity of wildlife species are infected in natural settings, where exposure to SARS-CoV-2 is likely to be indirect and at a lower exposure dose.

Since the emergence of SARS-CoV-2 in 2019, numerous variants have been detected in humans and animals. Many variants that have become dominant have mutations that increase their infectivity in humans16, and may also impact the virus’s ability to infect new wildlife species. SARS-CoV-2 collected from white-tailed deer have included lineages circulating in humans, caused by human-to-deer transmission5, but have also included lineages with unique mutations suggestive of deer-to-deer transmission17. This implies that only minimal adaptation may be needed for transmission to occur among deer following initial human-to-animal transmission events18. Other human peridomestic species, such as deer mice (Peromyscus maniculatus)12,13 and skunks (Mephitis mephitis)14 have been shown to be capable of viral shedding in laboratory settings11. Collectively, these studies raise important questions about the extent of human-to-wildlife transmission and the ability of other wildlife species to sustain transmission.

Establishment of SARS-CoV-2 infections in wildlife communities could result in novel mutations that increase virulence, transmissibility, or confer immune escape, negatively impacting both human and wildlife populations. Furthermore, as SARS-CoV-2 adapts to not only human hosts, but potentially a wide diversity of wildlife species, SARS-CoV-2 evolution may become more unpredictable19. This could present several challenges for human health, including concerns related to vaccine development targeting human-specific lineages, and novel impacts on pathogenicity and transmissibility of the virus.

Here, we examine how widespread SARS-CoV-2 exposure has been in wildlife communities between May 2022 and September 2023. We used quantitative reverse transcription polymerase chain reaction (RT-qPCR) to examine 789 nasopharyngeal/oropharyngeal samples from 23 species sampled across Virginia and Washington D.C., USA and documented the presence of SARS-CoV-2 RNA in six of these species. In addition, we analyzed 126 serum samples from six species collected before and after the arrival of SARS-CoV-2 and detected neutralizing antibody titers in five of the six species. Finally, we detected an effect of urbanization and human use on seropositivity in animals, and examined genomic data associated with positive samples.

Results
SARS-CoV-2 RNA detections in multiple species
We amplified SARS-CoV-2 RNA extracted from nasopharyngeal/oropharyngeal samples by RT-qPCR from six of the 23 species examined (26.1% of species sampled) and had a total of 23 unique individual animals that were positive (2.9% of samples tested; n=789; Fig. 1a, b, Supplementary Tables 1–3, Supplementary Data 1). This included eight deer mice (P. maniculatus; 4.7%, N=172), four Virginia opossums (Didelphis virginiana; 2.9%, N=140), four raccoons (Procyon lotor; 4.8%, N=84), three Eastern cottontail rabbits (Sylvilagus floridanus; 2.5%, N=118), three groundhogs (Marmota monax; 9.7%, N=31), and one Eastern red bat (Lasiurus borealis; 8.3%, N=12) (Fig. 1b). We had slightly higher positivity rates in field collected samples compared to samples collected at wildlife rehabilitation centers (4.04% compared to 2.24%), which may reflect repeated sampling of wildlife at a site during a SARS-CoV-2 outbreak.

SARS-CoV-2 neutralizing antibodies in mammal communities
We detected SARS-CoV-2 antibodies in five of the six species we collected serum samples from in 2022 (60% neutralization cutoff, Fig. 1c, d and Supplementary Data 2), including the Virginia opossum (37.5%, N=8), raccoon (36.4%, N=11), Eastern gray squirrel (Sciurus carolinensis; 57.1%, N=7), white-footed mouse (Peromyscus leucopus; 16.7%, N=6) and the deer mouse (7.1%, N=14). Percent neutralization values from samples prior to the arrival of SARS-CoV-2 (pre-2020) were significantly lower than those collected after SARS-CoV-2 arrival (t=−10.774, p<0.001, Fig. 1d). Furthermore, four samples (two raccoons, one opossum, and one white-footed mouse) had a percent neutralization above a more conservative 80% cut-off, further supporting that previous SARS-CoV-2 exposure is likely in these species (Fig. 1d). In 2022, an opossum with a positive detection was trapped one month later and collected serum revealed a percent neutralization value of 51.1% which was below our 60% cutoff, suggesting this cutoff may be too conservative for some species.

We found a positive relationship between urbanization (imperviousness) and wildlife seroprevalence (intercept: −1.665±0.48 SE, urbanization slope 0.039±0.02 SE, p=0.031, Fig. 2a, b). However, antibody detections were highest (80%, N=5) at one of the least urban sites, which is a highly visited state park (average imperviousness 1.5%), and more closely matched seroprevalence at our more urbanized sites (50%, N=10 and 33%, N=12). Human visitation at the state park was similar to human activity in urbanized sites (Fig. 2c), and we found a positive relationship between human presence and seroprevalence of SARS-CoV-2 (univariate linear mixed model with species as a random effect; intercept: -1.132±0.36 SE, human presence coeff: 0.705±0.35 SE, p=0.044; Fig. 2c, Supplementary Fig. 1). This relationship was maintained across different neutralization cutoffs (40–65%; Supplementary Tables 4 and 5).

Whole genome sequencing from nasal and oropharyngeal swabs
SARS-CoV-2 sequences were obtained for 12 of the 23 RNA positive samples (Fig. 3, Supplementary Figs. 2–9, Supplementary Tables 2 and 3) and Pango lineages were determined for nine total individuals (39.1%; Supplementary Table 2). The SARS-CoV-2 sequence from an opossum trapped in 2022 was assigned to BA.2.10.1 (opossum, N=1) and shared defining mutations in ORF1a/b, S, E, and M genes found in the BA.2 Omicron lineage (Fig. 3b and Supplementary Fig. 2 and 11). All eight sequences collected in 2023 were assigned to the XBB* Pango lineages (Supplementary Figs. 3–9, Supplementary Table 2). These lineages were circulating among humans in Virginia during the time of collection (Fig. 3a). The lineages of the sequences obtained from wildlife included XBB (deer mouse, N=1; Supplementary Fig. 3), XBB.1.5 (raccoon, N=1; Supplementary Fig. 4), XBB.1.5.10 (opossum N=1; Supplementary Fig. 5), XBB.1.16 (Eastern cottontail, N=1; Supplementary Fig. 6), XBB.1.5.45 (groundhog, N=1; Supplementary Fig. 7), EG.5.1.1 (deer mouse, N=2; Supplementary Fig. 8), and JD.1 (deer mouse, N=1; Supplementary Fig. 9). The remaining three sequences included two from Eastern red bats and one from a deer mouse, which were generated using amplicon sequencing of the S gene. These matched the SARS-CoV-2 reference sequence (Supplementary Table 3) with 95–100% identity but were at insufficient length for phylogenetic analysis and lineage assignment.

The relative similarity of the sequences obtained from wild animals compared to the closest sequences obtained from humans (Supplementary Table 2) suggest recent introductions of these SARS-CoV-2 lineages into wildlife with at least seven independent introductions of SARS-CoV-2 over a several month period. Interestingly, two sequences from the same lineage (EG.5.1.1) were collected from deer mice on the same day and location (site=PP; Fig. 3 and Supplementary Fig. 8). These two sequences clustered together along with a sequence collected from an infected human in Vienna, Virginia (~330 km distance) three days prior (Fig. 3 and Supplementary Fig. 8). This clustering suggests the deer mice might have been exposed to the same source of the virus or could represent a mouse-to-mouse transmission event (Fig. 3c). Sequence from a third mouse collected from the same day and site (site=PP) aligned with another lineage (JD.1, Fig. 3b and Supplementary Table 2), and suggests a separate human-to-mouse exposure (Fig. 3c). At a nearby site (site=NRT; located ~52 km away) sequences produced for two additional individuals (one deer mouse and one raccoon, Fig. 3b and Supplementary Table 2) were also from separate lineages (XBB and XBB.1.5, respectively), further suggesting independent human-to-animal introductions (Fig. 3c).

All sequences from wildlife were aligned to their closest related human sequence based on phylogeny, and then assessed for any unique amino acid substitutions. Whole genome sequencing (WGS) of the SARS-CoV-2 variant collected from a positive opossum in July 2022 revealed several mutations shared within the Omicron clade as well as a unique amino acid substitution not previously identified in other SARS-CoV-2 virus samples collected from humans at the time of sampling (Table 1).

The Omicron sub-lineage, BA.2.10.1 (Fig. 3), includes the G798D mutation found in the Spike (S) protein. The nearest neighbor (EPI_ISL_14334179), which was collected in New York state 11 days prior to our sampling, contains all the nucleotide and amino acid mismatches found in the opossum sequence, compared to the reference sequence from Wuhan, China, except for A22974T (S:E471V) which appears to be unique to the opossum. Only one mutation was found among the eight other SARS-CoV-2 sequences from 2023 as compared to closely related human sequences. However, this mutation (S:H146Q), in a positive detection from a groundhog, is a common mutation found in other Pango XBB* lineages circulating in 202320.

Unique SARS-CoV-2 mutations collected from wildlife
We investigated the unique E471V mutation detected in the opossum sequence from 2022, which was located in the receptor-binding motif (RBM, residues 437 to 508) within the receptor-binding domain (RBD, residues 319–541) of the S protein. The RBD is necessary for binding the viral S protein to ACE2 in both human and animal cells and is therapeutically targeted by neutralizing monoclonal antibodies21. We conducted in silico studies to examine the local impact of the E471V and G798D mutations on the S protein. We employed molecular modeling and molecular mechanics/generalized born surface area (MM/GBSA) free energy calculations to predict the free energy of binding interaction (ΔG binding). Our focus was on examining the local impact of E471V and G798D on the favorability of S protein-hACE2 interaction and the modifications in the S2 site as a hypothesis for the experimentally observed traits found in SARS-CoV-2 wildlife reservoir. Here, the open conformation structure of the S protein trimer from the glycosylated BA.2 Omicron variant was used for residue variant mapping [PDB 7XO8]22. Because the crystal structure shows the trimer bound to the human ACE2 (hACE2) receptor, it is possible to investigate the local impact of mutations in the RBD on interactions between SARS-CoV-2 S protein and the human ACE2, as a prediction for mechanistic impact on S protein-ACE2 interaction favorability as previously reported and validated experimentally23 (Fig. 4). Our free energy studies predict that ΔG binding between the S:E471V-carrying mutation and hACE2 was improved relative to that of BA.2 with hACE2 (-76.5 and -60.2 kcal/mol for E471V vs. BA.2, respectively). Large-scale low-mode conformational sampling was performed on the RBD of BA.2 and E471V Spike (residues 329-531) with hACE2 to probe ΔG binding on multiple potential conformations24,25. Consistent with the results observed in full-length Spike-hACE2 interactions, the mean ΔG binding for the top five most favorable E471V conformations (-151.7±19.2 kcal/mol) was more favorable than the mean for BA.2 conformations (-125.4±4.1 kcal/mol, p≤0:05). The sampled conformations of E471V exhibited a larger range of ΔG binding values relative to BA.2, suggesting a minor variation in the nature of interactions in the RBD due to the E471V mutation. Residue 471 is predicted to be part of a flexible, hydrophobic loop that interacts with the N-terminal domain of ACE2 (Fig. 4). Several point mutations (A475V, S477G, V483A, F490L) within this loop region have been experimentally shown to improve S protein–hACE2 binding affinity and increased resistance to neutralizing antibodies26,27. Our analysis using normalized Kyte-Doolittle hydropathy scores and surface maps show a more hydrophobic interaction surface at not just position 471, but also among flanking residues 469-474 (0.8) compared to BA.2 [-2.6]28. The increased favorability of the E471V mutation for hACE2 interaction aligns with previously reported trends of residues in the RBD influencing stabilization and enhanced affinity29. Any differences may be attributed to entropic effects, emphasizing the potential enthalpic contribution of this mutation to improve hACE2 interaction at the RBD interface.

We mapped the second missense mutation, G798D, in the S2 subunit of the S protein (residues 686 to 1,273) and within the fusion peptide domain (residues 788–806) near the N801 glycosylation site. When glycosylated, the N801 site significantly enhances viral entry, with very low mutation rates observed in adjacent residues (e.g., within 2–3 residues from position 801)30. Interestingly, D798 creates a small, charged pocket on this solvent-accessible S2 loop, decreasing the probability of glycosylation at N801 (0.48) compared to BA.2 [0.61]31. It is hypothesized that G798D could impact structural stability and membrane interaction32 (Fig. 4 and Supplementary Fig. 10) and have a negative impact on Spike maturation and infectivity. While this approach limits the investigation of these mutations and their role in the dynamic conformational states of the S protein, it provides a foundation for understanding the mechanistic impact and connection of structure to experimental observables.

Methods
Study Sites
We collected nasopharyngeal or oropharyngeal samples from wildlife across 43 counties in Virginia and Washington D.C., U.S.A. We collected samples from three wildlife rehabilitation centers in Boyce, Roanoke, and Waynesboro, VA. In addition, we actively captured wildlife from eight sites between 2022 and 2023 that spanned a rural to urban gradient in Giles, Montgomery, Roanoke, and Wythe counties (Fig. 2a, Supplementary Table 8, and see Supplemental Methods for site descriptions.). At six sites we actively trapped wildlife using live-traps from 9-May through 1-July 2022 and at four sites (two sites were the same as 2022) we trapped from 27-June through 11-September 2023.

Trapping and processing
We trapped each of six sites for a 2–4 day session between 9 May and 1 July 2022 (two sites were trapped for two sessions with at least 27 days between them). Additionally, we trapped four sites for a 4 day session between 27-June and 11-September 2023 (Fig. 2a). We used three different sized live traps to capture animals: Tomahawk non-folding traps (Tomahawk Live Trap Co., Tomahawk, WI, U.S.A; 46×15x15cm, model 103.5), Tomahawk folding traps (66×23x23cm, model 205 or 81×25x31, model 207), and Sherman folding traps for smaller animals (Sherman Traps, Tallahassee, FL, U.S.A; 8×9x23cm, model LFG). We set and baited traps in the evening just before dark and checked them the following morning before 10:00 AM to ensure animals did not get overheated.

We processed all animals in the area where they were trapped. We anesthetized larger animals in either a bucket chamber49 or in a 50 L clear box chamber (58×38x56cm, Supplementary Fig. 12). We used a tec-4 funnel fill with a cage mount manifold vaporizer (Ohmeda, West Yorkshire, UK) connected to a portable medical oxygen cylinder (size E) with a built-in regulator (Walk-02-Bout+TM, Airgas Healthcare, Radnor, PA). We used a 2 L/min oxygen mix to 3–4% isoflurane in the chambers. Smaller animals were placed in a small plastic canister and anesthetized with 0.5–1.2 mL isoflurane (dosage depended on species and body size), placed onto a cotton ball held within a perforated canister to prevent direct contact with the animal (Supplementary Fig. 12). Once anesthetized, we moved animals from the chamber, and then masked animals (except for mice; Supplementary Fig. 12) with a dose of 2 L/min oxygen and 2–3% isoflurane.

We marked individuals with aluminum ear tags in each ear (National Band and Tag Company, model 1005-1 or 1005-3 depending on body size). We collected morphometric, sex and reproductive status from each animal. For RT-qPCR detections, we swabbed each individual with a polyester swab (25-8001 PD swab for larger mammals and 25–1000 1PD swab for smaller animals e.g., juveniles, skunks, squirrels, bats, shrews, and mice; Puritan Medical Products, Guilford, ME, U.S.A). Given the small swab size compared to humans and larger mammals and the reduction in material that could be collected from these swabs, we collected two swabs per animal after July 2023. We used an oropharyngeal swab for animals that were too small for nasopharyngeal swabbing (mice, shrews, bats, etc.). All collection swabs were placed in tubes with transport media43 and kept on ice until they were transported back to the laboratory following fieldwork each day. Transport media consisted of a 1:1 volume of Dulbecco’s Modified Eagle’s Medium (DMEM) containing low glucose (1g/l), sodium pyruvate and L-glutamine (Corning) media and 2x DNA/RNA ShieldTM (Zymo Research), which can preserve samples up to 10 days at room temperature43. Field collected samples were all stored in the refrigerator at 4 °C and processed within five days following field collection.

We collected blood for antibody screening from the submandibular vein for mice and from the nail quick of one of the rear toenails for all other species. Blood was stored in microvette 500 z-gel tubes for later processing. We put all blood samples into a cooler immediately following field collection and stored at 4 °C. Three wildlife rehabilitation centers in Virginia provided wildlife oral and nasal swabs following methods described above. We analyzed all samples from rehabilitation centers within 45 days of collection. Protocols were approved under Virginia Tech IACUC protocol #22-061.

Personal protective equipment and contamination control
All personnel collecting and working with the samples wore fit-tested N95 respirators and gloves. We sterilized all equipment between each animal. We also sterilized all traps and equipment between each trapping session, and individuals collecting the samples were tested using the assay described below 1–2 times per week to confirm they were negative throughout the duration of sampling. During fall of 2023, when most positive samples were collected, all personnel were tested daily.

RNA Extraction and detection of SARS-CoV-2 viral RNA
The RT-qPCR based test we used for the identification of SARS-CoV-2 RNA in wildlife has been thoroughly described in Ref43. This SARS-CoV-2 assay, authorized for use in humans by the FDA during the pandemic and adhering to all regulatory standards, detects three targets including the nucleocapsid (N), envelope (E), and spike (S) genes of the SARS-CoV-2 virus, has high analytical specificity, and across all three targets has a 100% probability of detection at 10.0 copies per 10 µL.

We processed swab-containing media shortly after collection and total RNA was purified using 96-well spin columns, assessed for quality control, and subjected to synthesis and amplification using the Power SYBRTM Green RNA-to-CTTM 1-Step kit (Applied Biosystems). The cytochrome c (Cytochrome C Oxidase Subunit 1; Cox1) housekeeping gene, for which species-specific primers were designed (Supplementary Data 3), and HRPP30 housekeeping gene for human samples was used as a control for whether sufficient sample was collected from each individual. RT-qPCR reactions were performed in a CFX384 Touch Real-Time PCR detection system (Bio-Rad), no template controls, and standard curves on each plate as previously described43. Samples were reported conclusively positive when two or more SARS-CoV-2 targets (N, E, and S genes) and the housekeeping gene (Supplementary Data 1) amplified below the threshold established by corresponding standard curves43, which is similar to the criteria established for the Applied Biosystems TaqPath COVID-19 used for confirmation in other wildlife studies7. The use of multiple targets is particularly relevant when considering genetic variations of SARS-CoV-2 among circulating variants for which the sole amplification of a single gene could result in a false negative result if the mutation was in a region of the genome assessed by the test. Thus, a molecular diagnostic test developed to detect multiple genetic targets of SARS-CoV-2 is likely less susceptible to novel genetic variation. For the positive samples from the opossum collected in 2022, the remaining sample was sent to the USDA National Veterinary Services Laboratories (NVSL) in Ames, Iowa for confirmation of results. However, while insufficient sample remained for NVSL to conduct confirmatory testing, raw sequence data from this positive was independently analyzed and confirmed to be from an opossum.

The Virginia Division of Consolidated Lab Services (DCLS) confirmed the performance of the developed assay using a comparator assay that followed the CDC EUA IFU (CDC DOC 006-0099 rev.03) “CDC 2019-Novel Coronavirus (2019-nCoV) Real-Time RT-PCR Diagnostic Panel” that includes the same probes and N primers. Results of all blindly tested samples were incorporated into the Emergency Use Authorization which was approved by the Food and Drug Administration. To date, this assay has been used for sample analysis by the Virginia Department of Health in over 230,000 tests in Southwest Virginia. Lastly, the Molecular Diagnostics Laboratory that tested the samples is CLIA certified and participates in a proficiency testing program from the American Proficiency Institute.

Serology data
We obtained additional samples (Supplementary Data 2) to compare seropositivity prior to the emergence of SARS-CoV-2 to those samples we collected in 202250. The samples were transferred to gold-top microtainer tubes (BD, Franklin Lakes, NJ, USA) and serum was separated via centrifugation at 5000xg for 10 min and transferred to 1.7 mL tubes. Serum samples were then heat inactivated at 55 °C for 1 h. Plaque reduction neutralization assays (PRNTs) were performed similarly to prior reports51,52. Briefly, samples were diluted 1:10 in RPMI-1640 with 10 mM HEPES and 2% FetalPure bovine serum (Genesee Scientific 25-525H). Diluted serum samples were mixed in a 1:1 v/v ratio with a solution containing 1300 plaque-forming units per mL (PFU/mL) of SARS-CoV-2 Delta virus strain USA-GNL-1205/2021 (a generous gift from the World Reference Center for Emerging Viruses and Arboviruses), or ~40 PFU’s per well. Our serum samples were collected in June and July of 2022, following the Delta wave in fall of 2021, and given the life expectancy of the animal species (Supplementary Table 9) we were sampling, we assumed Delta would have been the most likely recent variant they were exposed to53. Each sample was tested in triplicate. The mixture was then incubated at 37 °C for 1 h, and then 50 µL of the virus-serum mixture was used to inoculate wells in a confluent 24-well plate of Vero E6 hACE2-TMPRSS2 cells and again incubated at 37 °C for 1 h. After the 1-h adsorption period, we added 500 µL of overlay media to each well and incubated the plates at 37 °C, as previously described54. After 2 days, the plates were fixed with 10% formalin and stained with crystal violet. The level of neutralization is based on the reduction of plaque counts compared to virus mixed with diluent. Based on previously published sero-surveillance studies55–58, we used the following testing criteria:

1. Greater than or equal to 90% Neutralization at a 1:20 Dilution—Strong Positive
2. Greater than or equal to 60% Neutralization at a 1:20 Dilution—Weak Positive
3. Less than 60% Neutralization at a 1:20 Dilution—Negative

Any individual with a 60% neutralization or higher was counted as positive in analyses analyzing the effects of urbanization and human presence on SARS-CoV-2 exposure.

Urbanization and human presence variables
We obtained estimates of urban imperviousness from the National Land Cover Database 201959,60. We obtained estimates of population density from 2020 census at the 1 km spatial scale61. We created buffers of variable sizes based on species home range size to calculate the mean estimates of imperviousness and population density for each individual capture location (Supplementary Table 10).

We obtained estimates of human presence at the five sites where we collected serology data (Supplementary Table 8), and all visitation data was standardized to number of visits per month. New River Trails State Park and Roanoke Parks and Recreation provided monthly human visitation estimates between March and August 2022. Mountain Lake Biological Research Station provided estimates for how many people used the facilities at the station between March and August 2022; they averaged about 50 people at the station per month with each person staying 7 days on average. Brush Mountain recreational area was closed to the public from March through August 2022, however, there was some trail maintenance work and the occasional trespasser, so we assumed ~10 people used the area per month. For the town of Blacksburg, we used population estimates from the 2020 census (TIGER/Line, U.S. Census Bureau, 2020 Blocks). We created a minimum convex polygon (MCP) in the adehabitatHR package62 in program R v 4.1.163 around the different trapping locations. We then added a 1 km buffer area around each census assuming both animals and residents move in and out of the trapping mcp area. We used ArcGIS Pro (ArcGIS Pro 2.9.0, Esri, Redlands, CA, USA) to calculate the population within each mcp area. We used the intersect tool field calculator, and summary statistics to calculate the population for each census block by multiplying the percent of the block within the mcp by the population count. We then summed the population within the mcp to estimate the total population. We then multiplied the population count by 30 (approximately days in each month) to compare the daily estimates of occupancy to the monthly estimates from the other sites.

Clinical specimen analysis
Human (13,221) and wildlife (789) clinical samples were analyzed for the presence of SARS-CoV-2 by RT-qPCR. We performed whole genome and amplicon sequencing for positive samples. We limited the description of variants in the human population to the window at which wild sample collection took place (May 1, 2022, to September 8, 2023, Fig. 3a and Supplementary Fig. 13).

Human clinical samples were collected as part of ongoing SARS-CoV-2 surveillance efforts in the region and collected under approval through Virginia Tech Institutional Review Board (IRB 20-852) and the Virginia Department of Health (VDH IRB 70046) and consent was obtained for the use of all human samples. The analysis of human clinical specimens identified 4,123 as positive, 8,350 as negative, 162 as inconclusive and 586 as invalid (no housekeeping gene detected), respectively. (Supplementary Fig. 13a). A few samples failed sequencing or were not assigned due to their quality/quantity and low sequence depth. As shown in Supplementary Fig. 13b, to further examine the sequence data from the positive opossum from 2022, various SARS-CoV-2 lineages expanded the window of our study with BA.1 (0.17%), BA.2 (15.85%), BA.2.12.1 (25.80%), BA.4 (14.83%), BA.5 (37.91%), BE (1.57%), BF (1.45%), BG.2 (0.05%), AY.103 (0.02%), not assigned (N/A, 2.24%) being the most prevalent among humans in the region.

Whole genome sequencing and Sanger amplification
Sanger sequencing was carried out using ARTIC primers (v3/v4.1) *. Samples were treated with ExoSAP-IT™ Express (Applied Biosystems) and submitted to Eurofins Genomics (Louisville, KY) for Sanger sequencing. To amplify the entire SARS-CoV-2 genome of the 2022 Virginia opossum, positive samples were reverse transcribed and amplified using SuperScript IV Reverse Transcriptase (Thermofisher) and ARTIC nCoV-2019 Amplicon Panel V4.1 primers (Integrated DNA technologies). Purified PCR products were barcoded using the plexWell™ 384 Library Preparation Kit (seqWell, MA) and the pooled library sequenced using a MiSeq System (Illumina, CA) following manufacturer’s instructions. SARS-CoV-2 sequencing analysis was accomplished using an established pipeline and optimized workflow for validation that combines tools to assess quality, sequence alignment, variation calling, and variant assignment. We processed lineage assignments using Pangolin software suite (https://cov-lineages.org/)64.

All remaining positive SARS-CoV-2 detections in 2023 by qRT-PCR were whole genome sequenced using a GridION X5 Sequencer (Oxford Nanopore Technologies, ONT) separately from the 2022 Virginia opossum. Amplifications were performed using the Lunascript RT Supermix (New England BioLabs) and Midnight RT PCR Expansion (Oxford Nanopore Technologies). Barcoding was performed using the Rapid Barcoding Kit 96 (Oxford Nanopore Technologies) before the cDNA pool was purified. The pure cDNA pool was quantified using either a Nanodrop 2000/2000c or Qubit Flex (Thermo Fisher Scientific) before loading into a flow cell (R.9.4.0 or R.9.4.1 Flow Cell Versions; Oxford Nanopore Technologies). Sequencing was performed using a GridION X5 Sequencer. Raw sequence data was base called and demultiplexed using ONT’s Guppy software, and analysis performed using ONT’s EPI2ME Labs software package running the wf-artic analysis protocol using scheme Midnight-ONT/V3 (Oxford Nanopore Technologies and ARTIC Network). The wf-artic analysis protocol utilizes Medaka, Nextflow, Nextclade, and Pangolin to produce a comprehensive analysis report that assesses sequence quality and alignment, as well as variant calling and assignment. Whole genome sequences of all human samples were deposited in the GISAID database (https://www.gisaid.org/) and NCBI SARS-CoV-2 sequence repositories.

Phylogenetic analysis
Consensus FASTA files containing available sequences from shotgun whole genome sequencing and targeted amplicon sequencing were assembled using DECIPHER software package65. These sequences were uploaded into Ultrafast Sample placement on Exiting tRee [UShER; https://genome.ucsc.edu/cgi-bin/hgPhyloPlace]66,67 for global phylogenetic tree analysis using 12,397,869 SARS-CoV-2 genomes obtained from GISAID, GenBank, COG-UK and CNCB that were available as of 2022-10-03 for the original opossum sequence and 16,284,433 SARS-CoV-2 genomes that were available as of 2023-10-12 for the other eight SARS-CoV-2 sequences obtained from wild animals. Sequences were also uploaded into NextClade [https://clades.nextstrain.org/]68 for clade assignment and mutation calling.

For visualization of these nine consensus whole genome sequences in comparison to Pango lineages, a phylogenetic tree was generated using a randomized subset of 90 SARS-CoV-2 sequences from NCBI Virus69 that represent 15 different Pango lineages (6 sequences per lineage). The sequence alignment and a Maximum Likelihood (ML) tree was generated using DECIPHER’s TreeLine function65 and viewed using iTOL v670. Each of these SARS-CoV-2 sequences collected from wild life were further analyzed for similarity to published SARS-CoV-2 sequences obtained from humans within our geographical region. We assembled separate phylogenetic trees for each Pango lineage using human SARS-CoV-2 sequences collected from Virginia and the four states that neighbor Southwest Virginia, including West Virginia, Kentucky, Tennessee, and North Carolina (Supplementary Tables 11 and Supplementary Data 4). To limit the size of some of the trees, we only included sequences collected around the same time as the dates for the wildlife sample collection. For lineages XBB.1.4.45 and JD.1, an insufficient number of sequences were available from these five states so all available sequences from North America were included for generating the phylogenetic trees. These human SARS-CoV-2 sequences were downloaded from GISAID20 and NCBI Virus. Duplicate sequences were removed prior to tree assembly. The trees were rooted using the NCBI Reference Sequence (NC_045512.2) for SARS-CoV-2. Sequence alignments and ML trees were generated as described above.

We used Nextclade’s mutation calling feature to identify novel amino acid substitutions that were unique to our wildlife samples as they were not also found in their closest neighboring sequence.

Molecular modeling
We carried out computational studies of the structure of the glycosylated Omicron BA.2 variant bound to the human ACE2 receptor [Protein Data Bank 7XO8]22. We incorporated point mutations in BA.2 using PyMOL (The PyMOL Molecular Graphics System v. 2.4.0, Schrödinger, LLC) to generate a new structure containing both E471V and G798D mutations. We then performed energy minimization using Schrödinger-Maestro (v. 2020.4) and OPLS3e force field71,72. We used Schrödinger Bioluminate to carry out molecular mechanics/generalized born surface area (MM/GBSA) energy calculations73,74 as well as surface mapping [v. 2020.4]28. In brief, the S protein and hACE complex underwent a second round of minimization using the local optimization feature, the OPLS3e force field, and adopting the Variable Dielectric Surface Generalized Born (VSGB) continuum solvation model. Prime MM-GBSA calculates the energy of each individual protein and the complex as a unit. The total predicted binding free energy for the complexes was then calculated using:

ΔG binding = G complex − (G Sprotein + G hACE2) (1)

ΔG binding = ΔE MM + ΔG GB + ΔG SA (2)

ΔE MM = ΔE Electrostatic + ΔE internal + ΔE VdW (3)

These calculations involve various energy terms, including molecular mechanics (MM), electrostatic, van der Waals, and solvation energies. ΔE MM represents the total gas phase energy in OPLS3e. ΔE internal includes bond, angle, and dihedral terms. These calculations incorporated contributions from gas phase energy, electrostatic interactions, van der Waals forces, and solvation effects via the generalized Born (GB) method and a nonpolar contribution. We conducted large-scale low-mode conformational sampling using MacroModel24,25,75. We employed the OPLS476 force field, a convergence threshold of 1.0, and an energy window for saving structures of 21.0 kJ/mol, along with a maximum atom deviation cutoff of 0.5 Å. We used the NetNGlyc 1.0 server77 to calculate the glycosylation propensity. The structure files and MM/GBSA energy components, such as Coulomb, vDW, solvent, complex, are included in the CSV files on our Open Science Framework Page (https://osf.io/82n73/) and on Zenodo (https://zenodo.org/records/11404190).

Statistical analyses
We used the “sf”78 and “exactextractr”79 packages in program R (4.1.1) to calculate the mean of the urbanization variables within each of our trap buffers. We calculated prevalence and 95% confidence intervals using the “PropCIs” package80 with the Agresti-Coull method due to small sample sizes. We used the “glmmTMB” package81 in program R to run univariate generalized linear mixed models to assess whether seroprevalence rates were impacted by the following covariates: urbanization (imperviousness), human presence, and resident population and included species as a random effect. All P-values presented represent a 2-tailed test.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Goldberg et al. evaluate the viral and serological evidence for SARS-CoV-2 infection in a total of 23 wild animal species that were caught in Virginia in 2022 and 2023. Virus was detected in oropharyngeal or nasopharyngeal swabs by RT-qPCR against 3 viral targets, providing convincing data of viral replication in six of 23 species examined, often with multiple positive individuals within each species. Further suggestion (although considered not conclusive) of wider infection among the species sampled was suggested by RT-qPCR positive results against just one of the PCR target sequences, with individuals from all six confirmed species and seven additional species showing a positive outcome. Sera was also collected from 6 species and neutralization titers against virus confirmed in individuals from 5 of these species, which suggests viral replication in two animal species that were not identified in the 2-target confirmed results of the RT-qPCR results. Analysis of the relationships between the presence of humans, termed imperviousness, and human presence, showed a positive correspondence between human presence and detection of virus in wild animals. And, further to this was sequence confirmation from about one half of the PCR positive samples confirmed multiple introductions of virus strains that were currently circulating in humans, suggesting incidence in these wild species was primarily related to cross-species transmission from humans to animals.\n\nThe data are very interesting and provide confirmation of the spread of SARS-CoV-2 into several wild animal species that were not previously known to be susceptible. This has important implications for identifying potential new reservoirs of the virus and for the potential contribute to viral evolution through adaptation in new species with potential for spread in the wild and reverse zoonosis. Specific comments regarding the manuscript are listed below.\n\nMajor comments:\n1. Seroprevalance is an important criterion to use to establish prior infection in animals, and it is probably more likely that prior infection could be detected by the presence of virus-specific antibodies than by actual detection of ongoing or very recent infection through the detection of viral RNA using RT-qPCR. For this reason is it is somewhat surprising that a neutralization assay was used to detect virus specific antibodies as it would require higher titers of antibodies specific to spike for detection than if an ELISA based assay was used. Of course species specific secondary antibodies for Indirect ELISA are unlikely to be available for most animal species studied but sensitive methods like Double Antigen Sandwich ELISA, which are commercially available, have been used successfully in similar sorts of studies on seroprevalence, and can be used to target not only spike specific antibodies but also antibodies to nucleoprotein, which are typically very abundant, to provide a more sensitive indicator of prior infection with SARS-CoV-2. Some discussion/rationale of the selection of method used to identify virus-specific should be provided.\n2. Why was no virus culture of RT-qPCR positive samples attempted? It would be very interesting to know whether there was detectable shedding of live virus since RNA levels tend to be much higher than titers of live virus and do not necessarily indicate significant live titers, at least in experimental animals. Shedding of live virus would suggest whether animals species were capable of animal to animal or animal to human spread and contribute to further understanding the risk that wild species infections pose.\n\nMinor comments\n1. References 12 are 19 are the same.\n2. Figure 1d and Methods lines 491-6. It is not clear what is meant by the levels of neutralization used to define strong or weak positive vs negative. Is this based on reduction in plaque counts with complete absence of 90% (strong +) or 60% of plaques (weak +) or is this defined by reduction in plaque sizes. This should be defined more clearly.\n3. In figure 1c, the meaning of the “Average Dilution” legend is unclear. Do the numbers 25, 50 and 75 mean 1:25, 1:50 and 1:75 dilutions to get at least 60% neutralization?\n4. Line 92-3 (Figure 1 legend) – black is mentioned twice.\n5. Figure 3. It is unclear why two different designations (Pango like XBB and NextStrain Clade like 23A) are being used together to refer to the same (at least in some instances) viruses. It would be good to at least provide some indication of how they overlap (i.e. XBB is 22F, XBB1.5 is 23A) to help with the relationship between them.\n6. Line 482. 1300PFU mixed with serum dilutions is quite a large amount of virus when attempting to design an assay that is just sensitive enough to detect neutralization. For example mixing 50-100PFU with serum dilutions and plating all of it could be more sensitive to the amount of neutralizing activity. If the assay is based on a standardized or referenced protocol this should be stated.\n7. Line 510-12, it is stated “Brush Mountain recreational area was closed to the public from March through August 2022 and thus, we assumed approximately 10 people used the area per month”. It is unclear how this estimate is derived. The estimates for other sites span the period of March to August and since Brush Mountain is closed in the same time frame, it appears 0 people would have used the area between March and August. This should be explained better.\n8. Line 520-21. How is the multiplier of 30 times the population derived in this estimation of number of people that used the area? Is this to suggest that the flow of people in and out the area is 30 times the number that live there continuously. The basis for this assumption is not clear and should be explained.", "reviewer_2": "The study sought to investigate SARS-CoV-2 spillovers in to wildlife species. The investigators used molecular detection and serological techniques to answer the questions. The detection of SARS-CoV-2 in free ranging host species other than white-tailed deer may be important to future pandemic prevention efforts. There are several studies looking a SARS-CoV-2 in many of these species, which did not result in viral detection (for which the reasons could be many). Overall, there are some substantial limitations.\n\nThe authors state that study \"greatly expanded the known host range of SARS-CoV-2,\" but this statement must be questioned. As the authors also note, there is a very real possibility that some these SARS-CoV-2 detections do not reflect active infections. All of the detections occur with rather high Ct values on their PCR testing, and only 1/3 of their samples produced sequence coverage >80% of the SARS-CoV-2 genome. Undoubtedly there was SARS-CoV-2 RNA in 23 samples, but the claim of infection is weak. The claim would be much stronger with the recovery of viable virus from these host species, but it appears as though the investigators made that impossible with their selection of viral transport media (major study design flaw). This reviewer questions why the team would go to the effort to sample all of these animals, but not use VTM suitable for viral recovery. Several places in the manuscript the authors use the terms \"isolate\" and \"isolated\" when not only was virus isolation not attempted, it would have been impossible due to their VTM.\n\nThe apparent high identity to SARS-CoV-2 circulating in humans concurrent to sample collection indicates little to no animal-to-animal transmission in these species, which means these may be dead end hosts (limits the significance of these findings).\n\nSerology: Neutralization was performed using a delta lineage virus. Given the expected life span of these species and the sample collection window, why not use omicron lineage?\n\nThe materials and methods mention testing 13,221 human clinical samples, but there seem to be no details about the origins of those samples. There is also no mention of IRB approval or exemption for use of these samples in research.\n\nI am not qualified to evaluate the phylogenetic methods.", "reviewer_3": "The article by Goldberg et al. regarding the widespread of the SARS-CoV-2 in certain communities of wildlife. In order to quantify this study the authors have explore experimentally the collection of samples between 2022-2023 and quantify the presence of SARS-CoV-2 variants in different species as well as employed molecular modelling of the RBD:ACE2 receptor to validates their findings. The adaptation of variants in different wildlife communities certainly can pose a new pathway for the evolution of SARS-CoV-2 beside the currently known host (humans). In this regard the article seems to tackle an interesting topic. Certainly I find an exhaustive work in term o the RT-qPCR detections in years (2022-2023) as well as the comparison with the same species sampled both prior to SARS-CoV-2 arrival in 2020. A correlation seems to reported.\n\nRegarding the in silico study starting in page 10, I am afraid the methodology employed and the analysis is not considering possible structural chagnes at the interface of the RBD:ACE2 and in S2 domain. Both responsible for the well-know stability of the RBD (see enhanced affinity doi.org/10.1038/s41467-021-27325-1 and additional mechanostability doi.org/10.1039/D0NR03969A). The authors have employed a common MD toolkit (i.e. Schordinger) to mutate the S protein with energy minization and not conducted equilibration via MD simulation (i.e. NVT, NPT and long microsecond equilibration). I mention this because in previous studies has been reported the crucial effect of mutation in the dynamics of proteins. Authors should present a full equilibrated MD trajectory of the S protein with E471V and G798D mutations. After that, the MM/GBSA energy calculation which is a single point technique designed for small protein set has to be validated in several snapshots/frames and provide an error in the measure of energies. At this moment there is not statistical significance in the energy values reported.\n\nIn addition, the author should deposite the structures (several frames or trajectory) of the equilibrated models (e.g. zenodo repository) with mutations and a better description of the protocol of MM/GBSA calculation for the whole S protein in present of glycans and mutations.\n\nThe authors considered the glycosilation state of the S protein in the energy calculation and they do not provide a clear role.\n\nI consider the need for a major revision with a better description of the dynamics of the S:ACE2 with those key mutations, Such study can provide a more in depth understanding.", "reviewer_4": "This study detected SARS-CoV-2 RNA and neutralizing antibodies in multiple species of wildlife in southwest Virginia. RNA was detected in deer mouse, Virginia opossum, raccoon, groundhog, Eastern cottontail rabbit, and Eastern red bat. They were able to sequence 12 SARS-CoV-2 genomes, allowing them to perform a more detailed phylogenetic analysis to examine the human-animal interface. All sequenced viruses were found to be different lineages of Omicron and appeared to be recent introductions from humans (possible mouse-to-mouse transmission). The detection of SARS-CoV-2 in so many new wildlife species is notable and should encourage more surveillance efforts targeting a wider range of species. The phylogenetic analysis, however, was not well presented and difficult to interpret.\n\n1. The use of different PANGO and NextStrain nomenclatures, even in side-by-side panels (e.g., Figure 3), is problematic. PANGO lineages are used more commonly than NextStrain in the scientific literature, so better to go with that one.\n2. The phylogenetic trees are unreadable. What do all the numbers mean? Please include strain names. Panel 3c is a nice summary, but where do those inferences come from? I’m confused by the question mark for the source of the BA.2.10.1 introduction into opossum. Why is it difficult to resolve the host species source for this introduction? What evidence is there that it is not from humans? Can you show the most closely related strains in a phylogenetic tree? Were the wildlife SARS-CoV2 sequences always most closely related to human SARS-CoV2 sequences from Virginia? Or sometimes another state? And label them by host species? And please include strain names in these trees. How were these datasets selected? Are these the most closely related viruses from GISAID to each wildlife strain as identified in UsHER? What is the outgroup? (It is also common practice to provide the raw tree files in a repository like GitHub.)\n3. What methods were used to generate these trees? Maximum likelihood? What substitution models were used? Was a bootstrap analysis performed? NextStrain is a great exploratory tool for getting a snapshot look at data, but it's not designed to produce publication quality trees. IQTree can generate fast ML trees; or time-scale trees in BEAST can also provide estimates of when human-animal transferred occurred.\n3. 12 SARS-CoV2 sequences were generated from wildlife for this study, but the supplementary table listing accession numbers only includes 9 viruses.\nIsolate name Animal species Accession numbers\nV0632866 Didelphis virginiana OR866905\nV0654090 Didelphis virginiana OR878666\nV0654124 Marmota monax OR866349\nV0654196 Peromyscus maniculatus OR866382\nV0654535 Peromyscus maniculatus OR866443\nV0654612 Peromyscus maniculatus OR878668\nV0655023 Peromyscus maniculatus OR866910\nV0654185 Procyon lotor OR878667\nV0654414 Sylvilagus floridanus OR866437\n4. For the serology, could you perform a sensitivity analysis for different thresholds? How many samples would be considered positive if the threshold was lowered to 50%?\n5. It would be nice to have a summary table with a row for each wildlife species tested and columns for (a) total number of animals tested and whether any animals were positive for SARS-CoV-2 by (b) serology; (c) PCR; (d) WGS; (e) partial sequence.\n6. Is it possible to tell from the serology whether the animal was infected with omicron versus an earlier variant? If animals were sampled at different time points, could you infer anything about the duration of antibody responses?\n7. Was the rate of SARS-CoV-2 detection higher at wildlife rehabilitation centers? How soon after capture was testing performed?\n8. References 7 and 45 are the same."}, "author_response": "Response: We appreciate this suggestion. While neutralization assays are possibly less sensitive than an ELISA-based approach, they are also much more specific since non-neutralizing antibodies are highly cross-reactive between coronaviruses (https://www.thelancet.com/journals/ebiom/article/PIIS2352-3964(21)00494-1/fulltext). Surprisingly, false positives have also been observed in binding-based assays in patients with other infections (https://www.nature.com/articles/s41598-022-26709-7) or even auto-immune diseases (https://www.frontiersin.org/articles/10.3389/fimmu.2021.666114/full). Thus, while there are potential benefits to using ELISA, neutralization tests still provide the most specific means of assessing prior exposure to SARS-CoV-2. We have added this rationale into the discussion.\n\n2. Why was no virus culture of RT-qPCR positive samples attempted? It would be very interesting to know whether there was detectable shedding of live virus since RNA levels tend to be much higher than titers of live virus and do not necessarily indicate significant live titers, at least in experimental animals. Shedding of live virus would suggest whether animal species were capable of animal to animal or animal to human spread and contribute to further understanding the risk that wild species infections pose.\n\nResponse: We did consider collecting additional samples to try and isolate live virus from any animals that tested positive. However, we ended up using a viral transport media that rendered any virus non-viable because of the added health and safety challenges for collecting and transporting potentially live virus and the subsequent time and money investment in processing of these samples. Given that samples were also collected from field locations, there were also additional considerations for preservation, transport, and maintenance of samples from some sites to the lab. Since we initially expected quite low SARS-CoV-2 prevalence, we selected an approach that optimized along these axes (e.g. health, safety, logistics, and costs). However, we agree with the reviewer that this would be helpful in the future and are working on a way we can collect a subset of samples that we could collect live virus from going forward.\n\nMinor comments\n1. References 12 are 19 are the same.\n\nResponse: Thank you very much for catching this. We have corrected the mistake.\n\n2. Figure 1d and Methods lines 491-6. It is not clear what is meant by the levels of neutralization used to define strong or weak positive vs negative. Is this based on reduction in plaque counts with complete absence of 90% (strong +) or 60% of plaques (weak +) or is this defined by reduction in plaque sizes. This should be defined more clearly.\n\nResponse: The level of neutralization is based on the reduction of plaque counts and not the plaque size. We have updated the methods to clarify our analysis.\n\n3. In figure 1c, the meaning of the “Average Dilution” legend is unclear. Do the numbers 25, 50 and 75 mean 1:25, 1:50 and 1:75 dilutions to get at least 60% neutralization?\n\nResponse: We apologize for any confusion. This is referring to the percent neutralization of a sample, which is indicated by the size of the circle. We have revised the figure text to more accurately reflect this and have added clarification to the figure legend.\n\n4. Line 92-3 (Figure 1 legend) – black is mentioned twice.\n\nResponse: Thank you very much for catching this. We have removed the duplicate word.\n\n5. Figure 3. It is unclear why two different designations (Pango like XBB and NextStrain Clade like 23A) are being used together to refer to the same (at least in some instances) viruses. It would be good to at least provide some indication of how they overlap (i.e. XBB is 22F, XBB1.5 is 23A) to help with the relationship between them.\n\nResponse: Thank you very much for pointing out that area of confusion. We agree and have removed NextStrain clades and now just use Pango lineages to make it easier to relate each of the panels within Figure 3 to one another and to address another reviewer’s comments.\n\n6. Line 482. 1300PFU mixed with serum dilutions is quite a large amount of virus when attempting to design an assay that is just sensitive enough to detect neutralization. For example, mixing 50-100PFU with serum dilutions and plating all of it could be more sensitive to the amount of neutralizing activity. If the assay is based on a standardized or referenced protocol this should be stated.\n\nResponse: We generated a solution of SARS-CoV-2 at a concentration of 1300 PFU/mL, which resulted in us adding ~40 PFUs per well, just below the range you suggested. We have added additional text to clarify this in the methods. This amount is consistent with various other reports (https://www.nature.com/articles/s41467-020-17892-0 and https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001284)\n\n7. Line 510-12, it is stated “Brush Mountain recreational area was closed to the public from March through August 2022 and thus, we assumed approximately 10 people used the area per month”. It is unclear how this estimate is derived. The estimates for other sites span the period of March to August and since Brush Mountain is closed in the same time frame, it appears 0 people would have used the area between March and August. This should be explained better.\n\nResponse: We have added some clarifying text to the methods to better explain how we came up with this estimate. All the estimates were on the same time scale: # of persons/month averaged over March – August 2022. We have added more detail about how we derived the estimate of 10 people/month for Brush Mountain. There was one person utilizing the area occasionally to remove trees for future trail work but they were only on the property a few days a month (personal communication with the property manager). We then assumed a few other visits by locals who may occasionally walk through the area (observed one person while working at the site) for recreation despite the fact that it was technically closed. Hence, we estimated about 10 visits per month on average between March and August 2022.\n\n8. Line 520-21. How is the multiplier of 30 times the population derived in this estimation of number of people that used the area? Is this to suggest that the flow of people in and out the area is 30 times the number that live there continuously. The basis for this assumption is not clear and should be explained.\n\nResponse: We apologize for the confusion. We are estimating the number of people using an area per month. As a result, we are assuming that the population is relatively constant each day, so we multiplied it by 30 to represent 30 days of people occupying an area. We agree there likely are some days where more people are in town (weekend) or fewer (weekdays), but we are working under the assumption that the average monthly usage is about 30 times the population. We have added additional text to the methods to make this more clear.\n\nReviewer #2 (Remarks to the Author):\n\n1. The authors state that study \"greatly expanded the known host range of SARS-CoV-2,\" but this statement must be questioned. As the authors also note, there is a very real possibility that some these SARS-CoV-2 detections do not reflect active infections. All of the detections occur with rather high Ct values on their PCR testing, and only 1/3 of their samples produced sequence coverage >80% of the SARS-CoV-2 genome. Undoubtedly there was SARS-CoV-2 RNA in 23 samples, but the claim of infection is weak. The claim would be much stronger with the recovery of viable virus from these host species, but it appears as though the investigators made that impossible with their selection of viral transport media (major study design flaw). This reviewer questions why the team would go to the effort to sample all of these animals, but not use VTM suitable for viral recovery. Several places in the manuscript the authors use the terms \"isolate\" and \"isolated\" when not only was virus isolation not attempted, it would have been impossible due to their VTM.\n\nResponse: We apologize for the misuse of “isolate/isolated” in the manuscript. We have removed the 4 instances where this was not applicable and have replaced them with more appropriate terminology (collected).\n\nWe thank the reviewer for the comments concerning live virus collection and agree that there are limitations to the interpretation of qPCR positive detections, which we have included throughout the manuscript. However, given the specificity of our qPCR, we can definitively conclude that SARS-CoV-2 is being detected in a larger range of hosts than previously reported, with important implications for the ecology of this disease and ongoing global surveillance efforts. While live virus was not isolated, we want to highlight that this would still be steps away from demonstrating animal-to-animal transmission, which would require experimental infections and subsequent transmission experiments, which was outside of the scope of our current project. We did discuss and consider collecting additional samples to try and isolate live virus from any animals that tested positive (see response #2 to reviewer 1). Ultimately, we decided not to go this route because of the added health and safety restrictions for collecting and transporting potentially live viruses from wild animals, the subsequent investment in processing of these samples. In addition, we would face limitations in interpreting the genomic data from virus that was isolated from animals and propagated in culture, as we would be unable to determine whether mutations emerged as result of that procedure. We think our findings still represent a significant result, given that many species are being exposed to SARS-CoV-2 in the wild. It demonstrates that human-to-animal transmission routes exist for a wide array of species, and there is high probability these species will eventually be exposed to a lineage that is capable of infecting them. This could result in viral shedding and subsequent transmission, if it has not already (see response below). Ultimately, we agree with the reviewer and are working on a way we can collect even a subset of samples so we can isolate live virus going forward.\n\n2. The apparent high identity to SARS-CoV-2 circulating in humans concurrent to sample collection indicates little to no animal-to-animal transmission in these species, which means these may be dead end hosts (limits the significance of these findings).\n\nResponse: We agree that our findings indicate a large number of human-to-animal transmission events but have no evidence to suggest these species are dead-end hosts, particularly because several of them have been shown to be capable of transmission relevant shedding rates in the lab (e.g. deer mice; Table S9). In addition, while some hosts might not be less important for transmission, wildlife can still suffer disease impacts, as has been observed in mink, and thus documenting SARS-CoV-2 exposure in wild mammals has conservation relevance. Lastly, as mentioned above, if exposure continues at the rate we observed, there is a high probability that some of these species will eventually be exposed to a lineage that can cause animal-to-animal transmission. We also did find two instances where the phylogeny is unresolved (an opossum and deer mouse), which could indicate animal-to-animal transmission. We tried to be abundantly conservative in our interpretation of these data but given the information we have and experimental work on these species, animal-to-animal transmission may be the most parsimonious explanation.\n\n3. Serology: Neutralization was performed using a delta lineage virus. Given the expected life span of these species and the sample collection window, why not use omicron lineage?\n\nResponse: We collected our serum samples in 2022, which was just after the end of the Delta wave so we assumed many of the neutralizing antibodies in these captured animals would likely reflect prior exposure to the Delta variant. Given the life expectancy of the species we were sampling, they would have extended well into the period that the Delta variant was circulating. We have added the table below to the supplementary materials (Supplementary Table 12), and additional justification into the methods.\n\nSpecies Average lifespan in the wild (years); max is in parenthesis\nProcyon lotor 5a (20b)\nDidelphis virginiana 1.5-2a\nMephitis mephitis <1a (6b)\nSciurus carolinensis (12.5b)\nPeromyscus maniculatus <1a\nPeromyscus leucopus 1a\nBlarina brevicauda (2.5b)\nMarmota monax 4-6a\nLasiurus borealis\nSylvilagus floridanus <3a (5b)\nTamias striatus <2a (3b)\nVulpes vulpes 3a (7b)\na Myers, P., R. Espinosa, C. S. Parr, T. Jones, G. S. Hammond, and T. A. Dewey. 2024. The Animal Diversity Web (online). Accessed at https://animaldiversity.org\nbCarey, J. R., & Judge, D. S. 2002. Longevity records: Life spans of mammals, birds, amphibians, reptiles, and fish. Monographs on population aging, (8).\n\n4. The materials and methods mention testing 13,221 human clinical samples, but there seem to be no details about the origins of those samples. There is also no mention of IRB approval or exemption for use of these samples in research.\n\nResponse: We have added additional details into the methods about the human clinical samples and the IRB that this sample collection was conducted under. Samples were collected by the Virginia Department of Health and submitted to the Molecular Diagnostics Laboratory as part of a genomic surveillance agreement.\n\nReviewer #3 (Remarks to the Author):\n\nThe authors appreciate the reviewer’s comments and recognition of the work presented in this article by our group.\n\n1. Regarding the in-silico study starting in page 10, I am afraid the methodology employed, and the analysis is not considering possible structural changes at the interface of the RBD:ACE2 and in S2 domain. Both responsible for the well-known stability of the RBD (see enhanced affinity doi.org/10.1038/s41467-021-27325-1 and additional mechanostability doi.org/10.1039/D0NR03969A).\n\nResponse: The authors acknowledge the importance of considering possible structural changes and dynamics at the interface of RBD:ACE2 and in the S2 domain, as they are relevant to the stability of RBD. In response to the reviewer’s suggestion, we have now included contextualization and discussion of the relevant literature regarding the dynamics of the RBD in relation to ACE2 binding.\n\nThe new text now reads as follow:\n“While this approach limits the investigation of these mutations and their role on dynamic conformational states of S protein, it lays the foundation for understanding the mechanistic impact and connection of structure to experimental observables. Future work will involve molecular dynamics (MD) simulations to assess the impact of these mutations on structural stability, membrane interaction, and conformational states. This research will contribute to our understanding of evolutionarily pressures and cross-species transmission.“\n\n2. The authors have employed a common MD toolkit (i.e. Schordinger) to mutate the S protein with energy minization and not conducted equilibration via MD simulation (i.e. NVT, NPT and long microsecond equilibration). I mention this because in previous studies has been reported the crucial effect of mutation in the dynamics of proteins. Authors should present a full equilibrated MD trajectory of the S protein with E471V and G798D mutations.\n\nResponse: We appreciate your suggestion regarding the equilibration process and thank the reviewer for thoughtful comments. While we utilize the Schrodinger toolkit for mutation and energy minimization, we acknowledge the importance of presenting a fully equilibrated MD trajectory for a comprehensive understanding of the dynamics, especially in light of previous studies emphasizing the effect of mutations. However, performing MD simulations for over 1 microsecond on multiple mutations within the complete Spike:ACE2 complex to fully describe the dynamic events is not feasible within a reasonable time frame (it would have taken roughly 6 to 10 months using all the computing resources at Virginia Tech based on our calculations) due to the scale of the entire system (1 million+ atoms) and the scope of this study. Moreover, it is not suitable for this work to only simulate the clipped RBD + S protein, as including the S2 region is essential for clarity and comprehensiveness in a thorough MD investigation. Further MD simulation studies would be better suited as standalone papers to address the combined impact of mutation and provide a comprehensive discussion of the dynamic state of Spike bound to ACE2. In this work, our goal was to offer context and initial rationalization of the location and interaction impact between the Spike RBD and ACE2, focusing on the E471V and local structural impact of the G798D mutations. Indeed, for E471V, the local point-based mutations had an impact on the interaction energy, highlighting a potential enthalpic contribution of this mutation on the RBD interface interacting with ACE2. We have incorporated this additional context into the manuscript and will carefully consider your recommendation in a follow-up work.\n\nThe text now reads as follow:\n“We conducted in silico studies to probe the local impact of the E471V and G798D mutations on the S protein. We employed molecular modeling and molecular mechanics/generalized borne surface area (MM/GBSA) free energy calculations to predict the free energy of binding interaction (DG ). It is important to highlight that our focus in this study is on probing the local impact of E471V and G798D on the favorability of S protein-hACE2 interaction and the modifications in the S2 site as hypothesis for the experimentally observed traits found in SARS-CoV-2 wildlife reservoir.”\n“The increased favorability of the E471V mutation for hACE2 interaction aligns with the previously reported trend of residues in the RBD influencing stabilization and enhanced affinity (Koehler et al. 2021 Nat Comm). Any differences may be attributed to entropic effects, emphasizing the potential enthalpic contribution of this mutation to improve hACE2 interaction at the RBD interface.”\n“Future work will involve molecular dynamics (MD) simulations to evaluate the impact of these mutations on structural stability, membrane interaction, and conformational states. This research will contribute to our understanding of evolutionarily pressures and cross-species transmission.”\n\n3. After that, the MM/GBSA energy calculation which is a single point technique designed for small protein set has to be validated in several snapshots/frames and provide an error in the measure of energies. At this moment there is not statistical significance in the energy values reported.\n\nResponse: To generate multiple conformations for statistical analysis, we conducted large-scale low-mode conformational sampling on the RBD (residues 329-531) of both BA.2 and E471V Spike in complex with hACE2 (Watts et al. 2014 J Chem Info). To reduce computational expense, Spike was truncated to include only the RBD. The top five sampled conformations of E471V-hACE2, based on energy, exhibited significantly more favorable binding interaction energy than BA.2-hACE2. This information has been included in the computational section for its statistical value.\n\nThe text now reads as follows:\n“Large-scale low-mode conformational sampling was performed on the RBD of BA.2 and E471V Spike (residues 329-531) with hACE2 to probe ΔG across multiple potential conformations (Watts et al. 2014 J Chem Info, Kolossvarv 1996 J Amer Chem Soc). Consistent with the results observed in full-length Spike-hACE2 interactions, the mean ΔG for the top five most favorable E471V conformations (-151.7 ± 19.2 kcal/mol) was more favorable than the mean for BA.2 conformations (125.4 ± 4.1 kcal/mol, p0.05). The sampled conformations of E471V exhibited a larger range of ΔG values relative to BA.2 , reflecting not only increased interaction flexibility but also increased flexibility in the RBD due to the E471V mutation.”\n\nThe methods section has been updated to reflect these changes:\n“We conducted large-scale low-mode conformational sampling using MacroModel. The OPLS4 force field, a convergence threshold of 1.0, and an energy window for saving structures of 21.0 kJ/mol, along with a maximum atom deviation cutoff of 0.5 Å, were employed.”\n\n4. In addition, the author should deposit the structures (several frames or trajectory) of the equilibrated models (e.g. zenodo repository) with mutations and a better description of the protocol of MM/GBSA calculation for the whole S protein in present of glycans and mutations.\n\nResponse: At the reviewer’s request, we have deposited the structures and data files, along with detailed information on the MM/GBSA calculation protocol for the entire S protein in the presence of glycans and mutations. These files are accessible on our Open Science Framework repository at https://osf.io/82n73/. Additionally, we have expanded the methods section to provide further clarity on the conducted MM/GBSA calculation procedure.\n\nThe text now reads as follows:\n“In brief, the S protein and hACE2 complex underwent a second round of minimization using the local optimization feature, the OPLS3e force field, and adopting the Variable Dielectric Surface Generalized Born (VSGB) continuum solvation model. Prime MM-GBSA calculates the energy of each individual protein and the complex as a unit. The total predicted binding free energy for the complexes was then calculated using:\nΔG = G –(G +G )\nbinding complex S protein hACE2\nΔG = ΔE + ΔG + ΔG\nbinding MM GB SA\nΔE = ΔE + ΔE + ΔE\nMM Electrostatic internal VdW\nThese calculations involve various energy terms, including molecular mechanics (MM), electrostatic, van der Waals, and solvation energies. ΔE represents the total gas phase MM energy in OPLS3e. ΔE includes bond, angle, and dihedral terms. These calculations incorporated contributions from gas phase energy, electrostatic interactions, van der Waals forces, and solvation effects via the generalized Born (GB) method and a nonpolar contribution.“\n\n5. The authors considered the glycosilation state of the S protein in the energy calculation and they do not provide a clear role.\n\nResponse: We appreciate the reviewer’s attention to the glycosylation state of the S protein in our energy calculation. We have now extended our explanation regarding the role of glycosylation, particularly at the S2 site. The inclusion of glycosylation sites in our structural investigation aligns with previous work (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8260498/) and provides an important context for understanding structure-based modifications. We trust this additional clarification enhances the comprehensiveness of our study.\n\nThe text now reads as follows:\n“When glycosylated, the N801 glycosylation site significantly influences the reduction in viral entry, with very low mutation rates observed (e.g., within 2-3 residues from position 801). Interestingly, D798 creates a small, charged pocket on this solvent-accessible S2 loop, decreasing the probability of glycosylation at N801 (0.48) compared to BA.2 [0.61]29. It is hypothesized that G798D could impact structural stability and membrane interaction30 (Fig. 4 and Fig. S10) and have a negative impact on Spike maturation and infectivity. While this approach limits the investigation of these mutations and their role in the dynamic conformational states of the S protein, it provides a foundation for understanding the mechanistic impact and connection of structure to experimental observables.”\n“Glycosylation propensity was calculated using the NetNGlyc 1.0 server68. The structure files and MM/GBSA energy components, such as Coulomb, vDW, solvent, complex, are included in the CSV files on our Open Science Framework Page (https://osf.io/82n73/).”\n\n6. I consider the need for a major revision with a better description of the dynamics of the S:ACE2 with those key mutations, such study can provide a more in depth understanding.\n\nResponse: We thank the reviewer for their thoughtful and thorough review of the manuscript and their reconsideration for major revisions.\n\nReviewer #4 (Remarks to the Author):\n\n1. The use of different PANGO and NextStrain nomenclatures, even in side-by-side panels (e.g., Figure 3), is problematic. PANGO lineages are used more commonly than NextStrain in the scientific literature, so better to go with that one.\n\nResponse: We appreciate this comment and we have changed Figure 3 to only present the Pango lineages and reduce any confusion.\n\n2a. The phylogenetic trees are unreadable. What do all the numbers mean? Please include strain names.\n\nResponse: We apologize that the trees were too small to read. We have saved them in a taller format so they can be stretched to fit a page and hopefully make them easier to read. We have also uploaded the main phylogenetic figure (Fig. 3) as a separate manuscript file, and included strain names, state of collection, and year when available for easier interpretation.\n\n2b. Panel 3c is a nice summary, but where do those inferences come from? I’m confused by the question mark for the source of the BA.2.10.1 introduction into opossum. Why is it difficult to resolve the host species source for this introduction? What evidence is there that it is not from humans?\n\nResponse: For Figure 3c, the inferences about transmission were determined by phylogenetic analysis using local, related human host sequences and any of our sequences collected from wild animals as shown in Supplementary Figures 2-9 and in Supplementary Table 3. We also considered any unique mutations as those mutations which may have come from the originating host’s virus.\n\nFor the two events where we included question marks, these were based on two different scenarios. For the opossum (BA.2.10.1), we reported a mutation in the sample collected that has not been identified in any human sample collected. This of course could be from a human that was never sampled, arisen in the opossum itself, or occurred in another animal. Given the uncertainty, we designed this as a question mark. For mice, we had two samples collected from mice (Fig. S8) that clustered together. Again, we cautiously interpreted this as either both mice were exposed to the same infected human, or the mouse was exposed to the other mouse. Given that deer mice have been shown experimentally to be capable of transmitting, we felt either of these scenarios could be likely. We have added additional text to the Fig. 3 legend to clarify how we determined whether the source was likely human or unknown.\n\n2c. Can you show the most closely related strains in a phylogenetic tree?\n\nResponse: We have included all the most closely related strains in individual phylogenetic trees in the supplemental information (Fig. S2-S9). The two mice described above are included in the same tree. All of the SARS-CoV-2 sequences used for phylogenetic analysis were retrieved from human hosts except for our sequences collected from wild animals.\n\n2d. Were the wildlife SARS-CoV2 sequences always most closely related to human SARS-CoV2 sequences from Virginia? Or sometimes another state?\n\nResponse: In some cases, yes. For example, the most closely matched human sample to the two deer mice mentioned above was from Vienna, VA. However, the opossum with the unique mutation has the closest match to a human sample collected from NY. We have now added the state abbreviation to all the branch tip names to aid in seeing the origin location. We also have the closest human sample and state of origin in Table S3.\n\n2e. And label them by host species? And please include strain names in these trees.\n\nResponse: Phylogenetic trees have been modified to include strain names and host species are now indicated with color tips.\n\n2f. How were these datasets selected? Are these the most closely related viruses from GISAID to each wildlife strain as identified in UsHER? What is the outgroup? (It is also common practice to provide the raw tree files in a repository like GitHub.)\n\nResponse: All SARS-CoV-2 sequences used for phylogenetic analysis were retrieved from human hosts except for our sequences collected from wild animals. Selection of the datasets used for our phylogenetic trees are described in the Methods section lines 596-600 and 605-615. UShER was not used to select the sequences only to identify which PANGO lineages to use for which virus. All of our trees used the Reference Sequences for SARS-CoV-2 from NCBI (NC_045512.2). We appreciate the suggestion to submit the sequences used for phylogenetic analysis to GitHub for public availability. However, due to the use of GISAID sequences which are only available for registered users, we’ve listed the sequences for our analyses in Supplementary Table 11.\n\n3. What methods were used to generate these trees? Maximum likelihood? What substitution models were used? Was a bootstrap analysis performed? NextStrain is a great exploratory tool for getting a snapshot look at data, but it's not designed to produce publication quality trees. IQTree can generate fast ML trees; or time-scale trees in BEAST can also provide estimates of when human-animal transferred occurred.\n\nResponse: Our phylogenetic trees are Maximum Likelihood (ML) trees generated using TreeLine (Methods line 601) in R. TreeLine selected the substitution model based on Akaike information criterion for each tree. Usually, GTR+G4 was selected as the substitution model. We have added “maximum likelihood” and the substitution model used to each phylogenetic tree’s figure legend and to Supplementary Table 10. Boot strap analysis was not performed. We agree with the reviewer that NexStrain is a great tool for exploring data initially and so it wasn’t used for tree generation in this manuscript. We used TreeLine to generate the tree and iTOL to view the trees (Methods line 601-602). We will explore the suggested programs (IQTree and BEAST) for our future phylogenetic needs to compare them to TreeLine.\n\n3. 12 SARS-CoV2 sequences were generated from wildlife for this study, but the supplementary table listing accession numbers only includes 9 viruses.\n\nResponse: We apologize for any confusion the remaining 3 virus samples are in a separate Table (S4) since we were only able to get partial sequence. Accession numbers have been added to both Supplementary Tables 3 and 4.\n\n4. For the serology, could you perform a sensitivity analysis for different thresholds? How many samples would be considered positive if the threshold was lowered to 50%?\n\nResponse: We re-ran the generalized linear mixed model using different cutoffs for percent neutralization (40, 50, 65, and 70), which is now available in Supplementary Table 6. We find statistical support for this relationship at all cut-offs except for the most conservative (70% neutralization). However, even at this higher cutoff, we continue to find a negative relationship between these two variables. Additionally, we have included the differences in percent positivity using different threshold cutoffs, and the sensitivity of our analysis examining the relationship between imperviousness and seroprevalence using these cutoffs in the table below.\n\n% Neutralization Intercept Imperviousness slope\ncutoff chosen β SE p β SE p\n40% -0.054 0.372 0.884 0.053 0.024 0.026*\n50% -0.608 0.380 0.109 0.042 0.019 0.024*\n60% -1.655 0.477 0.001 0.038 0.018 0.032*\n65% -1.892 0.514 0.000 0.044 0.019 0.017*\n70% -1.885 0.561 0.001 0.005 0.179 0.858\n\n5. It would be nice to have a summary table with a row for each wildlife species tested and columns for (a) total number of animals tested and whether any animals were positive for SARS-CoV-2 by (b) serology; (c) PCR; (d) WGS; (e) partial sequence.\n\nResponse: We agree and have added a table (Supplemental Table 7) that summarizes all of this information. Thank you for this suggestion.\n\n6. Is it possible to tell from the serology whether the animal was infected with omicron versus an earlier variant? If animals were sampled at different time points, could you infer anything about the duration of antibody responses?\n\nResponse: This is a great question. Previous published work has suggested that certain variants can be distinguished based on the strength of neutralization, and this is something we are eager to understand in the future. We are hoping to obtain enough recaptured animals where we know what variant they were previously exposed to, so we can look at this in known individuals. Unfortunately, the serology data in this study was collected over a very brief period of time ~6 weeks, so we were not able to examine this question in this manuscript. We do hope to be able to address this question in future manuscripts when we will have completed regular sampling over a longer period of time with recaptured individuals.\n\n7. Was the rate of SARS-CoV-2 detection higher at wildlife rehabilitation centers? How soon after capture was testing performed?\n\nResponse: Thank you for that interesting question. The rehabilitation centers had 2.24% positivity across all species while we had 4.04% in our field collected animals. We will note that species like mice are not seen by wildlife rehabilitation centers which is the primary contributing factor to the differences. If we examine positivity rate by species with the greatest number of samples collected in both datasets: we find the following:\n\nSpecies Rehabilitation Field data\nVirginia opossum 2.83% 2.94%\nRaccoon 3.39% 8.00%\n\nIt does suggest we have slightly higher positivity rates in the field compared to the rehabilitation clinics for raccoons. We would have hypothesized the opposite – animals are typically being brought to clinics because they are sick or injured. However, our field sampling is generally more intensive in a given area whereas submissions to rehab centers are a bit more diffuse. As you can see if Figure 3a, most of the positive samples were collected during fall 2023, and we had multiple positives at a site during a four-day trapping period. This means if there was a change in human infection rates or activity that results in more exposure in wildlife, we are more likely to get multiple positives at that time. We agree with the reviewer that these are important points for considering surveillance going forward and have added some text about this in the manuscript. For the timing of our sampling, animals were trapped overnight and then processed in the morning, for the rehabilitation centers, animals are brought to the center (typically within 24 hr. following encounters) and then processed upon arrival.\n\n8. References 7 and 45 are the same.\n\nResponse: We have corrected this error."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have assessed the responses of the authors to my original review of the manuscript. The authors had adequately addressed all comments and I have no further comments for consideration.", "reviewer_2": "My comments have been addressed to the extent possible given the study design.", "reviewer_3": "In the revised article I find some improvements but additional in silico work has to be done to improve the message of the MS.\n\nMajor comments;\nHaving snapshot of the complex system that have been only energy minimized are not sufficient to prove what they author suggest in the MS. I have suggested to perform common MD protocols with NVT and NPT equilibration runs, typical each of few nanoseconds. Regarding the dynamics of the whole S:ACE2 complex with mutations, it is important to consider time scale of microseconds. The authors described the difficulties in terms of computational cost. However, I believe there was a misunderstanding regarding the time scale, yes microsecond but one can consider several parallel runs each of 200 ns. In this way, one can have 5 independent replicas, better statistics and obtain some dynamics of the possible conformations induced by those mutations. As the author mentioned, it will take 6-10 months to runs 1 microsecond in Virgitina Tech HPC, well, in practice it may take less than a 1-2 month considering replica study. The suggested parallel study is common in MD.\n\nI consider this computational study crucial, as in the section regarding free energy calculations. the authors now report values of free energy after perform large-scale low-mode conformational sampling which is another form of the simple NMA for large structures. The problem I can see is to have applied to a complex and not a single chain. For the case of the full S:ACE2 complex considering the mutations. The values are very contradictory as now they posses error bars and the absolute values are quite different from each other: E471V-RBD conformations (-151.7 ± 19.2 kcal/mol) and for BA.2-RBD conformations (125.4 ± 4.1 kcal/mol). Here the problem is reflected in the free energy difference, which is ~-300 kcal/mol between both variants, this is totally off. I believe this is the results of having conformations from NMA which are not consistent with the true conformational landscape of the protein complex.\n\nIf the authors can provide some biophysical insight into the conformations via MD dynamics, it will simply enrich the whole message. I totally undestand the need to use MD, but it has to be carefully used.", "reviewer_4": "TreeLine is not a standard program for phylogenetic (only 41 citations compared to >6500 for IQTree2). I suspect because the trees are not so well resolved, with many polytomies like in Figure 3 and not able to run bootstraps. But this is not a phylogenetics paper and the tree in Figure 3 is not necessarily wrong. The addition of strain names is a great improvement."}, "author_response": "Reviewer 3, Comment #1: “In the revised article I find some improvements but additional in silico work has to be done to improve the message of the MS.”\n\nAuthors’ response: The authors thank the reviewer for taking the time to review our manuscript and are pleased to learn that improvements were found.\n\nReviewer 3, Comment #2: “Having snapshot of the complex system that have been only energy minimized are not sufficient to prove what they author suggest in the MS. I have suggested to perform common MD protocols with NVT and NPT equilibration runs, typical each of few nanoseconds. Regarding the dynamics of the whole S:ACE2 complex with mutations, it is important to consider time scale of microseconds. The authors described the difficulties in terms of computational cost. However, I believe there was a misunderstanding regarding the time scale, yes microsecond but one can consider several parallel runs each of 200 ns. In this way, one can have 5 independent replicas, better statistics and obtain some dynamics of the possible conformations induced by those mutations. As the author mentioned, it will take 6-10 months to runs 1 microsecond in Virginia Tech HPC, well, in practice it may take less than a 1-2 month considering replica study. The suggested parallel study is common in MD.”\n\nAuthors’ response: The authors understand the reviewer’s concern regarding the importance of conformational changes in protein-protein binding. The new suggestion from the reviewer, to conduct replicates of shorter time, differs from the initial request for continuous microsecond trajectories (Reviewer’s #3 comment from initial submission). Capturing large-scale conformational changes in such a large system would likely require continuous simulation for over a microsecond, as originally indicated by the reviewer. Substituting shorter simulations to achieve the same aggregate time would not allow for the observation of the same effects. The suggested approach with replicates would allow for local conformational changes, but most of the time would likely be spent reaching equilibrium rather than true equilibrium sampling.\nPrevious literature on the ensemble-based simulation of the partial S protein trimer complex has been performed for 7.5 us microseconds, a task that demanded considerable effort and resources (80 P100 GPUs) (Casalino L, et al. AI-driven multiscale simulations illuminate mechanisms of SARS-CoV-2 spike dynamics. Int J High Perform Comput Appl. 2021 Sep;35(5):432-451. doi: 10.1177/10943420211006452). Substantial simulation time is necessary to observe S protein RBD conformational dynamics, particularly concerning the impact of mutations at the interface. Simply simulating a single full chain of the S protein would overstate the observed flexibility in the trimer and further would not offer new insights into the relevant binding interface, as suggested by the simulations.\n\nFinally, we consulted with an independent MD simulations expert, Dr. Justin Lemkul, who currently serves as a developer for CHARMM, one of the most popular MD simulation engines. Although, Dr. Lemkul did not contribute directly to this study, he acknowledged that our assertions are justified and fall within the interpretative resolution of the presented data. While the reviewer’s suggestion emphasizes thoroughness, it’s crucial to consider whether such simulations align with the primary focus of this project and whether extensive MD simulations would significantly affect the outcome. We have acknowledged Dr. Lemkul’s expertise in the Acknowledgments section of the manuscript.\n\nReviewer 3, Comment #3: “the authors now report values of free energy after perform large-scale low-mode conformational sampling which is another form of the simple NMA for large structures. The problem I can see is to have applied to a complex and not a single chain. For the case of the full S:ACE2 complex considering the mutations. The values are very contradictory as now they posses error bars and the absolute values are quite different from each other: E471V-RBD conformations (-151.7 ± 19.2 kcal/mol) and for BA.2-RBD conformations (125.4 ± 4.1 kcal/mol). Here the problem is reflected in the free energy difference, which is ~-300 kcal/mol between both variants, this is totally off. I believe this is the results of having conformations from NMA which are not consistent with the true conformational landscape of the protein complex.”\n\nAuthors’ response: The authors thank the reviewer for identifying this error and sorry that it caused confusion. There was a typographical mistake in the paragraph (Line 243) discussing the large-scale low-mode conformational sampling (please refer to the snapshot of the text below with the typo highlighted). The difference between WT and mutation is only ~25 kcal/mol, consistent with the other numerical values. This error has been corrected in the revised version of the manuscript.\n\nReviewer 4, Comment #1: “TreeLine is not a standard program for phylogenetic (only 41 citations compared to >6500 for IQTree2). I suspect because the trees are not so well resolved, with many polytomies like in Figure 3 and not able to run bootstraps. But this is not a phylogenetics paper and the tree in Figure 3 is not necessarily wrong.”\n\nAuthors’ response: The authors thank the reviewer for their positive feedback. We acknowledge the reviewer’s suggestion to try IQTree2 and explore the comparison between tree topologies.\n\nFirst, the TreeLine function in the DECIPHER package is a relatively new tool for generating phylogenetic trees (Published in April 2022), so it is not expected to have as many publications as more established programs like IQ-Tree (Published in January 2015). We compared the two trees generated using TreeLine (previously in the paper) and IQ-Tree of the same alignment of the sequences used in Figure 3 were remarkably similar to each other including the many polytomies (see output images included at the end of this document). For the IQ-tree the line thickness indicates the bootstrap value. Thicker lines have higher boot strap values. Circled in red are the 3 isolates (highlighted in red boxes below) of 100 samples that were found in different locations between the two trees and each had low bootstrap values. Thus, the trees are of comparable quality. The authors believe that the likely many polytomies seen are due to the sequences lacking enough information to figure out how those they are related. Our wild animal sample (V0654185) that \"moved\" is a sequence with ~25% coverage so it lacks sufficient coverage for optimal phylogenetic analysis."}, {"round": 3, "reviewer_comments": {"reviewer_3": "I see the authors corrected the confusing values of the free energy difference. I would simply advice the authors to provide the information of results available via zenodo.org. The current link https://osf.io/82n73/ connected to the website of the authors does not allow for proper verification of the MD data."}, "author_response": "We thank the reviewer for their suggestion. We have uploaded the results to zenodo.org as requested and can be found at the link (https://zenodo.org/records/11404190) and DOI 10.5281/zenodo.11404189. These have been added to the data availability statement as well."}]

Decision letter