Phylogenomics and bioinformatics of SARS-CoV.

Can phylogeneticists and bioinformaticians help virologists to tackle SARS-CoV (severe acute respiratory syndrome coronavirus)? Phylogenetic methodology has progressed almost beyond recognition in the past decade and the study of phylogenetic relationships among species is now a valuable source of information in a variety of biological fields. The widespread use of a reliable statistical formalisation in phylogenetic and bioinformatic studies is necessary to extract the maximum information from sequence data [1].

Models of virus evolution

Recent work, although reporting impressive insights into the mechanisms of pathogenesis of SARS-CoV, has assessed the phylogeny of proteins of this virus using overly simple algorithms 2, 3. The distance methods used to assess SARS-CoV phylogeny present several disadvantages. First, by converting a sequence alignment to pairwise distances we necessarily lose some of the evolutionary information contained within the analysed sequences [1]. Second, distance methods are known to compromise the accuracy of estimates of evolutionary divergence, which are fundamental to understanding the rate and mode of viral evolution. Third, there are no known methods to test evolutionary models and estimated trees produced using the pairwise distance methodology [1].

Here we use state-of-the-art phylogenetic methods to analyse all the available coronaviridae and SARs-CoV sequence datasets to gain an insight into the origin and evolution of SARS-CoV and to narrow down the list of potential regions of its genome that might be interesting targets for drug design.

The statistically most robust method that can be used to achieve these aims is to consider the phylogenetic inference problem in a likelihood framework, using a valid model of evolution for viral genomes [1]. The choice of such a model for single-stranded RNA virus genomes is difficult. A parametric model, based on chemical or biological properties of RNA, might underestimate or completely miss important unknown constraints; for example, packaging of single-stranded RNA genomes that requires interaction with coat proteins [4]. An alternative approach is to use empirical models that are generated through comparisons of observed sequences; for example, simply counting apparent replacements between closely related sequences. Given that sequence databases are biased towards mammalian and bacterial DNA sequences, there are relatively few coding single-stranded RNA sequences to be aligned, and non-coding RNA sequences, such as rRNA or tRNA, might be subjected to different selective (e.g. structural) constraints. Because the proteins encoded by RNA genes might be subjected to functional constraints in a similar manner to non-viral proteins, it might be better to use empirical amino acid substitution models that describe the probability of fixation of amino acid changes rather than RNA models. Furthermore, relative to primary structure, the secondary structure of homologous proteins persists long after any statistically significant sequence similarity has vanished; sequences with 25% amino acid identity probably have the same secondary structure [5]. Amino acid models of evolution that incorporate protein structural information perform better than simple amino acid models [6]. Moreover, selection pressures act on protein function, which in turn is closely related to structure. Therefore, incorporating structure information into evolutionary analysis can assist in incorporating selective constraints. Here, the programme Passml-TM, which implements protein secondary structure-based models of evolution, has been used for analyses 7, 8. The first undertaking is the determination and rooting of the coronavirus phylogeny, that is, the putative origin of the sequences of interest.

Rooting the Coronaviridae and SARS-CoV phylogenies

