Full-Genome Deep Sequencing and Phylogenetic Analysis of Novel Human Betacoronavirus

A novel betacoronavirus associated with lethal respiratory and renal complications was recently identified in patients from several countries in the Middle East. We report the deep genome sequencing of the virus directly from a patient’s sputum sample. Our high-throughput sequencing yielded a substantial depth of genome sequence assembly and showed the minority viral variants in the specimen. Detailed phylogenetic analysis of the virus genome (England/Qatar/2012) revealed its close relationship to European bat coronaviruses circulating among the bat species of the Vespertilionidae family. Molecular clock analysis showed that the 2 human infections of this betacoronavirus in June 2012 (EMC/2012) and September 2012 (England/Qatar/2012) share a common virus ancestor most likely considerably before early 2012, suggesting the human diversity is the result of multiple zoonotic events.

A novel betacoronavirus associated with lethal respiratory and renal complications was recently identified in patients from several countries in the Middle East. We report the deep genome sequencing of the virus directly from a patient's sputum sample. Our high-throughput sequencing yielded a substantial depth of genome sequence assembly and showed the minority viral variants in the specimen. Detailed phylogenetic analysis of the virus genome (England/ Qatar/2012) revealed its close relationship to European bat coronaviruses circulating among the bat species of the Vespertilionidae family. Molecular clock analysis showed that the 2 human infections of this betacoronavirus in June 2012 (EMC/2012) and September 2012 (England/Qatar/2012) share a common virus ancestor most likely considerably before early 2012, suggesting the human diversity is the result of multiple zoonotic events.
T he ability of coronaviruses (CoVs) to infect multiple species and to rapidly change through recombination presents a continuing human health threat. The epidemic of severe acute respiratory syndrome (SARS) during 2003-2004 during which a CoV transmitted from bats to civet cats and then to humans demonstrated this potential (reviewed in [1,2]). Recently, a novel human betacoronavirus (betaCoV) was found to be associated with at least 13 human infections, 7 of which were fatal (3)(4)(5)(6)(7). One of the viruses, EMC/2012, has been sequenced, and its sequence similarity to several bat CoVs suggested an animal origin, but a definitive bat species of origin has not yet been identified (3). Additional genome sequences from this virus are needed to aid diagnostics, monitor population dynamics, identify the animal source, and characterize mechanisms of pathogenesis. The large size (30,000 nt) and high variability of CoV RNA genomes present a challenge for sequencing.
We describe a strategy for rapidly designing the primers necessary for reverse transcription and cDNA amplification of such diverse RNA viruses and report the full-genome determination of the novel CoV directly from patient sputum using next-generation short-read sequencing. Full genomes from 2 epidemiologically unlinked novel CoV infections separated in time by >2 months were analyzed to gain clues to 2 major questions: what are the precursors of the virus, and how long has the virus been circulating in its current form?

Primer Design
Fifteen betaCoV full genomes with the closest homology to EMC/2012 (GenBank accession no. JX869059) were analyzed to identify all possible primer-like sequences (melting temperature 58°C-60°C, guanine plus cytosine content 35%-60%) yielding 357,198 potential primers. The 30-kb JX869059 CoV genome was divided into fifteen 2.5-kb overlapping amplicons, and the 2 highest frequency sequences mapping in the 5′ and 3′ 250 bp of each amplicon were selected. Reverse complements of the 3′ mapping sequences were prepared, resulting in a set of 60 primers for the 15 amplicons (4 primers/amplicon). Three additional primers were added for the extreme termini of the genome. The algorithm also prepared a virtual PCR map of the predicted binding and amplicon size for the primer set on CoV genomes. A map of the primer mapping positions Full-Genome Deep Sequencing and Phylogenetic Analysis of Novel Human Betacoronavirus and the predicted PCR products using EMC/2012 as the target is shown in Figure 1, panel A. The primer sequences and details of their use are listed in the online Technical Appendix Table (wwwnc.cdc.gov/EID/article/19/5/13-0057-Techapp1.pdf).

