Discovery and Evolution of a Divergent Coronavirus in the Plateau Pika From China That Extends the Host Range of Alphacoronaviruses

Samples

In July 2019, a total of 157 plateau pikas were captured from four sites (3,890–3,970 meters above sea level) in mountainous regions of Yushu Tibetan Autonomous Prefecture in Qinghai Province of China. The distance between sample collection sites is about 50 kilometers. All captured animals were clinically healthy and identified through morphological examination. Pikas were trapped and euthanized, and their respiratory tracts and intestinal contents were removed and preserved in 5-ml tubes directly. All sampling work was approved by the ethics committee of the National Institute for Communicable Disease Control and Prevention of China CDC (ICDC-2019012), and conducted by Yushu Prefecture Center for Disease Control and Prevention as part of plague surveillance.

RNA Extraction and Coronavirus Screening

Viral RNA of each fecal sample was extracted using QIAamp Viral RNA Mini kit (QIAGEN) and resuspended in 50μl of DNase-free, RNase-free water. Conserved primers were used to amplify the 440-bp fragment of RNA-dependent RNA polymerase (RdRp) of coronaviruses (Woo et al., 2005). Reverse transcription-polymerase chain reaction (RT-PCR) was performed by PrimeScript™ One Step RT-PCR Kit Ver.2 (Takara) with 50°C for 30min and 33cycles of 94°C for 30s, 48°C for 30s, 72°C for 30s, and a final extension at 72°C for 10min. The PCR products were gel-purified before conducting Sanger sequencing.

Complete Genome Sequencing

The total RNA of a coronavirus positive sample was firstly rRNA-removed using the Ribo-Zero Gold rRNA Removal Kit. The TruSeq stranded total RNA library prep gold kit was used for library construction according to instructions. Paired reads (2×150bp) were obtained using the Illumina HiSeq2000 platform, which were assembled into contigs using Trinity v2.8 (Grabherr et al., 2011). The assembled contigs were annotated in NCBI non-redundant protein database using Diamond (Buchfink et al., 2015). To confirm the assembled genome, clean data were mapped back to the genome using Bowtie 2 (Langmead and Salzberg, 2012) and assembled using SPAdes (Bankevich et al., 2012). Moreover, several primers were designed to amplify the NS2, S, and N genes followed by Sanger sequencing (Supplementary Table S2).

Genome and Phylogenetic Analyses

The deduced amino acid and nucleotide sequences of the open reading frames (ORFs) were found using ORF finder with default parameters and compared to those of other coronaviruses (Gao et al., 2003). The proteinase cleavage sites of CoVs were predicted using an online tool (Gao et al., 2003). The amino acid and nucleotide identities were calculated using BioAider version 1.334 (Zhou et al., 2020). Protein families and homology searches were performed by InterPro (Blum et al., 2021) and HMMER (Potter et al., 2018), respectively. Transmembrane helices in proteins were predicted by TMHMM with a hidden Markov model (Krogh et al., 2001). N-linked glycosylation and O-GalNAc (mucin type) glycosylation sites were predicated using web server.¹ The amino acid sequences were aligned using the MAFFT program (Katoh and Standley, 2013). Aligned regions were polished using Gblocks (Talavera and Castresana, 2007). Phylogenetic analyses were constructed using the maximum-likelihood method by PhyML 3.0 with IBV as the outgroup and bootstrap values of 1,000 (Guindon et al., 2010). Based on AIC criteria, the substitution models were selected using ModelFinder (Kalyaanamoorthy et al., 2017).

Estimation of Divergence Dates

The complete RdRp and M genes of PPCoV and related coronaviruses were acquired and aligned using the MAFFT program (Katoh and Standley, 2013). Sampling times were collected, and the divergence was estimated using a Bayesian Markov chain Monte Carlo (MCMC) method implemented in BEAST v 1.10 (Lau et al., 2012a). Analyses were performed using the SRD06 substitution model with uncorrelated relaxed clock. The MCMC run was 1×10⁸ steps long, with sampling every 1,000 steps. The mean time of the most recent common ancestor (tMRCA) and the highest posterior density regions at 95% (HPDs) were calculated. All the ESS values of statistic parameters were greater than 200 and visualized using Tracer v1.7.1 (Suchard et al., 2001). The trees were summarized in a target tree (a 10% burn-in) by the Tree Annotator program included in the BEAST package.

Virus-Host Co-divergence and Recombination Detection

The host species were collected and used to infer phylogenetic trees of alphacoronavirus mammal hosts using TimeTree (Kumar et al., 2017). A phylogenetic virus tree was built based on polyproteins as described above. A virus-host co-divergence tree was constructed using Jane software package (Conow et al., 2010) with generations set at 100 and population size set at 100. Event costs were set at 0 for co-divergence, 1 for duplications, 1 for host-switching, and 1 for loss. Recombination analysis was conducted using both Simplot version 3.5.1 with setting parameters of the F84 model, window size 1,000bp and step 200bp (Lau et al., 2012b) and Recombination Detection Program (RDP) version 4.97 with default parameters (Martin and Rybicki, 2000).

