Recombination Patterns between SARS-CoV-2 Relatives

To reconstruct a reliable phylogeny for a set of viruses, sufficient information needs to be present in the underlying sequence alignment. Thus, even though a whole-genome alignment can be split into shorter subalignments with the aim of getting rid of all independent recombination events, it is unlikely that all subalignments can produce reliable phylogenies. To overcome this trade-off, we performed a secondary, more conservative, recombination analysis using GARD and identified the locations of 21 recombination breakpoints that strongly impact the inferred phylogenetic relationships of the analyzed sequences when mosaic patterns are ignored (supplementary table S3, Supplementary Material online). In contrast to the RDP5 method used above for assessing breakpoint clustering, the GARD method focuses on extracting recombination signal for the entire alignment, and so is better suited for producing putatively nonrecombinant phylogenies. We then determined the phylogenetic relationships of the viral sequences in each of the 22 putatively nonrecombinant genome regions bounded by each identifiable breakpoint (fig. 3A). The 20 nCoV viruses identified in the nonrecombinant whole-genome phylogeny above (fig. 1A) were used to inform the clade annotation for the 22 new nonrecombinant phylogenies.

Open in a separate window

Fig. 3

Recombination analysis and geographic distribution of Sarbecoviruses. Maximum clade credibility (MCC) dated phylogeny of RBP region 5 of 78 Sarbecoviruses (A). All tips are annotated with the geographic region the viruses have been sampled in and notable viruses are annotated with genome schematics separated into the 22 inferred RBP regions, each colored based on phylogenetic distance from SARS-CoV-2 (see scale and Materials and Methods). RBP region 21 has been removed from the schematic due to limited phylogenetic information in the alignment. The GX cluster annotated with an asterisk contains the five pangolin coronaviruses collected in Guangxi. Map of East Asia with geographic regions (provinces within China, countries outside China) colored based on Sarbecoviruses sampling (B): blue for regions with only non-nCoV clade samples, pink for regions where nCoV viruses have been sampled. Shading in the nCoV regions corresponds to phylogenetic distance from SARS-CoV-2 (see scale). Notable nCoV viruses and pangolin trafficking routes (adapted from Xu et al. [2016]) are annotated onto the map.

The two genetically closest relatives of SARS-CoV-2 that were identified shortly after its emergence were the bat sarbecoviruses, RaTG13 and subsequently RmYN02, both from samples collected in Yunnan (Zhou, Chen, et al. 2020; Zhou, Yang, et al. 2020). We find RmYN02 shares a common ancestor with SARS-CoV-2 about 40years ago and RaTG13—about 50years ago (fig. 4A) consistent with previous estimates (Boni et al. 2020; MacLean et al. 2021; Wang et al. 2021). Although SARS-CoV-2 is most similar to RmYN02 across most of its genome, the region corresponding to the first half of the RmYN02 Spike ORF appears to have been derived through recombination from a parental sequence residing outside the nCoV clade (fig. 1A). Two more viruses very recently identified in Yunnan, RpYN06 and PrC31, are most closely related to RmYN02 for part of their genomes (Li et al. 2021; Zhou et al. 2021). In the portion of the genome corresponding to recombination breakpoint partitioned (RBP) regions 2–5, the three Yunnan viruses (RmYN02, RpYN06, and PrC31) cluster with strong support in a sister clade to SARS-CoV-2 (fig. 2A and supplementary fig. S1, Supplementary Material online). This pattern suggests that bat sampling efforts in Yunnan have uncovered a related viral population that has relatively recently shared a common ancestor with SARS-CoV-2’s proximal ancestor. Molecular dating of the RBP region 5 phylogeny (corresponding to the C-terminal part of nsp3; fig. 4A) indicates that this “Yunnan cluster” shared a common ancestor with SARS-CoV-2 around 1982 (95% HPD: 1970–1994). This analysis further allows us to date the node between PrC31 and RmYN02 to 2005 (95% HPD: 1998–2010), which is one of the most recent nodes in the phylogeny (fig. 4A).

Open in a separate window

Fig. 2

