A comparative recombination analysis of human coronaviruses and implications for the SARS-CoV-2 pandemic

The SARS-CoV-2 pandemic prompts evaluation of recombination in human coronavirus (hCoV) evolution. We undertook recombination analyses of 158,118 public seasonal hCoV, SARS-CoV-1, SARS-CoV-2 and MERS-CoV genome sequences using the RDP4 software. We found moderate evidence for 8 SARS-CoV-2 recombination events, two of which involved the spike gene, and low evidence for one SARS-CoV-1 recombination event. Within MERS-CoV, 229E, OC43, NL63 and HKU1 datasets, we noted 7, 1, 9, 14, and 1 high-confidence recombination events, respectively. There was propensity for recombination breakpoints in the non-ORF1 region of the genome containing structural genes, and recombination severely skewed the temporal structure of these data, especially for NL63 and OC43. Bayesian time-scaled analyses on recombinant-free data indicated the sampled diversity of seasonal CoVs emerged in the last 70 years, with 229E displaying continuous lineage replacements. These findings emphasize the importance of genomic based surveillance to detect recombination in SARS-CoV-2, particularly if recombination may lead to immune evasion.

The emergence of SARS-CoV-2 has generated interest in role of recombination in the evolution of this and other human coronaviruses (hCoV). Recombination has been observed in many RNA viruses and is noted to occur at a higher frequency in positive-sense RNA viruses, a category that includes SARS-CoV-2 and other medically important coronaviruses [1][2][3] . From an evolutionary biology perspective, it remains unclear why recombination occurs in RNA viruses. Hypotheses include recombination being an incidental outcome of RNA polymerase function, or an evolutionary favorable purge of deleterious genotypes and/or generation of advantageous genotypes 1 .
RNA virus recombination has been associated with changes in host range, host response and virulence 1 . Identifying the presence of recombination, or predicting the risk of recombination, in viral populations of SARS-CoV-2 is critical for several reasons. First, circulating recombinants may complicate molecular diagnostic performance. Second, recombinants may cause rapid escape from naturally acquired immunity, as has been observed in the norovirus genus, which has caused pandemics due to the rapid emergence of new genotypes generated by recombination of structural genes 20 . For SARS-CoV-2, such an event may have major implications, especially if circulating recombinant results in escape from both natural and vaccine induced immunity 21 . Finally, genomic epidemiology has increasingly been shown to be an important public health tool for SARS-CoV-2 and failing to accommodate for recombinant data may lead to incorrect epidemiological inference because of possible phylogenetic incongruence 22,23 .
We undertook a comparative recombination analysis of all published SARS-CoV-1, SARS-CoV-2, MERS-CoV and seasonal hCoV (OC43, 229E, NL63 and HKU1) genomes. Specifically, we aimed to identify the frequency Recombination is relatively frequent in seasonal endemic coronaviruses and has a propensity for the non-ORF1 genes. Within the 229E, OC43, NL63 and HKU1 datasets, we noted 1, 9, 14, and 1 high confidence recombination events, respectively (Table 1). These recombination events were found in 13.6%, 20.3%, 89.2%, and 35.1% of the analyzed sequences, respectively. In the OC43 and NL63 datasets, we noted significantly more breakpoints in the non-ORF1 region of the genome containing structural genes than in the ORF1 (p = 0.0004, p = 0.0032 respectively) ( Fig. 1).
Recombinants were noted across entire clades of HKU1, NL63 and OC43 viruses, with clusters of genomes sharing identical recombination patterns, indicating spread of a recombinant virus following its emergence through a recombination event. The OC43 and NL63 datasets also contained genomes with unique recombination patterns (singletons) (Figs. 2, 3). Furthermore, we noticed some singletons falling within already recombinant clades, indicating presence of successful second generation recombination (Figs. 2, 3). A 10th OC43 event with moderate recombination evidence based on RDP4 results showed topological incongruence in subgenomic phylogenetic trees and was therefore removed for further time-scaled analyses.