A large fraction (∼70%) of the SARS-CoV genome encodes a replicase polyprotein, which has significant sequence homology to all of the replicase polyproteins that have been sequenced to date from the order Nidovirales (comprising the Coronavirus, Torovirus, Okavirus and Arterivirus genera). Therefore, this protein is a good choice for investigating the phylogenetic relationships among the family Coronaviridae. Viral sequences have high mutation rates and, consequently, alignments are usually difficult to prepare. Here, ClustalW is applied using standard parameters [9], followed by careful refinement of the alignments both by eye and by using the protein secondary structural information for each nidovirus sequence as predicted by PHD (http://cubic.bioc.columbia.edu/predictprotein/) [10] and PSI-PRED (http://bioinf.cs.ucl.ac.uk/psipred) [11]. Figure 1 shows the maximum likelihood tree produced using a set of homologous replicases from five SARS-CoV strains, 12 other coronaviruses representing both groups 1 and 2 of the genus 2, 3, one torovirus (Breda virus) and one okavirus [yellow head (YH) virus], which were determined to most closely represent the consensus coronavirus sequence by a PSI-Blast search [12]. The coronavirus sequences allow the determination of the root of the SARS-CoV strains, whereas the torovirus and okavirus sequences provide insights into the rooting of the family Coronaviridae phylogeny, which is closer to group 1 than group 2 of the Coronavirus genus.

Figure 1

Maximum likelihood tree produced using Passml-TM [7] based on the replicase proteins from members of of the order Nidovirales, comprising 17 coronaviruses [from group 1: human coronavirus 229E (HCoV-229E), porcine epidemic diarrhea virus (PEDV), porcine transmissible gastroenteritis virus (TGEV), canine coronavirus (CCoV) and feline coronavirus (FCoV); from group 2: bovine coronavirus (BCoV), rat coronavirus (RtCoV), murine hepatitis virus (MHV), human coronavirus OC43 (HCoV-OC43) and porcine hemagglutinating encephalomyelitis virus (HEV); intermediate between groups 1 and 2: avian infectious bronchitis virus (IBV) and turkey coronavirus (TCV); and SARS-CoVs from strains Tor2, bj01, bj02, bj03 and CUHKW1], one torovirus (Breda) and one okavirus (yellow head virus, YH). The scale bar indicates evolutionary divergence corresponding to a mean of 0.1 amino acid replacements per site.

As found previously by Marra et al. [2] and Rota et al. [3], the root of SARS-CoV is closer to coronavirus group 2. All SARS-CoV strains are almost completely identical in sequence (∼99% DNA sequence homology) and therefore it is not possible to get any meaningful phylogeny within the SARS-CoV group. The avian infectious bronchitis virus (IBV), which causes respiratory disease in chickens, and the turkey virus (TCV), which causes enteric disease, are clustered together; their ancestor divides groups 1 and 2 of the mammal-infecting coronaviruses. The close clustering of the chicken and turkey viruses suggests that the difference between enteric and respiratory tropisms might require only a few amino acid changes. Experiments have shown that just two point mutations in the spike (S) glycoprotein can change porcine transmissible gastroenteritis coronavirus (TGEV), a mostly enteric virus that can kill piglets, into a non-deadly virus that excels at the respiratory route but replicates poorly in the gut 13, 14, 15. To infect enteric tract cells with TGEV, two different domains of the S protein of TGEV, mapping to between amino acids 522 and 744 and close to amino acid 219, are involved [13]. The first domain binds to aminopeptidase N (pAPN); many viruses use co-receptors, and it is probable that the second domain maps a co-receptor essential for the enteric tropism of TGEV 14, 15.

The clustering of murine hepatitis virus (MHV) and rat coronavirus (RtCoV) might reflect the relatively close proximity in which the hosts reside and perhaps the similarity of murine and rat target receptors. Note that MHV receptors, including CEACAM1, have recently been identified 16, 17. The clustering of the human OC43 (HCoV-OC43), bovine (BCoV) and porcine haemagglutinating encephalomyelitis (HEV) coronaviruses might reflect conditions contributing to cross-infection in farming. The Breda torovirus is enteric, whereas the YH okavirus infects gill tissue in prawns. This indicates that the switch between enteric and respiratory tropism is a general characteristic of the order Nidovirales.

Analysis of the spike protein

Coronaviruses attach to host cells through the S glycoprotein 15, 18. This protein is translated as a large polypeptide that is subsequently cleaved into a receptor-binding peripheral subunit (S1) that remains non-covalently associated with a fusion-inducing membrane-spanning (S2) fragment [18].

Studies have shown that the entry of the porcine coronavirus TGEV into cells is mediated by the interactions of S1 with pAPN, an ectoenzyme abundantly expressed at the apical membrane of enterocytes covering the villi of the small intestine [19]. The fact that the S protein mediates the first interaction of the virus with human cells suggests that it might represent an excellent target for effective anti-SARS-CoV drugs; unfortunately, structural information is only available for the 3CLpro proteinase, which is part of the coronavirus replication complex (PDB accessions 1Q2W, 1P9U, 1P9S and 1LVO) [20]. The S protein shows relatively high sequence homology within two of the major coronavirus groups (>60% within group 1; >38% within group 2), whereas the homology between groups is lower (15–21% between groups 1 and 2; 15–21% between SARS-CoV and group 1 or 2). Comparative sequence analysis between the SARS-CoV sequences and sequences from groups 1 and 2 of the coronavirus genus reveals three regions of varying sequence conservation: from amino acid positions 1 to ∼400, 401–680 and 681–1255. Interestingly, the first two regions correspond to the S1 fragment (we subsequently refer to these regions as S1A and S2A) and the third to the S2 fragment. Figure 2a shows a cartoon of the gene for the S protein. The region S1A is poorly conserved between groups 1 and 2 and SARS-CoV, whereas S1B is flanked by two very well conserved motifs making it easier to align the internal sequences. Notably, this region is homologous to a TGEV spike region that is reported to contain determinants for tissue tropism 13, 14, 15. The region S2 is more conserved than S1A and S1B and consequently the alignment is easily determined; Passml-TM, which uses the evolutionary relationships of the sequences analysed to improve its predictions of secondary structure and accessibility, predicts a transmembrane helix at location 1196–1218 of the S protein. The maximum likelihood trees of S1B and S2 are shown in Figures 2b and 2c, respectively. The length of the trees and the distances between sequences reflect the sequence homologies according to the model of evolution, explaining why the S2 tree is much shorter than S1B tree.

Figure 2

(a) Schematic diagram of the spike protein gene showing S1A, S1B (faster evolving) and S2 (conserved) regions. The inferred transmembrane region is represented in black. The heptad repeats (r) and cysteine-rich domain (c) are also shown. Maximum likelihood trees of the (b) S1B and (c) S2 regions of the spike protein. The scale bars indicate the mean numbers of amino acid replacements per site. Species used include 17 coronaviruses [from group 1: human coronavirus 229E (HCoV-229E), porcine epidemic diarrhea virus (PEDV), porcine transmissible gastroenteritis virus (TGEV), canine coronavirus (CCoV) and feline coronavirus (FCoV); from group 2: bovine coronavirus (BCoV), rat coronavirus (RtCoV), murine hepatitis virus (MHV), human coronavirus OC43 (HCoV-OC43) and porcine hemagglutinating encephalomyelitis virus (HEV); intermediate between groups 1 and 2: avian infectious bronchitis virus (IBV) and turkey coronavirus (TCV); and SARS-CoVs from strains Tor2, bj01, bj02, bj03 and CUHKW1], one torovirus (Breda) and one okavirus (yellow head virus, YH). Coloured areas indicate coronavirus groups 1 (blue) and 2 (yellow). (d) Predicted structure and mutation rate for the SARS-CoV S1B region. Secondary structure has been predicted (row Pred) in three classes: helix (H), sheet (E) and coil (C). Accessibility has been predicted (row Acc) in three classes: buried (b), exposed (e) and intermediate (i). Mutation rates (row Mut) are partitioned into eight classes (classes 1–8 have relative rates of evolution 0.10, 0.27, 0.44, 0.64, 0.88, 1.19, 1.65 and 2.83, respectively), inferred using the empirical Bayes method [37]. The sequence homologous to the region 522–744 of the TGEV spike protein is shown in red. Maximum likelihood trees and mutation rate analyses were computed using Passml-TM [7].

The phylogenetic tree produced from the analysis of the spike S2 fragment from several coronaviruses indicates that SARS-CoV is closer to group 2 of the coronavirus genus than to group 1. Phylogenetic analysis of the S1 fragment from several coronaviruses indicates that SARS-CoV has an even closer relationship to group 2 viruses than the previous analysis suggests. Phylogenetic analysis of a 300 amino acid region at the terminus of S1, which we denote S1B, indicates that SARS-CoV belongs to group 2 and is closely related to MHV and RtCoV. The differences between these results might be a result of recombination events involving SARS-CoV or convergent evolution, or they might also be caused simply by chance; however, it is apparent that SARS-CoV is most closely related to group 2 of the coronavirus genus. Because the S1 fragment of MHV binds to CEACAM1, we suggest that the S1B region of the SARS-CoV spike might bind to a human CEACAM1 receptor instead of pAPN. Holmes and collaborators, who have identified receptors for MHV (murine CEACAM1a), HCoV-229E (human pAPN) and feline coronaviruses (feline pAPN) 16, 17, are currently investigating this hypothesis (K. Holmes, pers. commun.). In support of the proposal that S1B is involved in receptor-binding is the fact that this region is homologous to a domain of the TGEV S protein, which is located between amino acids 522 and 744 and is also involved in receptor-binding 13, 16.

Notably, TGEV mutants that lack sialic acid-binding activity contain single point mutations in the S protein (Cys155Phe, Met195Val, Arg196Ser, Asp208Asn or Leu209Pro) 21, 22. Sialic acid-binding activity might help TGEV to resist detergent-like substances encountered during gastrointestinal passage and therefore facilitate infection of the intestinal epithelium [23]. We found that only Cys155 is conserved in SARS-CoV; this is in agreement with clinical findings that show that 20–50% of SARS patients present gastrointestinal symptoms [24]. The low conservation of the S1A region among coronavirus sequences suggests that once more strains of SARS-CoV or other closely related species are available it will become possible to use innovative comparative sequence analyses to examine positive selection that acts in this region [1].

Several important functional determinants have been discovered in fragment S2. It contains a cytoplasmic tail enriched in cysteine residues (1217–1236; Figure 2a); this is a common feature among coronaviruses and appears to be related to membrane fusion [25]. Several authors have discovered that SARS-CoV S2 contains two conserved regions of heptad repeats (913–1000 and 1151–1185; Figure 2a) 26, 27 (see also a press release by W.R. Gallaher and R.F. Garry http://www.virology.net/sars/s2model.html). These heptads suggest that SARS-CoV uses mechanisms to gain entry to a cell that are similar to those used by human immunodeficiency virus (HIV), the virus that causes AIDS (acquired immunodeficiency syndrome), orthomyxoviruses and paramyxoviruses, and also Ebola. It is known that peptides derived from these repeat regions in HIV and the paramyxoviruses can specifically inhibit virus entry and subsequently viral replication [26]. Currently, SARS treatment is modeled after the drug known as T20 (http://www.hivmedicine.com/textbook/drugs/t20.htm). This drug is a complex peptide that is difficult to manufacture, has a short half-life in the human body and must be injected. This suggests that other regions of the genome, such as S1B, which might be of some interest for drug design, should be described. We report in Figure 2d the consensus protein structure estimate of S1B obtained using Passml-TM, PHD [10] and PSI-PRED [11]. Passml-TM also estimates the distribution of mutation rates along the protein and the sitewise mutation rate (Figure 2d).

Mutational spectra of SARS-CoV genomes

The SARS-CoV genome is AT rich (59%). Asymmetries in strand composition can reveal mutation bias (for example, cytosine deamination) or selection [28]. We found that the GC [(G−C)/(G+C)=0.02] and AT skew [(A−T)/(A+T)=0.037] in the SARS-CoV genome are smaller than those of the HIV genome (GC skew is 0.15 and AT skew is 0.23), which has a double-stranded RNA genome with a similar AT content. This suggests the existence of some selection on G and C distribution along the sequence to control the types of RNA secondary structure that form [29]. CG is the only dinucleotide statistically under-represented {f(CG)/(f(C)f(G)=0.46, where f(CG) is the frequency of CG dinucleotides and significance is assessed according to Refs. 30, 31}. Because this depletion is also found in the HIV genome but not in that of Tobacco mosaic virus, it might occur as a result of mutational bias in vertebrate cells.

Comparison of the mutation patterns in the SARS-CoV genome sequences from 16 patients shows that a large number (38/84) of the base substitutions detected at 84 sites occur within or near to single base and dinucleotide repeat stretches. Despite the absence of a pairing rule, the ratio of rates of transition and transversion mutations is ∼2, as is often found in double-stranded DNA. The low GC and AT skews and the low number of mutations suggest that evolvability of SARS-CoV might be restricted by selective constraints acting on the RNA structure and packaging of the genome, and therefore it might also be restricted by the low fitness of its mutational neighbours. Sequence features that form stems and loops that are potentially involved in coronavirus genome packaging have been described 32, 33. Mutational neighbours with different fitness might explain why, although some RNA viruses evolve at high rates, some RNA viruses are highly stable 34, 35.

Interestingly, the ease of tropism switching as exemplified by the closeness between turkey and chicken coronaviruses (Figure 1) is favoured by the large number of viral particles in each host, by their mutation rates, by the large populations of hosts (birds and other species) and by the aerial mode of viral spread (for instance, through sneezing and faeces). These factors suggest that birds might act as powerful engines for virus evolution.

Conclusion

In this review we have highlighted that the S1 and S2 fragments of the SARS-CoV S protein have different mutational patterns. On the basis of phylogenetic evidence (Figure 2b) and the homology with the TGEV 522–744 region, it is suggested that a short region of the S1 fragment, which we denote S1B, located at positions 400–680, might be of particular interest to virologists, structural biologists and biotechnologists. The sitewise secondary structure, solvent accessibility and mutation rate of this region have been estimated; our work in progress includes further structural characterization and fold family determination. At the moment, SARS appears to be under control despite doctors having neither drugs nor a vaccine to protect against it. Because it could reappear in the future, research should proceed and hopefully our findings might assist in maintaining a feed-forward loop on SARS-CoV research between bioinformatics analysis and experimental work from microbiologists and virologists. As a final comment, it is notable that, in all the phylogenies, human coronaviruses HCoV-229 and HCoV-OC43 always cluster with porcine coronaviruses. Because Ericsson and collaborators [36] reported the identification of two homologous human proteins that act as receptors for porcine endogenous retrovirus, the benefits and risks of porcine–human xenotransplantation should be carefully balanced.