A Case for the Ancient Origin of Coronaviruses

Representative coronavirus genomes (n = 43; Table 1) from the four genera were selected for this study. For several viral taxa, genomic sequences from multiple sampling years were included, though these dates were not used explicitly in our analyses. To avoid highly recombinant sequences and/or sequence misalignments, we specifically selected partial viral sequences from five relatively conserved genomic regions (Fig. 1 and Table 1): Nsp15-16 (1,320 nucleotides [nt]), the matrix protein (640 nt), papain-like protease 2 (PLP2; 620 nt), the RdRp (1,860 nt), and the Y domain (400 nt) (18). These nucleotide sequences were aligned at the amino acid level by MUSCLE (19) as described previously (20).

Each of the five target sequences was screened for recombination with a genetic algorithm tool, GARD (21). Maximum-likelihood phylogenies were constructed for each nonrecombinant region with PHYML 3.0 (22), implemented in Seaview 4.0 (23), with a subtree-pruning and regrafting search algorithm and the general time-reversible substitution model with a four-bin gamma rate distribution (GTR + Γ₄).

With these topologies, branch lengths were re-estimated in HyPhy (24) under GTR + Γ₄ and a branch site random effects likelihood (BS-REL) model (25). The BS-REL model was implemented to account for the effects of variable selection pressure across codon sites and phylogenetic lineages by inferring three unconstrained ω classes (dN/dS: nonsynonymous substitution rate/synonymous substitution rate) for each branch and estimating the proportion of sites in the alignment that evolved under each ω class.

To estimate the variance in branch length estimates produced by the BS-REL model, we used a modified Latin hypercube sampling (LHS) importance resampling scheme (M = 500 samples), described in detail previously (26). LHS is a standard technique that can be efficiently parallelized, an important consideration here. It is used to assess the variability of high-dimensional distributions by discretizing the volume of parameter space around the maximum-likelihood estimate into N bins (10,000 in our case) of approximately equal probability along each coordinate and sampling the bin for each parameter in a way that no two parameters share the same bin index (this technique maximizes space coverage). Samples are reweighted according to their likelihood, and resamples (M = 500) are drawn from this distribution. A modified version of LHS (26), which we use here, has also been shown to compare favorably to other techniques for the assessment of parameter uncertainty.

Phylogenies with interpretable branch lengths can be inferred only when analyzing nonrecombinant regions. Therefore, we screened each of the five highly conserved coronavirus target regions with GARD (21). Within the Nsp15-16 region, two recombination breakpoints were detected, and phylogenetic incongruity was confirmed by a Kishino-Hasegawa test (P < 0.01) (27, 28). The first nonrecombinant section encompassed Nsp15 and 80 nt of Nsp16, which we refer to as Nsp15⁺; the second section fell entirely within Nsp16. Because of its short length (300 nt), the third section was not included in later analyses. Within the RdRp region, a single recombination breakpoint was detected with GARD, but it is unlikely that this breakpoint represents a true recombination event since the two topologies were not significantly different according to the Kishino-Hasegawa test. This putative breakpoint detected in the RdRp region was likely due to elevated rate variation across the sequence (21). No recombination breakpoints were detected in the other three loci.

The maximum-likelihood phylogenies inferred from the six nonrecombinant loci were almost all significantly different from each other (Shimodaira-Hasegawa test, P < 0.05) (29). This finding suggests that the different regions of the coronavirus genome have distinct evolutionary histories and should not be treated as a single region in a phylogenetic analysis. The lone exception was the comparison of the topologies from RdRp and Nsp15⁺ (P = 0.131). Nevertheless, we chose be conservative and analyze each of the six loci independently.

We employed two models of nucleotide sequence evolution (GTR + Γ₄ and BS-REL) to re-estimate the branch lengths of the inferred maximum-likelihood phylogenies. GTR + Γ₄ is a standard nucleotide substitution model that has been commonly in many evolutionary virology studies. In contrast, the BS-REL model is a codon model that explicitly accounts for variation in selection pressures across sites and lineages (25). Both evolutionary models produced comparable length estimates for shorter branches (≤0.05 substitution per site in Fig. 2). However, compared to results inferred from BS-REL, long branches in the coronavirus phylogenies were underestimated by GTR + Γ₄ (Fig. 2). Notably, there were a substantial number of long branches in which the number of substitutions per site approached saturation in the BS-REL model, ranging from 8 branches in Nsp16 to 16 branches in the Y domain (Table 2). The lengths of these saturated branches cannot be reliably estimated, although a meaningful lower bound on their lengths can be obtained (e.g., with LHS; see below).

