PLoS Pathogens
image
The coronavirus proofreading exoribonuclease mediates extensive viral recombination
DOI 10.1371/journal.ppat.1009226 , Volume: 17 , Issue: 1
Article Type: research-article, Article History
Abstract

Recombination is an essential part of normal coronavirus replication, required for the generation of the sub-genomic mRNAs as well as defective viral genome (DVGs) and is also implicated in novel strain emergence. However, the molecular mechanisms and determinants of RNA recombination in CoVs are unknown. Here, we compare recombination in 3 divergent beta-coronaviruses; murine hepatitis virus (MHV), MERS-CoV, and SARS-CoV-2. We show that they have striking similarities in the populations of RNA produced and in the sequences surrounding recombination junctions. Further, we demonstrate that the coronavirus proofreading exoribonuclease (nsp14-ExoN) is required to maintain the rates and loci of recombination generated during infection. These data suggest that recombination and the coronavirus exoribonuclease are conserved and important determinants of replication that may be targeted for inhibition and attenuation to control the ongoing pandemic of SARS-CoV-2 and prevent future outbreaks of novel coronaviruses.

Gribble, Stevens, Agostini, Anderson-Daniels, Chappell, Lu, Pruijssers, Routh, Denison, and Pekosz: The coronavirus proofreading exoribonuclease mediates extensive viral recombination

Introduction

The ongoing severe global pandemic of SARS-CoV-2, the etiological agent of coronavirus disease 2019 (COVID-19) underlines the importance of defining the determinants of coronavirus (CoV) evolution and emergence into human populations [1]. Studies comparing CoV strains that are closely related to SARS-CoV-2 have proposed that SARS-CoV-2 acquired the ability to infect human cells through recombination within the spike protein sequence [24]. Further, a study of genetic variation in patient SARS-CoV-2 samples has suggested that recombination may be occurring during infections in humans [5]. Recombination is also implicated in the emergence of severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV) [610]. Together, these data support the hypothesis that generation of novel CoVs, cross-species movement, and adaptation may be driven by recombination events in nature. CoV recombination has been reported to be associated with increased spread and severe disease, and has resulted in vaccine failure of multiple livestock CoVs [11,12]. Thus, targeting the ability of the virus to recombine is a critical consideration for vaccine development in the ongoing SARS-CoV-2 pandemic as well as future animal and zoonotic CoVs.

Coronaviruses are a family of positive-sense, single-stranded RNA viruses with genomes ranging in size between 26 and 32 kb (S1A Fig). During normal replication, the putative CoV replication-transcription complex (RTC), formed by multiple nonstructural proteins (nsp) encoded in ORF1ab, drives RNA synthesis and encompasses many enzymatic functions [1316]. Previous reports indicate that CoVs readily perform both inter-molecular recombination between 2 distinct molecules and intra-molecular recombination within the same molecule (S1B Fig). Co-infection with related strains of the model β-CoV murine hepatitis virus (MHV) results in chimeric viral genomes that are generated by inter-molecular recombination [17,18]. The CoV RTC performs intra-molecular recombination at virus-specific transcription regulatory sequences (TRSs) to generate a set of subgenomic mRNAs (sgmRNAs) with common 5’ and 3’ ends (S1A and S1B Fig) [19,20]. sgmRNAs are subsequently translated into structural and accessory proteins [19]. CoVs also generate defective viral genomes (DVGs) that contain multiple deletions of genomic sequence while retaining intact 5’ and 3’ genomic untranslated regions (5’ and 3’ UTRs). DVGs are amplified by RTC machinery supplied by co-infecting full-length helper CoVs [2124]. DVGs in respiratory viruses can act as pathogen-associated molecular patterns (PAMPs) and stimulate the innate immune system [25,26]. The role of DVGs in CoV biology is largely unknown, although some DVGs interfere with viral replication [27,28]. Therefore, CoVs perform recombination as a normal part of their replication, producing complex populations of recombined RNA molecules. Prior to the advent of Next Generation Sequencing (NGS), direct analysis of recombined CoV RNAs was not possible and the determinants of recombination could not be identified.

In other RNA virus families including picornaviruses and alphaviruses, regulation of recombination has been mapped to replication fidelity determinants in the viral RNA-dependent RNA polymerase (RdRp) [2932]. In contrast to these viruses, CoV replication fidelity is primarily determined by the 3’-to-5’ exoribonuclease encoded in nonstructural protein 14 (nsp14-ExoN) that proofreads RNA during replication through excision of mismatched incorporated nucleotides [3338]. Viral exonucleases are essential for recombination in DNA viruses, including vaccinia virus and herpes simplex virus 1 [39,40]. In contrast, a role of the nsp14-ExoN in CoV RNA recombination had not previously been defined. In our lab, viral mutants of MHV with engineered inactivation of nsp14-ExoN (ExoN(-)) resulted in reduced abundance of sgmRNA2. In another program, rescue of viable ExoN(-) human CoV 229E (HCoV-229E) was unsuccessful, but limited replication was associated with decreased detection of sgmRNAs [34,41]. Although these reports did not study recombination or molecular mechanisms, they support the hypothesis that CoV nsp14-ExoN activity RNA synthesis and possibly recombination, in addition to the known functions of nsp14-ExoN in CoV replication fidelity, viral fitness, in vivo virulence, resistance to nucleoside analogues, and immune antagonism [36,42,43].

In this study, we sought to define the frequency and patterns of recombination of divergent β-CoVs SARS-CoV-2, MERS-CoV, and MHV, and to test the role of nsp14-ExoN in recombination. We used both short-read Illumina RNA-sequencing (RNA-seq) and long-read direct RNA Nanopore sequencing for all three viruses to show that they perform extensive recombination during replication in vitro with broadly similar patterns of recombination, and generate diverse yet similar populations of recombined molecules. We further demonstrate that genetic inactivation of MHV nsp14-ExoN results in a significant decrease in recombination frequency, altered recombination junction patterns across the genome, and altered junction site selection. These defects and alterations result in a marked change in MHV-ExoN(-) recombined RNA populations, including defective viral genomes (DVGs). These results support future studies aimed at illuminating the role of SARS-CoV-2 nsp14-ExoN activity in RNA recombination, the regulation of sgmRNA expression, and its contribution to novel CoV zoonotic emergence. Combined with the multiple critical integrated functions of nsp14-ExoN, the role in recombination further defines nsp14-ExoN as a conserved, vulnerable, and highly specific target for inhibition by antiviral treatments and viral attenuation.

Results

SARS-CoV-2 and MERS-CoV generated extensive populations of recombination junctions

We first sought to quantify recombination frequency and identify recombination patterns in zoonotic CoVs by sequencing both MERS-CoV and SARS-CoV-2 RNA. In three independent experiments for each virus, Vero cell cultures were infected with either MERS-CoV or SARS-CoV-2 until the monolayer displayed >70% virus-induced cytopathic effect (CPE). Total RNA from infected cells was isolated and poly(A)-selected to capture all viral RNA containing poly-A tails, including genomic, subgenomic, and defective viral genome (DVG) RNA molecules. Equal amounts of total cell RNA from each of the three independent experiments for each virus was sequenced by short-read Illumina RNA-sequencing (RNA-seq), and by long-read direct RNA Nanopore sequencing. The depth and low error rate of RNA-seq facilitated the detection and quantification of both high- and low-abundance unique junctions. Long-read direct RNA sequencing on the Oxford Nanopore Technologies MinION platform was used to sequence complete RNA molecules, to define the organization of junctions in the context of intact RNA molecules. By comparing short- and long-read RNA sequencing, we accomplished high-confidence detection and quantification of recombination junctions as well as description of the genetic architectures of molecules formed by the junctions.

For RNA-seq, reads were aligned to the respective viral genomes (S1A Fig) using a recombination-aware mapper, ViReMa (Virus Recombination Ma pper) [44]. ViReMa detected recombination events that generated deletions greater than 5 base-pairs and that were flanked by a 25 base-pair alignment both upstream and downstream of the junction site. ViReMa-detected junctions may be formed from either inter-molecular or intra-molecular recombination during replication. ViReMa aligned both recombined and non-recombined reads in the library and reported the total number of nucleotides aligned to the genome and all detected recombination junctions.

Alignment of MERS-CoV and SARS-CoV-2 with ViReMa demonstrated nearly identical read coverages for MERS-CoV (1118) and SARS-CoV-2 (1122) (S2A and S2B Fig). Further, 82.95% of MERS-CoV RNA-seq reads and 77.48% of SARS-CoV-2 reads mapped to the viral genome, demonstrating RNA-seq libraries in both viruses had a similar proportion of viral RNA (S1 Table). To quantify recombination, recombination junction frequency (Jfreq ) was calculated for MERS-CoV and SARS-CoV-2 (Fig 1A). Jfreq refers to the number of nucleotides in all detected junctions normalized to viral RNA amount in a sample (total mapped nucleotides). Thus, Jfreq was not biased by the number of virus-mapping reads. Jfreq was multiplied by 104 to scale for library size and was reported as the number of junctions per 104 mapped nucleotides. MERS-CoV had a mean Jfreq of 37.80 junctions detected per 104 mapped nucleotides. SARS-CoV-2 had a mean Jfreq of 475.7 junctions per 104 mapped nucleotides (Fig 1A). This was a surprising difference in Jfreq between the two viruses that were infected at similar multiplicity of infections (MOIs), were collected when the cells displayed similar levels of CPE, and had similar viral abundance in sequenced RNA. We considered the possibility that the observed >10-fold difference between Jfreq of each virus could be due to the replication capacity of the parental virus. We compared the number of unique junctions generated by each virus to remove any potential viral replication bias. SARS-CoV-2 generated an average of 56,082 unique junctions per experiment, while MERS-CoV generated an average of 19,367 unique junctions per experiment (S2C Fig). Thus, both the number of recombination junctions and Jfreq were similarly higher in SARS-CoV-2 compared to MERS-CoV, suggesting that these differences are not solely due to an increased replication capacity or viral amplification of recombined species. This will be an important area for future study to determine if SARS-CoV-2 is associated with increased recombination in other cell types, in vivo models, or clinical samples. In any case, quantification of both recombination junction frequency and the number of unique recombination junctions in MERS-CoV and SARS-CoV-2 showed that both viruses produce abundant recombination junctions during replication in culture.