Sample History
The patient, originally from Qatar, visited Saudi Arabia, where he became unwell with severe respiratory disease and renal failure and was transferred to a hospital in London, UK (4). Viral RNA was extracted from a sputum sample collected on September 19, 2012, sixteen days after symptom onset. Real-time PCR was performed by using an assay targeting the upE gene (6)

Reverse Transcription and PCR Amplification
Reverse transcription was performed at 50°C for 60 min by using the amplicon reverse primers. PCRs were then performed with forward and reverse primers for each amplicon (15 reactions) for 35 cycles (98°C, 10 sec; 54.5°C, 30 sec; 72°C, 2.5 min) with a final 10-min 72°C extension. A 3-µL aliquot of each reaction was analyzed by electrophoresis, and each reaction showed the expected 2-2.5-kb products ( Figure 1, panel B). A detailed protocol is available on request from the authors.

Sequencing
The PCR products of amplicons were combined into 4 pools of approximately equal molarity and processed into standard Illumina multiplex libraries (Illumina, San Diego, CA, USA) with each pool bearing a unique bar code sequence. Libraries were sequenced by Illumina MiSeq generating 150-bp paired end reads. The resulting Illumina read sets were processed to remove adaptor and primer sequences and quality controlled to ensure a median read Phred quality score of 30 by using QUASR (8). The 5′and 3′-most primers were derived directly from the first 27 nt or last 31 nt of the EMC/2012 genome. These primer sequences were not trimmed from the reads, which might mask sequence difference between EMC/2012 and England/Qatar/2012 at these positions.

Genome Assembly
Reference-based mapping and de novo assembly methods were applied to the data for assembly into viral genomes. Reference-based mapping was performed against the EMC/2012 genome by using the Burrows-Wheeler Aligner software package [9]). For de novo assembly, maximum contig lengths were obtained by using subsets of 30,000-60,000 reads. Therefore, random subsets of reads were extracted from the readset and assembled by using Velvet version 1.2.07 (10,11) and VelvetOptimiser (12). The de novo assembled sequences were used to confirm the validity of the referencebased sequence and showed that the de novo assembly and the reference-based mapping produced identical sequences, apart from several small gaps near the termini of the de novo assembled sequence (results not shown). The complete genome sequenced here is named England/Qatar/2012 and is available in GenBank (accession no. KC667074).

Sequence Alignment
The complete England/Qatar/2012 genome was combined with 46 previously published complete genomes of  the α-(group 1), β-(group 2), γ-(group 3), and δ-(group 4) CoVs. The betaCoV complete genomes encompassed the a, b, c, and d lineages, as well as the 2 genomes of the novel human betaCoV, EMC/2012 (GenBank accession no. JX869059) and England1 (GenBank accession no. KC164505), the latter from the same patient as studied here. The sequences were aligned by using MUSCLE (13) within MEGA5 (14). Poorly aligned regions and the genes absent in some CoV genotypes were excluded from the alignment, resulting in a virtual concatenation of the open reading frame (ORF) 1ab, S, E, M, and N genes. A second alignment using a 396-bp region of the RdRp was generated similarly. This shorter alignment included strains that are genetically close to England/Qatar/2012 but lack complete genome sequences in GenBank.

Phylogenetic Methods
Maximum-likelihood (ML) phylogenetic trees were inferred from the sequence alignments by using PhyML version 3.0 (15). Phylogenies inferred from nucleotide and amino acid sequences used the general-time reversible and Whelan-and-Goldman substitution models, respectively. A discrete-γ distribution of 4 rate categories (Γ 4 ) was used to model among-site heterogeneity. The robustness of the tree topology was assessed by bootstrap analysis of 1,000 pseudo-replicates of the sequences.