MERS-CoV but not SARS-CoV-1 is characterized by frequent recombination with a propensity for non-ORF1 genes.
Within the MERS-CoV dataset, we detected 7 recombination events with a high level of confidence (detected by 3 or more methods and not explained by another evolutionary process) ( Table 1). These recombination events were found in 14.5% of all analyzed MERS-CoV genomes. Of these, 6 were found across clades, suggesting recombinants were sufficiently fit for onward transmission (Fig. 4). We noted recombinant clades defined by camel hosts, as well as camel and human hosts (Fig. 4), suggestive of interhost recombinant spread. Moreover, we noted significantly more breakpoints in the non-ORF1 region of the genome containing structural genes compared to ORF1 (p < 0.001) (Fig. 1). We noted only a single low evidence recombination event involving the structural region of the SARS-CoV1 genome (4503 -25,998 nt), but this was not confirmed by multiple methods (Table 1).
Recombination in seasonal coronaviruses substantially alters estimates of temporal structure. Our approach to identification of recombinants and their subsequent removal from the datasets led to major improvements in estimated temporal structure of hCoV-OC43, hCoV-NL63 and hCoV-HKU1 (Table 2, Fig S7-S10). These changes involved both the regression coefficient and the regression intercept, which serve as crude estimates of evolutionary rates and TMRCA, and are often used to assess the clock-likeness (linear Table 1. Frequency of recombination events detected in 229E, NL63, OC43, HKU1, MERS-CoV, SARS-CoV1 and SARS-CoV-2, stratified by level of evidence. a Randomly subsampled into 100 × n = 300 independent datasets. www.nature.com/scientificreports/ relationship of genetic distance across sequence sampling times) of the data for further analyses ( Table 2, Fig  S7-S10). In contrast, removal of the single 229E recombinant did not cause substantive change in estimated temporal structure as estimated by regression coefficient and TMRCA ( Table 2).
The current global diversity of seasonal hCoVs arose across the last 70 years. We estimated a TMRCA date (year A.D) of 1989, 1970, 1964, and 1951 for the 229E, OC43, NL63 and HKU1 lineages, respectively, indicating relatively recent emergence of the current seasonal hCoV lineages (Table 3, Supplementary Figures S1-S4). Importantly, these do not necessarily represent de novo emergence of these viruses from animal origins but rather may represent the divergence from older OC43, NL63 and HKU1 lineages, respectively. Indeed, the TMRCA for 229E is preceded by historical descriptions of the circulation of this virus in humans, which was discovered in 1962 24 . We therefore extended our 229E full genome analysis and inferred TMRCA estimates from partial genome datasets of the complete N gene (N = 101, length = 1167nt), the complete S gene (Supplementary Figure

Discussion
We performed a comprehensive recombination analysis across all medically important human coronaviruses, with an overarching aim to identify the current and future risk of recombinant emergence in the SARS-CoV-2 pandemic. We note moderate evidence for SARS-CoV-2 recombination during the first year of the COVID-19 pandemic. In other hCoV species, we note that recombination has a predilection for the non-ORF1 region of the genome containing structural genes, and is relatively frequent in most medically significant hCoV over a relatively short evolutionary timescale. These findings are timely with the announcement of a more recent unpublished recombinant SARS-CoV-2 strain 25 , as well as the recent focus on how insights can be learned from the functional implications of antigenic evolution now demonstrated in seasonal coronaviruses 26 . From more than 100,000 SARS-CoV-2 genomes, we noted 8 instances of recombination with moderate confidence. Two of these events were noted in the spike gene; none involved the furin cleavage site. In our analysis the SARS-CoV-2 estimated breakpoints were found in regions with typically low but variable GC content (median = 40%, IQR = 30-60%, range 0-100%), and breakpoints occurred at locations with variable predicted RNA secondary structure with icSHAPE scores ranging from 0 to 1.
However, all SARS-CoV-2 recombination events were flagged by the RDP4 software as being possibly driven by other processes despite support by three or more recombination detection methods. This may reflect the relatively lower viral diversity across the first year of the pandemic. Similarly, we did not detect a high-confidence recombination signal in SARS-CoV-1, a virus with a limited temporal distribution. In contrast, we show that recombination was relatively frequent in seasonal coronavirus and MERS-CoV datasets comprising a longer period of sampling, including recombinants sufficiently fit for onward transmission. Furthermore, our analyses www.nature.com/scientificreports/ show that three of the coronaviruses (MERS, OC43 and NL63) had breakpoint propensity for the non-ORF1 region of the genome containing structural genes. Several reasons for this may exist, such as inherent structural similarities of the coronaviruses in this region causing enhanced enzyme slippage, or positive selection pressure. The latter may be correlated with evasion of the human immune system, although recombination in RNA viruses is not generally thought as an evolutionary process which is driven by natural selection to favor advantageous genotypes 1 .
Our endemic coronavirus analyses also highlighted that recombination affected the estimated temporal structure of coronavirus sequence datasets, particularly in the OC43 and NL63 types. This serves as a caution for genomic epidemiology studies which do not identify and account for recombinants in SARS-CoV-2 and other coronavirus analyses. Indeed, recombination has long been known to be a cause of phylogenetic incongruence for other viruses 23 . Our time-scaled evolutionary analyses, adjusted for recombination to restore a molecular clock signal, yielded insights into the recent epidemiology of seasonal coronaviruses. We noted that the current sampled diversity of seasonal coronaviruses has emerged within a 70 year period, punctuated by new lineage emergence at intervals ranging from 5 to 20 years. For certain seasonal coronaviruses the uncertainty interval of these TMCRA estimates did overlap with the first reported cases. However, historical epidemiological data on the time-scales of human coronavirus emergence also have uncertainty. For instance, while OC43 was first isolated from a human case in 1967, potentially cross reactive sera have been identified as far back as 1965 27,28 . Moreover, caution is required in inferring that these viruses spilled over into humans at these timepoints, however, as this may reflect the divergence of new lineages from prior, unsampled older viral populations circulating in humans. In the case of 229E hCoV, we noted that the full genome TMRCA estimates were preceded by the clinical reports of 229E infection in the 1960s 24 , and a spike gene analysis incorporating older partial genome sequences from the earlier twentieth century showed lineage extinction and replacement. This is consistent with previous analyses suggesting that 229E evolution is characterized by prior lineage extinction and new lineage emergence 29 .
It is important to note that our RNA genomic recombination detection in sequence data remains a statistical estimation only, as previously discussed in detail for dengue viruses 30 . In addition, the size of our SARS-CoV-2  www.nature.com/scientificreports/ dataset was computationally prohibitive to perform recombination detection across all data, which might have resulted in missing additional recombination events. Also, unsampled data are a pervasive technical risk for recombination detection, as demonstrated by the variable finding of recombinant events across subsampled SARS-CoV2 data. Finally, our strict criteria for identifying recombination events with high confidence may have resulted in the removal of some true recombinant events. Still, these findings provide critical insights into the possible projected evolution of SARS-CoV2. Recombination in other RNA viruses has been associated with changes in host tropism, virulence or epidemiology 1 . Ongoing genomic-based COVID-19 surveillance has recently been highlighted as a critical public health tool to detect novel SARS-CoV-2 variants, such as the B.1.1.7, B.1.351, and P.1 variants 31,32 . Robust genomic surveillance will be essential for the timely detection of recombination in SARS-CoV-2, which may have implications for the diagnosis of and immunity to this pandemic pathogen.
All available full SARS-CoV-2 genomes (n = 157,439) up to October 23, 2020 were downloaded from the GISAID database 37 . These data were curated by removing (i) any partial genomes (< 90% full genome length), (ii) any genomes with > 100 continuous ambiguous base calls (Ns), and (iii) removing bat and pangolin non-human genomes. A single genome (Italy/CAM-IZSM-45946/2020) was removed due to the presence of non-IUPAC nucleotide codes. The remaining 100,296 genomes were aligned to the MN908947.3 reference using mafftc (7.471) in the MAFFT software 35 . Alignment ends (first 265 and last 229 nt) were trimmed using the trimal command 38 . This yielded a final alignment of length 29,409 bp. Datasets 46 . Following the detection of a 'recombination signal' with these methods, the approximate breakpoint positions were determined using a hidden Markov model, BURT, and the recombinant sequence identified using the PHYLPRO 47 , and VISRD 48 methods.

Recombination detection and determination of breakpoints.
For SARS-CoV2, due to the prohibitive computational demand of recombination detection analysis in a dataset of this size, we randomly subsampled, with replacement, n = 300 sequences with 100 × iterations (n = 30,000 full genomes) to perform recombination detection.
As individual recombination detection methods may have limited specificity, we developed a customized framework to ascertain the level of evidence for those recombination events detected in these datasets. Recombination events identified by only one or two methods in RDP4 were assigned a 'low' level of confidence, and those identified by at least three methods in RDP4 were assigned at least a 'moderate' level of confidence. We further assigned a 'high' level of confidence for those recombination events identified by at least three RDP4 methods with no other identified process which may have explained the recombination signal 39 . The evolutionary processes that might be misinterpreted as recombination are typically caused by a combination of mutation rate variation between sites superimposed on mutation rate variation along lineages (M. Darren, personal communication, Dec 4, 2020). For those recombination events with a high level of evidence, the location of breakpoints were plotted across the whole genome and a χ 2 test used to compare frequency of breakpoints at non-ORF1 region of the genome containing structural genes, versus non-structural ORF1. Breakpoints that could not be placed with certainty due to overlapping recombination or other reasons were not included in the analyses. We determined both the number of recombination events by hCOV type, in addition to whether such events were found in multiple genomes for that hCoV type. A single recombination event can contain several genomes, meaning that the recombination occurred within the ancestor of these genomes, followed by its successful spread. Therefore, the number of recombinant genomes may be higher than the number of observed recombination events.
Estimation of temporal structure of hCoV with and without recombinant strains. Datasets for each hCoV species underwent nucleotide model substitution selection using JModelTest2 49 , with model selection as follows: 229E = SYM + I + G, HKU1 = GTR + G + I, NL63 = GTR + I, OC43 = GTR + G + I. A maximum likelihood phylogeny was inferred using the PhyML software 50 , with aLRT for node support and tips labeled by date of collection. Root-to-tip regression was performed using the Temp-Est tool 51 , with slope coefficient and  www.nature.com/scientificreports/ intercept values used as preliminary estimates of evolutionary rate and time-to-most-common recent ancestor (TMRCA), respectively. Confirmed recombinant sequences were annotated on these phylogenies to identify recombinant clades and singletons. Phylogenetic analyses were repeated with the recombinants removed to estimate the impact of recombination on time-scales and temporal structure of hCoV evolution. For NL63, removal of recombinant genomes resulted in a small dataset and a small remaining fraction of the initial phylogenetic tree. Therefore, NL63 genomes were screened for a common region without any recombination signal. A region of 7105 nts (region 13093-20198 of the alignment) was found in which most genomes (n = 56) had no recombination signal, and this region was used for subsequent time-scaled phylogenetic analyses, with the 9 genomes with recombination breakpoints in this region removed from the dataset.
Estimating time-scales of emergence of the currently circulating seasonal hCoV lineages. We leveraged recombination-free endemic seasonal hCoV datasets, in addition to removing other root-to-tip regression outlier genomes (n = 3 + 4 identical genomes for OC43, n = 1 for NL63, n = 0 for 229e, n = 2 for HKU1), to reconstruct the time-scale of emergence of currently circulating seasonal hCoV lineages across the 229E, OC43, NL63 and HKU1 types. We focused on these coronaviruses because they are well established viruses within human populations and may serve as a model for the projected evolutionary future of SARS-CoV-2. This is in contrast to the now extinct SARS-CoV-1 virus, and MERS-CoV, the latter which continues to be defined by more sporadic and discrete spillover events 52 . Time-scaled genealogies of these viruses were inferred using the BEAST software 1.8.4 53 . To minimize the risk of model misspecification, we inferred maximum clade credibility phylogenies with combinations of demographic models (constant, exponential and skyline population models) and molecular clock models (strict versus relaxed) (Supplementary Table S4). For each hCoV dataset, the optimal combination of demographic and molecular clock model was selected by logarithmic marginal likelihoods inferred by the path-sampling/ stepping-stone method 54 .
Predicted features of underlying RNA template near putative SARS-CoV-2 recombination breakpoints. GC content around putative SARS-CoV-2 breakpoints was determined by calculating the percentage of G (guanine) and C (cytosine) bases in 5-base windows via the UCSC SARS-CoV-2 genome browser tool using 5 nt windows 55,56 . Secondary RNA structure at putative SARS-CoV-2 breakpoints was predicted using selective 2-hydroxyl acylation and profiling experiment (SHAPE) reactivity 57 , SHAPE Shannon entropy (com-SuperFold) 57 , and icSHAPE in vivo scores via the UCSC SARS-CoV-2 genome browser 58 . Disclaimer. Material has been reviewed by the Walter Reed Army Institute of Research. There is no objection to its presentation and/or publication. The View(s) expressed are those of the authors and do not necessarily reflect the official views of the Uniformed Services University of the Health Sciences, Henry M. Jackson Founda- Table 2. Root-to-tip regression coefficient and intercept of seasonal hCoV phylogenies with and without recombinants removed. a Approximates evolutionary rate (substitutions/site/year). b Approximates TMRCA. *Non-recombinant region was used, with genomes containing breakpoints in this region removed (N = 9). www.nature.com/scientificreports/ tion for the Advancement of Military Medicine, Inc., Department of Health and Human Services, the National Institutes of Health, the Departments of the Army, Navy or Air Force, the Department of Defense, or the U.S. Government.

Data availability
All genomes are available in public sequence repositories (GISAID and ViPR) and their accession numbers are reported in Supplementary Information. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. This is a U.S. Government work and not under copyright protection in the US; foreign copyright protection may apply 2021