Genome-wide recombination generates populations of diverse RNA molecules in MERS-CoV and SARS-CoV-2.
Fig 1
MERS-CoV total cell lysates (black) and SARS-CoV-2 infected cell monolayers (violet) were sequenced by RNA-seq. (A) Junction frequency (Jfreq) was calculated by normalizing number of nucleotides in ViReMa-detected junctions to viral RNA (total mapped nucleotides) and multiplying by 10,000 to express Jfreq as the number of junctions per 104 mapped nucleotides. Error bars represent standard errors of the mean (SEM) for three independent sequencing libraries (N = 3). (B) Recombination junctions are mapped according to their genomic position (5’ junction site, Start Position; 3’ junction site, Stop Position) and colored according to their frequency in the population of all junctions in MERS-CoV and SARS-CoV-2. The highest frequency junctions are magenta and completely opaque. The lowest frequency junctions are red and the most transparent. Dashed boxes represent clusters of junctions: (i) 5’ ➔ 3’; (ii) mid-genome ➔ 3’ UTR; (iii) 3’ ➔ 3’; (iv) local deletions; (v) 5’ UTR ➔ rest of genome. (C) The Jfreq of DVGs, canonical sgmRNAs, and alternative sgmRNAs was calculated and compared in MERS-CoV (black) and SARS-CoV-2 (violet). Error bars represent SEM for 3 independent sequencing libraries (N = 3) of each virus. 2-way ANOVA with multiple comparisons corrected by statistical hypothesis testing (Sidak test). *** p < 0.001, **** p < 0.0001. Mean recombination frequency is quantified at each position across the MERS-CoV (D) and SARS-CoV-2 (E) genomes (N = 3). Recombination frequency was calculated by dividing the number of nucleotides in detected junctions at that position (start and stop sites) by the total number of mapped nucleotides at the position. See also S2 Fig and S1 Table. (F) The percent adenosine (A), cytosine (C), guanine (G), and uracil (U) at each position in a 30-base pair region flanking DVG junction start and stop sites in MERS-CoV (black) and SARS-CoV-2 (violet). Each point represents a mean (N = 3) and error bars represent SEM. The junction site is denoted as a carat (^) and with a solid red line. Positions upstream from the junction are labelled -30 to -1 and positions downstream are labelled +1 to +30. The expected nucleotide percentage based on the composition of the viral genome is marked as a dashed line (black = MERS-CoV, violet = SARS-CoV-2). (G) Distribution of sequence microhomology in MERS-CoV (black) and SARS-CoV-2 (violet) compared to an expected probability distribution (gray). The frequency of each nucleotide overlap length is displayed as a mean (N = 3) and error bars represent SEM.Genome-wide recombination generates populations of diverse RNA molecules in MERS-CoV and SARS-CoV-2.

To define the patterns of the detected recombination junctions, we mapped forward (5’ ➔ 3’) recombination junctions according to their genomic position (Figs 1B and S2C and S2D). Both MERS-CoV and SARS-CoV-2 displayed clusters of junctions in multiple conserved patterns: 1) between the 5’ and 3’ ends of the genome; 2) between intermediate genomic positions and the 3’ end of the genome; 3) within the 3’ end of the genome; 4) representing local deletions across the genome; and 5) between the 5’ untranslated region (UTR) and the rest of the genome. (Fig 1B). SARS-CoV-2 also had many low-frequency junctions distributed across the genome and horizontal clusters of low-frequency junctions between common start sites at position ~2000 and ~8000 and the rest of the genome (Fig 1B). Overall, these data demonstrate that extensive RNA recombination during replication of both MERS-CoV and SARS-CoV-2 generates diverse populations of junctions with similar high-abundance clusters.

MERS-CoV and SARS-CoV-2 recombination generated defective viral genomes and subgenomic mRNAs

We next sought to define and quantify the populations of recombined RNA molecules produced in both MERS-CoV and SARS-CoV-2. SARS-CoV-2 sgmRNAs were identified by the location of recombination junctions within previously defined 65 base-pair regions containing the transcription regulatory sequence (TRSs) of each sgmRNA [45]. Similarly, 65 base-pair windows were defined encompassing the MERS-CoV TRS core sequences for each sgmRNA. Junctions between the 5’ TRS-L and sgmRNA-specific TRS were filtered. The most abundant sgmRNAs were designated as “canonical”, and other sgmRNA species were designated “alternative sgmRNAs”. Recombination junctions outside of the TRS-L and the sgmRNA-specific TRSs were designated as DVG junctions.

For each virus, the frequencies of DVGs, canonical sgmRNAs, and alternative sgmRNAs were normalized to total virus RNA. For both MERS-CoV and SARS-CoV-2, canonical and alternative junctions were detected for all sgmRNAs (Figs 1C and S2E and S2F). MERS-CoV and SARS-CoV-2 alternative sgmRNA was detected at similar frequencies (Fig 1C). In contrast, SARS-CoV-2 generated significantly higher frequencies of DVGs and canonical sgmRNAs than MERS-CoV (Fig 1C).

We next calculated the mean recombination frequency at each genomic position by comparing the number of nucleotides in detected junctions (both start and stop sites) at that position, and normalized to nucleotide depth at that position. Further, we determined genomic positions with a mean recombination frequency greater than 50% (Fig 1D and 1E). In MERS-CoV, there were 5 positions >50%; 4 of these mapped to TRS positions and 1 position was located in ORF5 (Fig 1D). In SARS-CoV-2, there were 26 positions with >50% recombination frequency, with13 mapping to TRS positions. SARS-CoV-2 also had high recombination frequency at positions in the nsp2 coding sequence, the S gene, M gene, and N gene (Fig 1E). In summary, the genomic positions with the highest frequency for both MERS-CoV and SARS-CoV-2 mapped to TRSs that form sgmRNA leader-body junctions. However, positions with high recombination frequency were identified at other locations across the genomes and relatively more in SARS-CoV-2 than MERS-CoV.

MERS-CoV and SARS-CoV-2 defective viral genomes demonstrated distinct nucleotide compositions in the sequences flanking junctions

For both SARS-CoV-2 and MERS-CoV, the nucleotide composition of the start and stop sequences resulting in junctions forming DVGs in MERS-CoV and SARS-CoV-2 was determined and compared to the expected nucleotide percentage based on the parental viral genomes (Fig 1F). Sequences upstream (-30 to -1) and downstream (+1 to +30) of both the genomic start and stop sites of DVG junctions were analyzed. DVGs formed by junctions would contain sequences upstream of the start site (-30 to -1) and downstream of the stop site (+1 to +30) (S1C Fig). For both MERS-CoV and SARS-CoV-2, start and stop sequences upstream of the junction were enriched for uracil (U) and depleted for adenosine (A) and guanine (G). Downstream of the junction in both start and stop sites, both viruses were enriched for guanine (G) and adenosine (A) and depleted for uracil (U). MERS-CoV demonstrated a preference for U(U/C)^(G/A/C)(A/C)C in DVG start sites and UU^(G/C/A)C(G/C) in DVG stop sites. SARS-CoV-2 DVG sequences favored AUUU^(G/A)AAA in the start site sequences and ACUU^G(C/A)(C/A) in the stop site sequences. The nucleotide composition of MERS-CoV and SARS-CoV-2 differ from TRS-like sequences of MERS-CoV (AACGAA) [46] and SARS-CoV-2 (ACGAAC) [47], and therefore represent a selection of separate sequences for DVG formation.

MERS-CoV and SARS-CoV-2 exhibited sequence microhomology at recombination junctions

We next tested whether MERS-CoV and SARS-CoV-2 junction sites favored regions of sequence microhomology at recombination junctions, defined as 2–20 nt regions of identical overlap [48]. The distribution of frequencies of 0–10 overlapping nucleotides at the start and stop sites of detected recombination junctions in both MERS-CoV and SARS-CoV-2 were compared to an expected probability distribution. Both MERS-CoV and SARS-CoV-2 junction sites demonstrated increased frequencies of overlaps of 2–7 nt (Fig 1G). Thus, MERS-CoV and SARS-CoV-2 favor the formation of recombined RNAs at junction sites exhibiting sequence microhomology.

Direct RNA Nanopore sequencing of MERS-CoV and SARS-CoV-2 defined the architecture of full-length genome, sgmRNAs, and DVGs

We performed direct RNA Nanopore sequencing on the same RNA used for short-read RNA-seq. We analyzed three independent experiments for each virus and sequenced 178,658 MERS-CoV RNA molecules and 1,725,862 SARS-CoV-2 RNA molecules that had 85.6% and 82.2% identity to the parental genome, respectively (S2 Table). To remove prematurely truncated sequences, we computationally selected only Nanopore reads containing both genomic termini. We obtained 3 full-length direct RNA sequences of the SARS-CoV-2 genome containing over 29,850 consecutive nucleotides that aligned to the SARS-CoV-2 genome (S3 Table). In MERS-CoV RNA, we detected 451 full-length molecules containing genomic termini and 473 unique junctions (Fig 2A and S2 and S4 Tables). SARS-CoV-2 RNA generated 172,191 complete molecules and 181,770 unique junctions (Fig 2B and S2 and S4 Tables). To confirm junctions in detected by direct RNA sequencing, we compared unique junctions detected in filtered complete RNA molecules with 20 bp windows at both the start and stop sites to unique junctions detected in short-read Illumina RNA-seq datasets reported in Figs 1 and S2. 89.29% of MERS-CoV and 97.97% of SARS-CoV-2 Nanopore junctions were also detected in RNA-seq datasets S2 Table).