Molecular Clock Dating
The evolutionary time scale of 2 novel human beta-CoVs, EMC/2012 and England/Qatar/2012, was estimated by using a strict clock model under the Bayesian Markov Chain Monte Carlo framework in BEAST version 1.7.4 (16). The Hasegawa-Kishino-Yano nucleotide substitution model was used with a discrete-γ distribution of 4 rate categories (Γ 4 ) to enable among-site heterogeneity. A simple constant coalescent prior on the age of the divergence was used. With only 2 samples, estimating a rate of evolution was not possible, so we conditioned our date estimates on a range of fixed rates from 1.0 × 10 -4 to 5.0 × 10 -3 substitutions/site/year to enable comparison with plausible values. For each fixed rate value, a total of 10 7 Bayesian Markov Chain Monte Carlo states were computed, with the first 10% discarded as burn-in.

Comparison of England1, England/Qatar/2012, and EMC/2012
The index genome for this virus, EMC/2012, was originally obtained from a Saudi Arabian patient in July 2012 after 6 passages in cell culture (3). Our England/Qatar/2012 sequence was derived directly from the sputum of a Qatari patient receiving care in London in September 2012. A consensus genome from the same Qatari patient, but from a lower respiratory tract sample obtained later in the infection, is also available (England1, GenBank accession no. KC164505). The 3 genomes were aligned to determine sequence differences.
The England/Qatar/2012 and England1 genomes show only nucleotide differences at the genome termini, with the England1 genome lacking 46 nt at the 5′ end and 42 nt at the 3′ end ( Figure 2, panel A, Appendix, wwwnc.cdc. gov/EID/article/19/5/13-0057-F2.htm). The England1 genome was generated by using 80 PCR products sequenced by Sanger dideoxy methods (17), and the high level of sequence identity between the England/Qatar/2012 and England1 sequences strongly validates the novel rapid sequencing methods described here.
The England/Qatar/2012 consensus genome has 30,021 (99.67%) of 30,119 nucleotides identical to the EMC/2012 genome. In addition to the single nucleotide differences, EMC/2012 shows an insertion of 6 nt starting at position 29639 and a single A insertion at 30661 relative to the England/Qatar/2012 genome. The sequence differences are evenly spread across the genome (Figure 2, panel A, Appendix; Table 1).
Deep sequencing data enable the generation of a consensus sequence from the majority nucleotide at each genome position and the identification of nonconsensus nucleotides at each position. The England/Qatar/2012 sample was sequenced to a high level of coverage (mean coverage = 4,444) compared with the previously reported England1 genome (2-11-fold), enabling insights into the variation consequent on within-host evolution. Setting a conservative estimate for total sequencing errors at 1% enables nucleotide variants present at >1% frequency to be considered true variants in the virus genome ( Figure 2, panel B, Appendix). Variation is clearly detected across the virus genome with certain regions, such as nonstructural protein (NSP) NSP3, NSP5, NSP6, NSP12, and NSP14 showing increased levels of nonconsensus nucleotides ( Table 2).
The predicted ORFs in England/Qatar/2012 are shown in Figure 2 An assessment of specific ORFs of the 2 viral genomes was performed because evolution of specific proteins is a key determinant of host range and defines cross-species transmission events that may be recent for this virus. As expected from the overall close homology of the 2 viruses, several ORFs show little or no change between the 2 genomes (Table 1). However, several ORFs show higher amino acid differences, including the nucleocapsid ORF N, and ORFs 3, 4a, and 8b (Table 1).
Of the 16 predicted nonstructural proteins encoded by ORF 1a and ORF 1b, 11 show no change at the amino acid level (Table 1), whereas 5 (NSP2, NSP3, NSP4, NSP13, and NSP15) show >0.3% difference ( Table 1). As more sequences become available from this virus, such comparisons will yield clues about the adaptation to humans.

