Introduction
Coronaviruses (CoVs) infect and cause disease in a wide variety of species, including bats, birds, cats, dogs, pigs, mice, horses, whales, and humans (
1,
2). Recent studies suggest that bats act as a natural reservoir for coronaviruses (
3–8). Coronaviruses may cause respiratory, enteric, hepatic, or neurological diseases with highly variable severity in their hosts. Until 2003, only two coronaviruses were known to infect humans. Human coronaviruses (HCoVs) HCoV-229E and HCoV-OC43 were identified in the 1960s as the causative agents of—generally mild—respiratory illnesses (
9,
10). In 2002 to 2003, a previously unknown coronavirus—severe acute respiratory syndrome coronavirus (SARS-CoV)—caused a widespread outbreak of respiratory disease in humans, resulting in approximately 800 deaths and affecting around 30 countries (
11–14). As a consequence of the renewed interest in coronaviruses after the SARS outbreak, two additional human coronaviruses were discovered after 2003: HCoV-NL63 in 2004 (
15,
16) and HCoV-HKU1 in 2005 (
17). A recent analysis of a large collection of human nasopharyngeal specimens using a
Coronaviridae-wide primer set suggested that HCoV-229E, -OC43, -NL63, and -HKU1 are the only coronaviruses circulating in the human population (
18).
Coronaviruses are enveloped single-stranded positive-sense RNA viruses with genomes of 25 to 32 kb, and the group includes the largest known genomes among the RNA viruses (
1,
19). The coronaviruses form a subfamily (
Coronavirinae) within the family
Coronaviridae of the order
Nidovirales. The International Committee on Taxonomy of Viruses (ICTV) has recognized four genera within the
Coronavirinae subfamily:
Alphacoronavirus,
Betacoronavirus, and
Gammacoronavirus, which were previously referred to as coronavirus groups 1, 2, and 3, and
Deltacoronavirus (
20). Coronaviruses are assigned to a genus on the basis of rooted phylogeny and calculation of pairwise evolutionary distances for seven highly conserved domains in the replicase polyprotein (
1,
21) (C. Lauber and A. E. Gorbalenya, unpublished data). HCoV-229E and HCoV-NL63 are viruses belonging to the genus
Alphacoronavirus (
1). Four monophyletic lineages (A through D) with no formal taxonomic standing, some of them encompassing multiple virus species, are commonly recognized within the genus
Betacoronavirus. Lineage A includes HCoV-OC43 and HCoV-HKU1 and lineage B SARS-CoV, all of which belong to different species. Lineages C and D include viruses detected only in bats, such as
Rousettus bat coronavirus HKU9 (BtCoV-HKU9) (lineage D),
Tylonycteris bat coronavirus HKU4 (BtCoV-HKU4), and
Pipistrellus bat coronavirus HKU5 (BtCoV-HKU5) (both lineage C) (
1). The genetic diversity of coronaviruses is likely facilitated by a high frequency of RNA recombination and the ability of their unusually large RNA genomes to both gain and lose domains (
1,
22,
23). These factors are believed to have promoted the emergence of viruses with novel traits that are able to adapt to new hosts and ecological niches, sometimes causing zoonotic events.
For the present study, we report and analyze the complete genome sequence of the recently identified HCoV-EMC/2012, which was isolated from the sputum of a 60-year-old man who died in a hospital in Jeddah, Saudi Arabia, after developing acute respiratory distress syndrome (ARDS) and multiple organ dysfunction syndrome (MODS) in June 2012 (
24). This virus appears to be closely related to the HCoV detected in a second patient who was transported from a hospital in Qatar to a hospital in London, 3 months after hospitalization of the first patient (
25). These two cases of human infection with very similar or identical coronaviruses alarmed health authorities globally, as it was a reminder of the potential threat of coronaviruses to human health that was first highlighted by the SARS outbreak of 2003 (
25). The sequence analysis of a small reverse transcription-PCR (RT-PCR) fragment that was first amplified from the HCoV-EMC/2012 genome revealed the highest similarity to two betacoronaviruses circulating in bats, BtCoV-HKU4 and -HKU5 (
24). Here we present the complete genome sequence of the newly isolated HCoV-EMC/2012, accompanied by a detailed annotation of its genome organization and expression strategy. Furthermore, comparative genomic analysis and state-of-the-art classification and phylogenetic analyses were applied to determine the position of the novel agent with respect to previously characterized coronaviruses. We conclude that the HCoV-EMC/2012 genome organization and expression indeed most closely resemble those of BtCoV-HKU4 and -HKU5. However, based on our analysis and in line with the ICTV guidelines for the demarcation of coronavirus species, HCoV-EMC/2012 clearly qualifies to be recognized as the prototype of a novel species, which would thus constitute the first human coronavirus in lineage C of the genus
Betacoronavirus.
DISCUSSION
Coronaviruses have been known for quite some time as viruses that cause a variety of diseases in humans and animals (
32,
33). The discovery of a coronavirus as the causative agent of SARS revived the interest in coronaviruses and resulted in a rapid increase of the number of identified coronaviruses, as well as of the number of full coronavirus genome sequences. Until this study, lineage C of the genus
Betacoronavirus (formerly known as subgroup 2c) included virus isolates from bats. Here, we determined and analyzed the complete genome sequence of a previously unknown lineage C betacoronavirus that was isolated from the sputum of a 60-year-old male suffering from acute pneumonia and renal failure in the Kingdom of Saudi Arabia whose death was probably a consequence of this infection (
24).
The sequencing of the full HCoV-EMC/2012 genome was greatly facilitated by the advent of high-throughput techniques. Using an optimized random amplification deep-sequencing approach, approximately 90% of the virus genome was covered with high accuracy in a single run. Using the data from this first run, primers could be designed to perform conventional Sanger sequencing for confirmation. This combination of techniques allowed the determination of the complete virus genome within a few days, without a requirement for prior knowledge of the virus genome under investigation. The error rate in 454 deep sequencing was generally higher than in Sanger sequencing, but the high coverage across the HCoV-EMC/2012 virus genome (up to 5,697 reads per nucleotide position) corrected for most of the incorrect base callings. The sequence obtained using the 454 platform aligned almost perfectly with that obtained by Sanger sequencing, with the exception of two nucleotide positions. The deep-sequencing data revealed variation at position 11623 (U or G), with G occurring in 44% of the reads, suggesting that ORF1a-encoded residue 3782 can be either valine (codon GUC) or glycine (codon GGC). The valine codon was the more abundant codon at this position in HCoV-EMC/2012, and valine is also present in most other betacoronaviruses. At position 27162, both G and A were detected in different runs, with an A in 45% of the reads. This G-to-A substitution introduces a premature stop codon (UGG to UAG) in ORF5. The virus stock that we sequenced was derived from passage of the virus from a sputum specimen six times in Vero cell culture. Hence, the observed sequence variants may reflect either natural heterogeneity or emerging genetic changes that occurred during virus passage in cell culture. Additional HCoV-EMC/2012 virus isolates or patient materials are currently not available to verify these genome sequence ambiguities at positions 11623 and 27162.
Adaptation to cell culture leading to a loss of functionality of genes, and in particular in relation to the so-called “accessory protein genes,” has previously been described for a variety of coronaviruses, including SARS-CoV (
2,
34). These genes, like ORF3 through ORF5 of HCoV-EMC/2012, are dispersed between the structural protein genes (
35) and in some cases may even overlap such a gene, as in the case of the ORF overlapping the N protein gene in betacoronaviruses (
Fig. 1A) (23, 36). The origin of most accessory protein genes remains unclear, although for some, acquisition by recombination with cellular or heterologous viral sequences seems plausible (
37,
38). Accessory gene functions have been probed by reverse genetics (knockout mutants) for a variety of coronaviruses, including SARS coronavirus (
39), which established that they are not essential for replication in cell culture systems. In animal models, on the other hand, profound effects on pathogenesis after the inactivation (or transfer to a heterologous coronavirus) of accessory protein genes have been previously described (
40–42). In some cases, accessory gene products have been implicated in immune evasion, e.g., by interfering with cellular innate immune signaling (
43).
The apparent absence of selection pressure on coronavirus accessory protein genes during cell culture passage may explain the relatively high frequency with which loss of functionality appears to occur. The detection of an internal termination codon in part of the HCoV-EMC/2012 ORF5 sequences (45% of the reads) may constitute another example of such an event, which would lead to the truncation of the ORF5 protein after 107 amino acids. This would resemble a 29-nt deletion that occurred in the SARS-CoV genome, which resulted in the truncation of ORF8 (
34,
44), and a 45-nt in-frame deletion in ORF7b of the same virus that emerged upon cell culture passage (
23).
Our analysis identified a potential ORF underlying the N protein gene (ORF8a), which is a common feature in betacoronaviruses. This ORF was not previously described for BtCoV-HKU4 and BtCoV-HKU5 (
8) but is conserved in the genome sequences of both viruses (see
Fig. S1J in the supplemental material). Remarkably, in HCoV-EMC/2012, both the 5′ and 3′ parts of the ORF appear to have been truncated. In BtCoV-HKU4 and BtCoV-HKU5, the ORF8b AUG codon would be the second AUG on sg mRNA8, making leaky ribosomal scanning a likely mechanism for translation initiation. In HCoV-EMC/2012, however, this AUG codon (positions 28606 to 28608) seems to have been mutated to AUA. Conservation of the sequence immediately downstream of this position, which is now formally upstream of ORF8b in HCoV-EMC/2012, was observed with BtCoV-HKU4 and BtCoV-HKU5, suggesting that the putative loss of this AUG codon may also have been a relatively recent event. In the 3′ part of ORF8b, sequence alignment of HCoV-EMC/2012 with BtCoV-HKU4 and BtCoV-HKU5 suggests that the former acquired a premature termination codon at positions 29099 to 29101 (UAA). Although we cannot at present assess the timing of these events in HCoV-EMC/2012 evolution, due to the lack of alternative samples for this species, the presumed loss of ORF8b functionality may also be a consequence of virus passage in cell culture.
To classify newly identified coronaviruses as the prototype of a novel virus species, it is required that the amino acid sequence identity in the conserved replicase domains in all intervirus pairwise comparisons is below the 90% threshold (
1). Here, we propose HCoV-EMC/2012 to represent a novel species of the betacoronavirus genus, since the amino acid sequence identities between HCoV-EMC/2012 and its closest relatives BtCoV-HKU4 and BtCoV-HKU5 in the seven conserved domains of ORF1ab were 75% and 77%, respectively. These viruses were originally detected in Asia in lesser bamboo bats (
Tylonycteris pachypus) and Japanese house bats (
Pipistrellus abramus), respectively (
8). This proposed classification will remain provisional until approved by ICTV.
The ICTV guidelines for coronavirus species demarcation require the availability of a (nearly) complete genome sequence prior to virus classification. However, there is considerable correlation between the results based on full-genome sequence analysis and those determined using the most conserved part of the ORF1b-encoded RdRp domain, which is commonly used in screening for new coronaviruses. In 2010, this partial sequence was reported for a betacoronavirus (VM314/2008) that was isolated 2 years earlier from a
Pipistrellus pipistrellus bat in The Netherlands. This virus was provisionally classified a betacoronavirus based on a 332-nt fragment from the RdRp-encoding domain of ORF1b (
31), which shares 88% nucleotide sequence identity with HCoV-EMC/2012, the highest identity with any coronavirus sequence available in the public domain. Although this high similarity is not sufficient to resolve the taxonomic relation between HCoV-EMC/2012 and isolate VM314/2008, it suggests that they may both belong to the same coronavirus species. Establishing the genome sequence of VM314/2008, or closely related viruses, is urgently required to verify this hypothesis. Based on the genetic relation between HCoV-EMC/2012 and bat coronaviruses, it is tempting to speculate that HCoV-EMC/2012 emerged from bats—either directly or via an intermediate animal host, possibly
Pipistrellus bats. This bat species is known to be present in the Kingdom of Saudi Arabia and neighboring countries.
Although most infections of human coronaviruses are relatively mild, the infection by HCoV-EMC/2012 with fatal outcome, and a similar severe case of an infection with a closely related coronavirus in London (
25), is a reminder that certain coronaviruses may cause severe and sometimes fatal infections in humans. It is important to develop an animal model that can be used to fulfill Koch’s postulates for the novel virus, by demonstrating that the isolated virus can indeed cause the observed disease. The availability of the HCoV-EMC/2012 genome sequence will facilitate the development of a variety of diagnostic assays that can be used to study the prevalence and clinical impact of HCoV-EMC/2012 infections in humans. The first generation of assays for this purpose has recently been described (
45). We anticipate that the availability of this full-length virus genome sequence will be valuable for the development of additional applied and fundamental research.
MATERIALS AND METHODS
Virus propagation.
Patient material had been subjected to passage in Vero cells four times in the Dr. Soliman Fakeeh Hospital, Jeddah, Saudi Arabia. Subsequently, in the Erasmus Medical Center, Rotterdam, The Netherlands, LLC-MK2 cells were inoculated with HCoV-EMC/2012 in minimal essential medium (MEM-Eagle) with Earle’s salts (BioWhittaker, Verviers, Belgium), supplemented with 2% serum, 100 U/ml penicillin, 100 mg/ml streptomycin, and 2 mM glutamine. Vero cells were inoculated with virus in Dulbecco’s modified Eagle medium (BioWhittaker) supplemented with 1% serum, 100 U/ml penicillin, 100 mg/ml streptomycin, and 2 mM glutamine. After inoculation, the cultures were incubated at 37°C in a CO2 incubator and checked daily for cytopathic changes. Three days after inoculation, supernatant from Vero cells was collected and used for virus genome characterization.
Arbitrarily primed PCR and virus genome sequencing.
To characterize the viral genome, we used a random amplification deep-sequencing approach as the first step. Virus-containing supernatant was centrifuged for 10 min at 3,000 rpm to remove cellular debris. This supernatant was then filtered through a 0.45-µm-pore-size centrifugal filter unit (Millipore, Amsterdam, The Netherlands) to minimize bacterial contamination. Omnicleave endonuclease (Epicenter Biotechnologies, Madison, WI) was used to remove any free DNA and RNA, according to the manufacturer’s protocol. Subsequently, viral RNA was extracted from the purified, infected cell culture supernatant using a High Pure RNA isolation kit (Roche Diagnostics, Almere, The Netherlands). To remove contaminating mammalian rRNA, a Ribo-Zero RZH110424 rRNA removal kit (Epicenter Biotechnologies, Madison, WI) was used according to the manufacturer’s protocol. RNA was reverse transcribed using circular permuted primers (
46) that were extended with random hexamer sequences, namely,
CCCACCACCAGAGAGAAAN(6),
ACCAGAGAGAAACCCACCN(6),
GAGAAACCCACCACCAGAN(6),
GGAGGCAAGCGAACGCAAN(6),
AAGCGAACGCAAGGAGGCN(6), and
ACGCAAGGAGGCAAGCGAN(6). Per reaction, reverse transcription mixtures contained 6 µl RNA, 1 µl primer (20 pmol), 0.5 µl (20 U) RNase inhibitor (Promega, Leiden, The Netherlands), 1 µl (10 mM each) deoxynucleoside triphosphates (Roche), and 5 µl water. After a 5 min incubation at 65°C for optimal primer hybridization to the template, 4 µl (10×) first-strand buffer, 1 µl (200 U/µl) SuperScript III reverse transcriptase (Invitrogen, Bleiswijk, The Netherlands), 1 µl (0.1 M) dithiothreitol (DTT), and 0.5 µl (20 U) RNase inhibitor (Promega) were added to the mixture in a 20-µl volume. To obtain cDNA, the reverse transcription mixture was sequentially incubated at 25°C for 5 min and at 42°C for 1 h. After 3 min at 95°C and 2 min on ice, 1 µl Klenow DNA polymerase (5 U) (New England BioLabs Inc., Ipswich, MA) was added and the mixture was sequentially incubated at 25°C for 5 min, 37°C for 1 h, and 75°C for 20 min to obtain double-stranded cDNA. The cDNA was purified using a MinElute PCR purification kit (Qiagen, Venlo, The Netherlands) according to the instructions of the manufacturer. To amplify the purified cDNA, a PCR with the individual circular permuted primers without the random hexamer was performed. The PCR mixture contained 2 µl primer (40 pmol), 2 µl purified cDNA, 1.25 µl (10 mM each) deoxynucleoside triphosphate (Roche), 5 µl (10×) PfuUltra II Rxn buffer, and 1 µl (2.5 U) PfuUltra II DNA polymerase (Stratagene, Amsterdam, The Netherlands). Water was added to reach a final volume of 50 µl. The PCR mixture was incubated at 95°C for 2 min and then for 40 cycles of 95°C for 20 s, 56°C for 1 min, and 72°C for 2 min, followed by a final extension at 72°C for 10 min. Fragments were purified using a MinElute PCR purification kit (Qiagen) according to the instructions of the manufacturer.
Amplified fragments were sequenced using a 454/Roche GS Junior sequencing platform. A fragment library was created according to the manufacturer’s protocol without DNA fragmentation (GS FLX Titanium rapid library preparation; Roche), selecting for fragments larger than 100 bp. The emulsion-based PCR (emPCR) (amplification method Lib-L) and GS Junior sequencing run were performed according to the instructions of the manufacturer (Roche). The sequence reads were trimmed at 30 nt from the 3′ and 5′ ends to remove all primer sequences. Sequence reads were assembled into contigs using CLC Genomics 5.5.1 software (CLC Bio, Aarhus, Denmark). Using this deep-sequencing approach, approximately 90% of the virus genome sequence was obtained.
As a second step, specific primers were designed to amplify overlapping fragments of approximately 800 bp by RT-PCR. These PCR products were purified from agarose gels and sequenced using a BigDye Terminator v3.1 cycle sequencing kit (Applied Biosystems, Nieuwerkerk a/d IJssel, The Netherlands) and a 3130XL genetic analyzer (Applied Biosystems), according to the instructions of the manufacturers. The genomic 5′- and 3′-terminal sequences were determined using a FirstChoice RLM-RACE kit (Ambion, Bleiswijk, The Netherlands).
Virus classification.
Newly identified members of the family
Coronaviridae are generally assigned by the ICTV to a subfamily and genus on the basis of rooted phylogeny and calculation of pairwise evolutionary distances for seven replicase polyprotein domains (
1): the ADP-ribose 1″-phosphatase (ADRP) in nsp3, the coronavirus 3C-like (3 CL) protease (3CLpro, or “main protease”) in nsp5, the RNA-dependent RNA polymerase (RdRp) in nsp12, the helicase (Hel) in nsp13, the exoribonuclease (ExoN) in nsp14, the nidoviral endoribonuclease specific for U (NendoU) in nsp15, and the ribose-2′-
O-methyltransferase (O-MT) in nsp16. Amino acid sequence alignments were generated for each of these domains using ClustalW within the BioEdit (version 7.0.5.3) (
47) program and concatenated, after which the sequence identity of HCoV-EMC/2012 with closely related strains was calculated. For this purpose, the full genomes of 9 strains, derived from 3 species, belonging to
Betacoronavirus lineage C were available.
To support virus classification, protein-based phylogenetic trees were generated. Multiple amino acid alignments, including sequences of HCoV-EMC/2012 and one representative of each of the 20 recognized species of the subfamily
Coronavirinae, were produced for the following proteins, using the Viralis platform (
48) followed by manual correction: ADRP, the N-terminal part of PLP2, TM1, Y domain, nsp4 to nsp16, and the C-terminal part of the spike (S) protein (S2), envelope (E) protein, membrane (M) protein, and nucleocapsid (N) protein. From each protein alignment, the most informative blocks (
49) were extracted using the BAGG program (
50), and only these strongly conserved alignment regions were used for further analyses. Two concatenated alignments were used. The first included replicase pp1ab protein regions (4,110 aa positions, gap content of 0.9%), and the second included regions in the C-terminal domain of the S protein (S2) and the E, M, and N proteins (1,127 aa positions, gap content of 3.9%). ProtTest version 3.2 (
51) was used to select the best-fitting model of protein evolution. For both datasets, the LG model with rate heterogeneity (4 categories) ranked top among 112 models tested, with a relative weight of 0.98 under the Bayesian information criterion (BIC) and 0.74 under the corrected Akaike information criterion (AICc) for the first data set and 0.97/0.75 (BIC/AICc) for the second data set. Hence, this model was applied for maximum likelihood phylogeny reconstruction using PhyML version 3.0 (
52).
Phylogenetic reconstruction.
Nucleotide sequences were aligned using the ClustalW software running within the BioEdit (version 7.0.5.3) (
47) program and MAFFT version 6 (
53). Maximum likelihood phylogenetic trees with 100 bootstrap replicates were estimated under the general time-reversible model (GTR) + I + Γ4 and the transversion model (TVM) + I + Γ4 (determined by ModelTest [54]), using PhyML 3.0 software (
52). For both the 332-nt ORF1ab alignment that included isolate VM314/2008 and the alignment of the complete ORF1ab, the GTR + I + Γ4 model ranked top among 65 models tested, with relative weights of 0.8185 and 1.000 under AIC, respectively.
Nucleotide sequence accession number.
The final HCoV-EMC/2012 consensus sequence was submitted to GenBank under accession number JX869059.