Direct RNA Nanopore sequencing of MERS-CoV and SARS-CoV-2 reveals accumulation of distinct recombined RNA populations.
Fig 2
Direct RNA Nanopore sequencing of poly-adenylated MERS-CoV and SARS-CoV-2 RNA. Three sequencing experiments were performed for each virus. Nanopore reads passing quality control were combined and mapped to the viral genome using minimap2 [70]. Genome coverage maps and Sashimi plots visualizing junctions (arcs) in full-length (A) MERS-CoV (black) and (B) SARS-CoV-2 (violet) RNA reads. (C) Distinct RNA molecules identified in MERS-CoV (black) with at least 3 supporting reads are visualized. The number of sequenced reads containing the junction is listed (Count). Genetic sequences of each RNA molecule are represented by filled boxes and deleted regions are noted (Deleted Region(s)) and represented by dashed lines. (D) The 15 most abundant SARS-CoV-2 (violet) recombined RNA molecules and 3 full-genome reads are visualized. See also S2 Table, S3 Table, S4 Table.Direct RNA Nanopore sequencing of MERS-CoV and SARS-CoV-2 reveals accumulation of distinct recombined RNA populations.

To define the architectures of detected molecules, we filtered for junctions with at least 3 supporting Nanopore reads. For both viruses, junctions were categorized as either a DVG or sgmRNA junction using the same criteria as with the RNA-seq data. In MERS-CoV, we defined 5 distinct species, including 3 sgmRNAs (6, 7, and 8) and 2 DVGs (Fig 2C). In SARS-CoV-2, there were 1166 species with a single junction and 227 containing 2 junctions. The 15 most abundant species in SARS-CoV-2 included 11 predicted sgmRNA transcripts and 4 DVGs (Fig 2D). We also identified potential alternative transcripts corresponding to the ORF6, ORF7a, ORF8, and the M genes (Fig 2D). In summary, direct RNA Nanopore sequencing defined a diverse set of recombined RNAs generated by both MERS-CoV and SARS-CoV-2 with most DVGs containing only a singular recombination event rather than extensive genomic rearrangement. Thus, both MERS-CoV and SARS-CoV-2 engaged in extensive RNA recombination during replication, producing diverse junctions across the viral genomes and many recombined RNA species.

Genetic inactivation of the MHV nsp14-exoribonuclease (ExoN) resulted in significantly decreased and altered RNA recombination

We previously have reported that the nsp14 exoribonuclease (nsp14-ExoN) activity is required for high-fidelity replication and proofreading for the β-CoVs murine hepatitis virus (MHV) and SARS-CoV [3336]. We sought to determine whether nsp14-ExoN activity also contributed to the extensive recombination observed in coronaviruses. Since no proofreading-deficient nsp14-ExoN catalytic mutant is available for MERS-CoV or SARS-CoV-2, we used the MHV nsp14-ExoN inactivation mutant (MHV-ExoN(-)) and wild-type virus (MHV-WT) to compare recombination [49]. Murine DBT cells were infected with MHV-WT or MHV-ExoN(-) in three independent experiments, and RNA was isolated from infected cell monolayers and viral supernatant when the cell monolayer was intact and 90% cytopathic effect (CPE) was observed. Poly(A)-selected RNA-seq libraries were aligned to the MHV genome using ViReMa (AY910861.1). In both infected cell monolayers and viral supernatants, MHV-WT and MHV-ExoN(-) had similar mean coverages ranging between 1100 and 1700 reads (S4A and S4B Fig).

Previous studies have shown that MHV-ExoN(-) has decreased genome replication compared to WT [34]. We accounted for decreased MHV-ExoN(-) viral RNA by normalizing the number of nucleotides participating in detected junctions to the amount of viral RNA (total mapped nucleotides), and Jfreq was calculated as described for Fig 1A. MHV-ExoN(-) had significantly decreased Jfreq relative to MHV-WT in both infected cells and viral supernatant (Fig 3A and 3C). To address any potential viral replication bias resulting from the differences between MHV-WT and MHV-ExoN(-) replication that have been previously reported, we quantified and compared the unique detected recombination junctions. In both infected cell monolayers and in viral supernatant, MHV-ExoN(-) had significantly decreased unique recombination junctions compared to MHV-WT (S3C and S4C Figs). Thus, MHV-ExoN(-) had decreased recombination junction frequency and number of unique junctions compared to MHV-WT, showing that loss of nsp14-ExoN activity resulted in significantly less recombination during infection.

Loss of nsp14-ExoN activity decreases recombination frequency and alters recombination junction patterns across the genome.
Fig 3
Infected monolayer and viral supernatant RNA poly(A) selected, sequenced by RNA-seq, and aligned to the MHV genome using ViReMa. Junction frequency (Jfreq) in infected monolayer RNA (A) and viral supernatant RNA (C) was calculated by normalizing the number of nucleotides in ViReMa-detected junctions to total viral RNA (total mapped nucleotides) and multiplying by 10,000, expressing Jfreq as number of junctions per 104 mapped nucleotides. Error bars represent standard error of the means (SEM) (N = 3). Statistical significance was determined by the unpaired student’s t-test. * p < 0.05, **** p < 0.0001. Unique forward (5’ ➔ 3’) recombination junctions detected in infected monolayers (C) and viral supernatant (E) were mapped in MHV-WT and MHV-ExoN(-) according to their genomic position. Junctions are colored according to their frequency in the population (high frequency = magenta; low frequency = red). Clusters are marked by dashed boxes: (i) 5’ ➔ 3’; (ii) mid-genome ➔ 3’; (iii) 3’ ➔ 3’; (iv) local deletions; (v) 5’ UTR ➔ rest of genome. See also S3 and S4 Figs.Loss of nsp14-ExoN activity decreases recombination frequency and alters recombination junction patterns across the genome.

Recombination junctions were plotted according to their start (5’) and stop (3’) sites in infected cells and viral supernatant (Figs 3B, 3D, S3C, S3D, S4C and S4D). MHV-WT displayed clusters of junctions that were similar to those demonstrated in MERS-CoV and SARS-CoV-2, specifically: 1) between the 5’ and 3’ ends of the genome; 2) between intermediate genomic positions and the 3’ end of the genome; 3) between the 5’ UTR and the rest of the genome; 4) in local deletions across the genome; and 5) within the 3’ end of the genome (Fig 3B and 3D). While both WT and MHV-ExoN(-) accumulated junction clusters between the 5’ and 3’ ends of the genome and within the 3’ end of the genome, MHV-ExoN(-) had fewer junctions between the 5’ UTR and the rest of the genome and fewer junctions forming local deletions (Fig 3B and 3D). Thus, loss of MHV nsp14-ExoN activity resulted in decreased recombination frequency and altered junction patterns across the genome.

MHV-ExoN(-) had altered recombination at distinct positions across the genome

We next calculated and compared mean recombination frequency at each genomic position in MHV-WT and MHV-ExoN(-) (Fig 4A–4B). Both MHV-WT and MHV-ExoN(-) had high recombination frequency at the 5’ and 3’ ends of the genome as well as at distinct sites across the genome. Positions with >50% recombination frequency were localized to the TRS regions (Fig 4A and 4B). MHV-ExoN(-) had significantly altered recombination frequency at 765 positions in infected cell RNA and 499 positions in viral supernatant RNA (Figs 4A and 4B and S5). These positions were distributed across the genome, including the 5’ TRS-Leader, non-structural protein coding sequences, TRSs, structural and accessory ORFs, and 3’ UTR (S5A–S5E Fig). Thus, genetic inactivation of nsp14-ExoN altered recombination frequency at multiple positions across the genome.

Loss of nsp14-ExoN alters recombination at multiple genomic loci and skews recombined RNA populations.
Fig 4
Mean recombination frequency at each position across the MHV genome was compared in MHV-WT (blue) and MHV-ExoN(-) (orange) infected monolayer (A) and viral supernatant RNA (B). 2-way ANOVA with multiple comparisons (N = 3). The junction frequencies (Jfreq) of DVGs, canonical sgmRNAs, and alternative sgmRNAs were compared in MHV-WT (blue) and MHV-ExoN(-) (orange) infected monolayers (C) and viral supernatant (D). Error bars represent standard errors of the mean (SEM) (N = 3) and statistical significance was determined by a 2-way ANOVA with multiple comparisons correct by statistical hypothesis testing (Sidak test), ** p <0.01, *** p < 0.001, **** p < 0.0001. The Jfreq of canonical sgmRNA junctions was compared in MHV-WT (blue) and MHV-ExoN(-) (orange) infected monolayers (E) and viral supernatant (F). Error bars represent SEM (N = 3). Statistical significance was determined by a 2-way ANOVA with multiple comparisons corrected by statistical hypothesis testing (Sidak test), *** p < 0.001, **** p < 0.0001. The Jfreq of alternative sgmRNA molecules was quantified for MHV-WT (blue) and MHV-ExoN(-) (orange) infected cell monolayers (G) and viral supernatant (H). Error bars represent SEM (N = 3). Statistical significance was determined by a 2-way ANOVA with multiple comparisons corrected by statistical hypothesis testing (Sidak test), * p < 0.05, **** p < 0.0001. The abundance of junctions in MHV-ExoN(-) was compared to MHV-WT in infected monolayers (I) and viral supernatant (J) by DESeq2 . Junctions with statistically significant altered abundance (p < 0.05, N = 3) in MHV-ExoN(-) are mapped across the genome and colored according to their fold-change (red squares = decreased abundance, green circles = increased abundance). See also S3S5 Figs and S5 and S6 Tables.Loss of nsp14-ExoN alters recombination at multiple genomic loci and skews recombined RNA populations.

MHV-ExoN(-) had decreased abundance and altered ratios of DVGs and sgmRNAs

Compared with WT, MHV-ExoN(-) had significantly decreased frequencies of DVGs and both canonical and alternative sgmRNAs (Fig 4C). MHV-ExoN(-) viral supernatant also demonstrated a significant decrease in canonical sgmRNAs (Fig 4D). In addition to frequencies of DVGs and sgmRNAs in MHV-ExoN(-), the ratios of DVGs and both canonical and alternative sgmRNAs were skewed. Compared to WT, MHV-ExoN(-) had a significantly increased proportion of DVGs and significantly decreased proportions of both canonical and alternative sgmRNAs (S3E and S4E Figs). MHV-ExoN(-) also displayed significantly skewed proportions of individual canonical and alternative sgmRNA species (S3F and S3G Fig and S4F and S4G Fig). Decreased frequencies and aberrant proportions of DVGs and both canonical and alternative sgmRNAs show that nsp14-ExoN activity is a key determinant in recombination producing distinct RNA populations.

