Summary
The emergence of SARS-CoV in 2002 and SARS-CoV-2 in 2019 led to increased sampling of sarbecoviruses circulating in horseshoe bats. Employing phylogenetic inference while accounting for recombination of bat sarbecoviruses, we find that the closest-inferred bat virus ancestors of SARS-CoV and SARS-CoV-2 existed less than a decade prior to their emergence in humans. Phylogeographic analyses show bat sarbecoviruses traveled at rates approximating their horseshoe bat hosts and circulated in Asia for millennia. We find that the direct ancestors of SARS-CoV and SARS-CoV-2 are unlikely to have reached their respective sites of emergence via dispersal in the bat reservoir alone, supporting interactions with intermediate hosts through wildlife trade playing a role in zoonotic spillover. These results can guide future sampling efforts and demonstrate that viral genomic regions extremely closely related to SARS-CoV and SARS-CoV-2 were circulating in horseshoe bats, confirming their importance as the reservoir species for SARS viruses.
In brief
Recombination-aware evolutionary analyses of the entire genomes of SARS-CoV-1-like and SARS-CoV-2-like viruses indicate that SARS-CoV-1 and SARS-CoV-2 descend from bat coronaviruses that circulated as recently as one to six years before their respective dates of human emergence. These inferred virus ancestors likely needed to travel faster than closely related bat viruses in order to reach their respective sites of emergence.
Graphical Abstract
Introduction
Horseshoe bats (Rhinolophus spp.) are the main hosts of the Sarbecovirus subgenus1 (genus Betacoronavirus, family Coronaviridae) and, likely via virus transmission through transient intermediate species, were the source of SARS-CoV (referred to hereafter as SARS-CoV-1, now extinct) and SARS-CoV-22–4. Since the first emergence of SARS-CoV-1 in Guangzhou, Guangdong Province in 2002 and SARS-CoV-2 (jointly referred to as the SARS-CoVs) in Wuhan, Hubei Province in 2019, the sampling and analyses of sarbecoviruses in horseshoe bats has contributed to our understanding of sarbecovirus diversity and their geographical spread. Although sarbecoviruses exhibit substantial geographic structuring linked to the ranges of their host horseshoe bat species2–6, where and when their ancestors circulated, both in recent and ancient history7, is poorly understood.
Genome-wide sequence identity is typically used to compare bat sarbecoviruses to the SARS-CoVs, but because coronaviruses frequently recombine5,6, whole-genome identity of these bat viruses and the SARS-CoVs does not adequately reflect their complex evolutionary histories (Fig. S1a). Previous attempts to date the most recent common ancestor (MRCA) of SARS-CoV-2 and a closely related bat virus genome using long genomic fragments resulted in estimates of 30–40 years before 20194,6,8. However, when analyzing all non-recombinant regions (NRRs), we will inevitably find some NRRs with older and some with younger MRCAs than either this estimate from a long fragment or an estimate based on the entire genome, which would effectively be a weighted average of NRR time of the MRCAs (tMRCAs). To find a bat virus genome that has not undergone recombination relative to either of the SARS-CoVs, it would likely have been necessary to identify the sarbecovirus involved in each cross-species transmission event leading to the SARS-CoVs (which we refer to as the “direct ancestor”; Fig. S1b) or a very closely related virus from a bat in the same ecological community infected around the time of the transmission event. What we have, instead, is a relatively small sample of bat sarbecovirus genome sequences, gathered across different times and places. Hence, it is necessary to identify all the NRRs of sarbecoviruses — genomic regions within which there is little to no detectable signal of recombination among the analyzed genomes and representing, to the extent possible, a single evolutionary history. These genomic regions provide the clearest insights into the evolutionary history of SARS-CoV-1 and SARS-CoV-2 and the origins of their respective outbreaks.
Here, we separately analyze the recombination patterns and respective evolutionary histories of SARS-CoV-1 and its closely related viruses (referred to as SARS-CoV-1-like viruses) and SARS-CoV-2 and its closely related viruses (referred to as SARS-CoV-2-like viruses). We estimate the separate sets of NRRs for the SARS-CoV-1-like and SARS-CoV-2-like viruses. For each NRR we utilize a molecular clock rate specific for that portion of the genome and determine the MRCA of the SARS-CoV and the most closely related bat virus, or set of viruses, among currently sampled sarbecoviruses for the NRR. We refer to this MRCA as the “closest-inferred ancestor” for that NRR because it is the parent node of the SARS-CoV in the phylogeny (Fig. S1b) and is limited by the current sample set of bat sarbecoviruses. Additionally, the dating of deeper nodes in the phylogeny is complicated by substitution saturation: the repeated occurrence of nucleotide substitutions at the same site resulting in biased estimation of older divergence in viral phylogenies9,10. We therefore account for substitution saturation in our NRR-specific inferences when examining the deeper evolutionary history of sarbecoviruses.
We show that the available sequence data provide evidence of the closest-inferred bat virus ancestors of both SARS-CoV-1 and SARS-CoV-2 circulating in bats in the years immediately preceding the emergence of each SARS-CoV. By performing phylogeographic analyses of these sarbecoviruses, we demonstrate that the sarbecoviruses typically spread at a rate that approximates their bat hosts. Based on this rate and current sampling, we find that both SARS-CoV-1 and SARS-CoV-2 likely emerged far from where their direct bat virus ancestors likely circulated and were therefore unlikely to reach their provinces of emergence via bat host movement alone.
Results
Inference of non-recombinant regions
The evolutionary history of sarbecoviruses can only be accurately inferred when focusing on non-recombinant regions rather than entire genomes or specific genes because these viruses, as well as their ancestors, have undergone extensive recombination 4,6,8. To identify these recombination patterns, we aligned the 250 full-genome sarbecovirus sequences that were available at the time of this analysis, including SARS-CoV-1 and SARS-CoV-2. We used GARD11 to independently identify recombination breakpoints on genomes closely related to SARS-CoV-1 and SARS-CoV-2 (Fig. S2). We then performed separate analyses to reduce the potentially overlapping recombination signal between the sets and increase our power for confidently detecting breakpoint positions unique to each set. The two sequence sets were separated based on previously published recombination-aware phylogenetic analyses of most genomes used here4,12,13 (see Methods). The two recombination breakpoint analyses identified 31 putative NRRs for the SARS-CoV-1-like viruses, with a median length of 870 nucleotides (nt) and a range of 309 to 3,414 nt, and 44 putative NRRs for the SARS-CoV-2-like viruses, with a median length of 595 nt and a range of 168 to 2,706 nt (Fig. 1a).
Figure 1. Recombination patterns and dating of the closest-inferred bat virus ancestor of SARS-CoV-1 and SARS-CoV-2 for each NRR.
a, Distribution of SARS-CoV-1-like and SARS-CoV-2-like recombination breakpoints across the sarbecovirus genome alignment (32,975 nt). Major open-reading frames and notable gene domains are displayed. Breakpoint pairs with significantly short distances (falling below the 10% tail of the observed distance distribution) are highlighted in yellow. Numbers indicate NRRs delineated from inferred recombination breakpoints. b, Time of the closest-inferred bat virus ancestor of SARS-CoV-1 across the 31 NRRs of the SARS-CoV-1-like viruses. The 20 years before the emergence of SARS-CoV-1 is highlighted in gray in the upper panel and shown in the lower panel. c, Time of the closest-inferred bat virus ancestor of SARS-CoV-2 across the 44 NRRs of the SARS-CoV-2-like viruses. The 20 years before the emergence of SARS-CoV-2 is highlighted in gray in the upper panel and shown in the lower panel. The dashed line in each panel is the date of the earliest sampled human case of each SARS-CoV (16 November 2002 for SARS-CoV-1; 24 December 2019 for SARS-CoV-2). The violin plots correspond to the 95% highest posterior density (HPD) intervals and the dot within each violin plot indicates the median. The location of the Spike (S) gene is indicated in the top panels of b-c.
To examine whether the inferred positions of recombination breakpoints of the SARS-CoV-1-like and SARS-CoV-2-like NRR analyses were closer to each other than expected, we simulated randomly distributed recombination breakpoint sets along the alignment for the SARS-CoV-1-like and SARS-CoV-2-like viruses (30 and 43 breakpoints, respectively) and compared the distribution of distances between the two breakpoint sets. The empirical breakpoint sets had a mean genomic distance of 203 nt (90% CIs: 138–299 nt) separating the pairs of closest SARS-CoV-1-like and SARS-CoV-2-like virus breakpoints, whereas the simulated distances had a mean of 396 nt (90% CIs: 259–606 nt; Fig. S3). Additionally, there were seven breakpoint pairs between the two virus sets that were in the smallest 10% of the observed distribution of distances between breakpoints (10 to 51 nt apart; Fig. 1a), indicating recombination hotspots shared between the two sarbecovirus groups (e.g., near nsp3 and between the Spike NTD and RBD).
Recent bat virus ancestors of human SARS-CoVs
Focusing on each NRR identified in the recombination analysis, we examined when the closest-inferred bat viruses ancestral to SARS-CoV-1 and SARS-CoV-2 circulated. These ancestors are represented by the MRCA of each SARS-CoV and the most closely related bat viruses, which is akin to the parent node of each SARS-CoV in its respective phylogeny (Fig. S1b).
There is insufficient temporal signal when calibrating a molecular clock using tip dating with sarbecoviruses sampled from bats and pangolins (Manis javanica)6, likely as a consequence of limited sampling across space and time. Therefore, we used SARS-CoV-1 genomes to identify a suitable rate prior for the molecular clock in our time-scaled phylogenetic analyses (Fig. S2; see Methods), and these inferred substitution rates served as NRR-specific rate priors for subsequent Bayesian phylogenetic inference of a SARS-CoV-1 genome from 2002 and 139 SARS-CoV-1-like viruses (Fig. S4a, Table S1). The median time of the closest-inferred bat virus ancestor of SARS-CoV-1 within each NRR varied by several decades (Fig. 1b), and the median of these medians is 1996. However, the NRR with the most recent closest-inferred ancestor of SARS-CoV-1 had an inferred date of 2001 (95% highest posterior density [HPD]: 1998–2002; NRR 14 [Orf1b]; 1,065 nt), only one year before the emergence of SARS-CoV-1 in 2002. The viruses most closely related to SARS-CoV-1 within this NRR were YNLF_31C and YNLF_34C, sampled in Yunnan, China in 2013. Our findings were consistent, albeit slightly earlier, when we included additional human and masked palm civet (Paguma larvata) SARS-CoV-1 genomes from both the 2002–2003 and 2003–2004 outbreaks (Fig. S4c, Table S1).
We similarly used the evolutionary rate of SARS-CoV-2 across the NRRs of the SARS-CoV-2-like viruses (Fig. S4b, Table S1), inferred from SARS-CoV-2 genomes sampled up to late 2020, as priors for the phylogenetic analysis of the NRRs of the SARS-CoV-2-like viruses, comprising SARS-CoV-2 and 37 SARS-CoV-2-like genomes. The time of the closest-inferred bat virus ancestor of SARS-CoV-2 ranged from several years to several decades preceding late 2019 — when SARS-CoV-2 was introduced into humans — across the 44 NRRs (Fig. 1c). Although the median time of the closest-inferred ancestor of SARS-CoV-2 was in 2000, the two NRRs with the most recent time of the closest-inferred ancestor both had an inferred date of 2014 (95% HPDs: 2002–2019 and 2007–2018; NRR 37 [Orf3a and E] and 43 [N]; 347 and 331 nt, respectively), only five years prior to its introduction into humans14. The genome most closely related to SARS-CoV-2 within NRR 37 is RacCS203, sampled in Thailand in 2020, whereas the genome most closely related to SARS-CoV-2 within NRR 43 is RpYN06, sampled in Yunnan in 2021. We note that although these two NRRs are each less than 350 nt, NRR 42 is 456 nt and has a similar time of the closest-inferred ancestor of 2013 (95% HPD: 1999–2019). Within NRR 42, the genomes most closely related to SARS-CoV-2 are BANAL-20–52 and BANAL-20–236, both sampled in Laos in 2020. When using a faster rate prior based on SARS-CoV-2 genomes sampled in early 2020 (Fig. S4b, Table S1), we obtain slightly more recent estimates of the closest-inferred ancestor (Fig. S4d, Table S1). Our inferred ages for the closest-inferred ancestors of each of the SARS-CoVs across the NRRs indicate that several of the published sarbecovirus genomes include non-recombinant fragments that are descendant from viruses that circulated likely only several years before the emergence of SARS-CoV-1 and SARS-CoV-2.
The times of closest-inferred ancestor of SARS-CoV-1 across the NRRs are nearer to the date of human emergence than those for SARS-CoV-2 (Fig. 1b,c): 21 NRRs of the SARS-CoV-1-like viruses have a median time of the closest-inferred ancestor of SARS-CoV-1 within 10 years of the emergence of SARS-CoV-1, whereas only 7 NRRs of the SARS-CoV-2-like viruses have a closest-inferred ancestor of SARS-CoV-2 within 10 years of the emergence of SARS-CoV-2. This pattern could be in part explained by the sampling of a substantially larger number of SARS-CoV-1-like genomes in bats compared to SARS-CoV-2. However, using a faster rate prior based on SARS-CoV-2 genomes sampled in early 2020 results in 17 NRRs having a closest-inferred ancestor of SARS-CoV-2 within 10 years of its emergence (Fig. S4d). Both SARS-CoV-1-like and SARS-CoV-2-like viruses consistently have older closest-inferred ancestors across the Spike gene (Fig. 1b,c). The oldest time of the closest-inferred bat virus ancestor within a single NRR of SARS-CoV-1 and SARS-CoV-2 was, respectively, 1963 (95% HPD: 1934–1985; NRR 19 [Orf1b]) and 1944 (95% HPD: 1798–2009; NRR 28 [Spike]). Notably, only 1 NRR among the SARS-CoV-1-like viruses (NRR 19) and 6 NRRs among the SARS-CoV-2-like viruses (NRRs 10, 27, 28, 31, 32, and 38) have a median time of the closest-inferred ancestor more than 50 years earlier than the most recently sampled SARS-CoV-1-like or SARS-CoV-2-like virus, respectively.
Effect of increased sampling on the recombinant common ancestors
To test the influence of sampling, we quantified the extent to which the publication of more viruses closely related to SARS-CoV-1 and SARS-CoV-2 has contributed toward our understanding of the ancestor of each of these human viruses. The respective recombinant common ancestor (“recCA”) of SARS-CoV-1 and SARS-CoV-2 can be understood as the aggregate of the closest-inferred ancestors of each of the SARS-CoVs across each of its NRRs14, with each closest-inferred ancestor existing at different points in time and with potentially different viruses most closely related to the SARS-CoV. The recCA therefore accounts for the closest relative(s) across all non-recombinant segments. We inferred the genome of the recCA of each SARS-CoV over time, progressively adding published genomes to the reconstruction based on their publication date. Then, we calculated the similarity of each recCA to its respective SARS-CoV as a function of the publication date to examine how additional closely related genomes can affect ancestral reconstruction and reflect our increasing understanding of sarbecoviruses over time (Fig. 2a).
Figure 2. The recCA (recombinant common ancestor) over time.
The genetic similarity of the complete recCA and the most closely related sarbecovirus to each SARS-CoV over time (a,b) and the genetic similarity of the NRRs of the recCA to each SARS-CoV over time (c,d). a, The similarity of the recCA of SARS-CoV-1 and the most closely related sarbecovirus in the SARS-CoV-1 viruses to SARS-CoV-1 as a function of time. b, The similarity of the recCA of SARS-CoV-2 and the most closely related sarbecovirus in the SARS-CoV-2 viruses to SARS-CoV-2 as a function of time. The right panels for a and b are zoomed-in panels of the dashed box in the left panels. The dashed vertical line in b is the sampling date of the earliest sampled genome of SARS-CoV-2. c, The similarity of each NRR of the recCA of SARS-CoV-1 to SARS-CoV-1 as a function of time. d, The similarity of each NRR of the recCA of SARS-CoV-2 to SARS-CoV-2 as a function of time. Each square in (c, d) represents the similarity of the NRR of the recCA to its relative SARS-CoV at one point in time and is colored by date of publication of the most recent non-human sarbecovirus in the reconstruction.
We observe that, as more genomes were published, the SARS-CoV-1 recCA and SARS-CoV-2 recCA became more similar to SARS-CoV-1 and SARS-CoV-2, respectively (Fig. 2a,b). The full bat virus genomes most similar to SARS-CoV-1 and SARS-CoV-2 are, respectively, YN2020E (96.2% genetic identity to SARS-CoV-1) and BANAL-20–52 (96.8% genetic identity to SARS-CoV-2), both published in 2021. By the end of 2022, the recCA of SARS-CoV-2 shares 98.8% genetic identity with SARS-CoV-2 and the recCA of SARS-CoV-1 shares 98.7% genetic identity with SARS-CoV-1, both expectedly exceeding the genetic identity of the bat virus genomes most genetically similar to the human SARS-CoVs. Notably, the recCA repeatedly increases in similarity to its respective SARS-CoV even when newly published genomes share less genetic identity with the recCA than the most similar genome at the time (Fig. 2a,b). This trend can be observed in detail at the level of individual NRRs: sequentially published genomes increase the genetic identity of some of the NRRs of the recCA to the matching regions of its SARS-CoV (Fig. 2c,d). The shared genetic identity of the recCA to its respective SARS-CoV is broadly consistent when incorporating additional SARS-CoV-1 genomes with the SARS-CoV-1-like viruses (Fig. S5a,c) and when using a faster rate prior in the Bayesian time-calibrated phylogenetic inference of the SARS-CoV-2-like viruses (Fig. S5b,d).
The patterns observed here are the result of newly added genomes containing fragments that descend from a more closely related ancestor despite not being the most genetically similar overall, and further highlight the importance of inferring NRRs in evolutionary analyses. Additionally, the similar genetic identity of the recCAs to their respective SARS-CoVs, both genome-wide and often at the level of individual NRRs, suggest that the greater number of sampled SARS-CoV-1-like viruses than SARS-CoV-2-like viruses is not solely responsible for the times of closest-inferred ancestor of SARS-CoV-1 being consistently nearer to the date of human emergence than those for SARS-CoV-2.
Phylogeography of bat sarbecoviruses
Sarbecovirus NRR phylogenies have an appreciable degree of geographic structuring linked to their hosts’ ranges6. To understand the spatiotemporal spread of the SARS-CoV-1-like and SARS-CoV-2-like viruses in bats, we simultaneously fitted a spatially-explicit phylogeographic model to each NRR. As the sampling locations of SARS-CoV-1, SARS-CoV-2, and the pangolin sarbecoviruses likely do not represent where their direct bat virus ancestors circulated15–18, we excluded their locations from the phylogeographic analyses to avoid the impact of dispersal of non-bat hosts (Fig. S2).
Due to substitution saturation, the deeper evolutionary history of rapidly evolving viruses, including sarbecoviruses, are underestimated using standard molecular clock approaches7,10,19–21. To account for the effects of substitution saturation on the SARS-CoV-1-like and SARS-CoV-2-like phylogenies when trying to understand the long-term dispersal patterns of sarbecoviruses, we applied the Prisoner of War (PoW) model21 — a molecular clock model that estimates divergence times while accounting for time-dependent evolutionary rates — to the posterior phylogenies with phylogeographic estimates scaled in units of genetic distance (Fig. S2; see Methods). We find that substitution saturation primarily affected the inferred time of divergence events more than 50 years before the most recently sampled virus in our dataset, with younger internal nodes having relatively unchanged ages and indicating that our inferences of the time of the ancestors of human SARS-CoVs and their closest bat sarbecoviruses are unbiased (Fig. 3a-b). The effects of substitution saturation can, however, be particularly more pronounced for older divergence estimates (Fig. 3a-b). When examining the entire phylogenies, the tMRCAs were often hundreds of years old when using a standard clock model, whereas the PoW-transformed phylogenies had tMRCAs that were typically thousands to tens of thousands of years older (Fig. 3c-d).
Figure 3. Time to the most recent common ancestor (tMRCA) after PoW transformation.
The number of years of each internal node before the most recently sampled taxon compared to the number of years of the same node after applying a PoW transformation for (a) 100 SARS-CoV-1-like posterior phylogenies for each NRR and (b) 100 SARS-CoV-2-like posterior phylogenies for each NRR. Red dashed line at 50 years before the most recently sampled taxon indicates a threshold beyond which saturation effects start to substantially influence tMRCA estimates. Diagonal dotted line depicts a line of equality, with a one-to-one correspondence between the two inferences on the two axes. Estimates of time to the most recent common ancestor for the (c) SARS-CoV-1-like viruses and (d) SARS-CoV-2-like viruses for each NRR when using a standard molecular clock or the PoW model. Each square indicates the median tMRCA and the lines indicate the 95% HPD. Instances where confidence intervals of the tMRCA overlap when using a standard molecular clock and a PoW model suggest that the age of the bat virus ancestors are young enough as to where the standard molecular clock accounts for most of the substitution saturation.
To investigate the dispersal history of the two sets of viruses, we mapped the PoW-transformed phylogeny of each NRR across Asia (Fig. 4g,h, Data S2). The phylogeographic reconstructions of both sets of viruses suggest a few lineages having dispersed over long distances, but the number of these dispersals is likely an underestimate due to undersampling of sarbecoviruses from various regions in Asia. The inferred dispersal history of viral lineages looks consistent among NRRs (Fig. 4g,h, Data S2): we observe a similar spatial distribution of the sets of viruses and their connections on the different NRR maps. Moreover, the virus lineages are spatially structured, with phylogenetically similar viruses being generally sampled in geographically similar regions with almost no back and forth travel of viral lineages across relatively large regions of the study area, which we confirmed by computing the isolation-by-distance (IBD) signal 22 for each NRR (Table S2). Our analyses indeed reveal a moderate to strong IBD signal when computing the correlation between the patristic distance — the sum of the branch lengths connecting two nodes on a tree — and great-circle geographic distance computed for each pair of bat virus samples: a Spearman correlation coefficient of 0.45 for SARS-CoV-1-like viruses (95% HPD: 0.36–0.60) and of 0.67 for SARS-CoV-2-like viruses (95% HPD: 0.29–0.84) (see Table S2 for a comparison with other viruses).
Figure 4. Estimated evolutionary history of SARS-CoV-1-like (left) and SARS-CoV-2-like (right) viruses.
a, Global map with the study area highlighted in dark gray. b, Horseshoe bat species richness with the scale bar showing the number of co-occurring horseshoe bat species per km2. c-d, time-scaled maximum clade credibility (MCC) phylogenies of individual NRRs (NRR14 for SARS-CoV-1-like viruses, and NRR37 for SARS-CoV-2-like viruses) obtained after correction for substitution saturation, shown on a log-transformed time scale. Horizontal segments displayed at each internal node indicate the 95% HPD interval of the estimated node age, with a light grey vertical bar indicating the median node age. The closest-inferred ancestors of SARS-CoV-1 and SARS-CoV-2 are indicated by red arrows. e-f, weighted diffusion coefficient estimates based on phylogenetic branches occurring less than 75 (plain light gray curves) and 50 (dashed curves) years ago. g-h, reconstruction of the dispersal history of viral lineages for two NRRs (NRR14 for SARS-CoV-1-like viruses, and NRR37 for SARS-CoV-2-like viruses) obtained through continuous phylogeographic inference (see Data S2 for all other NRRs). Nodes are coloured from yellow (oldest most recent common ancestor) to blue (most recent collection date) and are superimposed on 80% HPD polygons coloured according to the same log-transformed time scale, reflecting the uncertainty of the Bayesian phylogeographic inference. 80% HPD polygons were only computed and reported for the last 1,000 years. Dispersal direction (anti-clockwise) of viral lineages is indicated by the edge curvature. In panels c-d-g-h, tip nodes are displayed as squares and are associated with a number indicating the host species from which the virus was sampled; internal nodes are displayed as circles. The host species are numbered by decreasing order of sampling frequency. i-j, estimated position of the closest-inferred bat virus ancestors of SARS-CoV-1 (i) and SARS-CoV-2 (j) based on the continuous phylogeographic reconstruction of all NRRs. The 95%, 75% and 50% HPD regions are displayed with an increasing red color darkness. Note that we use the PoW-transformed trees to illustrate the overall time-scale of the evolutionary histories, but refer to Figure 1 and Table S1 for estimates of the time of the closest-inferred ancestor.
Location of the closest-inferred ancestors of the human SARS-CoVs
Determining where the recent ancestors of SARS-CoV-1 and SARS-CoV-2 circulated could inform future sampling efforts and elucidate how these virus lineages eventually traveled to their provinces of emergence. As the PoW model enforces contemporaneously sampled tips, it can produce unreliable results for the most recent inferred ages. However, the PoW transformation shows that the timing of divergence events in the last ~50 years inferred from standard molecular clock models are unbiased (Fig. 3a-b). Therefore, we performed a separate tip-dated phylogeographic analysis to infer where the closest-inferred ancestors of each SARS-CoV circulated in the past 50 years (Fig. S2).
We find that the closest-inferred ancestors of SARS-CoV-1 likely circulated in western China (Yunnan, Sichuan, or Guizhou) (Fig. 4i). The closest-inferred ancestors of SARS-CoV-2 likely circulated in Yunnan, China or northern Laos, overlapping with contiguous karst and cave landscapes extending through these regions23,24 (Fig. 4j). Our results indicate that both SARS-CoV-1 and SARS-CoV-2 emerged in humans over a thousand kilometers from where their closest-inferred bat virus ancestors likely circulated, with SARS-CoV-1 emerging in Guangzhou and SARS-CoV-2 in Wuhan. The inferred geographic location of the closest-inferred ancestors of each SARS-CoV overlapped with regions with moderate horseshoe bat species richness (Fig. 4b and Data S2), including in regions known to harbor R. affinis, R. pusillus and R. ferrumequinum (Data S2), and the closest-inferred ancestor of SARS-CoV-2 additionally overlapped with known locations of R. malayanus circulation (Data S2). Although sarbecoviruses have been found in bats from Guangdong and Hubei13,25,26, the viruses most genetically similar to SARS-CoV-1 and SARS-CoV-2, both at the genome level and among individual NRRs, have been predominantly sampled in and near Southwest China27,28. Therefore, our findings are not solely a result of localized sampling efforts but reflect a natural pattern in the geographic distribution of bat sarbecoviruses, though more sampling in areas proximal to the Laos and Vietnam borders may identify additional CoVs descendant from relatively recent ancestors circulating in bats.
Sarbecoviruses spread at rates that approximate bat host dispersal
To investigate how fast the SARS-CoV-1-like and SARS-CoV-2-like viruses were recently spreading across Asia, we inferred the weighted diffusion coefficient29, a measure of the diffusivity of viral transmission, for the 50 years before the most recently sampled genome in each dataset (Fig. S2). We found that the SARS-CoV-1-like viruses spread with a weighted diffusion coefficient of 1,666 km2/yr (95% HPD: 811–3,399; Fig. 4e), and the SARS-CoV-2-like viruses circulated with a weighted diffusion coefficient of 740 km2/yr (95% HPD: 116–2,587; Fig. 4f). The SARS-CoV-1-like viruses were likely inferred to spread, on average, slightly more quickly than the SARS-CoV-2-like viruses because of the recent dispersal events to South Korea. These results were similar to the diffusion coefficients calculated for the 75 years preceding the most recently sampled genomes (Table S2). Horseshoe bats, the primary hosts of sarbecoviruses, tend to forage within around 2–3 km of their roost30–33 and have a reported diffusion coefficient of approximately 2,000 km2/yr (ref.34). Therefore, both the SARS-CoV-1-like and SARS-CoV-2-like viruses (in the context of host infections) traveled at a roughly similar rate as their hosts, suggesting that these viruses are able to efficiently explore new niches, likely due to the density of horseshoe bat roosts within the karst landscapes of Southern China.
Rapid movement of lineages leading to human SARS-CoVs
To test whether the bat virus ancestors of each human SARS-CoV could have spread to where the SARS-CoVs emerged in humans via dispersal of SARS-related coronaviruses across bat populations, we determined the minimum distance from the location of the closest-inferred bat virus ancestor to the province of emergence (Fig. 4i-j), and given the time between emergence and this ancestor (Fig. 1b-c), we calculated the dispersal velocity this distance implies. We compared the implied dispersal velocity for the tip branch leading to SARS-CoV-2 to the dispersal velocity of bat virus tip branches to obtain a posterior rank distribution (Fig. S6). The tip branches leading to SARS-CoV-1 and SARS-CoV-2 had the highest mean posterior dispersal velocity rank (Fig. 5a), but there were other tip branches leading to SARS-CoV-1-like viruses with consistently high posterior ranks as well. Strikingly, the posterior rank distributions for the tip branches of SARS-CoV-1 and SARS-CoV-2 exhibit the most extreme skewness away from a rank of 0.5 and towards 1 (Fig. 5b; see Methods), with their posterior rank distributions rarely including values below 0.5 (Fig. S6). This pattern further indicates the rapid nature of the dispersal velocities associated with the tip branches leading to SARS-CoV-1 and SARS-CoV-2.
Figure 5. Dispersal velocity rank and skewness of SARS-CoV-1-like and SARS-CoV-2-like viruses.
a, Mean posterior dispersal velocity rank of each tip branch in the SARS-CoV-1-like and SARS-CoV-2-like phylogenies across all posterior samples and NRRs. Dashed colored lines indicate the mean of all tip branch posterior dispersal velocity ranks. b, Skewness away from a rank of 0.5 of the posterior dispersal velocity rank distribution (Fig. S6) for each tip branch. Dotted line at a rank skewness of 0 indicates expected average rank. c, Mean posterior dispersal velocity rank of each tip branch in the SARS-CoV-1-like (top) and SARS-CoV-2-like (bottom) phylogenies across all posterior samples grouped by NRR. A star indicates that the SARS-CoV has the greatest mean posterior dispersal velocity rank for that NRR. Additional SARS-CoV-1-like viruses with fast dispersal velocities are also colored and identified by their strain name. NRRs were excluded if the closest-inferred ancestor of SARS-CoV-1 or SARS-CoV-2 had a tMRCA more than 50 years before the most recently sampled virus for at least 50% of the posterior samples (see Methods).
We next explored the contributions of individual NRRs to the relatively fast dispersal velocities (and, therefore, high posterior ranks) required for the closest-inferred ancestors of SARS-CoV-1 and SARS-CoV-2 to reach their respective provinces of emergence. The tip branches leading to SARS-CoV-1 and SARS-CoV-2 have the greatest mean posterior rank in, respectively, 12 and 28 of the NRRs, indicating that it is not only NRRs with recent closest-inferred ancestors of the SARS-CoVs (Fig. 1) that result in relatively high dispersal velocities for the tip branches leading to the SARS-CoVs. Considering the large distance separating the closest-inferred bat virus ancestors from where SARS-CoV-1 and SARS-CoV-2 each emerged (Fig. 4i-j) and the high dispersal velocities necessary to traverse this distance (Fig. 5), it is unlikely that the lineages descending from the closest-inferred bat virus ancestors of SARS-CoV-1 and SARS-CoV-2 reached Guangdong and Hubei, respectively, solely via dispersal through their bat reservoirs.
Phylogeography robust to hypothetical unsampled diversity
We evaluated whether our findings could be affected by the hypothetical discovery of viruses closely related to the human SARS-CoVs nearer to their provinces of emergence than has been observed in the empirical data. We conducted tip-dated phylogeographic analyses that incorporated hypothetical unsampled viruses, which represent unsequenced and unobserved viral diversity. These unsampled taxa were positioned between the provinces of emergence of SARS-CoV-1 and SARS-CoV-2 and the locations where the taxa most genetically similar to them were sampled, while also enforcing monophyly between the unsampled taxa and the respective SARS-CoV (Fig. S7; see Methods), thereby ensuring that these viruses were more phylogenetically similar to SARS-CoV-1 and SARS-CoV-2 across the entire genome.
The inclusion of these unsampled taxa resulted in a shift in the inferred location of the closest-inferred ancestor toward the province of emergence of the respective SARS-CoV, although, notably, the 50% and 75% HPDs were consistent (Fig. S8). There were slight reductions in the dispersal velocities necessary for the closest-inferred ancestors of SARS-CoV-1 and SARS-CoV-2 to reach their provinces of emergence (Fig. S9), primarily due to the shift in the inferred location of the closest-inferred ancestor. However, the tip branches leading to SARS-CoV-1 and SARS-CoV-2 still had the highest mean posterior dispersal velocity rank on average and across many NRRs, as well as the most positive skewness (Fig. 6). Other measures to quantify differences in dispersal velocity also generally support the relatively high ranks of the tip branches leading to SARS-CoV-1 and SARS-CoV-2, but are more sensitive to potentially unsampled diversity, in particular for SARS-CoV-1 (Fig. S10).
Figure 6. Dispersal velocity rank and skewness of SARS-CoV-1-like and SARS-CoV-2-like viruses when including hypothetical unsampled taxa.
a, Mean posterior dispersal velocity rank of each tip branch in the SARS-CoV-1-like and SARS-CoV-2-like phylogenies across all posterior samples and NRRs. Dashed colored lines indicate the mean of all tip branch posterior dispersal velocity ranks. b, Skewness away from a rank of 0.5 of the posterior dispersal velocity rank distribution for each tip branch. c, Mean posterior dispersal velocity rank of each tip branch in the SARS-CoV-1-like (top) and SARS-CoV-2-like (bottom) phylogenies across all posterior samples grouped by NRR. A star indicates that the SARS-CoV has the greatest mean posterior dispersal velocity rank for that NRR. Additional SARS-CoV-1-like viruses with fast dispersal velocities are also colored and identified by their strain name. NRRs were excluded if the closest-inferred ancestor of SARS-CoV-1 or SARS-CoV-2 had a tMRCA more than 50 years before the most recently sampled virus for at least 50% of the posterior samples (see Methods).
Discussion
By reconstructing the evolutionary history of sarbecoviruses whilst accounting for recombination, we show that the ancestors of SARS-CoV-1 and SARS-CoV-2 likely circulated in Horseshoe bat populations hundreds to thousands of kilometers away from the sites of the emergence of these viruses in humans and as recently as one to six years prior to this emergence. Our findings indicate that there would not have been sufficient time for the direct bat virus ancestor to reach the locations of emergence of the human SARS-CoVs via normal dispersal through bat populations alone.
Sarbecoviruses undergo extensive recombination and tend to share breakpoints at specific genome regions, so-called recombination hotspots4,35. Comparisons of NNR subgenomic segments when less data was available previously suggested there may have been decades of separation between the closest-inferred bat virus ancestor and each of the SARS-CoVs6,36. However, dating the closest-inferred ancestor across non-recombinant subgenomic regions spanning the entire genome after calibrating a molecular clock on a per region level reveals much more recent ancestry between the SARS-CoVs and their bat-circulating relatives. The non-recombinant fragments, therefore, must be comprehensively analyzed to properly understand the emergence of the human SARS-CoVs5,6. We additionally show here that adding new genomes to the recombination-aware ancestral reconstruction brought the sequence identity of the recombinant common ancestor of each SARS-CoV to its respective human virus to greater than 98% genetic identity.
The phylogeographic approach implemented here provides important insights into the sarbecoviruses’ extensive historical movement across China and Southeast Asia. Our analysis confirms that both the SARS-CoV-1-like and SARS-CoV-2-like viruses have each been circulating in these regions for at least tens, if not hundreds, of thousands of years. Conversely, the bat species endemic to Japan known to harbor sarbecoviruses, R. cornutus, and its phylogenetically closest mainland species infected by SARS-CoV-2-like viruses, R. pusillus, are estimated to have shared a common ancestor 13.3 million years ago37. This much more ancient separation between the host species indicates that the viruses did not strictly co-diverge with these bat species. We find that, in more recent history, the viral lineages diffuse at a similar rate to horseshoe bats. These bats typically have small home ranges (2–3 km2 movement per night on average), are rarely migratory32,38, and have a limited dispersal capacity with scant evidence of longer migrations, with such evidence limited to species in temperate rather than tropical regions39,40. After accounting for recombination, the closest-inferred bat virus ancestors of SARS-CoV-1 and SARS-CoV-2 are inferred to have circulated in regions spanning from Northern Laos to Western China. Furthermore, our examination of the effects of hypothetical unsampled taxa on the phylogeographic reconstructions suggest that our inference of the location of the closest-inferred bat virus ancestors and their relative dispersal velocities are robust, rather than a byproduct of subsampling.
SARS-CoV-1 is thought to have been moved from bats in Yunnan province to Guangdong province via intermediate mammalian hosts15. The movement of palm civets (Paguma larvata) and raccoon dogs (Nyctereutes procyonoides) as part of the wild and farmed animal trade has been implicated in the emergence of SARS-CoV-141, suggesting that human-mediated transport of these infected animals can allow for viral lineages to move faster than their bat hosts. Given (i) the presumed emergence of SARS-CoV-1 through the animal trade, (ii) the recency and location of the closest-inferred ancestor relative to SARS-CoV-2, (iii) the consistently relatively high dispersal velocity associated with SARS-CoV-2, and (iv) the clear evidence that the epicenter of the SARS-CoV-2 pandemic at one of only four markets in Wuhan that sold live wildlife from plausible intermediate host mammal species42–44, either the closest-inferred ancestor or the direct ancestor of SARS-CoV-2 likely moved from an area in or around Yunnan province, to Hubei province, via the wild and farmed animal trade.
Our phylogeographic model assumes that bat and human sarbecoviruses evolve at similar rates. If bat viruses have slower rates of evolution than those inferred here, the closest-inferred bat virus ancestors likely circulated earlier. However, if this were true, there would be a concomitant decrease in their inferred rates of travel and, therefore, still not enough time for them to reach the provinces of emergence of their respective SARS-CoVs.
The roughly similar rates of diffusion among viral lineages and horseshoe bats suggest these viruses are able to rapidly and efficiently move among distinct horseshoe bat populations. If virus circulation was restricted to the location of each bat population, with minimal transmission of strains to other bats, we would expect a considerably lower diffusion coefficient for the viruses, but there is substantial overlap between the virus lineage and bat diffusion coefficients. Therefore, the sarbecoviruses are likely being transmitted among dense bat populations of various horseshoe bat species. This finding is further supported by the high levels of recombination in bat sarbecoviruses, requiring frequent co-infection and co-circulation of distinct virus strains in a given bat population. Viral movement is likely to take place through horseshoe bats regularly moving among roosts, especially as roost switching tends to happen during breeding season or if bats travel or migrate to hibernacula33,45. Furthermore, in hibernacula horseshoe bats will roost in very dense aggregations with different species46, providing the potential for spillover between bat species whilst bodily functions are downregulated and body temperatures can reach very low levels46,47. Heterothermy during hibernation and variable temperatures among species may also impact vulnerability to pathogens, as metabolic processes are downregulated during this time48,49.
Despite limited evidence for long-distance travel of horseshoe bats, there are several relatively long-distance viral lineage dispersals between nodes in the phylogenies. These geographically distant nodes therefore likely represent a lack of sarbecovirus sampling from areas between the distant locations rather than long-range bat movements, though it could also indicate dispersal from more migratory conspecific bat species, such as Chaerephon which can share roosts or Hipposideros which can disperse longer distances and are phylogenetically close to Rhinolophidae. For example, internal branches of the SARS-CoV-2-like tree represent movement from the east of China (Zhejiang) to the west (Yunnan), to Thailand and Cambodia (Fig. 4h). There are likely many more, still-undiscovered SARS-CoV-2-like viruses that would bridge the geographic gap between these regions (e.g., in northern Vietnam, Guizhou, and Guangxi) as the diversity of bats remains high across these regions, such that increased sampling in underexplored and undersampled regions would be helpful (Fig. 4b).
Whilst horseshoe bat ranges extend across most of China, species richness decreases at increasing latitudes27,50, and virus sampling in bats has largely been limited to higher biodiversity areas, including Yunnan province in China and Laos. Additionally, other large areas with high diversity, such as Northern Vietnam, remain unsampled and host a similar community of bat species. The SARS-CoV-1-like bat viruses are still better represented (n = 134) and therefore have more phylogeographic resolution than the SARS-CoV-2-like bat viruses (n = 31). Although there are viruses closely related to the SARS-CoVs across most of the NRRs, the viruses are relatively more divergent across the NRRs that overlap with the Spike gene (Fig. 1). This relatively greater divergence suggests that Spike genes that can bind non-bat entry receptors circulate in lower frequencies in wild horseshoe bat populations and, therefore, are less likely to be sampled.
Since the emergence of SARS-CoV-1 and SARS-CoV-2 in humans, the descendants of the bat viruses that gave rise to them have experienced years of evolution and almost certainly recombined with yet-unsampled lineages of the sarbecovirus tree. Hence, we should not expect to find the direct ancestor of either SARS-CoV-1 or SARS-CoV-2 circulating in wild bats in future sampling. Although the search for, and analyses of, sarbecoviruses related to SARS-CoV-1 and SARS-CoV-2 have often focused on partial genomes, typically the RdRp protein13,51–55, future sequencing efforts should routinely aim for whole genome sequences to detect and characterize all genomic fragments that descend from closely related ancestors. The non-recombinant segments of these mosaic genomes are pieces of a complex puzzle that are crucial for understanding the past and future emergence of SARS-CoVs into the human population.
Limitations of the study
Because there is insufficient temporal signal when calibrating a molecular clock using tip dating with sarbecoviruses sampled from bats and pangolins, we needed to calibrate the clock using SARS-CoV-1 and SARS-CoV-2 genomes from viruses sampled from humans and civets. The non-recombinant regions reported here are based on distinct inferred evolutionary histories across the sarbecovirus genomes, but due to the need for a detectable signal they likely represent a lower estimate of how much recombination has likely occurred among these viruses. Phylogeographic reconstructions across deep evolutionary histories can only capture general patterns and necessarily make abstraction of potentially complex population and habitat dynamics. Due the absence of bat viruses sampled closer to where SARS-CoV-1 and SARS-CoV-2 emerged, the inference of the location of the closest-inferred bat virus ancestors could be biased toward where published, closely related viruses were sampled. Similarly, a lack of sampling and publication of sarbecoviruses in nearby and potentially difficult to access regions limits our understanding of the broader ecology of these viruses. Even if sampling efforts would target more regions, host shedding of virus will vary over time because it is likely to increase during periods of increased metabolic stress, such as pregnancy or lactation. Therefore, detectability will vary across the year and there may be temporal sampling biases. However, the lack of sampled bat sarbecoviruses very closely related to SARS-CoV-1 and SARS-CoV-2 in and around their provinces of emergence despite increased sampling of bat viruses there13,25,26 suggests our findings are not solely the result of ascertainment bias. Further sampling of closely related sarbecoviruses, methodological improvements in disentangling complex recombination patterns, and improved molecular clock calibration can allow for the inference of more detailed spatiotemporal trends and recombination patterns, as well as more complete evolutionary histories of each NRR.
Resource availability
Lead contact
Further information and requests for resources and additional data should be directed to and will be fulfilled by the lead contact, Jonathan E. Pekar (jpekar@ed.ac.uk).
Materials availability
This study did not generate new unique reagents or specific biological material. The materials used and generated in this study are available from the lead contact. As a computational project, the input, results, output, and code of this study are publicly available and described below in the “data and code availability” section.
Data and Code Availability
All genome accessions are available in Data S1. The XML files, code, and newick trees are available at https://github.com/phylogeography/sarbecovirus_phylogeography.
STAR Methods
EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS
The sources of bioinformatic data analyzed in this study are available in the key resources table.
Key resources table
METHOD DETAILS
Dataset compilation
For the animal sarbecovirus dataset we compiled a set of 250 whole-genome sequences (>28,500 nt) of the subgenus Sarbecovirus sampled primarily in horseshoe bat hosts, including reference genomes of the human viruses SARS-CoV-1 and SARS-CoV-2 and six pangolin-infecting viruses (Data S1). This set contains all publicly available whole-genome sarbecoviruses available on 4th of October 2023 retrieved from the NCBI Genbank, GISAID and National Genomics Data Center of China databases. A single genome was used for SARS-CoV-1 (HSZ-Cc; accession AY394995.1), SARS-CoV-2 (Wuhan-Hu-1; accession MN908947.3) and the Guangdong pangolin sarbecovirus (MT121216.1) which has multiple published genome accessions associated with the same sample16,18,56,57. Collection date, sampling location and host metadata associated with the sequences were retrieved from the corresponding sequence databases and manually checked with the published papers describing the sequences. For sequences where exact city/cave sampling location information was provided a single longitude and latitude data point corresponding to that location was recorded, while provinces were recorded as sampling locations otherwise. The metadata for these 250 genomes are available in Data S1.
For the SARS-CoV-2 datasets, we used 863 genomes from the early pandemic44 and 1,000 genomes from late 202058. For the SARS-CoV-1 dataset, we downloaded the 98 publicly available whole genomes of SARS-CoV-1 from humans and civets from Genbank, with dates and locations matched by metadata or literature review: https://github.com/andersen-lab/SARS-CoV-1_Evolution/blob/main/metadata_extended.csv. A list of all accessions used in this study are available in Data S1. All GISAID accessions are also available through GISAID when using the identifier EPI_SET_241121yn.
Recombination analysis
The 250 whole-genome sequences were aligned using MAFFT (v7.453, localpair option)59 and manually curated using BioEdit60. To explore the recombination patterns of SARS-CoV-1-like and SARS-CoV-2-like sarbecoviruses separately we selected a number of representatives from each set of genomes on which to perform recombination analysis using the Genetic Algorithm for Recombination Detection (GARD)11. We analysed the two sarbecovirus clades separately to reduce overlapping recombination signal between the sets that may skew the position of identified breakpoints unique to each clade. Even though inter-clade recombinants exist in the dataset, the GARD algorithm does not require potential recombination parents to be in the dataset, but assesses recombination simply by comparing tree topology incongruencies. Thus, signals of SARS-CoV-1-like/SARS-CoV-2-like inter-clade recombinants can be detected in genomes of each set without the parental lineages of recombinant genomes having to be present in each dataset. The genome selection was made by first performing a preliminary phylogenetic analysis of the relatedness between these genomes using the recombination breakpoints presented in Lytras et al. (2022)4. Based on the preliminary analysis, clusters of genomes grouping with SARS-CoV-1 and SARS-CoV-2 for at least one of their genomic segments and showing the same overall recombination pattern were manually identified and one representative genome was selected from each cluster for the corresponding downstream GARD analyses. For the SARS-CoV-2-like set, 24 genomes including the SARS-CoV-2 reference Wuhan-Hu-1 (GenBank accession MN908947) were used in the GARD analysis. For the SARS-CoV-1-like set, 38 genomes including the SARS-CoV-1 HSZ-Cc genome (accession AY394995) were used in the analysis (Data S1). The representative sequences for each set were retrieved from the full alignment. The ends of genomes are usually more difficult to sequence, leading to shorter or extended 3’ and 5’ ends in the sequences. Due to the incomplete coverage of these short genomic regions they are unlikely to contain any informative sites for the recombination analysis and for this reason were trimmed from each resulting subset alignment. GARD was performed for each subset alignment using the HyPhy software suite v.2.5.2961 under a general time reversible (GTR) substitution model and a 3-class general discrete distribution (GDD) to account for site-to-site rate variation. The SARS-CoV-2-like GARD analysis finished with a total of 92 consecutive breakpoint models, while the SARS-CoV-1-like analysis finished with 56. Based on these models, we accepted positions as confident recombination breakpoints if they were present in at least 1/2 of the consecutive models for the SARS-CoV-2-like set and in at least 1/3 of the models for the SARS-CoV-1-like set (proportional to the total number of models). Breakpoints closer than 100 alignment positions to one another were merged. The resulting breakpoint positions were mapped back to the full whole-genome alignment which was subsequently split into the corresponding non-recombinant regions (NRR) for each set of breakpoints. Since some NRRs have relatively short lengths, incorrect alignment of sites may have a larger effect on the phylogenetic inference. To avoid alignment errors from biasing the phylogenetic inferences, all NRR alignments were processed using HMMCleaner62, which removes sequence segments that fit poorly in the rest of the alignment, before performing any downstream analyses.
To separate all genomes into SARS-CoV-1-like and SARS-CoV-2-like viruses we first inferred maximum likelihood phylogenies for every NRR alignment using IQ-TREE (v.1.4.4)63 under a GTR+I+F+G4 substitution model, assessing node support by performing 10,000 replicates of the ultrafast bootstrap approximation for each inference64. All trees were rooted by the outgroup clade of sarbecovirus genomes sampled in Europe and Africa (accessions: NC_014470, MZ190137, KY352407, MT726043, MT726044, MT726045, MZ190138, MW719567) and nodes with ultrafast bootstrap support below 70 were collapsed using the ‘ete3’ python package65.
Based on the resulting topologies, the viruses clustering monophyletically with SARS-CoV-2 but not with SARS-CoV-1 for at least one of the SARS-CoV-2-like GARD NRR trees were placed in the SARS-CoV-2-like viruses and vice versa for the SARS-CoV-1-like virus set. This led to two non-independent sets of 38 SARS-CoV-2-like and 140 SARS-CoV-1-like genomes (each including SARS-CoV-2 and SARS-CoV-1 respectively) used in the downstream dating and phylogeographic analyses; there were 19 genomes shared between these two datasets. The set to which each genome is assigned is listed in Data S1.
Comparison of breakpoint positions
Based on the GARD recombination analysis, 30 and 43 breakpoint positions were detected for the two virus subgroups. Only one SARS-CoV-1 genome was used in both GARD sets (Data S1), so overall breakpoint inferences are expected to represent independent recombination patterns of known SARS-CoV-1-like and SARS-CoV-2-like viruses respectively. The lengths of the resulting NRRs on the alignment with all 250 sarbecovirus genomes had a median of 870 nt and ranged from 309 to 3,414 nt for SARS-CoV-1-like viruses, and a median length of 595 nt and ranged from 168 to 2,706 nt for SARS-CoV-2-like viruses. It should be noted that the last NRR is the longest due to extended ends at the 3’ end of the genomes.
We calculated the distances between the closest position pairs between the two breakpoint sets, resulting in 30 pair distances (each SARS-CoV-1-like position matched with its closest SARS-CoV-2-like position). To compare whether the pairs were more or less closely located than one would expect by chance, we performed 100 simulations of position sets randomly sampled across the same alignment length without replacement (one set with 30 positions and one with 43 for each simulation to match the empirical data on the number of breakpoints). We then quantified the distances between closest pairs in the simulated sets as done with the observed breakpoints. Both empirical and simulated distances formed distributions with log-normal shapes, hence 90% confidence intervals and means were estimated using Cox’s method for inferring log-normal confidence intervals66. Similarly, the significantly closest observed breakpoint pairs were defined as those with distances below the one-tailed lower 90% confidence interval of the observed log-normal distance distribution.
Molecular clock calibration
There is insufficient data to inform a molecular clock when using tip dating with sarbecoviruses sampled from bats and pangolins6. To identify a suitable rate prior for each NRR identified from the recombination analysis, we conducted a phylogenetic analysis of SARS-CoVs sampled from humans.
Molecular clock calibration across the 44 inferred NRRs of the SARS-CoV-2-like viruses was conducted using a Bayesian approach with BEAST v1.10.567 and using the BEAGLE library v3 to increase computational performance68. We inferred the substitution rate of SARS-CoV-2 across each of these NRRs by using a subset of SARS-CoV-2 genomes (n=1,000) from our previously published Bayesian phylogenetic analysis of 3,959 genomes sampled up to late 2020 in Europe58 (Data S1), a non-parametric skygrid prior with 19 grid points and a cutoff of 1 year, and a Hasegawa-Kishino-Yano (HKY85) substitution model69 with four gamma categories. The topologies and substitution models were shared across all the NRRs, but assuming an independent and estimable strict clock rate under an uninformative continuous-time Markov chain (CTMC) rate reference prior70 for each NRR. We ran one Markov chain Monte Carlo (MCMC) chain of 500 million generations, subsampling every 50,000 iterations to a continuous parameter log file. The first 10% of the chains were discarded as burn-in. Convergence and mixing was assessed in Tracer v1.7.171 and all relevant effective sample size (ESS) values for clock rates were >100. We also perform an additional clock calibration without any NRR partitioning to infer a single clock rate prior for the phylogeographic analysis (see below). The rest of the parameterizations are identical to the partitioned clock calibration inference.
To calibrate the molecular clock for the SARS-CoV-1-like viruses analysis, we first performed Bayesian phylogenetic inference on 98 human and civet SARS-CoV-1 genomes (Data S1), using an HKY85 substitution model with a gamma site heterogeneity model and 4 gamma categories (G4), a strict molecular clock, and an exponential growth coalescent prior. We ran a single chain of 100 million generations, subsampling every 10 thousand iterations to a continuous parameter log and tree files. We next inferred the substitution rate of SARS-CoV-1 across the 31 inferred NRRs of the SARS-CoV-1-like viruses by using an empirical tree distribution (n=1,000) from these results, sharing the topologies and substitution models across all the NRRs, with independent and estimable clock rates under CTMC rate priors. We ran one chain of one million generations, subsampling every one hundred iterations to a continuous parameter log file. For each step of the inference, the first 15% of the chains were discarded as burn-in. Convergence and mixing was assessed in Tracer and all relevant effective sample size (ESS) values were >200.
Bayesian divergence time estimation
We used BEAST to perform phylogenetic inference and ancestral state reconstruction across the 44 NRRs of the SARS-CoV-2-like viruses. For the primary analysis (Fig. S2), we used all of the 44 SARS-CoV-2-like genomes (including SARS-CoV-2), an uncorrelated relaxed clock using calibrated region-specific normally-distributed substitution rate priors for each NRR based on SARS-CoV-2 sampled from humans, a GTR+G4 substitution model, a non-parametric skygrid prior72 with 49 grid points and a cutoff of 1,000 (which translates to 1022 CE), and Hamiltonian Monte Carlo gradient-based sampling for node ages and branch rates73,74. Specifically, we set the mean and standard deviation of the rate prior for each NNR to equal the posterior mean and standard deviation estimated during molecular clock calibration with SARS-CoV-2 genomes. We tracked the stem age of SARS-CoV-2 (i.e., where SARS-CoV-2 connects to the rest of the phylogeny, also referred to as the time of the closest-inferred ancestor, with the parent node of SARS-CoV-2 referred to as the closest-inferred ancestor). For each NRR, we first ran one chain of 20 million generations, subsampling every two thousand iterations to continuous parameter log and tree files. The first 10% of the chain was discarded as burn-in for each analysis. Convergence and mixing was assessed in Tracer. For analyses of NRRs where any relevant ESS values were <200, we ran a second chain with the same specifications as the first, discarded the first 5–25% of the chain as burn-in, and combined the two chains. All relevant ESS values for the one- and two-chain analyses were >200.
We performed the same analysis for the SARS-CoV-1-like viruses for its 31 NRRs, using all of the 140 SARS-CoV-1-like genomes (including SARS-CoV-1). The BEAST parameterization was identical to the analysis for the SARS-CoV-2-like viruses, except we (1) tracked the stem age of SARS-CoV-1 (i.e., the time of the closest-inferred ancestor of SARS-CoV-1), (2) used the relevant calibrated molecular clock priors for the 31 NRRs of the SARS-CoV-1-like viruses, and (3) ran one chain of 100 million generations, sub-sampling every 10 thousand iterations to continuous parameter log and tree files.
Bayesian divergence time estimation sensitivity analyses
We performed a series of sensitivity analyses for the Bayesian divergence time estimation: (i) using SARS-CoV-2 genomes from early 2020, rather than late 2020, when calibrating the molecular clock prior for the subsequent analyses of SARS-CoV-2-like viruses, and (ii) including two additional human SARS-CoV-1 genomes and two civet SARS-CoV-1 genomes in the SARS-CoV-1-like viruses analyses, representing both the 2002–2003 and 2003–2004 SARS-CoV-1 outbreaks. For (i), we performed molecular clock calibration using the same approach as with the SARS-CoV-2 genomes sampled from late 2020, except here using 863 genomes from the early pandemic (Data S1), sharing the topologies and substitution models across all the NRRs, but assuming an independent and estimable strict clock rate under an uninformative CTMC rate prior70 for each NRR. We ran one MCMC chain of one million generations, subsampling every one thousand iterations to a continuous parameter log file. The first 10–15% of the chains were discarded as burn-in. Convergence and mixing was assessed in Tracer and all relevant effective sample size (ESS) values were >200. There was no temporal overlap between the datasets of SARS-CoV-2 genomes from early 2020 and late 2020; the latest genome in early 2020 dataset was sampled on 14 February 2020 and the earliest genome in the late 2020 dataset was sampled on 22 February 2020.
For (ii), SZ3 (accession AY304486) was selected as representative of the earliest sampled detection in the intermediate animal reservoir (market-based civets). ZS-B (AY394996) was selected to represent the early stages of the human epidemic from late January 2003 onwards. As both of these are from the 2002–2003 SARS-CoV-1 outbreak, we also included the civet virus HC/SZ/266/03 (accession AY545916) and the human virus GZ0401 (AY568539) from the 2003–2004 SARS-CoV-1 outbreak.
Convergence and mixing were assessed in Tracer. For analyses of NRRs where any relevant ESS values were <200, we ran up to two more chains with the same specifications as the initial chain, discarded the first 10% of the chain as burn-in, and combined all chains. All relevant ESS values for the combined chains were >200.
Constructing the recombinant common ancestor
For each NRR of the SARS-CoV-2-like viruses, we reconstructed the most likely nucleotide sequence belonging to the parent node of SARS-CoV-2 based on the tree distribution from the primary Bayesian divergence time estimation (Fig. S2). The most likely nucleotide sequences for the 44 NRRs are then concatenated to form the recombinant common ancestor (recCA)14 of SARS-CoV-2.
To analyze how the recCA changed as more genomes were released, we constructed the recCA of SARS-CoV-2 as the genomes of closely related sarbecoviruses were chronologically published. We performed ancestral state reconstruction using the posterior trees with all 38 taxa. We next identified the inferred sequence of the most recent common ancestor of SARS-CoV-2 and its most closely related taxa based on the published genomes as new genomes were sequenced. We then examined the genetic similarity of the recCA, as well as the most closely related published genome, to SARS-CoV-2 as a function of the publication date. The same analysis was performed for SARS-CoV-1 and the SARS-CoV-1-like viruses and each phylogenetic sensitivity analysis (see above).
Bayesian phylogeographic inference
We inferred the evolutionary and geographic history of the 44 NRRs of the SARS-CoV-2-like viruses with BEAST. In order to optimally inform a flexible phylogeographic reconstruction approach, we jointly fit a bivariate continuous diffusion model using a shared estimable precision matrix across all independent random NRRs phylogenies, rescaling this matrix across branches in each tree to model a relaxed random walk (RRW)75. All branch-specific RRW rate scalars are drawn from the same log-normal distribution with a mean of 1 and an estimable standard deviation. In the continuous phylogeographic diffusion model, we exclude the sampling locations for the human and pangolin viruses and integrate over the uncertainty for locations that lack precision. As a tree prior, we specify a joint non-parametric skygrid prior with 49 grid points and a cutoff of 1,000 (which translates to 1022 CE) or a cutoff of 0.5 substitutions per site for trees scaled in units of genetic distance (see below). We further specify a joint HKY85 nucleotide substitution model across NRRs. These relatively simple sequence evolutionary parameterizations were chosen to meet the requirements for the Prisoner of War (PoW) model21 to rescale the evolutionary histories (see below). We perform two versions of this analysis (Fig. S2): (i) using a hierarchical strict molecular clock model that allows for the clock rate of each NRR to vary (but must average to the overall rate), with the single inferred clock rate prior from above for the overall clock rate, in order to estimate the trees in units of time and (ii) using independent strict molecular clock models with clock rates fixed to 1 and contemporaneous tips in order to estimate trees in substitution units for PoW transformation. For (i), we also included gamma-distributed among-site rate heterogeneity with a hierarchical prior over the NRR-specific shape parameters. We ran two chains of 500 million generations, subsampling every 10 thousand iterations to continuous parameter log and tree files; the first 10% of the chain was discarded as burn-in. Convergence and mixing was assessed in Tracer, and all relevant ESS values were >200.
We performed the same analysis for the SARS-CoV-1-like virus set for its 31 NRRs. The BEAST parameterization was identical to the analysis for the SARS-CoV-2-like viruses, but here we use the rate priors with standard deviations divided by 5 to assist in convergence. We ran multiple chains for >200 million generations, subsampling every 500 thousand iterations to continuous parameter log and tree files, and discarded the necessary burn-in prior to combining the independent samples. Convergence and mixing was assessed in Tracer, and all relevant ESS values were >100.
To account for the time-dependent decline in the evolutionary rate of sarbecoviruses, we applied the PoW model, a mechanistic evolutionary method for estimating virus substitution rates over long evolutionary timescales. For this, we used the post-burn-in substitution trees inferred from the previous step along with the posterior rate distributions of the SARS-CoV-1-like viruses and the SARS-CoV-2-like viruses to rescale substitution trees into time trees using the PoW transformation (see ref Ghafari et al.21 for more details on the method).
We used TreeAnnotator v1.1067 to obtain and annotate a maximum clade credibility (MCC) tree for each continuous phylogeographic reconstruction. Finally, we used the R package ‘seraphim’76 to extract the spatiotemporal information embedded within 100 trees sampled from each post-burn-in posterior distribution and to visualize the different phylogeographic reconstructions. For each NRR, we estimated the weighted diffusion coefficients29 with the ‘spreadStatistics’ function of the R package ‘seraphim’76.
Bayesian phylogeographic inference sensitivity analysis
To examine the consequences of potentially biased sampling on the phylogeographic inference, we included in the Bayesian inference two hypothetical unsampled taxa that were associated with a specific location and a sampling date of 2020 but no observed sequence data. For the analysis of the SARS-CoV-2-like viruses, one taxon was given a sampling location associated with a limestone karst in Guizhou expected to have high bat diversity, and the other taxon was given a sampling location associated with caves in Hunan also expected to have high bat diversity. Both of these locations are between Hubei, where SARS-CoV-2 emerged, and Yunnan, China and Northern Laos, where the bat viruses most closely related to SARS-CoV-2 were sampled. For the analysis of the SARS-CoV-1-like viruses, both taxa were assigned sampling locations in Guangxi, between Guangdong, where SARS-CoV-1 emerged, and Yunnan, where the bat viruses most closely related to SARS-CoV-1 were sampled. The unsampled taxa were constrained to branch off from the branch connecting the human SARS-CoV and closest-inferred ancestor.
Ranking the dispersal velocity of the lineage descending from the closest-inferred ancestor
Using 100 trees sampled from each post-burn-in posterior tree distribution from each NRR from the tip-dated phylogeography analyses, we calculated the dispersal velocity for each tip branch, excluding the branches leading to the human or pangolin viruses. We next computed the minimum distance between the closest-inferred ancestor and Guangdong province or Hubei province for SARS-CoV-1 or SARS-CoV-2, respectively. The necessary dispersal velocity for the closest-inferred ancestor to reach either Guangdong or Hubei province was calculated using this minimum distance and the branch leading to, respectively, SARS-CoV-1 or SARS-CoV-2 (i.e., the time between closest-inferred ancestor and the SARS-CoV). Tip branches associated with taxa that have parent nodes dated to more than 50 years before the most recently sampled virus were excluded from downstream analyses. If the tip branch was excluded from at least 50% of the samples for an NRR, the tip branch was removed from the NRR entirely. If the tip branch associated with a human SARS-CoV was excluded from at least 50% of the samples for an NRR, the NRR was removed entirely from downstream analyses. The 50-year threshold was chosen based on the observed substitution saturation inflection point (Fig. 3a), thereby also focusing primarily on the recent movement of sarbecoviruses.
We then computed the rank of the dispersal velocity for each remaining tip branch in each tree sample. These posterior dispersal velocity ranks were pooled together for each tip across all samples for an overall posterior dispersal velocity rank distribution, or across all samples for each NRR for posterior dispersal velocity rank distributions per NRR.
To calculate how much the posterior dispersal velocity rank distribution deviates from an expected average rank, we use an adjusted formula for skewness with a given “mean” of 0.5:
Modeling of bat distributions
Our bat distribution database comprised of the data used in recently published studies27, recent Global Diversity Information Facility (GBIF) data77,78, and the Darkcide database79. The Darkcide database included recent data collected by bioacoustic surveys in Chiang Mai and Kanchanaburi in Thailand by A.C. Hughes using an echometer Pro, and multiple published datasets27,77,78,80–88. The data was integrated as we attempted to keep the database complete, even if not all data is used within our particular study. All Rhinolophid bat species, as well as Chaerephon plicatus and Aselliscus stoliczkanus were extracted from the dataset, and then the dataset was clipped to Southeast Asia (India to the West, Mongolia to the North, Japan to the East and Malaysia to the South) to buffer the region of interest with distribution data. Once all duplicate records had been removed, the clipped data for Rhinolophids included 10,254 distribution points for all species across the region, in addition to 235 records for A. stoliczkanus and 108 for C. plicatus.
For environmental data, we based variable selection on those used in Zhou et al. (2021)27 and known from our previous work to be important to bats in the region89,90. All analysis was conducted in ArcMap 10.8 (ESRI) at a resolution of 0.008°, which is approximately 1km2, and all layers were clipped and subsampled to the same resolution using a mask. Climate variables were downloaded from Chelsa91 (available at https://chelsa-climate.org/downloads/) based on a selection from our previous analyses (the bioclimatic variables: Annual Mean Temperature, Mean Diurnal temperature Range, Isothermality, Temperature Seasonality, Maximum and Minimum Temperature, Annual Precipitation, Precipitation seasonality, Maximum and Minimum precipitation)89,90. From Envirem92 (https://envirem.github.io/) we used the Thornthwaite Aridity index, climate moisture index, continentality, embergers pluviothermic quotient, growing degree days above 0, and both annual mean potential evapotranspiration and seasonality of potential evapotranspiration. In addition we used actual evapotranspiration93, a vegetation canopy height layer94, forest canopy height95, tree density96, distance to bedrock97 (the bdticm and bdricm aggregated grids available at https://github.com/ISRICWorldSoil/SoilGrids250m/), and bulk carbon stock 97. These were all trimmed to the same region that bat data was clipped to and then species distributions were modelled using Maxent98 using the GUI following the same approach as Hughes et al. (2012)89. For each species, three replicates were run and the average was used. These probabilistic distribution maps were then reclassified using the 10th percentile cumulative logistic threshold to produce a binary presence-absence map with ArcMap 10.8, and all models had an AUC exceeding 0.85. The models were then clipped to the study region, and any Rhinolophids with no modelled suitable habitat in the clipped study area were removed from analysis. As the projections show suitable habitat but do not account for species biogeography we then cross-checked the 40 remaining Rhinolophid maps against the IUCN database (https://www.iucnredlist.org/), and species not listed were checked against relevant publications on the species. Based on the IUCN mapped ranges species endemic to Japan, or species not found in Japan or at similar latitudes (areas of the South Korean Peninsula are climatically suitable for a number of species, but lack continuous habitat to reach these regions if largely distributed in Southeast Asia) were then clipped to the biogeographically suitable range. The binary maps were then stacked using the Mosaic to New Raster function in ArcMap 10.8 to find the maximum number of species that may be co-distributed.
QUANTIFICATION AND DATA ANALYSIS
All computational and statistical analyses were performed with the open-source software tools referenced in the STAR Methods along with the described procedures. Posterior estimates of dispersal velocities from the phylogenetic inference were compared using their rank distributions and skewness as post-hoc measures.
Supplementary Material
Data S1: Sarbecovirus genomes metadata and genome accessions, related to Figures 1–6. Data S1.1: Temporal, geographic and acknowledgement information for the 250 sarbecovirus genomes used in the paper. Sequences used for each SARS-CoV-1-like and SARS-CoV-2-like GARD recombination analysis and the virus group of each genome (SARS-CoV-1-like, SARS-CoV-2-like or outgroup) are also indicated in the table. Data S1.2: Genome accessions for all sarbecovirus genomes used in the paper.
Figure S1. Schematics describing the closest-inferred ancestor, direct ancestor, and consequences of recombination, related to Figures 1–6. a, Phylogenies illustrating the relationship between a SARS-CoV, the direct ancestor, and the closest-inferred ancestor. The closest-inferred ancestor could be the tMRCA of a given SARS-CoV and one non-human sarbecovirus (top) or the tMRCA of the SARS-CoV and multiple non-human sarbecoviruses (bottom); the closest-inferred ancestor is the parent node of the SARS-CoV. The direct ancestor is the bat sarbecovirus that was involved in each cross-species transmission event leading to the SARS-CoV and could have existed at any point on the branch leading from the closest-inferred ancestor to the SARS-CoV. b, Phylogenies and example genomes illustrating the effects of recombination in the phylogenetic inference when a whole-genome alignment is used for the inference (top) and how accounting for recombination can reveal multiple evolutionary histories (bottom) for recombinant genomes (genome C) when the alignment is split by recombination breakpoints into non-recombinant regions (NRRs).
Figure S2. Workflow to perform phylogenetic and phylogeographic analyses, related to Figures 1–6. Starting with 250 sarbecovirus genomes, we align and perform recombination analysis to determine the number of SARS-CoV-1-like and SARS-CoV-2-like viruses and their corresponding NRRs for subsequent analyses. We then calibrate the molecular clock for SARS-CoV-1-like viruses and SARS-CoV-2-like viruses using SARS-CoV-1 genomes and SARS-CoV-2 genomes (these datasets are separate from the SARS-CoV-1-like viruses and SARS-CoV-2-like viruses), respectively. Finally, we perform Bayesian divergence time estimation, recCA inference, and phylogeography analyses, each of which comprised both primary and sensitivity analyses. Ovals indicate the number of sarbecovirus genomes included in subsequent steps. Rectangles without inner vertical lines indicate analyses performed; rectangles with inner vertical lines indicate primary versus sensitivity analyses, with non-bolded text describing the details of the analysis performed.
Figure S3. Breakpoint pair distance distributions of the SARS-CoV-1-like and SARS-CoV-2-like viruses, related to Figure 1. Mean and 90% confidence intervals of the log-normal breakpoint pair distance distributions for the observed SARS-CoV-1-like and SARS-CoV-2-like estimated breakpoints and the simulated randomly sampled breakpoint sets.
Figure S4. Substitution rates for the SARS-CoV-1-like and SARS-CoV-2-like viruses and time of the closest-inferred ancestor using additional genomes or different molecular clock priors, related to Figure 1. Substitution rates (substitutions/site/year) across the (a) 31 NRRs for the SARS-CoV-1-like viruses and (b) 44 NRRs for the SARS-CoV-2-like viruses. Prior distributions are shown in pink and yellow and the dashed lines indicate two standard deviations from the mean. Posterior distributions are shown in blue and green and solid lines indicate the 95% highest posterior density (HPD) surrounding the median. Refer to Table S1 for numerical values. c, Time of the closest-inferred ancestor of SARS-CoV-1 across the 31 NRRs of the SARS-CoV-1-like viruses with the calibrated molecular clock from the SARS-CoV-1 analysis and when adding in two SARS-CoV-1 genomes from civets and and two SARS-CoV-1 genomes from humans. d, Time of the closest-inferred ancestor of SARS-CoV-2 across the 44 NRRs of the SARS-CoV-2-like viruses with the calibrated molecular clock from the SARS-CoV-2 analysis with early 2020 sampling, and with the calibrated molecular clock from the SARS-CoV-2 analysis with late 2020 sampling. The dots indicate the median and the lines are the 95% HPD. Refer to Table S1 for numerical values.
Figure S5. The recCA (recombinant common ancestor) of SARS-CoV-1 over time when including additional SARS-CoV-1 genomes and SARS-CoV-2 over time when using rate priors based on SARS-CoV-2 genomes sampled from early 2020, related to Figure 2. a, The similarity of the recCA of SARS-CoV-1 and the most closely related sarbecovirus among the SARS-CoV-1-like viruses to SARS-CoV-1 as a function of time. b, The similarity of each NRR of the recCA of SARS-CoV-1 to SARS-CoV-1 as a function of time. c, The similarity of the recCA of SARS-CoV-2 and the most closely related sarbecovirus in the SARS-CoV-2 viruses to SARS-CoV-2 as a function of time. The dashed vertical line is the sampling date of the earliest sampled genome of SARS-CoV-2. d, The similarity of each NRR of the recCA of SARS-CoV-2 to SARS-CoV-2 as a function of time. The right panels in (a, c) are zoomed-in panels of the dashed box in the left panels. Each square in (b, d) represents the similarity of the NRR of the recCA to its relative SARS-CoV at one point in time and is colored by date of the most recent non-human sarbecovirus in the reconstruction.
Figure S6. Posterior rank distributions of the tip branches from the SARS-CoV-1-like and SARS-CoV-2-like virus tip-dated phylogeographies, related to Figure 5. Violins sorted by mean, with violin corresponding to (a) SARS-CoV-1 colored blue and (b) SARS-CoV-2 colored pink.
Figure S7. Phylogenies including hypothetical unsampled taxa, related to Figure 6. Time-scaled maximum clade credibility (MCC) phylogenies of a single NRR of (a) the SARS-CoV-1-like viruses (NRR 19) and (b) and the SARS-CoV-2-like viruses (NRR 21) when including unsampled taxa in the phylogeographic reconstruction. The phylogenies are shown on a log-transformed time-scale, and horizontal segments displayed at each internal node indicate the 95% HPD interval of the estimated node age. The shaded box highlights the respective SARS-CoV and the two unsampled taxa, which are indicated by the dashed lines. Note that the respective SARS-CoV and two unsampled taxa are constrained to be monophyletic. Refer to Fig. S8 for sampling locations of the unsampled taxa.
Figure S8. Inferred location of the closest-inferred ancestor comparison when including hypothetical unsampled taxa, related to Figures 4–6. The estimated position of the closest-inferred ancestors of SARS-CoV-1 (a, b) and SARS-CoV-2 (c, d), according to continuous phylogeographic reconstructions that did not include (a, c) or did include (b, d) unsampled taxa in the analyses. Crosses indicate the locations of the unsampled taxa. As in Figure 4g-h, the estimated position of the closest-inferred ancestor is based on the continuous phylogeographic reconstruction of all NRRs. We here display the 95%, 75% and 50% HPD regions with an increasing red color darkness. Refer to Fig. S7 for examples of the phylogenetic relationship of the unsampled taxa to SARS-CoV-1 and SARS-CoV-2.
Figure S9. Comparison of dispersal velocity when including hypothetical unsampled taxa, related to Figures 4–6. Dispersal velocities for tip branches associated with (a) SARS-CoV-1 or (b) SARS-CoV-2 when excluding and including unsampled taxa in the phylogeographic analyses.
Figure S10. Additional dispersal velocity rank statistics of SARS-CoV-1-like and SARS-CoV-2-like viruses, related to Figures 5–6. Probability of posterior rank of 1 (a, c) and the Kullback–Leibler (KL) divergence from a uniform distribution (b, d) for the dispersal velocity rank distribution of each tip branch when excluding (a, b) or including (c, d) unsampled taxa in the phylogeographic analyses. The y-axes for the KL divergence are inverted because of a skewed distribution, as observed with the rank distributions for the dispersal velocities of the tip branches leading to SARS-CoV-1 and SARS-CoV-2 (Fig. S6), will be close to 0. Note the different y-axis scales for the SARS-CoV-1-like and SARS-CoV-2-like viruses in (a, c). Additional SARS-CoV-1-like viruses with fast dispersal velocities are colored and identified by their strain name.
Table S1. Posterior clock rates and times of the closest-inferred ancestor for SARS-CoV-1 and SARS-CoV-2 across each NRR, related to Figure 1. Table S1.1: Posterior clock rates across each NRR from the clock calibration using SARS-CoV-1 genomes and Bayesian divergence time analyses of SARS-CoV-1-like viruses. Table S1.2: Posterior clock rates across each NRR from the clock calibration using SARS-CoV-2 genomes sampled in late 2020 and Bayesian divergence time analyses of SARS-CoV-2-like viruses. Table S1.3: Posterior clock rates across each NRR from the clock calibration using SARS-CoV-2 genomes sampled in early 2020 and Bayesian divergence time analyses of SARS-CoV-2-like viruses. For Tables S1.1–S1.3, the clock calibration posterior results are used as a prior for the Bayesian divergence time analyses. Table S1.4: Posterior times (median and 95% HPD) of the closest-inferred ancestor of SARS-CoV-1 across each NRR. The SARS-CoV-1 clock calibration posterior serves as a prior for the analyses, with the mean and standard deviation of the clock calibration results varying across the analyses and one analysis having two additional civet SARS-CoV-1 genomes and two additional human SARS-CoV-1 genomes. Table S1.5: Posterior times (median and 95% HPD) of the closest-inferred ancestor of SARS-CoV-2 across each NRR. The SARS-CoV-2 clock calibration posterior serves as a prior for the analyses, with the standard deviation of the clock calibration results varying across the analyses.
Highlights.
Fragments of the human SARS-CoVs share very recent common ancestors with bat viruses
SARS-CoV-1-like and SARS-CoV-2-like viruses have circulated in Asia for millennia
Recent ancestors of human SARS-CoVs likely circulated in Western China and Northern Laos
These ancestors traveled unexpectedly fast to reach sites of human emergence
Acknowledgements
S.D. acknowledges support from the Fonds National de la Recherche Scientifique (F.R.S.-FNRS, Belgium; grant n°F.4515.22) and from the Research Foundation - Flanders (Fonds voor Wetenschappelijk Onderzoek - Vlaanderen, FWO, Belgium; grant n°G098321N). S.D. and P.L. acknowledge support from the European Union Horizon 2020 project MOOD (grant agreement n°874850). P.L., A.M., A.R. and M.A.S. acknowledge support from the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 725422-ReservoirDOCS), from the Wellcome Trust through project 206298/Z/17/Z and from the National Institutes of Health grants R01 AI153044, R01 AI162611 and U19 AI135995. P.L. also acknowledges support from the Research Foundation - Flanders (Fonds voor Wetenschappelijk Onderzoek - Vlaanderen G0D5117N and G051322N). T.I.V. acknowledges support from the Branco Weiss Fellowship. D.L.R. and J.H. acknowledge funding from the UK Medical Research Council (MRC, MC_UU_12014/12, MC_UU_00034/5 and MR/V01157X/1). S.L. was funded by an MRC studentship. M.G. is supported by the Biotechnology and Biological Science Research Council (BBSRC) grant number BB/M011224/1. A.K. acknowledges funding from European Research Council grant 101001623-PALVIREVOL. J.E.P. acknowledges support from NIH (T15LM011271) and the UC San Diego Merkin Fellowship. This work was funded in part with federal funds from the National Institute of Allergy and Infectious Diseases, National Institutes of Health (NIH), Department of Health and Human Services (contract no. 75N93021C00015 to M.W.). J.O.W. acknowledges support from NIH (R01AI135992). We gratefully acknowledge the authors from the originating laboratories and the submitting laboratories, who generated and shared through GISAID the viral genomic sequences and metadata on which this research is based (Data S1). We thank Kristian Andersen, Alexander Crits-Christoph, Fabian Leendertz, Niema Moshiri, and Richard Neher for fruitful conversations and insight. We gratefully acknowledge support from Advanced Micro Devices, Inc. with the donation of parallel computing resources used for this research.
Footnotes
Declaration of interests S.L. has received consulting fees from EcoHealth Alliance. M.A.S. receives contracts and grants from the US Food and Drug Administration, the US Department of Veterans Affairs, and Janssen Research and Development unrelated to this research. J.O.W. has received funding from the CDC (ongoing) through contracts or agreements to his institution unrelated to this research. J.E.P., M.A.S., A.R., M.W., and J.O.W. have received consulting fees and/or provided compensated expert testimony on SARS-CoV-2 and the COVID-19 pandemic.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
References
- 1.Wang L-F, Shi Z, Zhang S, Field H, Daszak P, and Eaton BT (2006). Review of bats and SARS. Emerg. Infect. Dis. 12, 1834–1840. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 2.Ren W, Li W, Yu M, Hao P, Zhang Y, Zhou P, Zhang S, Zhao G, Zhong Y, Wang S, et al. (2006). Full-length genome sequences of two SARS-like coronaviruses in horseshoe bats and genetic variation analysis. J. Gen. Virol. 87, 3355–3359. [DOI] [PubMed] [Google Scholar]
- 3.Lau SKP, Woo PCY, Li KSM, Huang Y, Tsoi H-W, Wong BHL, Wong SSY, Leung S-Y, Chan K-H, and Yuen K-Y (2005). Severe acute respiratory syndrome coronavirus-like virus in Chinese horseshoe bats. Proc. Natl. Acad. Sci. U. S. A. 102, 14040–14045. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 4.Lytras S, Hughes J, Martin D, Swanepoel P, de Klerk A, Lourens R, Kosakovsky Pond SL, Xia W, Jiang X, and Robertson DL (2022). Exploring the Natural Origins of SARS-CoV-2 in the Light of Recombination. Genome Biol. Evol. 14. 10.1093/gbe/evac018. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 5.Lin Xian-Dan, Wang Wen, Hao Zong-Yu, Wang Zhao-Xiao, Guo Wen-Ping, Guan Xiao-Qing, Wang Miao-Ruo, Wang Hong-Wei, Zhou Run-Hong, Li Ming-Hui, Tang Guang-Peng, Wu Jun, Holmes Edward C., Zhang Yong-Zhen (2017). Extensive diversity of coronaviruses in bats from China. Virology 507, 1–10. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 6.Boni MF, Lemey P, Jiang X, Lam TT-Y, Perry BW, Castoe TA, Rambaut A, and Robertson DL (2020). Evolutionary origins of the SARS-CoV-2 sarbecovirus lineage responsible for the COVID-19 pandemic. Nat Microbiol 5, 1408–1417. [DOI] [PubMed] [Google Scholar]
- 7.Wertheim JO, Chu DKW, Peiris JSM, Kosakovsky Pond SL, and Poon LLM (2013). A case for the ancient origin of coronaviruses. J. Virol. 87, 7039–7045. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 8.Müller NF, Kistler KE, and Bedford T (2022). A Bayesian approach to infer recombination patterns in coronaviruses. Nat. Commun. 13, 4186. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 9.Worobey M, Telfer P, Souquière S, Hunter M, Coleman CA, Metzger MJ, Reed P, Makuwa M, Hearn G, Honarvar S, et al. (2010). Island biogeography reveals the deep history of SIV. Science 329, 1487. [DOI] [PubMed] [Google Scholar]
- 10.Wertheim JO, and Kosakovsky Pond SL (2011). Purifying selection can obscure the ancient age of viral lineages. Mol. Biol. Evol. 28, 3355–3365. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 11.Kosakovsky Pond SL, Posada D, Gravenor MB, Woelk CH, and Frost SDW (2006). GARD: a genetic algorithm for recombination detection. Bioinformatics 22, 3096–3098. [DOI] [PubMed] [Google Scholar]
- 12.Temmam S, Vongphayloth K, Baquero E, Munier S, Bonomi M, Regnault B, Douangboubpha B, Karami Y, Chrétien D, Sanamxay D, et al. (2022). Bat coronaviruses related to SARS-CoV-2 and infectious for human cells. Nature 604, 330–336. [DOI] [PubMed] [Google Scholar]
- 13.Wu Z, Han Y, Wang Y, Liu B, Zhao L, Zhang J, Su H, Zhao W, Liu L, Bai S, et al. (2022). A comprehensive survey of bat sarbecoviruses across China in relation to the origins of SARS-CoV and SARS-CoV-2. Natl Sci Rev, nwac213. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 14.Pekar JE, Magee A, Parker E, Moshiri N, Izhikevich K, Havens JL, Gangavarapu K, Malpica Serrano LM, Crits-Christoph A, Matteson NL, et al. (2022). The molecular epidemiology of multiple zoonotic origins of SARS-CoV-2. Science 377, 960–966. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 15.Chinese SARS Molecular Epidemiology Consortium (2004). Molecular evolution of the SARS coronavirus during the course of the SARS epidemic in China. Science 303, 1666–1669. [DOI] [PubMed] [Google Scholar]
- 16.Xiao K, Zhai J, Feng Y, Zhou N, Zhang X, Zou J-J, Li N, Guo Y, Li X, Shen X, et al. (2020). Isolation of SARS-CoV-2-related coronavirus from Malayan pangolins. Nature 583, 286–289. [DOI] [PubMed] [Google Scholar]
- 17.Lam TT-Y, Jia N, Zhang Y-W, Shum MH-H, Jiang J-F, Zhu H-C, Tong Y-G, Shi Y-X, Ni X-B, Liao Y-S, et al. (2020). Identifying SARS-CoV-2-related coronaviruses in Malayan pangolins. Nature 583, 282–285. [DOI] [PubMed] [Google Scholar]
- 18.Liu P, Chen W, and Chen J-P (2019). Viral metagenomics revealed Sendai virus and Coronavirus infection of Malayan Pangolins (Manis javanica). Viruses 11, 979. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 19.Duchêne S, Holmes EC, and Ho SYW (2014). Analyses of evolutionary dynamics in viruses are hindered by a time-dependent bias in rate estimates. Proc. Biol. Sci. 281. 10.1098/rspb.2014.0732. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 20.Aiewsakun P, and Katzourakis A (2016). Time-Dependent Rate Phenomenon in Viruses. J. Virol. 90, 7184–7195. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 21.Ghafari M, Simmonds P, Pybus OG, and Katzourakis A (2021). A mechanistic evolutionary model explains the time-dependent pattern of substitution rates in viruses. Curr. Biol. 31, 4689–4696.e5. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 22.Dellicour S, Bastide P, Rocu P, Fargette D, Hardy OJ, Suchard MA, Guindon S, and Lemey P (2024). How fast are viruses spreading in the wild? bioRxivorg. 10.1101/2024.04.10.588821. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 23.Clements R, Sodhi NS, Schilthuizen M, and Ng PKL (2006). Limestone Karsts of Southeast Asia: Imperiled Arks of Biodiversity. Bioscience 56, 733–742. [Google Scholar]
- 24.Hassanin A, Tu VT, Curaudeau M, and Csorba G (2021). Inferring the ecological niche of bat viruses closely related to SARS-CoV-2 using phylogeographic analyses of Rhinolophus species. Sci. Rep. 11, 14276. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 25.Fan Y, Zhao K, Shi Z-L, and Zhou P (2019). Bat Coronaviruses in China. Viruses 11. 10.3390/v11030210. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 26.Han Y, Du J, Su H, Zhang J, Zhu G, Zhang S, Wu Z, and Jin Q (2019). Identification of Diverse Bat Alphacoronaviruses and Betacoronaviruses in China Provides New Insights Into the Evolution and Origin of Coronavirus-Related Diseases. Front. Microbiol. 10, 1900. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 27.Zhou H, Ji J, Chen X, Bi Y, Li J, Wang Q, Hu T, Song H, Zhao R, Chen Y, et al. (2021). Identification of novel bat coronaviruses sheds light on the evolutionary origins of SARS-CoV-2 and related viruses. Cell 184, 4380–4391.e14. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 28.Zhou H, Chen X, Hu T, Li J, Song H, Liu Y, Wang P, Liu D, Yang J, Holmes EC, et al. (2020). A Novel Bat Coronavirus Closely Related to SARS-CoV-2 Contains Natural Insertions at the S1/S2 Cleavage Site of the Spike Protein. Curr. Biol. 30, 2196–2203.e3. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 29.Trovão NS, Suchard MA, Baele G, Gilbert M, and Lemey P (2015). Bayesian Inference Reveals Host-Specific Contributions to the Epidemic Expansion of Influenza A H5N1. Mol. Biol. Evol. 32, 3264–3275. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 30.Funakoshi K, and Maeda F (2003). Foraging activity and night-roost usage in the Japanese greater horseshoe bat, Rhinolophus ferrumequinum nippon. Mammal Study 28, 1–10. [Google Scholar]
- 31.Bontadina F, Schofield H, and Naef-Daenzer B (2002). Radio-tracking reveals that lesser horseshoe bats (Rhinolophus hipposideros) forage in woodland. J. Zool. 258, 281–290. [Google Scholar]
- 32.Goiti U, Aihartza JR, Almenar D, Salsamendi E, and Garin I (2006). Seasonal foraging by Rhinolophus euryale (Rhinolophidae) in an Atlantic rural landscape in northern Iberian Peninsula. Acta Chiropt. 8, 141–155. [Google Scholar]
- 33.Rossiter SJ, Jones G, Ransome RD, and Barratt EM (2000). Parentage, reproductive success and breeding behaviour in the greater horseshoe bat (Rhinolophus ferrumequinum). Proc. Biol. Sci. 267, 545–551. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 34.Henley L, Jones O, Mathews F, and Woolley TE (2021). Coupling Bayesian theory and static acoustic detector data to model bat motion and locate roosts. arXiv [q-bio.QM]. [Google Scholar]
- 35.de Klerk A, Swanepoel P, Lourens R, Zondo M, Abodunran I, Lytras S, MacLean OA, Robertson D, Kosakovsky Pond SL, Zehr JD, et al. (2022). Conserved recombination patterns across coronavirus subgenera. Virus Evol 8, veac054. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 36.MacLean OA, Lytras S, Weaver S, Singer JB, Boni MF, Lemey P, Kosakovsky Pond SL, and Robertson DL (2021). Natural selection in the evolution of SARS-CoV-2 in bats created a generalist virus and highly capable human pathogen. PLoS Biol. 19, e3001115. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 37.Shi JJ, and Rabosky DL (2015). Speciation dynamics during the global radiation of extant bats. Evolution 69, 1528–1545. [DOI] [PubMed] [Google Scholar]
- 38.Russo D, Jones G, and Migliozzi A (2002). Habitat selection by the Mediterranean horseshoe bat, Rhinolophus euryale (Chiroptera: Rhinolophidae) in a rural area of southern Italy and implications for conservation. Biol. Conserv. 107, 71–81. [Google Scholar]
- 39.Lau SKP, Li KSM, Huang Y, Shek C-T, Tse H, Wang M, Choi GKY, Xu H, Lam CSF, Guo R, et al. (2010). Ecoepidemiology and complete genome comparison of different strains of severe acute respiratory syndrome-related Rhinolophus bat coronavirus in China reveal bats as a reservoir for acute, self-limiting infection that allows recombination events. J. Virol. 84, 2808–2819. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 40.Bihari Z (2001). Characteristics of the northernmost population ofRhinolophus ferrumequinum in the Carpathian Basin. Acta Theriol. 46, 13–21. [Google Scholar]
- 41.Guan Y, Zheng BJ, He YQ, Liu XL, Zhuang ZX, Cheung CL, Luo SW, Li PH, Zhang LJ, Guan YJ, et al. (2003). Isolation and characterization of viruses related to the SARS coronavirus from animals in southern China. Science 302, 276–278. [DOI] [PubMed] [Google Scholar]
- 42.Worobey M (2021). Dissecting the early COVID-19 cases in Wuhan. Science 374, 1202–1204. [DOI] [PubMed] [Google Scholar]
- 43.Worobey M, Levy JI, Serrano LM, Crits-Christoph A, Pekar JE, Goldstein SA, Rasmussen AL, Kraemer MUG, Newman C, Koopmans MPG, et al. (2022). The Huanan Seafood Wholesale Market in Wuhan was the early epicenter of the COVID-19 pandemic. Science, abp8715. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 44.Crits-Christoph A, Levy JI, Pekar JE, Goldstein SA, Singh R, Hensel Z, Gangavarapu K, Rogers MB, Moshiri N, Garry RF, et al. (2024). Genetic tracing of market wildlife and viruses at the epicenter of the COVID-19 pandemic. Cell 187, 5468–5482.e11. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 45.Ward HL, Ransome RD, Jones G, and Rossiter SJ (2014). Determinants and patterns of reproductive success in the greater horseshoe bat during a population recovery. PLoS One 9, e87199. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 46.Dragu A Species structure of the bat community hibernating in Muierilor Cave (Southern Carpathians, Romania: ). https://adatbank.ro/denever/cikkek/Dragu_2009.pdf. [Google Scholar]
- 47.López-Roig M, Bourhy H, Lavenir R, and Serra-Cobo J (2014). Seroprevalence dynamics of European bat lyssavirus type 1 in a multispecies bat colony. Viruses 6, 3386–3399. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 48.Hayman DTS, Cryan PM, Fricker PD, and Dannemiller NG (2017). Long‐term video surveillance and automated analyses reveal arousal patterns in groups of hibernating bats. Methods Ecol. Evol. 8, 1813–1821. [Google Scholar]
- 49.Stawski C, Willis CKR, and Geiser F (2014). The importance of temporal heterothermy in bats. J. Zool. [Google Scholar]
- 50.Hughes AC, Satasook C, Bates PJJ, Bumrungsri S, and Jones G (2011). Explaining the causes of the zoogeographic transition around the Isthmus of Kra: using bats as a case study. J. Biogeogr. 38, 2362–2372. [Google Scholar]
- 51.Wells HL, Letko M, Lasso G, Ssebide B, Nziza J, Byarugaba DK, Navarrete-Macias I, Liang E, Cranfield M, Han BA, et al. (2021). The evolutionary history of ACE2 usage within the coronavirus subgenus Sarbecovirus. Virus Evol 7, veab007. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 52.Latinne A, Hu B, Olival KJ, Zhu G, Zhang L, Li H, Chmura AA, Field HE, Zambrana-Torrelio C, Epstein JH, et al. (2020). Origin and cross-species transmission of bat coronaviruses in China. Nat. Commun. 11, 4235. [DOI] [PMC free article] [PubMed] [Google Scholar] [Retracted]
- 53.Suzuki J, Sato R, Kobayashi T, Aoi T, and Harasawa R (2014). Group B betacoronavirus in rhinolophid bats, Japan. J. Vet. Med. Sci. 76, 1267–1269. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 54.Kasibhatla SM, Kinikar M, Limaye S, Kale MM, and Kulkarni-Kale U (2020). Understanding evolution of SARS-CoV-2: A perspective from analysis of genetic diversity of RdRp gene. J. Med. Virol. 92, 1932–1937. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 55.Luk HKH, Li X, Fung J, Lau SKP, and Woo PCY (2019). Molecular epidemiology, evolution and phylogeny of SARS coronavirus. Infect. Genet. Evol. 71, 21–30. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 56.Liu P, Jiang J-Z, Wan X-F, Hua Y, Li L, Zhou J, Wang X, Hou F, Chen J, Zou J, et al. (2020). Are pangolins the intermediate host of the 2019 novel coronavirus (SARS-CoV-2)? PLoS Pathog. 16, e1008421. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 57.Liu P, Jiang J-Z, Wan X-F, Hua Y, Li L, Zhou J, Wang X, Hou F, Chen J, Zou J, et al. (2021). Correction: Are pangolins the intermediate host of the 2019 novel coronavirus (SARS-CoV-2)? PLoS Pathog. 17, e1009664. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 58.Lemey P, Ruktanonchai N, Hong SL, Colizza V, Poletto C, Van den Broeck F, Gill MS, Ji X, Levasseur A, Oude Munnink BB, et al. (2021). Untangling introductions and persistence in COVID-19 resurgence in Europe. Nature 595, 713–717. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 59.Katoh K, and Standley DM (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 60.Hall TA. 1999. BioEdit a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. In: Nucleic Acids Symposium Series 41. p. 95–98. Oxford: Oxford University Press. [Google Scholar]
- 61.Kosakovsky Pond SL, Poon AFY, Velazquez R, Weaver S, Hepler NL, Murrell B, Shank SD, Magalis BR, Bouvier D, Nekrutenko A, et al. (2020). HyPhy 2.5-A Customizable Platform for Evolutionary Hypothesis Testing Using Phylogenies. Mol. Biol. Evol. 37, 295–299. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 62.Di Franco A, Poujol R, Baurain D, and Philippe H (2019). Evaluating the usefulness of alignment filtering methods to reduce the impact of errors on evolutionary inferences. BMC Evol. Biol. 19, 21. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 63.Nguyen L-T, Schmidt HA, von Haeseler A, and Minh BQ (2015). IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 64.Minh BQ, Nguyen MAT, and von Haeseler A (2013). Ultrafast approximation for phylogenetic bootstrap. Mol. Biol. Evol. 30, 1188–1195. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 65.Huerta-Cepas J, Serra F, and Bork P (2016). ETE 3: Reconstruction, Analysis, and Visualization of Phylogenomic Data. Mol. Biol. Evol. 33, 1635–1638. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 66.Zhou XH, and Gao S (1997). Confidence intervals for the log-normal mean. Stat. Med. 16, 783–790. [DOI] [PubMed] [Google Scholar]
- 67.Suchard MA, Lemey P, Baele G, Ayres DL, Drummond AJ, and Rambaut A (2018). Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol 4, vey016. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 68.Ayres DL, Cummings MP, Baele G, Darling AE, Lewis PO, Swofford DL, Huelsenbeck JP, Lemey P, Rambaut A, and Suchard MA (2019). BEAGLE 3: Improved Performance, Scaling, and Usability for a High-Performance Computing Library for Statistical Phylogenetics. Syst. Biol. 68, 1052–1061. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 69.Hasegawa M, Kishino H, and Yano T (1985). Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22, 160–174. [DOI] [PubMed] [Google Scholar]
- 70.Ferreira MAR, and Suchard MA (2008). Bayesian analysis of elapsed times in continuous-time Markov chains. The Canadian Journal of Statistics / La Revue Canadienne de Statistique 36, 355–368. [Google Scholar]
- 71.Rambaut A, Drummond AJ, Xie D, Baele G, and Suchard MA (2018). Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Syst. Biol. 67, 901–904. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 72.Gill MS, Lemey P, Faria NR, Rambaut A, Shapiro B, and Suchard MA (2013). Improving Bayesian population dynamics inference: a coalescent-based model for multiple loci. Mol. Biol. Evol. 30, 713–724. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 73.Ji X, Fisher AA, Su S, Thorne JL, Potter B, Lemey P, Baele G, and Suchard MA (2021). Scalable Bayesian divergence time estimation with ratio transformations. arXiv [q-bio.PE]. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 74.Ji X, Zhang Z, Holbrook A, Nishimura A, Baele G, Rambaut A, Lemey P, and Suchard MA (2020). Gradients Do Grow on Trees: A Linear-Time O(N)-Dimensional Gradient for Statistical Phylogenetics. Mol. Biol. Evol. 37, 3047–3060. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 75.Fisher AA, Ji X, Zhang Z, Lemey P, and Suchard MA (2021). Relaxed Random Walks at Scale. Syst. Biol. 70, 258–267. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 76.Dellicour S, Rose R, Faria NR, Lemey P, and Pybus OG (2016). SERAPHIM: studying environmental rasters and phylogenetically informed movements. Bioinformatics 32, 3204–3206. [DOI] [PubMed] [Google Scholar]
- 77.GBIF.org (06 May 2023) GBIF Occurrence Download. 10.15468/dl.xc2jvq. [DOI]
- 78.GBIF.org (08 May 2022) GBIF Occurrence Download. 10.15468/dl.82sgat. [DOI]
- 79.Tanalgo KC, Tabora JAG, de Oliveira HFM, Haelewaters D, Beranek CT, Otálora-Ardila A, Bernard E, Gonçalves F, Eriksson A, Donnelly M, et al. (2022). DarkCideS 1.0, a global database for bats in karsts and caves. Sci Data 9, 155. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 80.Alroy J (2019). Latitudinal gradients in the ecology of New World bats. Glob. Ecol. Biogeogr. 28, 784–792. [Google Scholar]
- 81.Becker DJ, Crowley DE, Washburne AD, and Plowright RK (2019). Temporal and spatial limitations in global surveillance for bat filoviruses and henipaviruses. Biol. Lett. 15, 20190423. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 82.Di Gregorio C, Iannella M, and Biondi M (2021). Revealing the role of past and current climate in shaping the distribution of two parapatric European bats, Myotis daubentonii and M. capaccinii. The European Zoological Journal 88, 669–683. [Google Scholar]
- 83.Jargalsaikhan A, Ariunbold B, Dalannast M, Purevee E, Paek WK, Tsagaan K, and Ganbold O (2022). Molecular phylogenetic analysis of bats in the family Vespertilionidae in Mongolia. Journal of Asia-Pacific Biodiversity 15, 329–335. [Google Scholar]
- 84.Kruskop SV, Artyushin IV, Yuzefovich AP, Undrakhbayar E, Speranskya AS, Lisenkova AA, Bannikova AA, and Lebedev VS (2020). Genetic Diversity of Mongolian Long-Eared Bats (Plecotus; Vespertilionidae; Chiroptera). acta 22, 243–255. [Google Scholar]
- 85.Kuo H-C, Chen S-F, Fang Y-P, Flanders J, and Rossiter SJ (2014). Comparative rangewide phylogeography of four endemic Taiwanese bat species. Mol. Ecol. 23, 3566–3586. [DOI] [PubMed] [Google Scholar]
- 86.Mokrani Y, Univ. of Oum El-Bouaghi, Algeria, Mimeche, F., Nouidjem, Y., Saheb, M., Univ. of M’sila, Algeria, Univ. of M’sila, Algeria, and Univ. of Oum El-Bouaghi, Algeria (2018). Rapid assessment of cave-dwelling bat diversity in the Chebket ES-Sellaoua Mountains (Eastern Algeria). Arx. Misc. Zool. 16, 112–120. [Google Scholar]
- 87.Raman S, Shameer TT, Pooja U, and Hughes AC (2022). Identifying priority areas for bat conservation in the Western Ghats mountain range, peninsular India. J. Mammal. 104, 49–61. [Google Scholar]
- 88.Scherrer D, Christe P, and Guisan A (2019). Modelling bat distributions and diversity in a mountain landscape using focal predictors in ensemble of small models. Divers. Distrib. 25, 770–782. [Google Scholar]
- 89.Hughes AC, Satasook C, Bates PJJ, Bumrungsri S, and Jones G (2012). The projected effects of climatic and vegetation changes on the distribution and diversity of Southeast Asian bats. Glob. Chang. Biol. 18, 1854–1865. [Google Scholar]
- 90.Hughes AC (2017). Mapping priorities for conservation in Southeast Asia. Biol. Conserv. 209, 395–405. [Google Scholar]
- 91.Brun P, Zimmermann NE, Hari C, Pellissier L, Karger D (2022): Data from: CHELSA-BIOCLIM+ A novel set of global climate-related predictors at kilometre-resolution. EnviDat. 10.16904/envidat.332. [DOI] [Google Scholar]
- 92.Title PO, and Bemmels JB (2018). ENVIREM: an expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modeling. Ecography (Cop.) 41, 291–307. [Google Scholar]
- 93.Trabucco A, Zomer RJ (2019): Global High-Resolution Soil-Water Balance. figshare. Dataset. 10.6084/m9.figshare.7707605.v3. [DOI] [Google Scholar]
- 94.Lang N, Jetz W, Schindler K, Wegner JD( 2019) A high-resolution canopy height model of the Earth. Available at https://langnico.github.io/globalcanopyheight/. [DOI] [PMC free article] [PubMed]
- 95.Potapov P, Li X, Hernandez-Serna A, Tyukavina A, Hansen MC, Kommareddy A, Pickens A, Turubanova S, Tang H, Silva CE, et al. (2021). Mapping global forest canopy height through integration of GEDI and Landsat data. Remote Sens. Environ. 253, 112165. [Google Scholar]
- 96.Crowther TW, Glick HB, Covey KR, Bettigole C, Maynard DS, Thomas SM, Smith JR, Hintler G, Duguid MC, Amatulli G, et al. (2015). Mapping tree density at a global scale. Nature 525, 201–205. [DOI] [PubMed] [Google Scholar]
- 97.Hengl T, Mendes de Jesus J, Heuvelink GBM, Ruiperez Gonzalez M, Kilibarda M, Blagotić A, Shangguan W, Wright MN, Geng X, Bauer-Marschallinger B, et al. (2017). SoilGrids250m: Global gridded soil information based on machine learning. PLoS One 12, e0169748. [DOI] [PMC free article] [PubMed] [Google Scholar]
- 98.Phillips SJ, Dudík M, and Schapire RE (2004). A maximum entropy approach to species distribution modeling. In Proceedings of the twenty-first international conference on Machine learning ICML ‘04. (Association for Computing Machinery; ), p. 83. [Google Scholar]
Associated Data
This section collects any data citations, data availability statements, or supplementary materials included in this article.
Supplementary Materials
Data S1: Sarbecovirus genomes metadata and genome accessions, related to Figures 1–6. Data S1.1: Temporal, geographic and acknowledgement information for the 250 sarbecovirus genomes used in the paper. Sequences used for each SARS-CoV-1-like and SARS-CoV-2-like GARD recombination analysis and the virus group of each genome (SARS-CoV-1-like, SARS-CoV-2-like or outgroup) are also indicated in the table. Data S1.2: Genome accessions for all sarbecovirus genomes used in the paper.
Figure S1. Schematics describing the closest-inferred ancestor, direct ancestor, and consequences of recombination, related to Figures 1–6. a, Phylogenies illustrating the relationship between a SARS-CoV, the direct ancestor, and the closest-inferred ancestor. The closest-inferred ancestor could be the tMRCA of a given SARS-CoV and one non-human sarbecovirus (top) or the tMRCA of the SARS-CoV and multiple non-human sarbecoviruses (bottom); the closest-inferred ancestor is the parent node of the SARS-CoV. The direct ancestor is the bat sarbecovirus that was involved in each cross-species transmission event leading to the SARS-CoV and could have existed at any point on the branch leading from the closest-inferred ancestor to the SARS-CoV. b, Phylogenies and example genomes illustrating the effects of recombination in the phylogenetic inference when a whole-genome alignment is used for the inference (top) and how accounting for recombination can reveal multiple evolutionary histories (bottom) for recombinant genomes (genome C) when the alignment is split by recombination breakpoints into non-recombinant regions (NRRs).
Figure S2. Workflow to perform phylogenetic and phylogeographic analyses, related to Figures 1–6. Starting with 250 sarbecovirus genomes, we align and perform recombination analysis to determine the number of SARS-CoV-1-like and SARS-CoV-2-like viruses and their corresponding NRRs for subsequent analyses. We then calibrate the molecular clock for SARS-CoV-1-like viruses and SARS-CoV-2-like viruses using SARS-CoV-1 genomes and SARS-CoV-2 genomes (these datasets are separate from the SARS-CoV-1-like viruses and SARS-CoV-2-like viruses), respectively. Finally, we perform Bayesian divergence time estimation, recCA inference, and phylogeography analyses, each of which comprised both primary and sensitivity analyses. Ovals indicate the number of sarbecovirus genomes included in subsequent steps. Rectangles without inner vertical lines indicate analyses performed; rectangles with inner vertical lines indicate primary versus sensitivity analyses, with non-bolded text describing the details of the analysis performed.
Figure S3. Breakpoint pair distance distributions of the SARS-CoV-1-like and SARS-CoV-2-like viruses, related to Figure 1. Mean and 90% confidence intervals of the log-normal breakpoint pair distance distributions for the observed SARS-CoV-1-like and SARS-CoV-2-like estimated breakpoints and the simulated randomly sampled breakpoint sets.
Figure S4. Substitution rates for the SARS-CoV-1-like and SARS-CoV-2-like viruses and time of the closest-inferred ancestor using additional genomes or different molecular clock priors, related to Figure 1. Substitution rates (substitutions/site/year) across the (a) 31 NRRs for the SARS-CoV-1-like viruses and (b) 44 NRRs for the SARS-CoV-2-like viruses. Prior distributions are shown in pink and yellow and the dashed lines indicate two standard deviations from the mean. Posterior distributions are shown in blue and green and solid lines indicate the 95% highest posterior density (HPD) surrounding the median. Refer to Table S1 for numerical values. c, Time of the closest-inferred ancestor of SARS-CoV-1 across the 31 NRRs of the SARS-CoV-1-like viruses with the calibrated molecular clock from the SARS-CoV-1 analysis and when adding in two SARS-CoV-1 genomes from civets and and two SARS-CoV-1 genomes from humans. d, Time of the closest-inferred ancestor of SARS-CoV-2 across the 44 NRRs of the SARS-CoV-2-like viruses with the calibrated molecular clock from the SARS-CoV-2 analysis with early 2020 sampling, and with the calibrated molecular clock from the SARS-CoV-2 analysis with late 2020 sampling. The dots indicate the median and the lines are the 95% HPD. Refer to Table S1 for numerical values.
Figure S5. The recCA (recombinant common ancestor) of SARS-CoV-1 over time when including additional SARS-CoV-1 genomes and SARS-CoV-2 over time when using rate priors based on SARS-CoV-2 genomes sampled from early 2020, related to Figure 2. a, The similarity of the recCA of SARS-CoV-1 and the most closely related sarbecovirus among the SARS-CoV-1-like viruses to SARS-CoV-1 as a function of time. b, The similarity of each NRR of the recCA of SARS-CoV-1 to SARS-CoV-1 as a function of time. c, The similarity of the recCA of SARS-CoV-2 and the most closely related sarbecovirus in the SARS-CoV-2 viruses to SARS-CoV-2 as a function of time. The dashed vertical line is the sampling date of the earliest sampled genome of SARS-CoV-2. d, The similarity of each NRR of the recCA of SARS-CoV-2 to SARS-CoV-2 as a function of time. The right panels in (a, c) are zoomed-in panels of the dashed box in the left panels. Each square in (b, d) represents the similarity of the NRR of the recCA to its relative SARS-CoV at one point in time and is colored by date of the most recent non-human sarbecovirus in the reconstruction.
Figure S6. Posterior rank distributions of the tip branches from the SARS-CoV-1-like and SARS-CoV-2-like virus tip-dated phylogeographies, related to Figure 5. Violins sorted by mean, with violin corresponding to (a) SARS-CoV-1 colored blue and (b) SARS-CoV-2 colored pink.
Figure S7. Phylogenies including hypothetical unsampled taxa, related to Figure 6. Time-scaled maximum clade credibility (MCC) phylogenies of a single NRR of (a) the SARS-CoV-1-like viruses (NRR 19) and (b) and the SARS-CoV-2-like viruses (NRR 21) when including unsampled taxa in the phylogeographic reconstruction. The phylogenies are shown on a log-transformed time-scale, and horizontal segments displayed at each internal node indicate the 95% HPD interval of the estimated node age. The shaded box highlights the respective SARS-CoV and the two unsampled taxa, which are indicated by the dashed lines. Note that the respective SARS-CoV and two unsampled taxa are constrained to be monophyletic. Refer to Fig. S8 for sampling locations of the unsampled taxa.
Figure S8. Inferred location of the closest-inferred ancestor comparison when including hypothetical unsampled taxa, related to Figures 4–6. The estimated position of the closest-inferred ancestors of SARS-CoV-1 (a, b) and SARS-CoV-2 (c, d), according to continuous phylogeographic reconstructions that did not include (a, c) or did include (b, d) unsampled taxa in the analyses. Crosses indicate the locations of the unsampled taxa. As in Figure 4g-h, the estimated position of the closest-inferred ancestor is based on the continuous phylogeographic reconstruction of all NRRs. We here display the 95%, 75% and 50% HPD regions with an increasing red color darkness. Refer to Fig. S7 for examples of the phylogenetic relationship of the unsampled taxa to SARS-CoV-1 and SARS-CoV-2.
Figure S9. Comparison of dispersal velocity when including hypothetical unsampled taxa, related to Figures 4–6. Dispersal velocities for tip branches associated with (a) SARS-CoV-1 or (b) SARS-CoV-2 when excluding and including unsampled taxa in the phylogeographic analyses.
Figure S10. Additional dispersal velocity rank statistics of SARS-CoV-1-like and SARS-CoV-2-like viruses, related to Figures 5–6. Probability of posterior rank of 1 (a, c) and the Kullback–Leibler (KL) divergence from a uniform distribution (b, d) for the dispersal velocity rank distribution of each tip branch when excluding (a, b) or including (c, d) unsampled taxa in the phylogeographic analyses. The y-axes for the KL divergence are inverted because of a skewed distribution, as observed with the rank distributions for the dispersal velocities of the tip branches leading to SARS-CoV-1 and SARS-CoV-2 (Fig. S6), will be close to 0. Note the different y-axis scales for the SARS-CoV-1-like and SARS-CoV-2-like viruses in (a, c). Additional SARS-CoV-1-like viruses with fast dispersal velocities are colored and identified by their strain name.
Table S1. Posterior clock rates and times of the closest-inferred ancestor for SARS-CoV-1 and SARS-CoV-2 across each NRR, related to Figure 1. Table S1.1: Posterior clock rates across each NRR from the clock calibration using SARS-CoV-1 genomes and Bayesian divergence time analyses of SARS-CoV-1-like viruses. Table S1.2: Posterior clock rates across each NRR from the clock calibration using SARS-CoV-2 genomes sampled in late 2020 and Bayesian divergence time analyses of SARS-CoV-2-like viruses. Table S1.3: Posterior clock rates across each NRR from the clock calibration using SARS-CoV-2 genomes sampled in early 2020 and Bayesian divergence time analyses of SARS-CoV-2-like viruses. For Tables S1.1–S1.3, the clock calibration posterior results are used as a prior for the Bayesian divergence time analyses. Table S1.4: Posterior times (median and 95% HPD) of the closest-inferred ancestor of SARS-CoV-1 across each NRR. The SARS-CoV-1 clock calibration posterior serves as a prior for the analyses, with the mean and standard deviation of the clock calibration results varying across the analyses and one analysis having two additional civet SARS-CoV-1 genomes and two additional human SARS-CoV-1 genomes. Table S1.5: Posterior times (median and 95% HPD) of the closest-inferred ancestor of SARS-CoV-2 across each NRR. The SARS-CoV-2 clock calibration posterior serves as a prior for the analyses, with the standard deviation of the clock calibration results varying across the analyses.
Data Availability Statement
All genome accessions are available in Data S1. The XML files, code, and newick trees are available at https://github.com/phylogeography/sarbecovirus_phylogeography.