Nonrecombinant topologies of SARS-CoV-2 relatives. Zoomed in regions of selected RBP region maximum likelihood phylogenies (A). Branches within the nCoV clade are colored in red and outside the nCoV clade in green. Genome schematics of close SARS-CoV-2 relatives with recombinant Spike regions (B). RBP regions 15 and 16 are highlighted and the non-nCoV subclades of the maximum likelihood phylogenies containing the relevant viruses are presented. The coloring of nonrecombinant segments indicates patristic distance to SARS-CoV-2 (see fig. 3 legend). Nodes with bootstrap confidence values below 80% have been collapsed.

Open in a separate window

Fig. 4

Molecular dating and Rhinolophus host geographic distributions. Tip-dated Bayesian phylogeny of RBP region 5 showing the nine closest relatives to SARS-CoV-2 (A). Tree nodes have been adjusted to the mean age estimates and posterior distributions are shown for each node with mean age estimate and 95% HPD confidence intervals presented to their left. Tips are annotated with the host species they were sampled in, bat silhouette colors correspond to panel (B). Geographic ranges of Rhinolophus species the SARS-CoV-2 closest relatives have been sampled in (B). Maps are restricted to East Asia and separated into province-level within China and country-level outside China.

The recombination analysis, however, reveals a much more complex evolutionary history for the rest of the PrC31 genome (Li et al. 2021). As seen in the consensus whole-genome phylogeny (fig. 1A), most of its genome clusters with viruses CoVZC45 and CoVZXC21 sampled in Zhejiang, a coastal province in East China (Lin et al. 2017; Hu et al. 2018). Across the majority of their genomes (excluding segments of Orf1ab and Spike) these viruses are members of the nCoV clade and share a common ancestor with SARS-CoV-2 that existed before 1934 (95% HPD: 1907–1957) according to molecular dating of RBP region 5 (fig. 4A). However, in RBP regions 8–12, the sequences of these viruses cluster outside the nCoV clade, and are most closely related to Zhejiang virus Longquan_140 and the HKU3 set of closely related bat sarbecoviruses sampled in Hong Kong (bordering Guangdong province) (fig. 2A and supplementary fig. S1, Supplementary Material online). The link between SARS-CoV-2’s closest relatives and viral populations in the southeast of South China becomes even more apparent in the phylogeny of RBP region 2 where Longquan_140 clusters within the nCoV clade along with CoVZC45 and CoVZXC21 (fig. 2A and supplementary fig. S1, Supplementary Material online, RBP region 2 tree). These relationships indicate ancestral movement of the nCoV viruses across large geographic ranges in China, spanning Yunnan in southwest China and Zhejiang on the east coast (fig. 3B).

As more countries initiate wildlife-infecting coronavirus sampling and sequencing efforts, the geographic range of the nCoV clade linked to bat host species will be further refined, evident from the recent reporting of bat sarbecoviruses closely related to SARS-CoV-2 from: 1) two samples collected in Cambodia from Rhinolophusshameli (RShSTT182 and RShSTT200) confirmed by whole-genome analysis (Delaune et al. 2021), and 2) five bat samples from Rhinolophusacuminatus collected in Thailand with one fully sequenced genome of virus RacCS203 (Wacharapluesadee et al. 2021). These viruses are, after the China sampled CoVs mentioned above, the next closest relatives to SARS-CoV-2 with common ancestor age estimates (using RBP region 5) around 1907 (95% HPD: 1873–1938) and 1883 (95% HPD: 1841–1921), respectively (fig. 4A). Similar to the other nCoV viruses, the recombination analysis uncovers more intricate phylogenetic relations for some parts of the genome. Notably, RShSTT182 and RShSTT200, despite being sampled in Cambodia, cluster with RaTG13 for RBP regions 8 and 9 (fig. 2A and supplementary fig. S1, Supplementary Material online), while in RBP region 4 of the genome RacCS203, from Thailand, clusters together with SARS-CoV-2 within the Yunnan clade (fig. 2A). This indicates that cocirculation and recombination between these viruses in the last few centuries is responsible for the observed patterns in their inferred evolutionary history, despite the current geographic range of at least 2,500km. This wide distribution of related viruses, including shared recombination breakpoints, highlights an important feature of bat species: Their frequently overlapping/sympatric ranges will provide ample opportunities for transmissions of viral variants from one bat species (or subspecies) to another.