MHV-ExoN(-) had altered junction site selection

We next identified junctions with altered abundances in MHV-ExoN(-) compared to MHV-WT using DESeq2 [50]. MHV-ExoN(-) generated recombination junctions with significantly increased or decreased abundance relative to MHV-WT (S5F and S5G Fig and S6 Table). Clusters of junctions with either increased or decreased abundance in MHV-ExoN(-) compared to WT were localized to distinct genomic regions. Recombination junctions significantly enriched in MHV-ExoN(-) were mainly found between the 5’ and 3’ ends of the genome (Fig 4I and 4J). Junctions with significantly decreased abundance in MHV-ExoN(-) clustered between the 5’ UTR and the rest of the genome and local deletions of 10–50 bp in length across the genome (Fig 4I and 4J). Thus, the populations of recombination junctions that were differentially abundant in MHV-ExoN(-) were not randomly distributed across the genome, suggesting specific changes to junction site selection.

MHV-ExoN(-) DVG junction-flanking sequences demonstrated altered nucleotide composition while retaining microhomology at junction sites

To test whether MHV-ExoN(-) has altered sequence composition at its recombination junctions, we filtered DVG junctions and quantified nucleotide composition of adenosine (A), cytosine (C), guanine (G), and uracil (U) in the start and stop sequences flanking junction sites. Both MHV-WT and MHV-ExoN(-) demonstrated similar patterns of depletion and enrichment of nucleotides in infected cell monolayers and viral supernatant (Figs 5A and S6A). Start site sequences favored sequences of UUU(U/A)(U/A)^GG and were depleted for C upstream of the junction. Stop site sequences were relatively enriched for the sequence AAA(U/A)(U/A)^AA(G/A). These patterns and sequence preferences were similar to the sequence composition of both MERS-CoV and SARS-CoV-2 DVG recombination junctions (Fig 1F). In all three viruses, a preference for UUG spanning junction start sites was defined. Further, the DVG junction sequence preference differed from sequence composition of TRS-like sequences for MHV (AAUCUAUAC) [51] and represented a different selection of sequences for DVG formation. Loss of nsp14-ExoN(-) activity resulted in significantly altered nucleotide composition at multiple positions for all nucleotides in both the start and stop sites (Figs 5A and S6A). For both MHV-WT and MHV-ExoN(-), junction sites encoded more and longer microhomology overlaps of up to 8bp than would be expected by chance (Figs 5B and S6B).Thus, while loss of nsp14-ExoN activity altered nucleotide composition at multiple positions surrounding DVG junction sites, the overall patterns of enrichment and depletion were maintained and microhomology at the junction sites remained unchanged.

MHV-ExoN(-) DVG junction sites display both WT-like patterns of sequence composition and multiple alterations in nucleotide frequency, revealing microhomology at junctions.
Fig 5
(A) Nucleotide composition was calculated as the percent adenosine (A), cytosine (C), guanine (G), and uracil (U) at each position in a 30-base pair region flanking DVG junction start and stop sites in MHV-WT (blue) and MHV-ExoN(-) (orange) infected monolayer RNA. The junction is labelled as a carat (^) and a solid red line with upstream positions numbered -30 to -1 and downstream positions +1 to +30. The expected nucleotide percentage was calculated based on the overall MHV genome and represented as a dashed black line. Each point represents a mean (N = 3) and error bars represent SEM. 2-way ANOVA with multiple comparisons corrected for false discovery rate (FDR) by the Benjamini-Hochberg method. * q < 0.05, ** q < 0.01, *** q < 0.001, **** q < 0.0001. (B) Distribution of microhomology overlaps in MHV-WT (blue) and MHV-ExoN(-) (orange) compared to an expected probability distribution (gray). The frequency of each overlap length is displayed as a mean (N = 3) and error bars represent SEM. See also S5 Fig.MHV-ExoN(-) DVG junction sites display both WT-like patterns of sequence composition and multiple alterations in nucleotide frequency, revealing microhomology at junctions.

Direct RNA Nanopore sequencing identified changes in MHV-ExoN(-) full-length recombined RNA populations

To test the alterations of recombined RNAs due to loss of nsp14-ExoN proofreading activity, we sequenced MHV-WT and MHV-ExoN(-) viral supernatant RNA by direct RNA Nanopore sequencing. When reads were mapped to the MHV genome using minimap2 , MHV-WT datasets contained 102,367 viral molecules and MHV-ExoN(-) contained 19,445 (Fig 6A and S2 Table). We validated MHV-WT and MHV-ExoN(-) Nanopore junctions by comparing to RNA-seq datasets. 96.00% of MHV-WT and 97.50% of MHV-ExoN(-) Nanopore junctions were also detected in RNA-seq datasets (S2 Table).

Direct RNA Nanopore sequencing of MHV full-length recombined RNA molecules.
Fig 6
Direct RNA Nanopore sequencing of MHV viral supernatant RNA. (A) Genome coverage maps of full-length MHV-WT (blue) and MHV-ExoN(-) (orange) Nanopore reads aligned to the MHV-A59 genome using minimap2 . (B) Sashimi plot visualizing junctions (arcs) in MHV-WT (blue) and MHV-ExoN(-) (orange). (C) RNA molecule genetic architectures with at least 3 supporting reads identified in both MHV-WT and MHV-ExoN(-) (yellow) and unique to MHV-WT (blue). Genetic sequences of the RNA molecule are represented by filled boxes. Deleted regions are reported (Deleted Region) and represented by dashed lined. The number of reads supporting each species are noted (Count). See also S2 Table, S3 Table, and S4 Table.Direct RNA Nanopore sequencing of MHV full-length recombined RNA molecules.

MHV-ExoN(-) had a global decrease in the number of junctions across the genome (Fig 6B and S2 and S4 Tables). We filtered MHV-WT and MHV-ExoN(-) datasets for RNA molecules containing both 5’ and 3’ genomic ends that were supported by at least three reads. Nine such architectures were identified in MHV-WT (Fig 6C). These populations contained both DVGs and sgmRNAs. The four most abundant species were also detected in MHV-ExoN(-) viral supernatant RNA, which corresponded to a DVG and sgmRNAs 4,6 and 7 (Fig 6C). We did not detect unique MHV-ExoN(-) variants with at least 3 supporting reads, potentially due to their low frequency in the population. These data demonstrate that loss of nsp14-ExoN activity drives the accumulation altered recombined RNA populations and skewed DVG species diversity.

Discussion

While CoV recombination has long been proposed as a driver of novel strain emergence and is known to be a constitutive aspect of CoV replication, the diversity of recombination products and sequence and protein determinants had not previously been defined. In this study, we show the diversity of the CoV recombination landscape in the β-coronaviruses SARS-CoV-2, MERS-CoV, and murine hepatitis virus (MHV), and we demonstrate that loss of the nsp14 exoribonuclease activity in MHV results in decreased recombination and altered site selection of recombination junctions. Our results support a model in which nsp14-ExoN activity is required for normal recombination. Thus, nsp14-ExoN is a key component of CoV recombination, adding another essential function to the repertoire of those already reported for nsp14-ExoN, specifically CoV high-fidelity replication, RNA synthesis, resistance to antiviral nucleoside analogues, fitness, immune antagonism, and virulence.

Divergent β-CoVs generate extensive and similar recombination networks yielding diverse populations of RNA species

We show that MHV, MERS-CoV, and SARS-CoV-2 perform extensive recombination and generate diverse populations of RNA molecules, demonstrated by independent short-read Illumina RNA-seq and long-read, direct RNA Nanopore sequencing. These divergent group 2a (MHV), 2b (SARS-CoV-2), and 2c (MERS-CoV) β-CoVs demonstrated many strong similarities in their patterns of recombination junctions across the genomes and in the types of recombined RNAs produced. Specifically, the similarities across all three viruses in the nucleotide composition of sequences flanking DVG junctions and the common increased junction sequence microhomology support the conclusion that recombination mechanisms have been conserved across different evolutionary trajectories and host species specificity.

There also were distinct recombination patterns for each virus that were confirmed across independent experiments and by agreement between RNA-seq and Nanopore datasets for all viruses. These differences most likely represent evolutionary divergence of recombination in distinct viruses or sub-genera represented by MHV, SARS-CoV-2 and MERS-CoV. However, it remains possible that observed differences could be impacted by the diversity of the original sample or replication in different cell types. SARS-CoV-2 stock virus was a low passage (P5) population from a clinical isolate that had been passaged in Vero cells, while MERS-CoV and MHV were low passage stocks generated from isogenic cDNA clones. It will be important for future studies to determine the role of the diversity of the viral population, cell environment, virus-specific RNA synthesis kinetics, and virus adaptation/evolution in viral recombination. The extent of the pandemic and availability of genetically diverse viruses will allow investigators to test whether patterns of SARS-CoV-2 recombination show alterations between early and later pandemic isolates, and if any identified differences correlate with or predict changes in other replication or pathogenesis.

Sequences containing microhomology are likely determinants of recombination resulting in CoV defective viral genome formation

High-resolution analysis of DVG junctions produced during replication by MERS-CoV, SARS-CoV-2, and MHV reveals that a significant preference for a UUG motif, suggesting a possible conserved core sequence for DVG synthesis that differs from sgmRNA transcriptional regulatory sequences. These results support a model across multiple divergent β-CoVs in which DVGs result from recombination junction selection by the RTC based on both broadly similar sequence identity and specific sequence microhomology of 2–10 bp (S1D Fig). This model would be most similar to microhomology-mediated end-joining, a mechanism of genomic repair in eukaryotic DNA systems that results in large sequence deletions [52,53]. The presence of sequence homology-driven recombination and DVG formation suggests an selection for specific DVG biogenesis, supporting the hypothesis that DVGs play specific roles in coronavirus replication, pathogenesis and evolution. The results of this study will form the basis for direct genetic studies of DVGs as well as ability to target templates for study of the viral replicase functions.