Among the six nonrecombinant regions, there was substantial variation in the total tree length expansion between GTR + Γ₄ and BS-REL, ranging over 3 orders of magnitude (Table 2). Differing selective pressures along lineages among loci appear to account for this variation, as we observed a correlation between the strength of purifying selection along a branch (i.e., lower mean ω) and the relative increase in branch length under BS-REL compared to GTR + Γ₄ (Spearman's nonparametric correlation, P ≤ 0.0001; branches in which either model inferred a length of zero substitutions per site were excluded from this analysis). Thus, lineages that bear the mark of stronger purifying selection generally experienced a greater increase in length under the BS-REL model. Moreover, the number of branches approaching saturation was correlated with a greater expansion in the total length under BS-REL (Spearman's nonparametric correlation, P < 0.05). Interestingly, simpler, site-based measurements of selection did not correlate with total tree length expansion (e.g., mean dN/dS [P = 0.87] and number of sites of pervasive purifying selection [P = 0.46]), though analyzing only six loci limits the statistical power of these tests. Nonetheless, the relationship between ω and branch length expansion supports the hypothesis that increased purifying selection is responsible for the underestimation of the lengths of long branches in RNA virus phylogenies.

The variance in branch length expansion from GTR + Γ₄ to BS-REL, as approximated with LHS, differed among coronavirus loci (Table 2). For the Nsp15⁺, Nsp16, matrix protein, and PLP2 loci, the upper and lower 95% confidence limits differed by orders of magnitude, suggesting highly imprecise estimates of branch length expansion in these loci. This lack of precision is not surprising, given the extent to which these loci have experienced saturation along their phylogeny. Once a site becomes saturated with repeated substitutions, it is very difficult for an evolutionary model to determine the exact number of substitutions that have taken place. Therefore, the maximum-likelihood estimate will be imprecise. Nevertheless, these results clearly suggested that there is substantial underestimation of the evolutionary history of coronaviruses.

To examine the effects of saturation and determine the limit of reliable divergence inference with BS-REL, we simulated sequence alignments in HyPhy (24) under the maximum-likelihood parameter values from the BS-REL model fitted to the RdRp locus. HyPhy scripts and input files needed to perform the simulations can be downloaded at https://github.com/veg/pubs/tree/master/SARS. The tree inferred by using this locus experienced a dramatic expansion in a previous BS-REL analysis. For example, a single internal branch at this locus accounted for >70% of the total inferred tree length, but it also experienced incredibly strong purifying selection (dN/dS = 0) at 97.5% of the sites and very strong diversifying selection (dN/dS > 100) at the remaining 2.5% of the sites. In the simulations (n = 100), we kept the relative lengths of branches fixed and scaled the entire tree length from 1 to 10,000 expected substitutions/site. We then fitted the GTR + Γ₄ and BS-REL models to the replicates to evaluate the saturation behavior of the models. There were striking differences between the models (Fig. 3). GTR + Γ₄ consistently, and progressively more badly, underestimated the simulated tree length, starting at the lowest simulated values; this behavior is due to a long internal branch predominantly subject to strong purifying selection, which poses a challenge for the nucleotide-based model (30). Note that this type of saturation does not lead to infinite branch lengths because a majority of the sites are maintained by strong purifying selection. In contrast, the BS-REL model did not saturate until the total length of the tree was on the order of 1,000 substitutions per site. Interestingly, BS-REL estimates are biased slightly upward, largely because of the saturation of nonsynonymous substitutions at positively selected sites and corresponding infinite values of dN/dS.

The current best estimate of the tMRCA of all four coronavirus genera of 10,141 years ago was obtained by Woo et al. (10) by using the RdRp locus; their tMRCA estimate is the only one that includes all four coronavirus genera. Notably, their tMRCA estimates for alpha- and betacoronaviruses are similar to the dates inferred by a similar methodology from a different data set (31). Our results suggest that, despite the comprehensive nature of these studies, there are intrinsic shortcomings in the standard evolutionary models used to estimate branch lengths and the tMRCA of ancient viral lineages like coronaviruses. Nevertheless, Woo et al. calibrated the coronavirus molecular clock with serially sampled sequences, which rely predominantly on the relationship between the branches near the tips of the phylogeny to estimate substitution rates. Notably, it is possible to infer consistent RNA virus substitution rates from these types of serially sampled sequences, even if the tMRCA of older lineages deep in the tree cannot be reliably estimated because of saturation along internal branches (e.g., simian immunodeficiency virus/human immunodeficiency virus [32, 33] and avian influenza virus [17, 34]).

When variation in selection pressure was taken into account using BS-REL, we observed a dramatic expansion of the total length of the RdRp phylogeny: 28,860.4-fold (95% confidence interval [CI], 18,695.7-fold to 48,256.0-fold) (Table 2). If one were to extrapolate the evolutionary ages estimated by Woo et al. to the BS-REL phylogeny, the adjusted tMRCA of coronaviruses would be 293 (95% CI, 190 to 489) million years ago, 4 orders of magnitude greater than that previously inferred. The robustness of this tMRCA estimate is highly dependent on both the rate of evolution estimated by Woo et al. (10), which is inferred primarily from short branches near the tips of the phylogeny, and the accuracy of the new branch lengths, which appear extremely variable. Nevertheless, even with the conservative lower bounds of the tMRCA of coronaviruses and branch length expansion (974 BCE [10] and 18,695.7-fold expansion [Table 2]), the adjusted coronavirus tMRCA of 55.8 million years ago would still be 3 orders of magnitude greater than the current estimate.