Consistent with the Spike S1 recombination hotspots revealed in the initial analysis (fig. 1B and C), closest relatives of SARS-CoV-2 presented here have non-nCoV derived recombinant sequences at the start of the Spike gene (fig. 2B). Despite one collected from Yunnan, China and the other from Thailand, viruses RmYN02 and RacCS203 share a closely related non-nCoV sequence in RBP regions 15 and 16 (encompassing the Spike NTD and RBD, respectively; fig. 2B) having a distinct RBD compared with that of SARS-CoV-2. On the other hand, viruses RpYN06, PrC31, CoVZC45, and CoVZXC21 cluster within the nCoV clade for region 15 but, similar to the RmYN02 and RacCS203, form a distinct cluster in the non-nCoV clade for region 16 (fig. 2B; Wells et al. 2021). We speculate that some of the apparent patterns of recombination-mediated exchange between nCoV and non-nCoV viruses can be partly explained by sequential recombination, that is, “overprinting” of recombination events involving different ancestral parental viruses. This will occur when an nCoV virus acquires a non-nCoV genomic sequence through ancestral recombination but its progenitors cocirculating with other nCoV viruses incur subsequent recombination events that overlap portions of the original non-nCoV recombinant sequence, producing the more complex “patchy” patterns we see in the currently sampled viruses. Note, overprinting of recombination regions will result in reduced confidence in the breakpoints at deeper nodes in the phylogeny.

The finding that Sunda (also known as Malayan) pangolins, Manis javanica, nonnative to China, are the other mammal species from which nCoV sarbecoviruses have been sampled in Guangxi and Guangdong provinces in South China (Lam et al. 2020; Xiao et al. 2020), indicates these animals are likely being infected in this part of the country (fig. 3B). Pangolins are one of the most frequently trafficked animals with multiple smuggling routes leading to southern China (Xu et al. 2016). The most common routes involve moving the animals from Southeast Asia (Myanmar, Malaysia, Laos, Indonesia, Vietnam) to Guangxi, Guangdong, and Yunnan. The most likely scenario that is consistent with both the reported respiratory distress that the sampled pangolins exhibited (Liu et al. 2019; Xiao et al. 2020) and the lack of confirmed CoV infections among Sunda pangolins in Malaysia (Lee et al. 2020), is that the viruses obtained from these animals infected them (presumably from bat sources) after they were trafficked into southern China. Still, serological data of trafficked Sunda pangolins could suggest potential circulation of sarbecoviruses in the animals’ wild populations (Wacharapluesadee et al. 2021).

Although the recombination patterns inferred in the pangolin-derived virus genomes seem to be less complex than those of the bat nCoV genomes, the Guangdong Pangolin-CoV has a Spike receptor binding domain that is most similar to that of SARS-CoV-2. This finding was highlighted by Li, Giorgi, et al. (2020) and attributed to recombination between the SARS-CoV-2 and Pangolin-CoV proximal ancestors. However, based on the nucleotide divergence between the two viruses in this short Spike segment, the most likely explanation is recombination in an ancestor of RaTG13, making it more divergent than Pangolin-CoV compared with SARS-CoV-2 (Boni et al. 2020) (reflected in region 17, fig. 2A). The susceptibility of pangolins to an apparently new human coronavirus is not surprising given the well-documented generalist nature of SARS-CoV-2 (Conceicao et al. 2020), which has been found to readily transmit to multiple mammals with similar ACE2 receptors, most notably, on mink farms (Oude Munnink et al. 2021).

Genome-Specific Recombination Analysis