MHV nsp14-ExoN determines the extent, diversity, and junction site selection of RNA recombination during infection

MHV-ExoN(-) mutants showed decreased recombination junction frequency and altered populations of sgmRNAs and DVGs, demonstrating a previously unknown role for nsp14-ExoN in CoV RNA recombination. There is no precedent in RNA viruses for the regulation of recombination by a virus encoded exoribonuclease. In contrast, in DNA viruses such as poxviruses and herpesviruses, virus-encoded exonuclease activity stimulates recombination by single-strand annealing through both exonuclease degradation of nucleic acids and interactions with other proteins [39,40]. In the single-stranded, positive-sense RNA virus families picornaviridae and alphaviridae that lack any exonuclease, low-fidelity mutant viruses have altered polymerase speed and processivity [54] and these properties contribute to recombination and the generation of DVGs [32,55,56]. Our results suggest that CoVs have evolved to regulate both proofreading and recombination by the nsp14-ExoN protein. Mutation of the active site of nsp14-ExoN alters both these functions, supporting a complex interaction with other proteins in the CoV RTC, including the nsp12 RNA-dependent RNA polymerase. In the low-fidelity picornavirus and alphavirus mutants, it has been proposed that impaired fidelity alters polymerase processivity and speed, resulting in decreased recombination. It is possible that CoV nsp14-ExoN mutations may similarly impair polymerase speed and processivity, resulting in altered patterns of DVGs and non-canonical sgmRNAs. The direct role of polymerase speed and processivity and the potential mechanisms by which these principles influence recombination remains to be determined, but possibilities include altered RTC stability through the changes to the complex protein-protein interactions or RTC-RNA interactions.

ExoN is a powerful tool for understanding CoV replication, and a novel and conserved target for inhibition and attenuation

The similarities between the patterns of recombination across divergent WT β-CoVs, along with the differences observed between recombination in MHV WT and ExoN(-) viruses, support the hypothesis that ExoN mutants will inform our understanding of the evolution of the unique CoV multi-protein polymerase complex. Specifically, the model of DVG synthesis defined in MHV, MERS-CoV, and SARS-CoV-2 will allow for the direct testing of the roles of DVGs in CoV replication. Further, the role of ExoN in CoV recombination, along with the previously defined roles of ExoN in RNA proofreading during replication, native resistance to nucleoside analogues, immune evasion, and virulence and pathogenesis, highlight nsp14-ExoN as conserved and vulnerable target for both antiviral inhibitors and virus attenuation. ExoN(-) viruses are profoundly more sensitive to a range of antiviral nucleoside analogues, including remdesivir, ribavirin, 5-fluorouracil, and β-d-N4 -hydroxycytidine (NHC, EIDD 1931/2801) [33,38,57]. Nucleoside analogues and exonuclease inhibitors that target nsp14-ExoN can be tested for an additional impact on recombination and illuminate antiviral mechanisms of action. Finally, recombination has driven the vaccine escape in multiple CoVs [11,12]. The finding that MHV-ExoN(-) has decreased recombination during viral replication may have important implications for any design of live-attenuated SARS-CoV-2 or other animal or zoonotic CoVs. Our previous studies have shown that the ExoN(-) substitutions in MHV and SARS-CoV are evolutionarily stable over long-term passage in culture and in mice, and that a SARS-CoV ExoN(-) mutant is attenuated in mice while producing a robust and protective immune response against WT SARS-CoV infection [38,42,58,59]. The results in this paper raise the intriguing possibility that any CoV encoding ExoN(-) would have less recombination potential for repair or escape.

Materials and methods

Cell lines

DBT-9 (delayed brain tumor, murine astrocytoma clone 9) cells were maintained at 37°C as described previously [60]. DBT-9 cells were originally obtained from Ralph Baric at University of North Carolina-Chapel Hill and were maintained within 50 passages of this progenitor stock. Cells were maintained in Dulbecco’s modified Eagle medium (DMEM) (Gibco) supplemented with 10% fetal clone serum (FCS) (Invitrogen), 100 U/mL penicillin and streptomycin (Gibco), and 0.25 μg/mL amphotericin B (Corning). Cercopithecus aethiops Vero CCL-81 cells maintained in Dulbecco’s modified Eagle medium (DMEM) (Gibco) supplemented to final concentrations of 10% fetal calf serum (Gibco), 100 IU/ml penicillin (Mediatech), 100 mg/ml streptomycin (Mediatech), and 0.25 mg/ml amphotericin B (Mediatech) were used for MERS-CoV-2 infection. Vero CCL-81 cells were obtained from ATCC. Vero E6 cells maintained in Dulbecco’s modified Eagle medium (DMEM) (Gibco) supplemented to final concentrations of 10% fetal calf serum (Gibco), 100 IU/ml penicillin (Mediatech), 100 mg/ml streptomycin (Mediatech), and 0.25 mg/ml amphotericin B (Mediatech) were used for SARS-CoV-2 infections. Vero E6 cells were obtained from ATCC.

Viruses

All MHV work was performed using the recombinant WT strain MHV-A59 (GenBank accession number AY910861.1 [61]) at passage 4 and an engineered ExoN(-) strain of MHV-A59 at passage 2. The recovery of MHV-ExoN(-) were previously described include the four-nucleotide substitution of motif I residues resulting in alanine substitution (DE ➔ AA) [34] Experiments involving MERS-CoV were conducted using the human EMC/2012 strain recovered from an infectious clone (GenBank accession number JX869059.2) [62]. Experiments involving SARS-CoV-2 were conducted with a passage 5 virus inoculum generated from a Seattle, WA, USA COVID-19 patient (GenBank accession number MT020881.1). All virus manipulations were performed under stringent BSL-3 laboratory conditions according to strict protocols designed for safe and controlled handling of MERS-CoV and SARS-CoV-2.

MHV isolation and viral supernatant purification

Subconfluent 150-cm2 flasks were infected with either MHV-A59 or MHV-ExoN(-) at an MOI of 0.01 PFU/cell. Supernatant was harvested at either 16 hours post infection (MHV-A59) or 24 hours post infection (MHV-ExoN(-)) when the monolayer was >95% fused and remained intact. Infection supernatant was clarified by centrifugation at 1500 x g for 5 minutes at 4°C. Viral supernatant was purified on a 30% sucrose cushion by ultracentrifugation at 25,000 RPM at 4°C for 16 hours. The viral pellet was resuspended in MSE buffer (10mM MOPS, pH 6.8; 150mM NaCl; 1 mM EDTA). Viral RNA was extracted using the TRIzol-LS reagent according to manufacturer’s protocols. RNA was quantified using the Qubit RNA HS assay. Supernatant data in this paper is the result of three experiments sequenced independently from the infected cell monolayer samples.

MHV isolation from infected monolayers

In three independent experiments, a subconfluent 150-cm2 flask of DBT-9 cells was infected with either MHV-WT or MHV-ExoN(-) at an MOI or 0.01 PFU/cell. Monolayer was harvested at either 16 hpi (MHV-WT) or 24 hpi (MHV-ExoN(-)) when the monolayer was >95% fused and >75% intact. RNA was extracted with TRIzol according to manufacturer’s protocols. Infected monolayer data in this paper is the result of three independent experiments sequenced independently.

MERS-CoV infection

In three independent experiments, a nearly confluent 25-cm2 flask of Vero CCL-81 cells was infected with MERS-CoV at an MOI of 0.3 pfu/cell. Total infected cell lysates were collected at 72 hpi with the monolayer >70% fused. RNA was extracted in TRIzol according to manufacturer’s protocols.

SARS-CoV-2 infection

In three independent experiments, a total of 5 subconfluent 25-cm2 flasks of Vero E6 cells were infected at an MOI = 0.45 pfu/cell and cellular monolayers were harvested 60 hpi when the monolayer was >90% fused. RNA was extracted in TRIzol according to manufacturer’s protocols.

Short-read Illumina RNA-sequencing of viral RNA

Next generation sequencing (NGS) libraries were generated using 2 μg of RNA of each sample. RNA was submitted to Genewiz for library preparation and sequencing. Briefly, after quality control, polyadenylated RNA was selected during library preparation. Isolated RNA was heat fragmented, RT-PCR amplified with equivalent number of cycles, size-selected, and libraries were prepared for 2 x 150 nucleotide paired-end sequencing performed (Illumina). Genewiz performed basecalling and read demultiplexing.

Direct RNA Nanopore sequencing

RNA from ultracentrifuge-purified viral supernatant was prepared for direct RNA Nanopore sequencing on the Oxford Nanopore Technologies MinION platform according to the manufacturer’s protocols. Libraries were sequenced on fresh MinION R9.4 flow-cells for 24 hours, or until the pore occupancy was under 20%. Viral supernatant RNA from three independent experiments was sequenced on three separate flow cells for both MHV-WT and MHV-ExoN(-). MERS-CoV RNA from three independent cultures was sequenced on three separate flow cells. SARS-CoV-2 RNA isolated from three independent infections was sequenced on three separate flow cells.

Illumina RNA-seq processing and alignment

Raw reads were processed by first removing the Illumina TruSeq adapter using Trimmomatic [63] default settings (command line parameters java -jar trimmomatic.jar PE sample_R1.fastq.gz sample_R2.fastq.gz output_paired_R1.fastq output_unpaired_R1.fastq output_paired_R2.fastq output_unpaired_R2_unpaired.fastq ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGIWINDOW:4:15 MINLEN:36). Reads shorter than 36 bp were removed and low-quality bases (Q score < 30) were trimmed from read ends. The raw FASTQ files were aligned to the MHV-A59 genome (AY910861.1), the MERS-CoV genome (JX869059.2), and the SARS-CoV-2 genome (MT020881.1) using the Python2 script ViReMa (Viral Recombination Ma pper, version 0.15) [44] using the command line parameters python2 ViReMa.py reference_index input.fastq output.sam—OuputDir sample_virema/—OutputTag sample_virema -BED—MicroIndelLength 5. The sequence alignment map (SAM) file was processed using the samtools [64] suite to calculate nucleotide depth at each position in a sorted binary alignment map (BAM) file (using command line parameters samtools depth -a -m 0 sample_virema.sorted.bam > sample_virema.coverage).