Close Phylogenetic Relationship with European Bat CoVs
A ML phylogenetic tree inferred from the whole genome alignment indicated that the 3 novel human beta-CoVs sequences (England1, England/Qatar/2012, and EMC/2012) clustered closely, forming a monophyletic lineage that falls into group 2c (Figure 3, panel A, Appendix, wwwnc.cdc.gov/EID/article/19/5/13-0057-F3.htm). This novel human betaCoV lineage shares common ancestries with other group 2c bat CoV variants, including HKU4 and HKU5 strains isolated in southern People's Republic of China ( Figure 3, panel A, Appendix); however, the genetic distances between them remain substantial (≈70% similarity in nucleotide level). The novel human betaCoVs have even less similarity with group 2a, 2b, and 2d CoVs (<60% nt). Phylogenies constructed from individual ORFs including ORF 1ab, S, E, M, and N demonstrated a largely consistent phylogenetic position of the novel human betaCoVs, except for a small discordance in the order of branching between novel human betaCoV, HKU4, and HKU5 lineages (Figure 3, panels B-F, Appendix). Such phylogenetic incongruence could result from several possible evolutionary features of the virus, including evolutionary rate variation, homoplasy, and recombination, as well as from uncertainty in the alignment of distant sequences. Additional related sequences are needed to clarify this issue.
Because more CoV strains are available in GenBank for partial genomic sequences, we performed extensive BLAST searches (www.ncbi.nlm.nih.gov/blast) and identified several additional European bat CoV sequences sharing higher nucleotide sequence similarity (82.0%-87.7%) with the novel human betaCoVs. Because only partial RdRp sequences were available for these European bat CoVs, another ML phylogeny was constructed from this short region to examine their evolutionary relationships (Figure 3, panel G, Appendix). As noted, P.pipi/VM314/2008/NLD (GenBank accession no. GQ259977) identified in a Pipistrellus pipistrellus bat from the Netherlands (18) shows the closest sequence similarity (87.7%) and phylogenetic relationship to the novel human betaCoVs. In addition, 2 CoV sequences, H.sav/J/Spain/2007 (GenBank accession no. HQ184059; from Hypsugo savii, also known as P. savii) and E.isa/M/Spain/2007 (HQ184062; from Eptesicus isabellinus), obtained from bats from Spain (19) are also closely related to the novel human betaCoVs (Figure 3, panel G, Appendix).