We first performed an analysis for detecting unique recombination events specific to individual genome sequences using the RDP (Martin and Rybicki 2000), GENECONV (Padidam et al. 1999), BOOTSCAN (Martin et al. 2005), MAXCHI (Smith 1992), CHIMAERA (Posada and Crandall 2001), SISCAN (Gibbs et al. 2000), and 3SEQ (Boni et al. 2007) methods implemented in the program RDP5 (Martin et al. 2021). Default settings were used throughout except: 1) only potential recombination events detected by three or more of the above methods, coupled with phylogenetic evidence of recombination were considered significant and 2) sequences were treated as linear. We required supporting evidence from three or more of the recombination signal detection methods because none of three methods query the same recombination signals and all have varying power to detect recombination in data sets with different degrees of sequence diversity (Posada and Crandall 2001; Posada 2002). The recombinant sequence identification, recombination breakpoint verification, and shared recombination event verification steps used are outlined in Martin et al. (2017). The approximate breakpoint positions and recombinant sequence(s) inferred for every potential recombination event, were manually checked and adjusted where necessary using the phylogenetic and recombination signal analysis features available in RDP5. Breakpoint positions were classified as undetermined if the 95% confidence interval on their location overlapped: 1) the 5′ and 3′ ends of the alignment; or 2) the position of a second detected breakpoint within the same sequence that had a lower associated P value (in such cases it could not be discounted that the actual breakpoint might not have simply been lost due to a more recent recombination event). All of the remaining breakpoint positions were manually checked and adjusted when necessary using the BURT method with the MAXCHI matrix and LARD two breakpoint scan methods (Holmes et al. 1999) used to resolve ties. A putatively nonrecombinant version of the original whole-genome alignment was reconstructed by excluding all minor parent sequence segments based on the supervised RDP5 analysis.

Recombination Hotspot Analysis

The distribution of 236 unambiguously detected breakpoint positions defining 160 unique recombination events based on the RDP5 analysis described above were analyzed for evidence of recombination hotspots and coldspots using the permutation-based RRT (Simon-Loriere et al. 2009) and BDT (Heath et al. 2006). The RRT accounts for site-to-site variations in the detectability of individual recombination events and examines the distribution of point estimates of pairs of breakpoint locations bounding each of the unique recombination events detected by RDP5. Rather than using point estimates of recombination breakpoint locations, the BDT accounts for underlying uncertainties in the estimation of individual breakpoint locations as determined from the state transition likelihoods yielded by the hidden Markov model-based recombination breakpoint detection method, BURT (described in the RDP5 program manual at http://web.cbio.uct.ac.za/~darren/rdp.html).

To verify whether the recombination breakpoint clusters detected with these tests were consistent with the presence of recombination hotspots, we simulated recombination with SANTA-SIM (Jariani et al. 2019). Four data sets of 100×10 kb long sequences that had: 1) approximately the same degree of genetic diversity as the analyzed sarbecovirus alignment and 2) approximately the same numbers of detectable recombination events and recombination breakpoints per nucleotide as those detected in the analyzed sarbecovirus alignment. The SANTA-SIM settings used were: population size=4,500, inoculum=all, mutation rate=3.5×10⁻⁵, rate bias matrix=(0.42, 2.49, 0.29, 1.73, 0.23, 4.73, 6.99, 9.20, 0.60, 1.02, 2.56, 0.88), dual infection probability=0.1, background recombination probability=0.06, and generation number=5,000. Simulated recombination events had a maximum of two breakpoints: a setting that required the use of a slightly modified version of SANTA-SIM that can be obtained from https://github.com/phillipswanepoel/santa-sim/tree/Recomb_and_align. Whereas one of the four data sets had no simulated recombination hotspots, the other three each had a single 100-nt-long hotspot between alignment positions 6000 and 6100 wherein recombination frequencies were 4×, 8×, or 16× higher than the background level.

All data sets were analyzed for recombination by RDP5 without any supervision, and RRT and BDT plots were produced for each data set (all with the same program settings used to analyze the actual sarbecovirus data set).

The true positive rate of the BDT was estimated as the proportion of 200-nt windows centered on nucleotides between positions 6000 and 6100, that is, within the simulated hotspot, that contained a number of breakpoints greater than the upper bound of the 99% confidence interval of the breakpoint clustering distribution expected under random recombination (e.g., indicated by the light gray areas of the plots in fig. 1C). Since a 200-nt sliding window was used for both breakpoint clustering tests, all windows overlapping with the hotspot (positions 5801 to positions 6299) were ignored when determining the BDT and RRT false positive rates. The false positive rate of BDT was calculated as the proportion (across all 100 simulated alignments of each of the four data sets) of the examined 200-nt windows centered on nucleotides outside region 5801–6299 that contained a number of breakpoints greater than the upper bound of the 99% confidence interval of the breakpoint clustering distribution expected under random recombination.