Novel Alphacoronavirus in Plateau Pikas

A total of 157 plateau pikas were obtained from four sites in Qinghai-Tibet plateau (Qinghai Province) in the northwest of China. RT-PCR targeting a 440bp fragment of the RdRp gene was conducted to screen potentially novel coronaviruses (CoVs). Viral RNA was positive in seven (4.5%) of 157 intestinal samples, with only positive in two respiratory samples of these seven pikas. All seven positive intestinal samples were collected from the same site, with overall detection rate of 13.7% at this site. Sequences of the screening PCR products revealed that they had <85.8% nucleotide identities to the corresponding sequences of known alphacoronaviruses, with 100% nucleotide identity to each other. Primer for multiple genome sites including NS2, S, and N genes were used to confirm the genome of all the positive samples, showing 100% nucleotide similarity among them. These results indicated that a new CoV (sharing <85.8% nucleotide identities to known alphacoronaviruses based on the 440bp fragment) circulates in plateau pikas of this site. No obvious diseases were observed in plateau pikas infected with the new CoV. To further confirm the pathogenesis, sections of respiratory tracts of plateau pikas that were positive for CoV were checked by staining with hematoxylin–eosin (HE). No inflammation was observed in respiratory samples (Supplementary Figure S1).

Complete Genome Characterization

Since partial RdRp sequences indicated that the CoV may represent a new member of the genus Alphacoronavirus, we selected one (designated P83) of the CoV positive samples for high-throughput sequencing. A total of 144,123,796 reads were obtained. To verify the assembled viral genome, a total of 154,419 reads were mapped back to the coronavirus genome. The genome was temporarily named plateau pika coronavirus (PPCoV) P83 (GenBank accession number MZ577265) with a size of 28,312 nucleotides (without 5' and 3' rapid amplification) and 35.8% G+C content. The genome organization of PPCoV was similar to that of other related alphacoronaviruses, with classical gene order of 5'-replicase ORF1ab-spike (S)-envelope (E)-membrane (M)-nucleocapsid (N)-3' (Figure 1; Table 1). A putative transcription regulatory sequence (TRS) motif of 5'-AACUAA-3' was located upstream in majority of the genes (Table 1), which was similar to that in other alphacoronaviruses. Some variations of TRS motif were observed in the PPCoV genome for NS2a (5'-TACUUUAA-3'), E (5'-UACUAA-3'), and NS7a (5'-AAGUAA-3') genes (Table 1), but no TRS motif was observed for NS9 gene. The putative nonstructural proteins (nsp) from ORF1ab with their corresponding cleavage sites are listed in Supplementary Table S1. The lengths of nsp2 and nsp4 in PPCoV differed from those of related alphacoronaviruses (Supplementary Table S1).

The genome of PPCoV showed 69.2 to 71.4% nucleotide sequence identities to Lucheng Rn rat coronaviruses (LRNVs), and <50.7% to other alphacoronaviruses (Table 2). The ORF1ab of PPCoV shared 75.1 to 76.2% amino acid (aa) identities with the ORF1ab of LRNVs, and <48.6% with that of other CoVs. The replicase domains of 3C-like protease (3CLpro), 3'-to-5' exonuclease (ExoN), nidoviral uridylate-specific endoribonuclease (NendoU) possessed <62.9%, <86.2%, and <85.0% aa identities, respectively, to the corresponding domains of other CoVs (Table 2).

NS2a gene between ORF1ab and S genes was minimally present in rat coronaviruses with exception of alphacoronavirus UKRn3, LRNV Lucheng-19, RtMruf-CoV-1/JL2014, and RtRl-CoV/FJ2015 (Table 2). An NS2a gene (20,246–21,079nt) was predicated in PPCoV with 277 aa (Table 1); it had 53.1 to 55.8% aa identities to related alphacoronaviruses (Table 2) and <41.0% aa identities to corresponding domains of betacoronaviruses. Multiple alignments of amino acid sequences of NS2a are shown in Supplementary Figure S2. Analyses indicated that NS2a of PPCoV contained no membrane helices, with N-glycosylated sites at 149 and 181, and no O-glycosylation sites. An InterProScan search showed that NS2a had homologs to the cyclic phosphodiesterase superfamily. However, the definite function of NS2a in CoVs remain unclear.

