Estimation of the potential recombination breakpoint location.

The data set was further analyzed using single-breakpoint estimation algorithms implemented in Genetic Algorithms for Recombination Detection (GARD) and Likelihood Analysis of Recombination in DNA (LARD). Based on the bootscan analysis, only the 2,000 nucleotides (nt) around the open reading frame 1b (ORF1b)/S junction (nt 20150 to 22202; all nucleotide numberings in this study are based on AY274119) were analyzed in order to increase the precision of recombination breakpoint estimation. Based on the RDP results, three selected taxa, Rp3, Tor2, and Rm1, were used in the analyses described below. Briefly, GARD uses a genetic algorithm to search for the best breakpoint locations (23). LARD uses a maximum likelihood (ML) method and a likelihood ratio test (LRT) to access the significance of the inferred breakpoint (15). To demonstrate that the detected recombination event is not likely to be a result of random chance (15), the likelihood ratio (LR) of our data set was evaluated against the null distributions of LRs of 1,000 simulated data sets, assuming no recombination, using Seq-Gen (42).

Investigation of the phylogenetic origin of the potential parents.

The genome regions 5′ upstream and 3′ downstream of the estimated breakpoint were designated major and minor parental regions, respectively. To investigate the phylogenetic origins of these potential parents, coding sequences of essential ORFs of the major (i.e., ORF1) and minor (i.e., S, E, M, and N genes) parental regions of selected CoV strains (n = 13) were aligned independently using ClustalX based on their codon sequences. The aligned ORFs of the two parental regions were degapped and concatenated separately, generating two alignments of 20,085 bp and 5,778 bp for the major and minor parental regions, respectively. For each of the parental regions, phylogenies were constructed using the Bayesian Markov chain Monte Carlo (BMCMC) method. The BMCMC analyses summarized the majority consensus trees produced by two sets of four tempered MCMC chains of 10⁷ states sampled every 1,000th generation, with the initial 10% of states discarded. The Bayesian phylogenetic analysis was performed with MRBAYES 3 (44) under the best-fit substitution model determined by MRMODELTEST 2 (http://people.scs.fsu.edu/~nylander/). According to the BMCMC phylogeny (see Fig. ​Fig.2A),2A), the major parental lineage of Rp3 is designated the human-bat SLCoV (HB-SLCoV) lineage based on its close phylogenetic relationship with the Hu-SCoV lineage.

Open in a separate window

FIG. 2.

Phylogenetic origins of the major and minor parental regions of Rp3. ML phylogenies were constructed from the concatenated sequences of the essential ORFs of the major (A) and minor (B) parental regions of selected CoVs. For the purposes of display, the phylogenies were midpoint rooted. The taxa were annotated according to their accession numbers and host species—civets (C), humans (H), or bats (B)—and strain names. The numbers on the left of the nodes refer to the BMCMC posterior probabilities. The percentages of support for all other internal nodes within the two lineages were omitted for simplicity. The recombinant strain Rp3, the most recent common ancestor of Hu-SCoVs (MRCA-Hu), and the divergence event between Hu-SCoVs and HB-SLCoVs (DIV-Hu/HB) are indicated. The scale bars are in units of nucleotide substitutions per site.

Estimation of the time of the divergence events.

To estimate the time of the most recent common ancestor (tMRCA) of Hu-SCoVs, as well as the time of divergence events (tDIV) between the Hu-SCoV and HB-SLCoV lineage (designated tMRCA-Hu and tDiv-Hu/HB, respectively) (see Fig. ​Fig.2),2), coding sequences of S1 (nt 21492 to 22784; n = 36) and ORF1 (nt 898 to 21479; n = 24) were analyzed under various molecular-clock models in both ML and Bayesian frameworks (details of the taxa in the two data sets are listed as supplementary material at http://evolution.hku.hk/SARS_dating.htm). The sampling times (i.e., the month and year) of the taxa were collected from the literature and used as calibration points in the clock models (41).

First, the strict molecular clock (i.e., a constant rate of evolution) of the two data sets was evaluated in an ML framework using PAML 3.15 as previously described (16, 41, 53). Briefly, the performances of the single-rate dated-tip (SRDT) (i.e., strict-clock) and the different-rate (DR) (i.e., no-clock) models in the data sets were compared using an LRT. Second, the two data sets were analyzed under the strict-clock model (CLOC), as well as the uncorrelated exponentially and lognormally distributed relaxed-clock models (UCED and UCLN) in a Bayesian framework. The CLOC model assumed a constant rate of evolution throughout the tree. The UCED and UCLN models assumed independent rates on different branches, which were drawn from an underlying exponential and lognormal distribution, respectively (6). These clock models are implemented in BEAST 1.4 (8). The MCMC chains were run for 5 × 10⁶ (S1 data set) or 1 × 10⁸ (ORF1 data set) states sampled every 1,000 generations with the initial 10% of burn-in samples discarded (7). For both data sets, the best-fit substitution model was the general time-reversible (GTR) model allowing four categories of gamma-distributed rate heterogeneity distribution and a proportion of invariant sites (GTR + Γ₄ + I), as determined by MODELTEST. Since the past population dynamics of the data sets were not the primary interest of our study, we assumed a constant coalescent tree prior for all analyses, with a Jeffreys prior on the constant population size hyperparameter (7). To investigate if this tree prior biased our date estimation, we also analyzed our data sets using a Yule tree prior, which assumes a constant speciation rate per lineage (6). All MCMC chains were independently run twice for the same analysis.

To use information from the S1 data set to improve our estimate of tDIV-Hu/HB from the ORF1 data set, an S1-derived prior distribution was specified on tMRCA-Hu, which is a divergence event shared by the phylogenies of both data sets. This prior distribution was based on the posterior distribution of tMRCA-Hu estimated from the S1 data set under the best-fit clock model. The mode and parameters of this distribution were estimated using distribution-fitting software, EasyFit 3.2 (MathWave Technologies). The MCMC chains for the ORF1 data set were rerun under the same configurations described above, except an S1-derived prior was specified on tMRCA-Hu. For all Bayesian analyses, median and the highest posterior density regions at 95% (HPD) of the parameters were summarized from two identical but independent MCMC chains using TRACER 1.3 (http://beast.bio.ed.ac.uk/). The adequacy of sampling was assessed via effective sample size, which was larger than 200 for all summary statistics investigated (all xml files for BEAST are available as supplementary material at http://evolution.hku.hk/SARS_dating.htm).

Phylogenetic origin of the putative parental strains.

To investigate the phylogenetic origin of the putative parents, two BMCMC phylogenies were constructed based on the major (Fig. ​(Fig.2A)2A) and minor (Fig. ​(Fig.2B)2B) parental regions, respectively. The minor parental region of Rp3 was clustered within the Bt-SLCoV lineage and shared monophyly with Rm1 and BtCoV/279/2005 (Fig. ​(Fig.2B).2B). This suggests that the potential minor parent of Rp3 is probably a Bt-SLCoV that shared a close phylogenetic relationship with Rm1 and BtCoV/279/2005. It has been suggested that there is species-specific host restriction of CoVs in bats, since most CoVs from a single bat species grouped together in phylogenetic analyses (48). Moreover, the S protein (which is located within the minor parental region) is the primary determinant of species specificity in CoVs (12, 36), and thus, we speculate that this minor parent may be a Bt-SLCoV residing in Rhinolophus pearsoni, i.e., the host species of Rp3.

On the other hand, the major parental region of Rp3 grouped with, but clustered outside of, the Hu-SCoV lineage (Fig. ​(Fig.2A).2A). Based on this observation, the potential major parent of Rp3 is possibly derived from an uncharacterized lineage that is phylogenetically closely related to Hu-SCoVs. The host species of this speculative parental lineage cannot be ascertained, as it was clustered within neither the Hu-SCoV nor the Bt-SLCoV lineage. Here, we outline three possibilities regarding the host species of this lineage. First, the lineage may originate from an unsampled group of phylogenetically distinct SCoVs residing in live-animal market mammals, like civets or racoon dogs. However, extensive surveillances of various mammalian species over a wide range of geographic locations have been performed, and only CoVs that are highly similar to SCoVs in humans were sampled (20). Thus, this possibility seems unlikely. Second, the lineage may originate from an unknown nonbat intermediate host species, which possibly acquired a SLCoV from bats and transmitted the virus to an amplifying host, such as civets, resulting in spillover in live-animal markets in southern China. However, one of the prerequisites for recombination is coinfection of parental strains within an individual. Therefore, recombination of parental strains residing in different species, i.e., bats and the unknown intermediate host in this case, may be rare due to the relatively strict tropism barrier of CoVs (12, 52). Third, the strain may originate from an unsampled SLCoV lineage residing in a bat species that is phylogenetically closer to Hu-SCoVs than all other currently sampled Bt-SLCoVs. Based on the relatively high genetic diversity among the currently sampled Bt-SLCoVs, the existence of an unsampled phylogenetically distinct lineage of Bt-SLCoV is highly likely, and therefore, the third hypothesis seems to be the most plausible. In the discussions below, this parental lineage is therefore referred to as the HB-SLCoV lineage, while the term “Bt-SLCoV lineage” refers to all other sampled Bt-SLCoVs (Fig. ​(Fig.2).2). This lineage is proposed to contain the major parent of Rp3 and other closely related strains, and we cannot exclude the possibility that the lineage may also contain the direct ancestor of Hu-SCoVs. To further investigate the time of this interspecies transmission event, tMRCA-Hu and tDIV-Hu/HB (Fig. ​(Fig.2)2) were estimated under various molecular-clock models in both ML and Bayesian frameworks.

Molecular clock-like behavior of the data sets and choice of Bayesian clock models.

For the ORF1 data set, under the ML framework, LRT analysis suggests that the SRDT model should be rejected in favor of the DR model (Table ​(Table1).1). Moreover, BF analysis suggests that the UCED model fits the ORF1 data set significantly better than the other two models (Table ​(Table2),2), implying that the rate variations among branches of the ORF1 phylogeny are significant and that a strict clock cannot be assumed. The Bt-SLCoV lineage may contribute to the rate variations in the ORF1 data set, since CoVs of different hosts (i.e., bats and humans or civets) may have different substitution rates.

For the S1 data set, LRT analysis suggests the SRDT model cannot be rejected (Table ​(Table1).1). Moreover, the performance of the CLOC model is not significantly worse than that of the UCLN model, implying that the rate variations may not be significant among branches in the S1 phylogeny (Table ​(Table2).2). However, the BF analysis also suggests the UCED model performed slightly better than the CLOC model. Nonetheless, the tMRCAs estimated under the relaxed- and strict-clock models are generally consistent (Fig. ​(Fig.3A),3A), suggesting that these rate variations did not have a significant impact on our estimates of tMRCA-Hu. Based on the marginal likelihoods and the BF analysis (Table ​(Table2),2), the estimates under the UCED model were taken as the final dating results of both data sets.

Open in a separate window

FIG. 3.

tMRCA-Hu estimated from the S1 data set. (A) tMRCA-Hu estimated from the S1 data set under various Bayesian clock models and the ML SRDT model. (B) Posterior MCMC samples (left y axis) of tMRCA-Hu estimated from the S1 data set under the UCED model and the lognormal distribution (right y axis) fitted using Easyfit. The values of the parameters for the lognormal distribution are as follows: σ = 0.56, μ = −1.00, and γ = 2.04.

tMRCA-Hu.

Based on the analysis of the S1 data set under the UCED model, tMRCA-Hu was estimated to be at a median of 2002.74 (HPD, 2002.18 to 2003.04). This time point refers to the emergence of the common ancestor of all Hu-SCoVs. Under modest assumptions, i.e., that the root had been sampled and the emergence was the result of a single cross-species infection of a single viral lineage, this time point can be considered an estimate of the theoretical onset of the 2003 SARS outbreak. Our estimation of tMRCA-Hu is consistent with previous estimations (5, 47, 55).

tDIV-Hu/HB.

A prior was specified on tMRCA-Hu as a lognormal distribution with parameters chosen to fit the posterior distribution estimated from the S1 data set (Fig. ​(Fig.3B).3B). Bayesian inference specifically provides for the incorporation of prior knowledge, and in this way, we were able to combine information from both data sets in the estimation of tDIV-Hu/HB. Under the UCED model, the medians of tMRCA-Hu estimated from the ORF1 data set with or without the S1-derived tMRCA-Hu prior were similar, and the posterior distribution of tMRCA-Hu was not solely dependent on its prior distribution (Fig. ​(Fig.4A),4A), suggesting that the ORF1 data set was providing additional information in the Bayesian inference. Moreover, tDIV-Hu/HB was consistently estimated at a median around the late 1990s with or without the S1-derived tMRCA-Hu prior (Fig. ​(Fig.4B).4B). It was noted that the specification of S1-derived tMRCA-Hu priors substantially narrowed the HDP of the tDIV-Hu/HB estimate by about 40%, i.e., it decreased from 12.8 to 7.7 years (Table ​(Table3).3). Similar results were observed under the UCLN model (the data are not shown for simplicity).

Open in a separate window

FIG. 4.

Specification of an S1-derived lognormal tMRCA-Hu prior in the analysis of the ORF1 data set under the UCED model. (A) Prior and posterior distributions of tMRCA-Hu. (B) Effects of the tMRCA-Hu prior on the posterior distribution of tDIV-Hu/HB.

The relatively low substitution rate of ORF1 and the rate variation among branches in the ORF1 phylogeny may limit the power of the molecular-clock analysis on the ORF1 data set. These factors add uncertainty to our analysis and may widen the credible intervals of our estimates. In contrast, the S1 coding region has been identified as the most variable region among SCoVs (40) and the S1 data set was found to be more clocklike than the ORF1 data set. Therefore, we specified an informative prior on tMRCA-Hu based on the S1 data set, allowing us to combine information across the data sets and to reduce the uncertainty of our divergence time estimates.

Assuming there was an interspecies transmission of HB-SLCoVs from bats to an amplifying host (e.g., civets), the upper and lower bounds of this event should be theoretically represented by tDIV-Hu/HB and tMRCA-Hu, respectively (Fig. ​(Fig.5).5). Therefore, the time period between these two events can be considered the most conservative estimation of the period between the cross-species event and the onset of the epidemic. The median and HPD of this period were summarized by sampling the length of a particular branch (i.e., branch A in Fig. ​Fig.5)5) of all time-scaled MCMC phylogenies under the UCED model. This period was estimated at a median of 4.08 years (HPD, 1.45 to 8.84 years) (Table ​(Table3).3). The estimated mean substitution rate of the ORF1 data set under the UCED model was 2.79 × 10⁻³ (HPD, 1.64 × 10⁻³ to 4.35 × 10⁻³) substitution per site per year. This estimate is comparable to a previous estimation for the whole genome of Hu-SCoV (i.e., 0.80 × 10⁻³ to 2.38 × 10⁻³) (55) and is at the same order of magnitude as in other RNA viruses (4, 9, 18, 37, 38, 50). In addition, the ORF1 data set was reanalyzed under the UCED model with a Yule tree prior assumption, and the estimate is generally consistent with the estimate under the constant coalescent tree prior assumption, suggesting our date estimation is robust for the choice of tree priors.

Open in a separate window

FIG. 5.

Estimation of the window period between the cross-species event and the onset of the 2003 SARS epidemic. This time-scaled phylogeny was summarized from all MCMC phylogenies of the ORF1 data set analyzed under the UCED model with the S1-derived tMRCA-Hu prior. The heights of the nodes are represented by the median of their estimates. The HPD of tMRCA-Hu and tDIV-Hu/HB are indicated by gray boxes at these nodes. The taxa were labeled in the same style as in Fig. ​Fig.2,2, except their sampling dates were annotated.

This work was supported by the Research Fund for the Control of Infectious Diseases (reference number 06060672) from the Hong Kong SAR government.

We thank Susanna K. P. Lau of the Department of Microbiology, Faculty of Medicine, University of Hong Kong, for her valuable comments on the manuscript.