Recombination junction analysis

Recombination junction frequency (Jfreq) was calculated by comparing the number of nucleotides in detected recombination junctions to the total number of mapped nucleotides in a library. Nucleotides in detected recombination junctions were quantified as a sum of nucleotide depth reported at each junction in the BED file generated by ViReMa. Total nucleotides mapped to the MHV-A59 genome were quantified as a sum of nucleotide depth at each position across the genome in the tab-delineated text file generated by the samtools. Jfreq was reported as junctions per 104 nucleotides sequenced. Mean Jfreq values were compared between MHV-WT and MHV-ExoN(-) and statistical significance determined by an unpaired student’s t-test. Junctions were mapped across the genome according to their start (5’) and stop (3’) positions. These junctions were first filtered in the forward (5’ ➔ 3’) direction using the dpylr package (RStudio). The frequency of each junction was calculated by comparing the depth of the unique junction to the total number of nucleotides in all detected junctions in a library. Junctions were plotted according to the genomic position and colored according to log10 of the frequency using ggplot2 in RStudio.

Recombination frequency was calculated at each genomic position by dividing the number of nucleotides in any junction mapping to the position divided by the total number of nucleotides sequenced at the position. Mean recombination frequencies were compared between MHV-WT and MHV-ExoN(-) for three independent sequencing experiments by a 2-way ANOVA statistical analysis with multiple comparisons corrected through statistical hypothesis testing using the Sidak test.

Identification of sgmRNA and DVG junctions

Forward recombination junctions were classified as either canonical sgmRNA junctions, alternative sgmRNA junctions or DVG junctions based on the position of their junction sites and filtered in Microsoft Excel. Briefly, junction start sites were filtered to those positioned within 30 nucleotides of the TRS-L for each virus. The stop sites were then filtered for those positioned within 30 nucleotides of each respective sgmRNA TRS. This window is supported by other reports defining the flexibility of the CoV transcriptome [45,65]. Canonical sgmRNAs were identified as the most abundant junction matching these criteria. Other, less abundant sgmRNA junctions were categorized as alternative sgmRNAs. The junction frequency (Jfreq) of each sgmRNA was calculated by dividing the number of nucleotides in a specific sgmRNA population by the total amount of viral RNA (total mapped nucleotides). This ratio is multiplied by 10,000 to scale for the number of nucleotides sequenced and is therefore expressed as the number of junctions per 104 mapped nucleotides. The filtered sgmRNA junctions were compiled and DVG junctions were filtered in RStudio by performing an exclusionary anti_join() using dplyr on forward junctions identified in each sample. DVG Jfreq was calculated by dividing the number of nucleotides in DVG junctions by the total amount of viral RNA in a sample (total mapped nucleotides). The ratio is multiplied by 10,000 to scale for number of nucleotides sequenced and the frequency is expressed as the number of junctions per 104 mapped nucleotides. The percentage of canonical and alternative sgmRNA and DVG junctions was calculated by comparing the depth of all filtered sgmRNA or DVG junctions to the sum of all detected forward junctions. Mean percent canonical and alternative sgmRNAs and DVG was compared between three independent sequencing experiments in viral supernatant RNA. Statistical significance was determined by a 2-way ANOVA test with multiple comparisons and corrected by statistical hypothesis testing using the Sidak test.

Differential abundance of junctions

To compare the abundance of junctions in MHV-A59 and MHV-ExoN(-), the ViReMa output list of junctions was analyzed by in-house scripts (https://github.com/DenisonLabVU) and the R package DESeq2 [50]. Junctions significantly up- or down-regulated in MHV-ExoN(-) were visualized using bioinfokit [66] and further mapped according to their genomic positions. Statistical significance was determined by the p-value of each junction calculated by the DESeq2 package in RStudio and junctions with a significant alteration of abundance in MHV-ExoN(-) compared to MHV-WT were visualized as either red or green in the graph generated by bioinfokit.

Nucleotide composition analysis

DVG junctions were filtered as described above and the nucleotide composition at each position was determined. To avoid bias of highly replicated DVGs and to more closely reflect the stochastic nature of RNA recombination, each unique detected junction was counted equally rather than weighting by read count [67]. Sequences were extracted from a sorted BED file listing the junctions using Rec_Site_Extraction.py with a 30-base pair window. Start site and stop site sequences were separated in Microsoft Excel and the nucleotide frequency at each position was calculated using the Biostrings [68] package in RStudio. The mean percentage of a nucleotide was compared between MHV-WT and MHV-ExoN(-) using a 2-way ANOVA test with multiple comparisons and were corrected for false-discovery rate (FDR) using the Benjamini-Hochberg method. Length of microhomology at junction sites were extracted from ViReMa SAM file using the Compiler_Module.py of ViReMa and -FuzzEntry—Defuzz 0 flags. The frequency of overlaps ranging from 0–10 bp was calculated and compared to an expected probability distribution using uHomology.py.

Direct RNA Nanopore alignment and analysis

Live basecalling was performed by Guppy in MinKNOW. Run statistics were generated from each sequencing experiment by NanoPlot [69]. Pass reads from all three experiments were concatenated for each virus and aligned to the genome using minimap2 [70] and FLAIR (Full Length Alternative Isoforms of R NA) [71] to generate alignment files and BED files listing deletions detected in each sequenced RNA molecule. Both BAM and BED files were filtered for full length molecules using samtools and Microsoft Excel, respectively. Full-length CoV molecules were defined as encoding coverage at in the 5’ UTR and 3’ UTR of the respective viruses. Nanopore junctions output in BED files were compared to junctions in ViReMa RNA-seq BED files to confirm its presence in both datasets. To account for noisiness in Nanopore datasets, a Nanopore junction was considered confirmed if at least 1 RNA-seq junction start and stop sites fell within 20 bp of the Nanopore start and stop sites, respectively. Filtering of Nanopore and RNA-seq datasets was performed in Microsoft Excel. BED files generated by the flair align module were parsed based on the number of junctions were identified. Nanopore reads containing only 1 junction were identified using Microsoft Excel and unique junctions were quantified in RStudio using base-R functions. Sequencing coverage maps were generated from samtools depth analysis of filtered BAM files. All junctions present in sequenced libraries were mapped in Sashimi plots generated by the Integrated Genome Viewer (IGV) [72]. Junctions present in full-length MHV RNA molecules with a single deletion were mapped according to their genomic positions as previously described. The genetic architectures of full-length RNA molecules sequenced by direct RNA Nanopore sequencing were visualized by filtering RNA molecules for at least 3 supporting reads. Low frequency variants were removed from this analysis.

Acknowledgements

We thank members of the Denison laboratory for insightful discussions regarding this study.

References

WuF, ZhaoS, YuB, ChenY-M, WangW, Song Z-G, et al A new coronavirus associated with human respiratory disease in China. Nature. 2020 2 3;15.

Patiño-Galindo, FilipI, AlQuraishiM, RabadanR. Recombination and lineage-specific mutations led to the emergence of SARS-CoV-2. bioRxiv. 2020 3 23;2020.02.10.942748. doi: 10.1101/2020.02.10.942748

HuangJ-M, JanSS, WeiX, WanY, OuyangS. Evidence of the Recombinant Origin and Ongoing Mutations in Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2). bioRxiv. 2020 3 19;2020.03.16.993816.

LiX, GiorgiEE, MarichannMH, FoleyB, XiaoC, KongX, et al Emergence of SARS-CoV-2 through Recombination and Strong Purifying Selection. bioRxiv. 2020 3 24;2020.03.20.000885.

YiH. 2019 novel coronavirus is undergoing active recombination. Clin Infect Dis [Internet]. 2020 [cited 2020 Mar 11]; Available from: https://academic.oup.com/cid/advance-article/doi/10.1093/cid/ciaa219/5781085

AnthonySJ, GilardiK, MenacheryVD, GoldsteinT, SsebideB, MbabaziR, et al Further Evidence for Bats as the Evolutionary Source of Middle East Respiratory Syndrome Coronavirus. mBio. 2017 4;8(2).

HonC-C, LamT-Y, ShiZ-L, DrummondAJ, YipC-W, ZengF, et al Evidence of the Recombinant Origin of a Bat Severe Acute Respiratory Syndrome (SARS)-Like Coronavirus and Its Implications on the Direct Ancestor of SARS Coronavirus. Journal of Virology. 2008 2 15;82(4):181926. doi: 10.1128/JVI.01926-07

LauSKP, FengY, ChenH, LukHKH, YangW-H, LiKSM, et al Severe Acute Respiratory Syndrome (SARS) Coronavirus ORF8 Protein Is Acquired from SARS-Related Coronavirus from Greater Horseshoe Bats through Recombination. PerlmanS, editor. J Virol. 2015 10 15;89(20):1053247. doi: 10.1128/JVI.01048-15

LiW, ShiZ, YuM, RenW, SmithC, EpsteinJH, et al Bats are natural reservoirs of SARS-like coronaviruses. Science (New York, NY). 2005 10;310(5748):676679. doi: 10.1126/science.1118391

10 

YusofMF, QueenK, EltahirYM, PadenCR, Al HammadiZMAH, TaoY, et al Diversity of Middle East respiratory syndrome coronaviruses in 109 dromedary camels based on full-genome sequencing, Abu Dhabi, United Arab Emirates. Emerging Microbes & Infections. 2017 Jan;6(1):110.

11 

ChenN, LiS, ZhouR, ZhuM, HeS, YeM, et al Two novel porcine epidemic diarrhea virus (PEDV) recombinants from a natural recombinant and distinct subtypes of PEDV variants. Virus Research. 2017 10 15;242:905. doi: 10.1016/j.virusres.2017.09.013

12 

FengKY, ChenT, ZhangX, ShaoGM, CaoY, ChenDK, et al Molecular characteristic and pathogenicity analysis of a virulent recombinant avain infectious bronchitis virus isolated in China. Poultry Science. 2018 10;97(10):351931. doi: 10.3382/ps/pey237