Although showing similar structures to those of related CoVs (Supplementary Figure S3), S protein of PPCoV shared only 69.5 to 73.0% aa identities with LRNVs, and <43.4% aa identities with other CoVs. A transmembrane domain (from residues 1,063 to 1,085) was predicated in the S protein, with most of the protein (residues 1 to 1,062) on the outside of the virus and a cytoplasmic tail (residues 1,086 to 1,120), which was similar to that of the Wencheng Sm shrew CoV and bat coronavirus HKU10 (Lau et al., 2012a; Wang et al., 2017). Meanwhile, N-glycosylated modifications of S protein of PPCoV were detected in 13 sites (93, 150, 190, 300, 326, 540, 657, 778, 919, 981, 997.1,010, and 1,050).

Other predicted proteins of NS3, E, M, NS7a and N proteins of PPCoV showed a difference (<90% aa identities) to corresponding proteins of other alpha CoVs, especially in NS8 (<40% aa identities). Notably, M protein of PPCoV was genetically close to that of RtMruf-CoV-1/JL2014 (87.6% aa identity). However, PPCoV shared the highest (70.7%) aa identity to Alpha-CoV UKRn3 and LRNV Lucheng-19 in N protein (Table 2). NS9 locating downstream from the N gene was also found; this is present in some genomes of alphacoronaviruses, such as LRNVs, transmissible gastroenteritis virus (TGEV), porcine respiratory coronavirus (PRCV), and Rhinolophus bat coronavirus HKU2.

Phylogenetic Relationships

In order to evaluate the evolutionary status of PPCoV, we constructed phylogenetic trees based on aa sequences of ORF1ab, S, M, and N from genera Alphacoronavirus and Betacoronavirus (Figures 2, 3 and Supplementary Figure S3). The phylogenetic tree of ORF1ab (Figure 2) showed that PPCoV formed a distinct lineage within the genus Alphacoronavirus, and was closely related to the cluster of rodent alphacoronaviruses (viruses obtained from rat species; Tsoleridis et al., 2019). A slightly different clustering pattern was observed in phylogenetic trees based on M and N protein sequences, showing that PPCoV was clustered within the rodent coronavirus clade (Supplementary Figure S4).

Evolutionary tree based on S protein sequences (Figure 3) revealed significantly different clustering patterns to phylogenetic trees based on ORF1ab, M, and N proteins. The tree of S indicated that PPCoV formed a divergent cluster with Shrew-CoV/Tibet2014, Wencheng Sm shrew CoV, SADS-CoV, bat coronavirus HKU2, and rodent alphacoronaviruses. However, the clade of the PPCoV as well as other phylogenetically related sub-genera all clustered together with the betacoronaviruses with high bootstrap values.

Recombination Analysis and Estimation of Divergence Times

All alpha CoVs clustered together in tree of ORF1ab, but some of alpha CoVs clustered together with the beta CoVs in tree of S. This discrepancy in clades between trees suggested a potential recombination for S gene. However, recombination analysis based on S genes of whole set of alpha- and betacoronaviruses indicated that no recombination event was detected for PPCoV (Supplementary Figure S5).

For tip-dating analysis, the RdRp and M genes of PPCoV and related alphacoronaviruses were selected and checked for potential recombination regions. Then, the entire sequence alignment (recombination region free) was used to evaluate the divergence time with relaxed clock model as previously described (Lau et al., 2012a,b). The most recent common ancestor (tMRCA) of PPCoV and rodent alphacoronaviruses was estimated to be in the year of 1935 (HPDs 1782 to 1987; approximately 84years ago; Figure 4). The molecular clock analysis based on the M genes also showed the tMRCA of PPCoV and rodent alphacoronaviruses was at 1935 (HPDs 1,420 to 2,008; Supplementary Figure S6). The mean substitution rate of RdRp and M genes was calculated to be 2.1×10⁻⁴ and 4.8×10⁻⁴ per site per year, respectively, which is similar to previous report (Woo et al., 2012).

Virus-Host Co-divergence Analysis

The evolutionary history of a virus can mirror that of its host over long-term evolutionary timescales (Shi et al., 2016). A tanglegram was obtained by comparing the tree topologies of viruses and their hosts, which revealed a significant codivergence for alphacoronaviruses (p<0.01). The viruses clearly clustered with their hosts at the order level (Figure 5). However, the topology of hosts did not fully match those of viruses with 14 host switching events, two loss events, and three duplication events. Until now, the host species of alphacoronaviruses were clustered into six orders, including orders Artiodactyla, Carnivora, Chiroptera, Eulipotyphla, Rodentia, and Primates (Figure 5). PPCoV was identified in Ochotona curzoniae belonging to order Lagomorpha, which is different from host taxon of other alphacoronaviruses. This extends the host range of alphacoronaviruses. Meanwhile, virus-host co-speciation analysis showed that plateau pikas appear to be infected with the ancestor of PPCoV through cross-species transmission (Figure 5), and indicated that the common ancestor for PPCoV and rodent alphacoronaviruses can now trace back to 11.6 (range 0.2–44.6) million years ago.