Timing the Origin of Zoonosis
The time to the most recent common ancestor (tMRCA) of EMC/2012 (isolated June 13, 2012) and England/Qatar/2012 (isolated September 19, 2012) was estimated with a strict molecular clock model. A plot of estimated tMRCA of the 2 human strains as a function of the fixed rate of molecular evolution shows that with the fastest measured rate of human CoVs (2 × 10 -3 substitutions/ site/year for SARS-CoV [20]) the tMRCA (i.e., October 2011; 95% highest posterior density: August-December 2011) is older than the start of 2012, but with slower rates it is calculated to be much older (Figure 4 [20,[21][22][23]). However, for evolutionary rates as high as those of human influenza A virus and HIV (e.g., 3-5 × 10 -3 substitutions/ site/year), the estimated tMRCA becomes compatible with the earliest disease report in Jordan in April 2012 (7).

Discussion
The reports of several human infections by similar strains of a novel betaCoV have raised global concern about a new SARS-like outbreak. Such worry is not misplaced, particularly when little is known about the origin and transmissibility of the virus and the natural history of the disease it causes. In this study, we presented a rapid deep sequencing method to obtain the complete genome sequence of the novel human betaCoV and its minor variants from a patient's sputum sample, which provides a useful tool to study the origin, evolution, and detection of this novel CoV. Our detailed phylogenetic analyses of the viral genomes also provide additional clues to the emergence and evolution of this virus.
Mammalian zoonoses are often identified as the origin of new human CoV infections including OC43, NL63, and SARS-CoVs (24)(25)(26)(27). In particular, bats, which maintain the largest diversity of CoVs (28), are thought to be a natural reservoir of the virus and thus a probable origin of the novel human betaCoV studied here. Indeed, the initial genetic study of the novel human betaCoV demonstrated its phylogenetic relatedness to the HKU5 and HKU4 bat CoVs in southern China (5). van Boheemen et al. found a bat CoV sequence from the Netherlands that clusters closely with EMC/2012 (3). In addition, most bat CoVs in the group 2c lineage, including the sequence from the Netherlands, were isolated from bat species of the Vespertilionidae family, whereas the SARS-like CoVs in group 2b lineage were mainly from bats from the Rhinolophidae family (29). In this study, we showed that 2 other CoVs from bats from Spain of Vespertilionidae again are also clustered near the novel human betaCoVs (Figure 3, panel G, Appendix), which further supports a European Vespertilionidae bat ancestry for this virus. However, such a genealogic link may be indirect and distant because a similar hypothesis for SARS-CoV has been under dispute: whether the CoV jumped from bats directly to humans or through civets or even through some other animals as intermediate hosts (30). Therefore, like SARS-CoV, the high dissimilarity between the bat CoVs and novel human betaCoVs makes it difficult to rule out the presence of other intermediate hosts transferring the virus from bats to humans and to confirm the geographic origin of the direct predecessor. With the current data, Europe and the Middle East would be plausible regions for more intensive pilot surveillance studies on bats and other animal reservoirs for this virus.
In the interest of public health, it is critical to determine whether these CoV infections in humans are the consequence of a single zoonotic event followed by ongoing human-to-human transmissions or whether the 3 geographic sites of infection (Jordan, Saudi Arabia, and Qatar) represent independent transmissions from a common nonhuman reservoir. The large genetic diversity of CoV maintained in animal reservoirs suggests that viruses that independently moved to humans from animals at different times and places are likely to be reasonably dissimilar in their genomes, possibly making the multiple transmission events model less likely. Further information is needed to confirm this point because the currently available data are limited. If we calibrate our molecular clock analysis using the evolutionary rate of Zhao et al. (20) estimated for SARS-CoV, we dated the tMRCA of EMC/2012 and England/Qatar/2012 viruses to early 2011. Therefore, if both sequenced viruses and the other cases descended from a single zoonotic event, then this tMRCA suggests that the novel virus has been circulating in human population for >1 year without detection and would suggest most infections were mild or asymptomatic. The rate would have to be considerably faster, of a magnitude observed for human influenza A virus, for the tMRCA to be compatible with the earliest known cases in April 2012. Perhaps more probable, therefore, is that the 13 known cases of this disease represent >1 independent zoonotic transmission from an unknown source. Viral sequence data from other patients infected with this novel human betaCoV will help to more accurately estimate the estimate a genomic evolutionary rate specific to this virus, which will then yield a tMRCA estimate closer to the actual time.
A major determinant of SARS pathology is the distribution of the host receptor that the virus has evolved to use for entry. A detailed analysis of EMC/2012 receptor use in cell culture (31) revealed that EMC/2012 was capable of infecting a range of mammalian cells, including human, pig, and bat cells, and the virus entry was independent of the ACE2 receptor. Comparing the S genes of England/Qatar/2012 predicts only 2 aa changes (in a 1,367-residue protein). Thus, this broad infection capability of EMC/2012 is likely to be valid for England/Qatar/2012 and would imply a higher possibility for the existence of other intermediate hosts transferring the virus from bats to humans.
Precise identification of the origin of this virus, defining its mode of evolution, and determining the mechanisms of viral pathogenesis will require full-genome sequences from all cases of human infection and substantially more sampling and sequencing from Vespertilionidae bats and other related animals. The sequencing method reported here markedly shortens the time required to process the clinical sample to genome assembly to 1 week and will provide a useful tool to study this novel virus.
Note added in proof. The current number of confirmed novel CoV cases is 17, including 11 deaths (www. who.int/csr/don/2013_02_21/en/index.html). An analysis of the 4 novel CoV genomes now available suggests an evolution rate of 1.77 × 10 -3 substitutions/site/year, which would give a tMRCA of mid-2011. The analysis is online (http://epidemic.bio.ed.ac.uk/ coronavirus_analysis).