13 

KirchdoerferRN, WardAB. Structure of the SARS-CoV nsp12 polymerase bound to nsp7 and nsp8 co-factors. Nat Commun. 2019 12;10(1):2342 doi: 10.1038/s41467-019-10280-3

14 

SmithEC, DenisonMR. Coronaviruses as DNA Wannabes: A New Model for the Regulation of RNA Virus Replication Fidelity. RacanielloV, editor. PLoS Pathog. 2013 12 5;9(12):e1003760 doi: 10.1371/journal.ppat.1003760

15 

SmithEC, SextonNR, DenisonMR. Thinking Outside the Triangle: Replication Fidelity of the Largest RNA Viruses. Annu Rev Virol. 2014 11 3;1(1):11132. doi: 10.1146/annurev-virology-031413-085507

16 

SubissiL, PosthumaCC, ColletA, Zevenhoven-DobbeJC, GorbalenyaAE, DecrolyE, et al One severe acute respiratory syndrome coronavirus protein complex integrates processive RNA polymerase and exonuclease activities. Proc Natl Acad Sci USA. 2014 9 16;111(37):E39009. doi: 10.1073/pnas.1323705111

17 

KeckJG, MatsushimaGK, MakinoS, FlemingJO, VannierDM, StohlmanSA, et al In vivo RNA-RNA recombination of coronavirus in mouse brain. Journal of Virology. 1988 5;62(5):18101813. doi: 10.1128/JVI.62.5.1810-1813.1988

18 

MakinoS, KeckJG, StohlmanSA, LaiMM. High-frequency RNA recombination of murine coronaviruses. Journal of Virology. 1986 3;57(3):729737. doi: 10.1128/JVI.57.3.729-737.1986

19 

DufourD, Mateos-GomezPA, EnjuanesL, GallegoJ, SolaI. Structure and functional relevance of a transcription-regulating sequence involved in coronavirus discontinuous RNA synthesis. Journal of Virology. 2011 5;85(10):49634973. doi: 10.1128/JVI.02317-10

20 

SolaI, AlmazánF, ZúñigaS, EnjuanesL. Continuous and Discontinuous RNA Synthesis in Coronaviruses. Annu Rev Virol. 2015 11;2(1):26588. doi: 10.1146/annurev-virology-100114-055218

21 

BrianDA, SpaanWJM. Recombination and coronavirus defective interfering RNAs. Seminars in Virology. 1997;8(2):10111. doi: 10.1006/smvy.1997.0109

22 

MakinoS, FujiokaN, FujiwaraK. Structure of the intracellular defective viral RNAs of defective interfering particles of mouse hepatitis virus. Journal of Virology. 1985 5;54(2):329336. doi: 10.1128/JVI.54.2.329-336.1985

23 

WuH-Y, BrianDA. Subgenomic messenger RNA amplification in coronaviruses. Proceedings of the National Academy of Sciences of the United States of America. 2010 7;107(27):1225712262. doi: 10.1073/pnas.1000378107

24 

SchaadMC, BaricRS. Genetics of mouse hepatitis virus transcription: evidence that subgenomic negative strands are functional templates. Journal of Virology. 1994 12;68(12):81698179. doi: 10.1128/JVI.68.12.8169-8179.1994

25 

SunY, JainD, Koziol-WhiteCJ, GenoyerE, GilbertM, TapiaK, et al Immunostimulatory Defective Viral Genomes from Respiratory Syncytial Virus Promote a Strong Innate Antiviral Response during Infection in Mice and Humans. PLoS Pathog [Internet]. 2015 9 3 [cited 2020 May 17];11(9). Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4559413/

26 

VasilijevicJ, ZamarreñoN, OliverosJC, Rodriguez-FrandsenA, GómezG, RodriguezG, et al Reduced accumulation of defective viral genomes contributes to severe outcome in influenza virus infected patients. PLoS Pathog [Internet]. 2017 10 12 [cited 2020 May 17];13(10). Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5638565/ doi: 10.1371/journal.ppat.1006650

27 

MéndezA, SmerdouC, IzetaA, GebauerF, EnjuanesL. Molecular characterization of transmissable gastroenteritis coronavirus defective interfering genomes: packaging and heterogeneity. Virology. 1996;217(2):495507. doi: 10.1006/viro.1996.0144

28 

PenzesZ, WroeC, BrownTDK, BrittonP, CavanaghD. Replication and Packaging of Coronavirus Infectious Bronchitis Virus Defective RNAs Lacking a Long Open Reading Frame. J VIROL. 1996;70:9 doi: 10.1128/JVI.70.12.8660-8668.1996

29 

LiC, WangH, ShiJ, YangD, ZhouG, ChangJ, et al Senecavirus-Specific Recombination Assays Reveal the Intimate Link between Polymerase Fidelity and RNA Recombination. PfeifferJK, editor. J Virol. 2019 4 17;93(13):e0057619, /jvi/93/13/JVI.00576-19.atom. doi: 10.1128/JVI.00576-19

30 

WoodmanA, LeeK-M, JanissenR, GongY-N, DekkerNH, ShihS-R, et al Predicting Intraserotypic Recombination in Enterovirus 71. PfeifferJK, editor. J Virol. 2018 11 28;93(4):e0205718, /jvi/93/4/JVI.02057-18.atom.

31 

KempfBJ, PeersenOB, BartonDJ. Poliovirus Polymerase Leu420 Facilitates RNA Recombination and Ribavirin Resistance. PerlmanS, editor. J Virol. 2016 10 1;90(19):841021. doi: 10.1128/JVI.00078-16

32 

PoirierEZ, MounceBC, Rozen-GagnonK, HooikaasPJ, StaplefordKA, MoratorioG, et al Low-Fidelity Polymerases of Alphaviruses Recombine at Higher Rates To Overproduce Defective Interfering Particles. J Virol. 2015 12 16;90(5):244654. doi: 10.1128/JVI.02921-15

33 

AgostiniML, AndresEL, SimsAC, GrahamRL, SheahanTP, LuX, et al Coronavirus Susceptibility to the Antiviral Remdesivir (GS-5734) Is Mediated by the Viral Polymerase and the Proofreading Exoribonuclease. SubbaraoK, editor. mBio. 2018 3 6;9(2):e0022118, /mbio/9/2/mBio.00221-18.atom. doi: 10.1128/mBio.00221-18

34 

EckerleLD, LuX, SperrySM, ChoiL, DenisonMR. High fidelity of murine hepatitis virus replication is decreased in nsp14 exoribonuclease mutants. Journal of Virology. 2007 11;81(22):1213512144. doi: 10.1128/JVI.01296-07

35 

EckerleLD, BeckerMM, HalpinRA, LiK, VenterE, LuX, et al Infidelity of SARS-CoV Nsp14-exonuclease mutant virus replication is revealed by complete genome sequencing. PLoS pathogens. 2010 5;6(5):e1000896 doi: 10.1371/journal.ppat.1000896

36 

FerronF, SubissiL, Silveira De MoraisAT, LeNTT, SevajolM, GluaisL, et al Structural and molecular basis of mismatch correction and ribavirin excision from coronavirus RNA. Proc Natl Acad Sci USA. 2018 1 9;115(2):E16271. doi: 10.1073/pnas.1718806115

37 

MaY, WuL, ShawN, GaoY, WangJ, SunY, et al Structural basis and functional analysis of the SARS coronavirus nsp14–nsp10 complex. Proc Natl Acad Sci USA. 2015 7 28;112(30):943641. doi: 10.1073/pnas.1508686112

38 

SmithEC, BlancH, VignuzziM, DenisonMR. Coronaviruses Lacking Exoribonuclease Activity Are Susceptible to Lethal Mutagenesis: Evidence for Proofreading and Potential Therapeutics. DiamondMS, editor. PLoS Pathog. 2013 8 15;9(8):e1003565 doi: 10.1371/journal.ppat.1003565

39 

GammonDB, EvansDH. The 3′-to-5′ Exonuclease Activity of Vaccinia Virus DNA Polymerase Is Essential and Plays a Role in Promoting Virus Genetic Recombination. J Virol. 2009 5;83(9):423650. doi: 10.1128/JVI.02255-08

40 

SchumacherAJ, MohniKN, KanY, HendricksonEA, StarkJM, WellerSK. The HSV-1 Exonuclease, UL12, Stimulates Recombination by a Single Strand Annealing Mechanism. PLOS Pathogens. 2012 8 9;8(8):e1002862 doi: 10.1371/journal.ppat.1002862

41 

MinskaiaE, HertzigT, GorbalenyaAE, CampanacciV, CambillauC, CanardB, et al Discovery of an RNA virus 3’->5’ exoribonuclease that is critically involved in coronavirus RNA synthesis. Proceedings of the National Academy of Sciences of the United States of America. 2006 3;103(13):51085113. doi: 10.1073/pnas.0508200103

42 

GraepelKW, AgostiniML, LuX, SextonNR, DenisonMR. Fitness Barriers Limit Reversion of a Proofreading-Deficient Coronavirus. GallagherT, editor. J Virol. 2019 7 24;93(20):e0071119, /jvi/93/20/JVI.00711-19.atom. doi: 10.1128/JVI.00711-19

43 

CaseJB, LiY, ElliottR, LuX, GraepelKW, SextonNR, et al Murine Hepatitis Virus nsp14 Exoribonuclease Activity Is Required for Resistance to Innate Immunity. GallagherT, editor. J Virol. 2017 10 18;92(1):e0153117. doi: 10.1128/JVI.01531-17

44 

RouthA, JohnsonJE. Discovery of functional genomic motifs in viruses with ViReMa–a Virus Recombination Mapper–for analysis of next-generation sequencing data. Nucleic Acids Res. 2014 1;42(2):e11 doi: 10.1093/nar/gkt916

45 

KimD, LeeJ-Y, YangJ-S, KimJW, KimVN, ChangH. The Architecture of SARS-CoV-2 Transcriptome. Cell. 2020 5 14;181(4):914921.e10. doi: 10.1016/j.cell.2020.04.011

46 