Similarly, the true positive rate of the RRT was estimated as the proportion, across all 100 simulated alignments in a data set, of 200-nt windows centered on nucleotides between positions 6000 and 6100, that is, within the simulated hotspot, that had associated breakpoint clustering permutation P values <0.01 (e.g., indicated by the upper bound of the light gray area of the plot in fig. 1B). The RRT false positive rate was calculated as the proportion, across all 100 simulated alignments in a data set, of the examined 200-nt windows centered on nucleotides outside region 5801–6339 that had associated permutation P values <0.01.

The true and false positive rates for BDT and RRT with respect to identifying the presence of the simulated recombination hotspots are indicated in supplementary table S2, Supplementary Material online. Note that, due to the nature of the simulations, it was not guaranteed that even with perfect recombination detection power and accuracy: 1) the recombination hotspot regions would contain any detectable excess of recombination breakpoints, and 2) the “normal” genome regions would contain no breakpoint clusters. What these simulations capture is the power of the two clustering tests to indirectly infer the locations of actual recombination hotspot regions that, due to chance during the simulation process, might not even contain any detectable recombination breakpoints. Nevertheless, as expected, the hotspot detection power of both BDT and RRT increases substantially with the intensity of the simulated recombination hotspots: from ~10% for both tests with a 4× increase in simulated breakpoint probabilities within the 100-nt hotspot region to ~60% for a 16× increase in breakpoint probabilities within the hotspot region. It is also noteworthy that the false positive rates for both tests are likely between 1.5 and 2× higher than the expected rate of 0.01 (which is expected given that the windows containing breakpoint clusters exceeding the 99% confidence interval were used to identify breakpoint hotspots). This false positive rate may not seem very high but, for a long alignment such as that examined for the sarbecoviruses that can be broken into ~150 non-overlapping 200-nt windows, it indicates that for such an alignment we might expect to find on average two to three significant clusters of breakpoints that are in fact not associated with any elevation in the underlying recombination rate.

Whole-Genome Alignment Recombination Analysis

Next, we sought to conservatively examine the entire genome alignment for the subset of recombination breakpoints that had the largest impacts on the inferred evolutionary relationships between the analyzed sarbecoviruses using the GARD method (Kosakovsky et al. 2006) implemented in Hyphy v2.5.29 (Kosakovsky Pond et al. 2020). Model goodness of fit was evaluated using the small sample Akaike Inference Criterion (c-AIC) (Akaike 1998). To improve computational efficiency and statistical efficiency (as GARD requires more statistical evidence of recombination for larger phylogenies, and the minimal length of detectable nonrecombinant fragments increases with the number of sequences) and focus on the closest relatives of SARS-CoV-2, 22 of the 78 viruses that are closest to SARS-CoV-2 or had preliminary evidence of clustering near detected interclade recombinants were included in the GARD analysis (supplementary table S1, Supplementary Material online). Only breakpoints present in more than 2/3 of the 64 GARD consecutive step-up models were retained to produce a final set of 21 likely breakpoints (positions corresponding to the SARS-CoV-2 reference genome Wuhan-Hu-1 in order: 1680, 3093, 3649, 4973, 8208, 11445, 12622, 14401, 15954, 16923, 19965, 20518, 21198, 21411, 22460, 23396, 24144, 24843, 26323, 27388, 27685). Based on these the whole-genome alignment was split into 22 RBP regions. The position of each region on the alignment and relative to the SARS-CoV-2 genome as well as their length is presented in supplementary table S3, Supplementary Material online.

We further used the GARD recombination analysis to validate the RDP5 recombination hotspot analysis. We performed a permutation test of breakpoint clustering by fixing the number of all inferred breakpoints (64) and the location of the 13,550 variable sites in the alignment. Then defined a sliding window so that each window would have an average of one breakpoint in it (alignment length/64) producing 474 windows. N=10,000 replicates were drawn where 64 variable sites were randomly chosen from one of the breakpoints. For each sliding window, we tabulated the distribution of randomly drawn breakpoints in the window. Two hotspots and 17 coldspot windows were identified, presented in supplementary figure S2, Supplementary Material online. This analysis is not expected to produce results identical to the RDP5-based hotspot analysis, since the GARD method does not distinguish between potential breakpoints in very near genomic proximity, so this post hoc test is unlikely to identify clustering of unique breakpoints that are very close to one another (in contrast to the RDP5 approach).

