Extensive Recombination-driven Coronavirus Diversification Expands the Pool of Potential Pandemic Pathogens

Abstract

The ongoing SARS-CoV-2 pandemic is the third zoonotic coronavirus identified in the last 20 years. Enzootic and epizootic coronaviruses of diverse lineages also pose a significant threat to livestock, as most recently observed for virulent strains of porcine epidemic diarrhea virus (PEDV) and swine acute diarrhea-associated coronavirus (SADS-CoV). Unique to RNA viruses, coronaviruses encode a proofreading exonuclease (ExoN) that lowers point mutation rates to increase the viability of large RNA virus genomes, which comes with the cost of limiting virus adaptation via point mutation. This limitation can be overcome by high rates of recombination that facilitate rapid increases in genetic diversification. To compare the dynamics of recombination between related sequences, we developed an open-source computational workflow (IDPlot) that bundles nucleotide identity, recombination, and phylogenetic analysis into a single pipeline. We analyzed recombination dynamics among three groups of coronaviruses with noteworthy impacts on human health and agriculture: SARSr-CoV, Betacoronavirus-1, and SADSr-CoV. We found that all three groups undergo recombination with highly diverged viruses from undersampled or unsampled lineages, including in typically highly conserved regions of the genome. In several cases, no parental origin of recombinant regions could be found in genetic databases, demonstrating our shallow characterization of coronavirus diversity and expanding the genetic pool that may contribute to future zoonotic events. Our results also illustrate the limitations of current sampling approaches for anticipating zoonotic threats to human and animal health.

Keywords: coronaviruses; evolution; virology.

Fig. 1.

AlphaCoV and BetaCoV phylogenetic relationships are genome region-dependent. (A) Basic coronavirus genome organization with the 5′ ∼20 kb comprising the replicase gene that is proteolytically processed into up to 16 individual proteins. The 3′ 10 kb comprises structural and genus-specific accessory genes. (B) Maximum-likelihood (ML) phylogenetic tree of alpha and betaCoVs full-length RNA-dependent RNA-polymerase encoding region of Orf1ab. (C) ML phylogenetic tree of full-length spike genes from viruses in the species Betacoronavirus 1 (red) rooted with the distantly related betacoronavirus mouse hepatitis virus. (D) ML phylogenetic tree of spike genes of SARSr-CoVs, with SARS-CoV-2-like viruses further analyzed in the paper highlighted in blue. (E) ML phylogenetic tree of spike genes from SADSr-CoVs (magenta) rooted with the distantly related alphacoronavirus RnCov/Lucheng-19.

Fig. 2.

IDPlot workflow. (A) Reference and query sequences are aligned using MAFFT. (B) Breakpoint detection is performed using GARD, capturing breakpoints across iterative refinements. (C) Phylogenetic trees based on breakpoints from each iteration and are created using FastTree 2. (D) Improvement in ΔAIC-c is plotted against the iteration. (E) Phylogenetic trees associated with the selected GARD iteration are displayed.

Fig. 3.

SARSr-CoV IDPlot analysis. (A) IDPlot analysis of SARS-CoV-2-like SARSr-CoVs with color-coded dashed lines defining divergent regions arising from recombination events with ancestral viruses. (B) ML tree of the RdRp-encoding region of SARS-2-like and other SARSr-CoVs showing close relationship between the SARS-CoV-2-like viruses. (C) ML tree of PangolinCoV/GD19 RR1 (which overlaps with BtCoV/RmYN02 RR1) showing different topology than the RdRp tree. (D) Schematic of spike proteins indicating divergent regions and nucleotide identity to the reference sequence and closest related sequence in GenBank. (E) ML tree of ORf8 showing that RmYN02 Orf8 is a divergent member of the SARS-CoV-like Orf8 branch.

Fig. 4.

Recombination analysis of Betacoronavirus-1. (A) Nucleotide identity plot and multiple sequence alignment of BetaCoV-1 viruses. Orange dashed lines indicate divergent regions of the ECoV-NC99 genome while black dashed lines are regions with high identity to the reference sequence bovine coronavirus (BCoV). (B) ML tree of nsp2-encoding region of Orf1ab, which falls within the divergent ECoV-NC99 Region 2. (C) ML tree of the RdRp-encoding region of Orf1ab. (D) Schematic depicting the spike gene diversity of BetaCoV1 demonstrating the divergence of ECoV-NC99 and PHEV. Top BLAST hits in bolded red indicate no GenBank entries with >80% nucleotide identity.

Fig. 5.

SADSr-CoV IDPlot analysis. (A) IDPlot nucleotide identity and multiple sequence alignment of eight SADSr-CoVs. Color-coded dashed lines indicate divergent regions in corresponding viruses owing to recombination events. (B) Schematic of spike genes of SADSr-CoVs along with nucleotide identity to the reference sequence and closest related sequences in GenBank for S1 and S2 domains. (C) Schematic of Orf7a diversity with nucleotide identity to the reference sequence and closest related sequences in GenBank. (D) Phylogenetic tree of SADSr-CoVs based on 3ClPro sequence illustrating the history of inferred recombination events indicated by arrowheads.