van BoheemenS, de GraafM, LauberC, BestebroerTM, RajVS, ZakiAM, et al Genomic characterization of a newly discovered coronavirus associated with acute respiratory distress syndrome in humans. mBio. 2012 11 20;3(6). doi: 10.1128/mBio.00473-12

47 

ChanJF-W, KokK-H, ZhuZ, ChuH, ToKK-W, YuanS, et al Genomic characterization of the 2019 novel human-pathogenic coronavirus isolated from a patient with atypical pneumonia after visiting Wuhan. Emerging Microbes & Infections. 2020 1 1;9(1):22136. doi: 10.1080/22221751.2020.1719902

48 

PeccoudJ, LequimeS, Moltini-ConcloisI, GiraudI, LambrechtsL, GilbertC. A Survey of Virus Recombination Uncovers Canonical Features of Artificial Chimeras Generated During Deep Sequencing Library Preparation. G3 (Bethesda). 2018 3 27;8(4):112938. doi: 10.1534/g3.117.300468

49 

OgandoNS, Zevenhoven-DobbeJC, MeerY van der, BredenbeekPJ, PosthumaCC, SnijderEJ. The enzymatic activity of the nsp14 exoribonuclease is critical for replication of MERS-CoV and SARS-CoV-2. Journal of Virology [Internet]. 2020 9 16 [cited 2020 Oct 2]; Available from: https://jvi.asm.org/content/early/2020/09/10/JVI.01246-20 doi: 10.1128/JVI.01246-20

50 

LoveMI, HuberW, AndersS. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014 12;15(12):550 doi: 10.1186/s13059-014-0550-8

51 

SawickiSG, SawickiDL. Coronavirus Transcription: A Perspective In: EnjuanesL, editor. Coronavirus Replication and Reverse Genetics [Internet]. Berlin, Heidelberg: Springer Berlin Heidelberg; 2005 [cited 2020 May 25]. p. 3155. (Compans RW, Cooper MD, Honjo T, Koprowski H, Melchers F, Oldstone MBA, et al., editors. Current Topics in Microbiology and Immunology; vol. 287). Available from: http://link.springer.com/10.1007/3-540-26765-4_2

52 

LeeK, LeeSE. Saccharomyces cerevisiae Sae2- and Tel1-Dependent Single-Strand DNA Formation at DNA Break Promotes Microhomology-Mediated End Joining. Genetics. 2007 8;176(4):200314. doi: 10.1534/genetics.107.076539

53 

MaJ-L, KimEM, HaberJE, LeeSE. Yeast Mre11 and Rad1 proteins define a Ku-independent mechanism to repair double-strand breaks lacking overlapping end sequences. Mol Cell Biol. 2003 12;23(23):88208. doi: 10.1128/mcb.23.23.8820-8828.2003

54 

CampagnolaG, McDonaldS, BeaucourtS, VignuzziM, PeersenOB. Structure-function relationships underlying the replication fidelity of viral RNA-dependent RNA polymerases. Journal of Virology. 2015 1;89(1):275286. doi: 10.1128/JVI.01574-14

55 

KimH, EllisVD, WoodmanA, ZhaoY, ArnoldJJ, CameronCE. RNA-Dependent RNA Polymerase Speed and Fidelity are not the Only Determinants of the Mechanism or Efficiency of Recombination. Genes (Basel) [Internet]. 2019 11 25 [cited 2020 Feb 18];10(12). Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6947342/

56 

LangsjoenR, MuruatoA, KunkelS, JaworskiE, RouthA. Differential alphavirus defective RNA diversity between intracellular and encapsidated compartments is driven by subgenomic recombination events [Internet]. Molecular Biology; 2020 3 [cited 2020 Apr 22]. Available from: http://biorxiv.org/lookup/doi/10.1101/2020.03.24.006353

57 

AgostiniML, PruijssersAJ, ChappellJD, GribbleJ, LuX, AndresEL, et al Small-Molecule Antiviral β-d-N4-Hydroxycytidine Inhibits a Proofreading-Intact Coronavirus with a High Genetic Barrier to Resistance. J Virol. 2019 15;93(24). doi: 10.1128/JVI.01348-19

58 

GrahamRL, BeckerMM, EckerleLD, BollesM, DenisonMR, BaricRS. A live, impaired-fidelity coronavirus vaccine protects in an aged, immunocompromised mouse model of lethal disease. Nat Med. 2012 12;18(12):18206. doi: 10.1038/nm.2972

59 

MenacheryVD, YountBL, JossetL, GralinskiLE, ScobeyT, AgnihothramS, et al Attenuation and restoration of severe acute respiratory syndrome coronavirus mutant lacking 2’-o-methyltransferase activity. Journal of Virology. 2014 4;88(8):42514264. doi: 10.1128/JVI.03571-13

60 

ChenW, BaricRS. Molecular anatomy of mouse hepatitis virus persistence: coevolution of increased host cell resistance and virus virulence. Journal of Virology. 1996 6 1;70(6):394760. doi: 10.1128/JVI.70.6.3947-3960.1996

61 

YountB, DenisonMR, WeissSR, BaricRS. Systematic Assembly of a Full-Length Infectious cDNA of Mouse Hepatitis Virus Strain A59. Journal of Virology. 2002 11 1;76(21):1106578. doi: 10.1128/jvi.76.21.11065-11078.2002

62 

ScobeyT, YountBL, SimsAC, DonaldsonEF, AgnihothramSS, MenacheryVD, et al Reverse genetics with a full-length infectious cDNA of the Middle East respiratory syndrome coronavirus. PNAS. 2013 10 1;110(40):1615762. doi: 10.1073/pnas.1311542110

63 

BolgerAM, LohseM, UsadelB. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014 8 1;30(15):211420. doi: 10.1093/bioinformatics/btu170

64 

LiH, HandsakerB, WysokerA, FennellT, RuanJ, HomerN, et al The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009 8 15;25(16):20789. doi: 10.1093/bioinformatics/btp352

65 

IrigoyenN, FirthAE, JonesJD, ChungBY-W, SiddellSG, BrierleyI. High-Resolution Analysis of Coronavirus Gene Expression by RNA Sequencing and Ribosome Profiling. PLoS Pathog [Internet]. 2016 2 26 [cited 2020 Feb 19];12(2). Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4769073/ doi: 10.1371/journal.ppat.1005473

66 

Bedre R. reneshbedre/bioinfokit [Internet]. 2020 [cited 2020 Feb 27]. Available from: https://github.com/reneshbedre/bioinfokit

67 

JaworskiE, RouthA. Parallel ClickSeq and Nanopore sequencing elucidates the rapid evolution of defective-interfering RNAs in Flock House virus. PLoS pathogens. 2017 5;13(5):e1006365 doi: 10.1371/journal.ppat.1006365

68 

PagèsH, AboyounP, GentlemanR, DebRoyS. Biostrings: Efficient manipulation of biological strings [Internet]. Bioconductor version: Release (3.10); 2020 [cited 2020 Apr 9]. Available from: https://bioconductor.org/packages/Biostrings/

69 

De CosterW, D’HertS, SchultzDT, CrutsM, Van BroeckhovenC. NanoPack: visualizing and processing long-read sequencing data. BergerB, editor. Bioinformatics. 2018 8 1;34(15):26669. doi: 10.1093/bioinformatics/bty149

70 

LiH. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018 9 15;34(18):3094100. doi: 10.1093/bioinformatics/bty191

71 

TangAD, SouletteCM, BarenMJ van, HartK, Hrabeta-RobinsonE, WuCJ, et al Full-length transcript characterization of SF3B1 mutation in chronic lymphocytic leukemia reveals downregulation of retained introns [Internet]. Genomics; 2018 9 [cited 2020 Jan 2]. Available from: http://biorxiv.org/lookup/doi/10.1101/410183

72 

RobinsonJT, ThorvaldsdóttirH, WincklerW, GuttmanM, LanderES, GetzG, et al Integrative Genomics Viewer. Nat Biotechnol. 2011 1;29(1):246. doi: 10.1038/nbt.1754
This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.1371/journal.ppat.1009226&title=The coronavirus proofreading exoribonuclease mediates extensive viral recombination&author=&keyword=&subject=Research Article,Biology and life sciences,Molecular biology,Molecular biology techniques,Sequencing techniques,RNA sequencing,Research and analysis methods,Molecular biology techniques,Sequencing techniques,RNA sequencing,Biology and life sciences,Organisms,Viruses,RNA viruses,Coronaviruses,SARS coronavirus,SARS CoV 2,Biology and life sciences,Microbiology,Medical microbiology,Microbial pathogens,Viral pathogens,Coronaviruses,SARS coronavirus,SARS CoV 2,Medicine and health sciences,Pathology and laboratory medicine,Pathogens,Microbial pathogens,Viral pathogens,Coronaviruses,SARS coronavirus,SARS CoV 2,Biology and life sciences,Organisms,Viruses,Viral pathogens,Coronaviruses,SARS coronavirus,SARS CoV 2,Biology and life sciences,Organisms,Viruses,RNA viruses,Coronaviruses,Biology and Life Sciences,Microbiology,Medical Microbiology,Microbial Pathogens,Viral Pathogens,Coronaviruses,Medicine and Health Sciences,Pathology and Laboratory Medicine,Pathogens,Microbial Pathogens,Viral Pathogens,Coronaviruses,Biology and Life Sciences,Organisms,Viruses,Viral Pathogens,Coronaviruses,Biology and Life Sciences,Genetics,Genomics,Biology and Life Sciences,Molecular Biology,Molecular Biology Techniques,Gene Mapping,Nucleotide Mapping,Research and Analysis Methods,Molecular Biology Techniques,Gene Mapping,Nucleotide Mapping,Biology and Life Sciences,Biochemistry,Proteins,Recombinant Proteins,Biology and Life Sciences,Molecular Biology,Molecular Biology Techniques,Sequencing Techniques,Nucleotide Sequencing,Research and Analysis Methods,Molecular Biology Techniques,Sequencing Techniques,Nucleotide Sequencing,Biology and Life Sciences,Microbiology,Virology,Viral Replication,