Phylogenetic Reconstruction

The phylogeny of each RBP alignment region based on the GARD analysis and the nonrecombinant whole-genome based on the RDP5 analysis were reconstructed using iqtree version 1.6.12 (Nguyen et al. 2015) under a general time reversible (GTR) substitution model assuming invariable sites and a four-category Γ distribution. Tree node confidence was determined using 10,000 ultrafast bootstrap replicates.

Based on the nonrecombinant whole-genome phylogeny, 20 viruses form a monophyletic nCoV clade (fig. 1A). To illustrate the distance of each virus from SARS-CoV-2 for each GARD determined genomic region, we defined the nCoV clade on each phylogeny as the subset of the aforementioned 20 nCoV viruses forming a monophyly with SARS-CoV-2 in each phylogeny. The rest of the viruses were classified as members of the non-nCoV clade for each RBP region. We then used an arbitrary tip distance scale normalized between all phylogenies so distances are comparable between regions. For each maximum likelihood tree, the patristic distance between each tip and SARS-CoV-2 is calculated using ETE 3 (Huerta-Cepas et al. 2016) as d1 for members of the nCoV clade and d2 for members of the non-nCoV clade. The distances are then normalized so that for nCoV clade members range between 0.1 and 1.1 (1.1 being SARS-CoV-2 itself and 0.1 being the most distant tip from SARS-CoV-2 within the nCoV clade) and between −0.1 and −1.1 for non-nCoV members (−0.1 being the closest non-nCoV virus to SARS-CoV-2 and −1.1 the most distant), as follows:

With d′₁ and d′₂ being the normalized values for each clade, variables denoted with “min” being the smallest distance and variables denoted with “max” being the largest distance in each given set.

Phylogenies were visualized using FigTree (http://tree.bio.ed.ac.uk/software/figtree/) and ETE 3 (Huerta-Cepas et al. 2016).

Molecular Dating

To provide temporal information to the phylogenetic history of the viruses, we performed a Bayesian phylogenetic analysis on RBP region 5, using BEAST v1.10.4 (Suchard et al. 2018). This region was selected due to its length, being one of the two longest nonrecombinant regions in the analysis (3,238bp), and because all 20 nCoV viruses form a monophyly in the respective tree. Based on the observation of an increased evolutionary rate specific to the deepest branch of the nCoV clade reported in MacLean et al. (2021) (MacLean et al. 2021), we adopted the same approach of fitting a separate local clock model to that branch from the rest of the phylogeny. A normal rate distribution with mean 5×10⁻⁴ and SD 2×10⁻⁴ was used as an informative prior on all other branches. The lineage containing the BtKY72 and BM48-31 bat viruses was constrained as the outgroup to maintain overall topology. Codon positions were partitioned and a GTR+Γ substitution model was specified independently for each partition. The maximum likelihood phylogeny reconstructed previously for RBP region 5 was used as a starting tree (rooted at the BtKY72 and BM48-31 clade). A constant size coalescent model was used for the tree prior and a lognormal prior with a mean of 6 and SD of 0.5 was specified on the population size. Two independent MCMC runs were performed for 500 million states for the data set. The two chains were inspected for convergence and combined using LogCombiner (Drummond and Rambaut 2007) using a 10% burn-in for each chain. The effective sample size for all estimated parameters was above 200.

We would like to thank all the authors who have kindly deposited and shared genome data on GISAID. Credit also needs to be given to the surveillance projects for generating the genome data that is available in GenBank and to the software developers for making the tools we have used freely available. A table with genome sequence acknowledgments can be found in the Supplementary Material online (supplementary table S1, Supplementary Material online). D.L.R. and J.H. are funded by the UK Medical Research Council (MRC, MC_UU_1201412) and D.L.R. by the Wellcome Trust (WT, 220977/Z/20/Z). S.L. is funded by an MRC studentship. D.M. is funded by the WT (222574/Z/21/Z). S.L.K.P. was supported in part by a grant from the National Institutes of Health (NIH, R01 AI134384 [NIH/NIAID]). X.J. is funded by the Jiangsu Province High-level Innovation and Entrepreneurship Talent Programme.