Identification of Diverse Bat Alphacoronaviruses and Betacoronaviruses in China Provides New Insights Into the Evolution and Origin of Coronavirus-Related Diseases

Outbreaks of severe acute respiratory syndrome (SARS) in 2002, Middle East respiratory syndrome in 2012 and fatal swine acute diarrhea syndrome in 2017 caused serious infectious diseases in humans and in livestock, resulting in serious public health threats and huge economic losses. All such coronaviruses (CoVs) were confirmed to originate from bats. To continuously monitor the epidemic-related CoVs in bats, virome analysis was used to classify CoVs from 831 bats of 15 species in Yunnan, Guangxi, and Sichuan Provinces between August 2016 and May 2017. We identified 11 CoV strains from 22 individual samples of four bat species. Identification of four alpha-CoVs from Scotophilus kuhlii in Guangxi, which was closely related to a previously reported bat CoV and porcine epidemic diarrhea virus (PEDV), revealed a bat-swine lineage under the genus Alphacoronavirus. A recombinant CoV showed that the PEDV probably originated from the CoV of S. kuhlii. Another alpha-CoV, α-YN2018, from Rhinolophus sinicus in Yunnan, suggested that this alpha-CoV lineage had multiple host origins, and α-YN2018 had recombined with CoVs of other bat species over time. We identified five SARS-related CoVs (SARSr-CoVs) in Rhinolophus bats from Sichuan and Yunnan and confirmed that angiotensin-converting enzyme 2 usable SARSr-CoVs were continuously circulating in Rhinolophus spp. in Yunnan. The other beta-CoV, strain β-GX2018, found in Cynopterus sphinx of Guangxi, represented an independently evolved lineage different from known CoVs of Rousettus and Eonycteris bats. The identification of diverse CoVs here provides new genetic data for understanding the distribution and source of pathogenic CoVs in China.

Introduction

Coronaviruses (CoVs) are a group of enveloped viruses with a large positive single-stranded RNA genome (∼26–32 kb in length) of the subfamily Coronavirinae under the family Coronaviridae. The complete genome of CoV contains five major open reading frames (ORFs) that encode replicase polyproteins (ORF1ab), spike glycoprotein (S), envelope protein (E), membrane protein (M), and nucleocapsid protein (N) flanked by a 5′- untranslated region (UTR) and a 3′- UTR. Currently, members of the subfamily Coronavirinae are classified into four genera, Alphacoronavirus, Betacoronavirus, Gammacoronavirus, and Deltacoronavirus (Fehr and Perlman, 2015; Su et al., 2016). CoVs can cause upper and lower respiratory diseases, gastroenteritis, and central nervous system infections in a wide variety of avian and mammalian hosts. Some CoVs are human pathogens that cause mild to severe disease, including NL63 and 229E of the genus Alphacoronavirus, and severe acute respiratory syndrome CoV (SARS-CoV), Middle East respiratory syndrome CoV (MERS-CoV), OC43, and HKU1 of the genus Betacoronavirus (Dijkman and van der Hoek, 2009; Desforges et al., 2014; van den Brand et al., 2015; Lim et al., 2016). The SARS pandemic originated in Guangdong Province in China between 2002 and 2003, and spread to 29 countries, resulting in nearly 8,000 cases and 800 deaths worldwide (Donnelly et al., 2004). Since being discovered in Middle Eastern countries in 2012, MERS-CoV has infected 2,260 people with a current fatality rate of 35.5% (Mackay and Arden, 2015; van den Brand et al., 2015; de Wit et al., 2016).

With the discovery of diverse wildlife-borne CoVs indifferent regions of the world, previous studies have indicated that bats are the main and original natural reservoirs of Alphacoronavirus and Betacoronavirus (Calisher et al., 2006; Woo et al., 2012). The special metabolic and immune systems allow bats to tolerate diverse viruses, and the social roosting behavior of many bat species lets them serve as ideal incubators for the occurrence of frequent co-infection, recombination, and intra-species transmission of CoVs (Zhang et al., 2013; O’Shea et al., 2014; Hu et al., 2015; Lu et al., 2015; Sola et al., 2015). Although 15 years have passed without a recurrence of the SARS outbreak, the consistent discovery of SARS-related CoVs (SARSr-CoVs) in Rhinolophus bats indicates an ongoing risk of SARS re-emergence (Li et al., 2005; Balboni et al., 2012; Yang et al., 2013; Wu et al., 2016b; Hu et al., 2017). Recently, a CoV of Rhinolophus bat origin, swine acute diarrhea syndrome CoV (SADS-CoV), was found to be responsible for the death of 24,693 piglets across four farms in Guangdong Province (Zhou et al., 2018). Besides similar clinical signs, the outbreak of SADS coincided with infection by another swine enteric CoV, porcine epidemic diarrhea virus (PEDV; Zhou et al., 2018). It is worth mentioning that a PEDV-related bat CoV, BtCoV/512, was previously identified in Scotophilus kuhlii bats (Tang et al., 2006; Wang et al., 2016; Gerdts and Zakhartchouk, 2017; Zhou et al., 2018). All of these findings emphasize the importance of continuous monitoring for ecological diversity of CoVs in bats to minimize the impact of potential CoV-related diseases on public health and economic growth.

As the original site of the SARS epidemic in humans and SADS outbreak in pigs (Heymann et al., 2013), Guangdong Province is also the primary region in China in which wildlife is being consumed. Since the supply of wildlife for consumption in Guangdong mainly comes from Yunnan and Guangxi Provinces (Wu et al., 2016b), understanding the presence and dynamic changes in SARS- or SADS- related bat CoVs within and between these provinces is important for predicting and tracing CoV-related diseases. Sichuan is next to Yunnan Province, and has high ecological diversity of animal resources, and it is the largest pig-farming province in China. Although PEDV-related pig diseases were reported previously (Cai et al., 2016), no systematic survey of bat-borne CoVs has been conducted in Sichuan to establish if there are SARS- or SADS- related CoVs.

In this study, a survey of bat-borne CoVs was performed in bats of Guangxi, Yunnan, and Sichuan provinces between 2016 and 2017. Samples from Rhinolophus bats in Sichuan and Yunnan Provinces, and S. kuhlii and Cynopterus sphinx in Guangxi Province were found to be CoV-positive, and a large number of sequencing reads classified into diverse alpha- and beta-CoVs were obtained using virome analysis. The characterization of 11 CoVs in bats from different regions provides clues for the evolutionary relationships of PEDV, SARS-CoV, and other bat-associated CoVs.

Materials and Methods

Bat Sampling

Bats were treated according to the guidelines of Regulations for the Administration of Laboratory Animals (Decree No. 2 of the State Science and Technology Commission of the People’s Republic of China, 1988). The sampling was approved by the Ethics Committee of Institute of Pathogen Biology, Chinese Academy of Medical Sciences & Peking Union Medical College (Approval number: IPB EC20100415). The bat species were initially determined morphologically and subsequently confirmed by sequence analysis of mitochondrial cytochrome b DNA (Tang et al., 2006; Du et al., 2016). Anal swab samples from captured bats were immersed in virus sampling tubes (Yocon, China) containing maintenance medium and temporarily stored at −20°C. The samples were then transported to the laboratory and stored at −80°C. The accurate sampling locations were recorded by place name, latitude, and longitude.

Viral DNA and RNA Library Construction and Next-Generation Sequencing

The samples of each species were pooled by adding 1 ml from each maintenance medium sample into one new sample tube. The pooled samples, classified by species, were processed with a virus-particle-protected, nucleic acid purification method as described in our previous studies (Wu et al., 2012, 2016a, 2018). The samples were homogenized and subsequently filtered through a 0.45 μm polyvinylidene difluoride filter (Millipore, Germany). The filtered samples were then centrifuged at 150,000 × g for 3 h at 4°C. To remove naked DNA and RNA, the pellet was digested in a cocktail of DNase and RNase enzymes. The viral DNA and RNA were simultaneously isolated using a QIAmp MinElute Virus Spin Kit (Qiagen, United States). First-strand viral cDNA was synthesized using the primer K-8N and a Superscript III system (Invitrogen, United States). The cDNA was converted into dsDNA by Klenow fragment (NEB, United States). Sequence-independent PCR amplification was conducted using primer K. The PCR products were analyzed by agarose gel electrophoresis. All DNA smears larger than 500 bp were extracted by a MinElute Gel Extraction Kit (Qiagen). The extracted nucleic acid libraries were then analyzed using an Illumina HiSeq2500 sequencer, for a single read of 100 bp. The sequence reads were filtered using previously described criteria (Wu et al., 2012, 2016a, 2018). Reads with no call sites, reads with similarity to the sequencing adaptor and the primer K sequence, duplicate reads, and low-complexity reads were removed.

Taxonomic Assignment

Sequence-similarity-based taxonomic assignments were conducted as described in our previous study (Wu et al., 2012). Valid sequence reads were aligned to sequences in the NCBI non-redundant nucleotide database and non-redundant protein database using BLASTn and BLASTx, respectively. The taxonomies of the aligned reads with the best BLAST scores (E score <10^–5) were parsed by the MEGAN 6 - MetaGenome Analyzer.

Viral Prevalence and Genome Sequencing

CoV screening of individual samples was performed by amplifying a 440-bp fragment of the RNA-dependent RNA polymerase (RdRp) gene of CoVs using conserved primers (5′-GGTTGGGACTATCCTAAGTGTGA-3′ and 5′-CCATC ATCAGATAGAATCATCATA-3′), as described previously (Du et al., 2016).

Sequence reads classified into the same virus genus were extracted. The accurate locations of the reads were determined based on the alignment results exported with MEGAN 6. Specific primers were designed from the located sequence reads. Fragments between reads were amplified with nested specific primers and then sequenced. The remaining genomic sequences were determined using 5′- and 3′- rapid amplification of cDNA ends.

Phylogenetic and Recombination Analysis

MEGA6.0 was used to align nucleotide sequences and deduced amino acid sequences using the MUSCLE package and default parameters. The best substitution model was then evaluated by the Model Selection package. Finally, we constructed a maximum-likelihood method using an appropriate model to process the phylogenetic analyses with 1,000 bootstrap replicates. Recombination among CoVs was detected with SimPlot software. All analyses were performed with a Kimura model, a window size of 1000 bp, and a step size of 100 bp.

Nucleotide Sequence Accession Numbers

All genome sequences were submitted to GenBank. The accession numbers for the five bat alpha-CoVs were MK211369 to MK211373. The accession numbers for the six bat beta-CoVs were MK211374 to MK211379.

Results

Virome and Prevalence Analysis of CoVs

Anal swabs from 831 bats of 15 species were collected from Yunnan, Guangxi, and Sichuan Provinces in China between August 2016 and May 2017. All anal swab samples were classified by species and then combined into 27 pools and subjected to virome analysis. CoV-related sequencing reads were detected in ten of the 27 sample pools. These 10 pools related to 22 Rhinolophus spp. samples from Sichuan, 176 Rhinolophus sinicus samples from Yunnan, 150 S. kuhlii samples from Guangxi, and 30 C. sphinx samples from Guangxi. Pan CoV screening revealed that 22 individual samples were CoV-positive (Table 1 and Figure 1). A total of 1,040,901 reads of 100 bp in length showed the best matches with Coronavirinae viral proteins in the NCBI non-redundant database, and these reads showed 40–100% amino acid (aa) identity with known alpha- or beta-CoVs.

TABLE 1

Table 1. Coronavirus distribution in different bat species and locations.

FIGURE 1

Figure 1. Numbers of bat samples from Sichuan, Yunnan, and Guangxi Provinces. The size of the pie chart is in proportion with the total number of specimens collected from each province. Red color represents samples positive for CoV, and blue color represents negative samples.

Eleven viruses of representative positive samples were selected for genomic sequencing as quasi-species. There were five viruses in the genus Alphacoronavirus; four from S. kuhlii in Guangxi Province were named as BtSk-AlphaCoV/GX2018A (α-GX2018A), BtSk-AlphaCoV/GX2018B (α-GX2018B), BtSk-AlphaCoV/GX2018C (α-GX2018C), and BtSk-AlphaCoV/GX2018D (α-GX2018D). The other virus from R. sinicus in Yunnan Province was named as BtRs-AlphaCoV/YN2018 (α-YN2018). There were six CoVs in the genus Betacoronavirus. The virus identified from Rhinolophus spp. in Sichuan Province was named as BtRl-BetaCoV/SC2018 (β-SC2018); four CoVs from R. sinicus in Yunnan Province were named as BtRs-BetaCoV/YN2018A (β-YN2018A), BtRs-BetaCoV/YN2018B (β-YN2018B), BtRs-BetaCoV/YN2018C (β-YN2018C), and BtRs-BetaCoV/YN2018D (β-YN2018D). One CoV from C. sphinx in Guangxi Province was named as BtCs-BetaCoV/GX2018 (β-GX2018).

Genome Organization

As shown in Table 2 and Figure 2, with the typical genome organization, the full genome sizes of 11 CoVs ranged from 28, 146 (α-GX2018C) to 30, 256 bases (β-YN2018B) and their G + C contents ranged from 0.38 to 0.41. There were several accessory ORFs in each CoV strain. These accessory ORFs were mainly present in the regions between S and E (named ORFs-e), M and N (named ORFm-n), or downstream of N (named ORFn). The accessory ORFs of α-GX2018 (A–D), α-YN2018, and β-GX2018 were mainly present in ORFs-e and ORFm-n regions; the small ORFs of β-SC2018 and β-YN2018 (A–D) were mainly present in ORFs-e and ORFm-n regions. We named these accessory ORFs according to previous studies.

TABLE 2

Table 2. Predicted ORFs in the genomes of bat CoVs.

FIGURE 2

Figure 2. Genome analysis of eleven novel CoVs. Color-coding for different genomic regions as follows. Green, non-structural polyproteins ORF1a and ORF1b; yellow, structural proteins S, E, M and N; blue, accessory proteins ORF3, ORF6, ORFx, ORF7 and ORF8; orange, untranslated regions.

As shown in Table 3, The core sequences of the leader transcriptional regulatory sequence (TRS; 5′-CUAAAC-3′) were identified in the 5′ untranslated sequences of α-GX2018 (A–D) and α-YN2018, which is unique to Alphacoronavirus (Madhugiri et al., 2014). The TRS motifs of S, ORF3, E, and M genes in α-GX2018 (A-D) differed from the core sequences of the leader TRS (Supplementary Table 1). The TRS motif of S was identified as 5′-CCAAAU-3′ or 5′-CCAAAC-3′; the TRS motif of ORF3 was identified as 5′-CAUUAC-3′; the TRS motif of E was identified as 5′-CUAGAC-3′; the TRS motif of M was identified as 5′-AUAAAC-3′. An alternative TRS motif (5′-GUAAAC-3′) of α-YN2018 was discovered preceding ORF3, which differed by 1 nucleotide (nt) from those of all other genes in α-YN2018. In addition to the E of β-GX2018, the TRS motifs of six betacoronaviruses all conformed to the consensus motif 5′-ACGAAC-3′, which is unique to Betacoronavirus (Woo et al., 2010). An alternative TRS motif (5′-UCGAAC-3′) was discovered preceding the E in β-GX2018.

TABLE 3

Table 3. Transcription regulatory sequences (TRS) for four representative bat CoVs^a,b.

Sequence Similarity Analysis

Four Alpha-CoVs of S. kuhlii in Guangxi Province

Four alpha-CoVs were identified from S. kuhlii. The full-length genomes of strains α-GX2018A, B, C, and D were sequenced. The genome sequence similarity among these four CoVs, PEDV, and BtCoV/512 was examined by Simplot analysis (Figure 3A) and pairwise alignment (Supplementary Table 2). The α-GX2018A, B, and C were closely related to a known CoV, BtCoV/512, identified in S. kuhlii of Hainan Province with 96.9–100% aa identity for ORF1b, ORF3, E, M, and N. The differences mainly occurred in S (93.4–96.1% aa identity) and ORFx (65.7–67.6% aa identity). There were also significant differences between α-GX2018A and B and BtCoV/512 at the front end of ORF1a (1000–2500 nt) with 81.3–82% nt identity. In addition, although the α-GX2018D was closely related to that of BtCoV/512 with 92.2–99.4% aa identity for ORF1a, ORF1b, ORF3, and M, it was unexpected that the aa identity of α-GX2018D and BtCoV/512 in S1 region was 39.7%, which was lower than that of α-GX2018D and PEDV (42.7% aa identity). Although the aa identity of α-GX2018D and PEDV in S1 region was only 42.7%, we have not found any known viruses with higher identity than this.

FIGURE 3

Figure 3. Similarity plot analysis. (A) Full-length genome sequence of BtCoV/512 was used as a query sequence and α-GX2018A, B, C, D, and PEDV as reference sequences. (B) Full-length genome sequence of α-YN2018 was used as a query sequence and Ro-BatCoV HKU10, Hi-batCov HKU10, Ca-batCov/Kenya, and α-HuB2013 as reference sequences. (C) Full-length genome sequence of SARS-CoV SZ3 was used as a query sequence and β-YN2018A, B, C, D and β-SC2018 as reference sequences. (D) Full length genome sequence of β-GX2018 was used as query sequence and BatCoV HKU9-4, BatCoV HKU9-10-1, Ro-BatCoV HKU9, and BatCoV Philippines/Diliman as reference sequences. These analyses were performed with the Kimura model, a window size of 1000 base pairs and a step size of 100 base pairs.

Novel Alpha-CoV of R. sinicus in Yunnan Province

A novel alpha-CoV (α-YN2018) was identified inR. sinicus of Yunnan Province. The genome sequencesimilarity among α-YN2018, BtRf-alphaCoV/HuB2013 (α-HuB2013), Hipposideros bat CoV HKU10 (Hi-batCov HKU10), Rousettus bat CoV HKU10 (Ro-BatCoV HKU10), and Cardioderma bat coronavirus/Kenya/KY43/2006 (Ca-batCov/Kenya) was examined by Simplot analysis (Figure 3B) and pairwise alignment (Supplementary Table 3). This virus showed low sequence similarity with any known alpha-CoVs. α-YN2018 showed the highest sequence similarity to α-HuB2013 of Rhinolophus ferrumequinum with 63.9–92.4% aa identity in ORF1a, ORF1b, and N; to Hi-batCov HKU10 with 50–85.7% aa identity in S, M, and ORF7b; to Ca-batCov/Kenya in ORF3a with 60.6% aa identity; and to Ro-BatCoV HKU10 in E with 74.3% aa identity. The sequences of ORF7a from α-YN2018 showed <25% aa identity with those of the afore-mentioned known CoVs. In addition, we did not find any sequences that were similar to ORF3b and ORF3c of α-YN2018 from known viruses. According to the ICTV criteria, CoVs that share <90% aa sequence identity in the conserved replicase domains are considered to belong to different species; thus, α-YN2018 could be considered a new species under the genus Alphacoronavirus.

Five Beta-CoVs of Rhinolophus spp. From Sichuan and Yunnan Provinces

We obtained five beta-CoVs from Rhinolophus spp. from Yunnan and Sichuan Provinces, which were named as β-YN2018A, β-YN2018B, β-YN2018C, β-YN2018D, and β-SC2018. All of these beta-CoVs had high sequence identity with SARSr-CoVs. These five SARSr-CoVs showed 92–97% nt identity with human SARS-CoV SZ3. The genome sequence similarity among the five SARSr-CoVs and SARS-CoV SZ3 strain was examined by Simplot analysis (Figure 3C). These five SARSr-CoVs were highly conserved and shared a uniformly high sequence similarity to SARS-CoV with 94.2–100% aa identity in ORF1a, ORF1b, E, M, and N (Supplementary Table 4). There were some differences between the SARSr-CoVs and SARS-CoV SZ3 in the coding regions of ORF3a, ORF3b, ORF7a, and ORF7b. After further analysis, the SARS-CoV SZ3 showed the highest sequence similarity to β-YN2018B with 92–96% aa identity in S, ORF3a, and ORF3b; and to β-YN2018A with 93.2–95.9% aa identity in ORF7a and ORF7b. In contrast, considerable genetic diversity was shown in the S (72.9–92.4% aa identity) and ORF8 (32.6–34.1% aa identity) among the five SARSr-CoVs and SARS-CoV SZ3 strain.

β-YN2018B shared 98.8% aa identity with SARSr-CoV WIV1, which was higher than that of the other four SARSr-CoVs (78.8–79.1% aa identity) (Figure 3C and Supplementary Table 4). Remarkably, WIV1 and WIV16 (or Rs4874) were identified between 2012 and 2013, and they could use the human angiotensin converting enzyme II (ACE2) receptor for entry, like SARS-CoV (Ge et al., 2013; Yang et al., 2015).

ORFx was also found in the genomes of β-YN2018B and β-YN2018D. ORFx of β-YN2018B and β-YN2018D both shared 98.8% aa identity with WIV1. Except for ORF3b (95.6% aa identity), β-YN2018B was highly conserved and shared a uniformly high sequence similarity to SARS-CoV WIV1 in every ORF (98.2–100% aa identity) (Supplementary Table 4).

Novel Beta-CoV of C. sphinx in Guangxi Province

A novel beta-CoV, β-GX2018, was identified from this bat species. The genome sequence similarity among β-GX2018 and four known CoVs, BatCoV HKU9-4, BatCoV HKU9-10-1, Ro-BatCoV HKU9, and BatCoV Philippines/Diliman was examined by Simplot analysis (Figure 3D) and pairwise alignment (Supplementary Table 5). The full-length genome sequence of this virus had sequence similarity to that of HKU9- related CoVs of Rousettus bats and GCCDC1-related CoVs of Eonycteris bats (Woo et al., 2007; Huang et al., 2016), but the sequence identity between them was low. The sequence similarity between β-GX2018 and HKU9- and GCCDC1-related CoV s were 51–90.8% aa identity in ORF1a, ORF1b, S, ORF3, E, M, and N. Some partially sequenced genome segments of a recently reported CoV, Philippines/Diliman1525G2/2008 (BatCoV Philippines/Diliman) (Watanabe et al., 2010), identified in C. brachyotis from the Philippines, showed the highest aa sequence identity with the β-GX2018 in M (96% aa identity), N (91% aa identity), ORF7a (80% aa identity), ORF7b (84.3% aa identity), and ORF7c (80.8% aa identity) (Watanabe et al., 2010). According to the ICTV criteria, β-GX2018 could be considered a new species under the genus Betacoronavirus.

Phylogenetic Analysis

To determine further the evolutionary relationships between these CoVs, phylogenetic trees were constructed using the aa sequences of RdRp, S, S1, E, M and N (Figure 4). As we have speculated before, α-GX2018A, α-GX2018B, α-GX2018C, α-GX2018D and α-YN2018 clustered with the genus Alphacoronavirus, and β-SC2018, β-YN2018A, β-YN2018B, β-YN2018C, β-YN2018D, and β-GX2018 clustered with the genus Betacoronavirus. The results of the phylogenetic analyses were consistent with those of the sequence identity analyses, and confirmed that the identified CoVs could be divided into four lineages.

FIGURE 4

Figure 4. Phylogenetic trees based on aa sequences of RdRp, S, S1, E, M, and N. The trees were constructed by the maximum likelihood method using appropriate models (WAG + G + I for RdRp and N; LG + G + I for S, S1, and E; rtREV + G for M) with bootstrap values determined by 1000 replicates. Only bootstraps >40% are shown. The scale bars represent 0.05 (RdRp), 0.5 (S and S1), and 0.2 (E, M, and N), respectively.

Lineage E

α-GX2018A, α-GX2018B and α-GX2018C clustered with BtCoV/512, and their branches were short, reflecting the high sequence similarities. In addition, phylogenetic analysis of the S protein supported our previous analysis that α-GX2018D represented an independent clade between PEDV and BtCoV/512. Further analysis of S1 revealed that α-GX2018D represented a separate evolution distant from all other members of lineage E. From the length of the branch, α-GX2018D had a closer relationship with PEDV, which means that α-GX2018D may have undergone recombination with the ancestors of PEDV and BtCoV/512.

Lineage F

α-YN2018 clustered with Hi-batCov HKU10, Ro-BatCoV HKU10, Ca-batCoV/Kenya, and α-HuB2013. However, α-YN2018 had closer relationships with α-HuB2013 of R. ferrumequinum in RdRp and N, and closer relationships with HKU10 of Hipposideros and Rousettus bats in S and M. This means that recombination may have occurred between α-YN2018 and these three related CoVs.

Lineage B

β-SC2018, β-YN2018A, β-YN2018B, β-YN2018C, and β-YN2018D clustered with SARS- and SARSr- CoVs under lineage B. In RdRp, E, M, and N trees, these novel SARSr-CoVs and SARS-CoVs had short branches, which indicated that they were closely related in these regions. Phylogenetic analysis of S and S1 revealed that β-YN2018B was closer to SARS-CoV SZ3 and SARS-CoV Tor2 as well as WIV1 and WIV16.

Lineage D

β-GX2018 of C. sphinx clustered with HKU9-related CoVs of Rousettus and GCCDC1-related CoVs of Eonycteris under lineage D; however, β-GX2018 always represented a separate lineage distinct from other members of lineage D. Although BatCoV Philippines/Diliman of C. brachyotis had a closer relationship with α-YN2018 in M and N, it was only partially sequenced, so we could not analyze its evolutionary relationship with α-YN2018 in RdRp, S, and E, which hindered further analysis.

Recombination Analysis

We found that recombinant events had occurred among α-GX2018D, other bat CoVs, and PEDV of lineage E. α-GX2018D showed the highest degree of similarity to BtCoV/512 in the ORF1ab region, and the highest degree of similarity to PEDV in the S1 region. Evolutionary analysis showed that the S of α-GX2018D may have acted as the common ancestors of BtCoV/512 and PEDV. We analyzed possible recombination events in lineage E using SimPlot software. Similarity plot and bootscan results demonstrated that the S-N region of the α-GX2018D had the highest degree of similarity to PEDV, and complicated recombination may have happened between bat CoV and PEDV in the ORF1ab region (Figure 5A).

FIGURE 5

Figure 5. Detection of potential recombination events by bootscan analysis. (A) Full-length genome sequence of α-GX2018D was used as a query sequence and BtCoV/512, α-GX2018A, and PEDV as reference sequences. (B) Full-length genome sequence of α-YN2018 was used as a query sequence and Ro-BatCoV HKU10, Hi-batCov HKU10, and α-HuB2013 as reference sequences. All analyses were performed with a Kimura model, a window size of 1000 base pairs, and a step size of 100 base pairs. The gene map of query genome sequences is used to position breakpoints.

For CoVs of lineage F, at least two recombination events had occurred (Figure 5B). Similarity plot and bootscan analysis showed that the S1 subunits of the α-YN2018 had the highest degree of similarity to Ro-BatCoV HKU10, while their S2 subunits had the highest degree of similarity to Hi-batCov HKU10 (Fig 5B). The complicated recombination history between these alpha-CoVs suggests frequent gene transfers, especially of S, among different CoVs, which may be the result of the cross-species transmission of these CoVs.

Discussion

The identification of SARSr-CoV in bats in 2005 first established a genetic relationship between bats and human SARS-CoVs (Lau et al., 2005; Li et al., 2005). After that, in the following 13 years, the constant discovery of bat SARSr-CoVs suggested the risk of re-emergence of SARS-CoV originating from bat species. Other than that, two CoV-related cases, human MERS in 2012 and SADS in 2017, which had a great impact on health and the economy were identified as being of bat origin (Omrani et al., 2015; Zhou et al., 2018). In our previous study, a large-scale virome analysis was conducted between 2010 and 2013 to understand the ecological diversity of bat viruses and the bat origin of emerging infectious diseases. Rhinolophus bats in Yunnan and Guangxi Provinces have provided genetic clues to the origin of SARS-CoV, and a SADS-CoV- and HKU2-related CoV, BtRf-alphaCoV/YN2012, has also been found in Yunnan Province (Wu et al., 2016a). Sichuan Province has diverse wildlife resources and the largest pig breeding in China, and a MERS-related CoV was reported in bats in this province in 2013 (Yang et al., 2014). It is necessary for us to conduct continuous surveillance for the prevalence, genetic diversity and geographical distribution of CoVs in some bat species in these potential high-risk regions. Using virome analysis, our study investigated the presence and genetic diversity of CoVs circulating in bat species from Sichuan, Yunnan, and Guangxi Provinces between 2016 and 2017.

The close relationship between five alpha-CoVs in S. kuhlii in Guangxi Province and a previously reported CoV in S. kuhlii in Hainan Province reveals that this bat species from different geographic locations contained the same CoV species, but with distinct S proteins. The S1 region of α-GX2018D is more closely related to that of PEDV than that of any other bat CoV, and is phylogenetically located at the root of the lineage of PEDV and S. kuhlii CoVs. This finding suggests the presence of a much closer common ancestor of PEDV in bats CoVs. The recombination analysis supports that the S region of PEDV may come from the ancestor of α-GX2018D. Considering that bat species such as S. kuhlii prefers to circulate around farmhouses and pigsties, and other PEDV-related sequences have been reported in bats in previous reports (Simas et al., 2015; Phan et al., 2018), the relationship between bat CoVs and PEDV should be further investigated.

The discovery of β-SC2018 provides evidence that Rhinolophus in Sichuan Province could harbor SARSr-CoVs, and reveals a broader geographical distribution of these viruses. The discovery of β-YN2018B here, combined with previously reported WIV1 and WIV16 (Ge et al., 2013; Yang et al., 2015), revealed that these ACE2- adaptable SARSr-CoVs were continuously circulating in R. sinicus of Yunnan Province at least from 2012 to 2016. These findings emphasize again the importance of continuous longitudinal monitoring for SARSr-CoVs in Rhinolophus over a wide area. Although no SADS and HKU2-related CoV were found in Rhinolophus bats in our study, the surveillance of this CoV type in bats is still needed in the future.

In our previous study, by the identification of correlated alpha-CoVs in different Miniopterus species from multiple regions, the characteristics of co-infection, recombination, and host-shifting for alpha-CoVs have been described. Here, the alpha-CoVs of lineage F found in more diverse bat species (including frugivorous and insectivorous bats) extends the known host range of these viruses. However, the recombination events found between α-YN2018, Ro-BatCoV HKU10, and Hi-batCov HKU10, and the multiple origins of the S of α-YN2018, further indicate that a more distant host-shifting in accordance with the recombination of the S region may have happened in the evolutionary history of these viruses among diverse host species.

When we collected S. kuhlii in Guangxi Province, we found that there were many C. sphinx at the sampling site. The Pteropodidae bats (C. sphinx bat belongs to Pteropodidae) are natural hosts of several emergent human pathogens such as Hendra virus, Nipah virus, and Ebola virus (Leroy et al., 2005; Halpin et al., 2011; Sendow et al., 2013; Mari Saez et al., 2015; Paez et al., 2017). An Ebola-related filovirus was found in the Chinese Pteropodidae bats recently (Yang et al., 2019), and thus, we collected samples of C. sphinx for virome analysis. The C. sphinx beta-CoV, β-GX2018, represents an evolved independent lineage different from HKU9 of Rousettus and GCCDC1 of Eonycteris under lineage D. Based on the partially sequenced region, CoV strain Philippines/Diliman1525G2/2008 previously identified in C. brachyotis from the Philippines, could also be clustered into this separate β-GX2018 clade. By providing full-length sequence data, this finding reveals a new host range of this virus in an unreported location.

Author Contributions

YH, QJ, and ZW conceived the experiments, analyzed the results, and wrote the manuscript. YH, JD, HS, and ZW conducted the experiments and analyzed the results. JZ, GZ, and SZ collected the specimens. All authors reviewed the manuscript.

Funding

This work was supported by the National S&T Major Project “China Mega-Project for Infectious Disease” (Grant No. 2018ZX10101001), the Non-profit Central Institute Fund of Chinese Academy of Medical Sciences (Grant No. 2018RC310018), the CAMS Innovation Fund for Medical Sciences (Grant Nos. 2016-I2M-1-014 and 2017-I2M-B&R-12), and the National Natural Science Foundation of China (Grant No. 81772228).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.01900/full#supplementary-material

References

Balboni, A., Battilani, M., and Prosperi, S. (2012). The SARS-like coronaviruses: the role of bats and evolutionary relationships with SARS coronavirus. New Microbiol. 35, 1–16.

PubMed Abstract | Google Scholar

Cai, Y., Yin, W., Zhou, Y., Li, B., Ai, L., Pan, M., et al. (2016). Molecular detection of Porcine astrovirus in Sichuan Province, China. Virol. J. 13:6. doi: 10.1186/s12985-015-0462-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Calisher, C. H., Childs, J. E., Field, H. E., Holmes, K. V., and Schountz, T. (2006). Bats: important reservoir hosts of emerging viruses. Clin. Microbiol. Rev. 19, 531–545. doi: 10.1128/cmr.00017-06

PubMed Abstract | CrossRef Full Text | Google Scholar

de Wit, E., Van Doremalen, N., Falzarano, D., and Munster, V. J. (2016). SARS and MERS: recent insights into emerging coronaviruses. Nat. Rev. Microbiol. 14, 523–534. doi: 10.1038/nrmicro.2016.81

PubMed Abstract | CrossRef Full Text | Google Scholar

Desforges, M., Le Coupanec, A., Stodola, J. K., Meessen-Pinard, M., and Talbot, P. J. (2014). Human coronaviruses: viral and cellular factors involved in neuroinvasiveness and neuropathogenesis. Virus Res. 194, 145–158. doi: 10.1016/j.virusres.2014.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Dijkman, R., and van der Hoek, L. (2009). Human coronaviruses 229E and NL63: close yet still so far. J. Formos. Med. Assoc. 108, 270–279. doi: 10.1016/S0929-6646(09)60066-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Donnelly, C. A., Fisher, M. C., Fraser, C., Ghani, A. C., Riley, S., Ferguson, N. M., et al. (2004). Epidemiological and genetic analysis of severe acute respiratory syndrome. Lancet Infect. Dis. 4, 672–683.

PubMed Abstract | Google Scholar

Du, J., Yang, L., Ren, X., Zhang, J., Dong, J., Sun, L., et al. (2016). Genetic diversity of coronaviruses in Miniopterus fuliginosus bats. Sci. China Life Sci. 59, 604–614. doi: 10.1007/s11427-016-5039-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fehr, A. R., and Perlman, S. (2015). Coronaviruses: an overview of their replication and pathogenesis. Methods Mol. Biol. 1282, 1–23. doi: 10.1007/978-1-4939-2438-7_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, X. Y., Li, J. L., Yang, X. L., Chmura, A. A., Zhu, G., Epstein, J. H., et al. (2013). Isolation and characterization of a bat SARS-like coronavirus that uses the ACE2 receptor. Nature 503, 535–538. doi: 10.1038/nature12711

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerdts, V., and Zakhartchouk, A. (2017). Vaccines for porcine epidemic diarrhea virus and other swine coronaviruses. Vet. Microbiol. 206, 45–51. doi: 10.1016/j.vetmic.2016.11.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Halpin, K., Hyatt, A. D., Fogarty, R., Middleton, D., Bingham, J., Epstein, J. H., et al. (2011). Pteropid bats are confirmed as the reservoir hosts of henipaviruses: a comprehensive experimental study of virus transmission. Am. J. Trop. Med. Hyg. 85, 946–951. doi: 10.4269/ajtmh.2011.10-0567

PubMed Abstract | CrossRef Full Text | Google Scholar

Heymann, D. L., Mackenzie, J. S., and Peiris, M. (2013). SARS legacy: outbreak reporting is expected and respected. Lancet 381, 779–781. doi: 10.1016/s0140-6736(13)60185-3

CrossRef Full Text | Google Scholar

Hu, B., Ge, X., Wang, L. F., and Shi, Z. (2015). Bat origin of human coronaviruses. Virol. J. 12:221. doi: 10.1186/s12985-015-0422-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, B., Zeng, L. P., Yang, X. L., Ge, X. Y., Zhang, W., Li, B., et al. (2017). Discovery of a rich gene pool of bat SARS-related coronaviruses provides new insights into the origin of SARS coronavirus. PLoS Pathog. 13:e1006698. doi: 10.1371/journal.ppat.1006698

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Liu, W. J., Xu, W., Jin, T., Zhao, Y., Song, J., et al. (2016). A Bat-derived putative cross-family recombinant coronavirus with a reovirus gene. PLoS Pathog. 12:e1005883. doi: 10.1371/journal.ppat.1005883

PubMed Abstract | CrossRef Full Text | Google Scholar

Lau, S. K., Woo, P. C., Li, K. S., Huang, Y., Tsoi, H. W., Wong, B. H., et al. (2005). Severe acute respiratory syndrome coronavirus-like virus in Chinese horseshoe bats. Proc. Natl. Acad. Sci. U.S.A. 102, 14040–14045. doi: 10.1073/pnas.0506735102

PubMed Abstract | CrossRef Full Text | Google Scholar

Leroy, E. M., Kumulungui, B., Pourrut, X., Rouquet, P., Hassanin, A., Yaba, P., et al. (2005). Fruit bats as reservoirs of Ebola virus. Nature 438, 575–576. doi: 10.1038/438575a

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Shi, Z., Yu, M., Ren, W., Smith, C., Epstein, J. H., et al. (2005). Bats are natural reservoirs of SARS-like coronaviruses. Science 310, 676–679. doi: 10.1126/science.1118391

PubMed Abstract | CrossRef Full Text | Google Scholar

Lim, Y. X., Ng, Y. L., Tam, J. P., and Liu, D. X. (2016). Human Coronaviruses: a review of virus-host interactions. Diseases 4:26. doi: 10.3390/diseases4030026

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, G., Wang, Q., and Gao, G. F. (2015). Bat-to-human: spike features determining ‘host jump’ of coronaviruses SARS-CoV, MERS-CoV, and beyond. Trends Microbiol. 23, 468–478. doi: 10.1016/j.tim.2015.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, I. M., and Arden, K. E. (2015). MERS coronavirus: diagnostics, epidemiology and transmission. Virol. J. 12:222. doi: 10.1186/s12985-015-0439-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Madhugiri, R., Fricke, M., Marz, M., and Ziebuhr, J. (2014). RNA structure analysis of alphacoronavirus terminal genome regions. Virus Res. 194, 76–89. doi: 10.1016/j.virusres.2014.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Mari Saez, A., Weiss, S., Nowak, K., Lapeyre, V., Zimmermann, F., Dux, A., et al. (2015). Investigating the zoonotic origin of the West African Ebola epidemic. EMBO Mol. Med. 7, 17–23. doi: 10.15252/emmm.201404792

PubMed Abstract | CrossRef Full Text | Google Scholar

Omrani, A. S., Al-Tawfiq, J. A., and Memish, Z. A. (2015). Middle East respiratory syndrome coronavirus (MERS-CoV): animal to human interaction. Pathog. Glob. Health 109, 354–362. doi: 10.1080/20477724.2015.1122852

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Shea, T. J., Cryan, P. M., Cunningham, A. A., Fooks, A. R., Hayman, D. T., Luis, A. D., et al. (2014). Bat flight and zoonotic viruses. Emerg. Infect. Dis. 20, 741–745. doi: 10.3201/eid2005.130539

PubMed Abstract | CrossRef Full Text | Google Scholar

Paez, D. J., Giles, J., Mccallum, H., Field, H., Jordan, D., Peel, A. J., et al. (2017). Conditions affecting the timing and magnitude of Hendra virus shedding across pteropodid bat populations in Australia. Epidemiol. Infect. 145, 3143–3153. doi: 10.1017/S0950268817002138

PubMed Abstract | CrossRef Full Text | Google Scholar

Phan, M. V. T., Ngo Tri, T., Hong Anh, P., Baker, S., Kellam, P., and Cotten, M. (2018). Identification and characterization of Coronaviridae genomes from Vietnamese bats and rats based on conserved protein domains. Virus. Evol. 4:vey035. doi: 10.1093/ve/vey035

PubMed Abstract | CrossRef Full Text | Google Scholar

Sendow, I., Ratnawati, A., Taylor, T., Adjid, R. M., Saepulloh, M., Barr, J., et al. (2013). Nipah virus in the fruit bat Pteropus vampyrus in Sumatera, Indonesia. PLoS One 8:e69544. doi: 10.1371/journal.pone.0069544

PubMed Abstract | CrossRef Full Text | Google Scholar

Simas, P. V., Barnabe, A. C., Duraes-Carvalho, R., Neto, D. F., Caserta, L. C., Artacho, L., et al. (2015). Bat coronavirus in Brazil related to appalachian ridge and Porcine epidemic diarrhea viruses. Emerg. Infect. Dis. 21, 729–731. doi: 10.3201/eid2104.141783

PubMed Abstract | CrossRef Full Text | Google Scholar

Sola, I., Almazan, F., Zuniga, S., and Enjuanes, L. (2015). Continuous and discontinuous RNA synthesis in coronaviruses. Annu. Rev. Virol. 2, 265–288. doi: 10.1146/annurev-virology-100114-055218

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, S., Wong, G., Shi, W., Liu, J., Lai, A. C. K., Zhou, J., et al. (2016). Epidemiology, genetic recombination, and pathogenesis of coronaviruses. Trends Microbiol. 24, 490–502. doi: 10.1016/j.tim.2016.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, X. C., Zhang, J. X., Zhang, S. Y., Wang, P., Fan, X. H., Li, L. F., et al. (2006). Prevalence and genetic diversity of coronaviruses in bats from China. J. Virol. 80, 7481–7490. doi: 10.1128/jvi.00697-06

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Brand, J. M., Smits, S. L., and Haagmans, B. L. (2015). Pathogenesis of Middle East respiratory syndrome coronavirus. J. Pathol. 235, 175–184. doi: 10.1002/path.4458

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D., Fang, L., and Xiao, S. (2016). Porcine epidemic diarrhea in China. Virus Res. 226, 7–13. doi: 10.1016/j.virusres.2016.05.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, S., Masangkay, J. S., Nagata, N., Morikawa, S., Mizutani, T., Fukushi, S., et al. (2010). Bat coronaviruses and experimental infection of bats, the Philippines. Emerg. Infect. Dis. 16, 1217–1223. doi: 10.3201/eid1608.100208

PubMed Abstract | CrossRef Full Text | Google Scholar

Woo, P. C., Lau, S. K., Lam, C. S., Lau, C. C., Tsang, A. K., Lau, J. H., et al. (2012). Discovery of seven novel Mammalian and avian coronaviruses in the genus deltacoronavirus supports bat coronaviruses as the gene source of alphacoronavirus and betacoronavirus and avian coronaviruses as the gene source of gammacoronavirus and deltacoronavirus. J. Virol. 86, 3995–4008. doi: 10.1128/JVI.06540-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Woo, P. C., Wang, M., Lau, S. K., Xu, H., Poon, R. W., Guo, R., et al. (2007). Comparative analysis of twelve genomes of three novel group 2c and group 2d coronaviruses reveals unique group and subgroup features. J. Virol. 81, 1574–1585. doi: 10.1128/jvi.02182-06

PubMed Abstract | CrossRef Full Text | Google Scholar

Woo, P. C. Y., Huang, Y., Lau, S. K. P., and Yuen, K. Y. (2010). Coronavirus genomics and bioinformatics analysis. Viruses 2, 1804–1820. doi: 10.3390/v2081803

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z., Li, Y., Ren, X., He, G., Zhang, J., Jian, Y., et al. (2016a). Deciphering the bat virome catalog to better understand the ecological diversity of bat viruses and the bat origin of emerging infectious diseases. ISME J. 10, 609–620. doi: 10.1038/ismej.2015.138

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z., Yang, L., Ren, X., Zhang, J., Yang, F., Zhang, S., et al. (2016b). ORF8-related genetic evidence for Chinese horseshoe bats as the source of human severe acute respiratory syndrome coronavirus. J. Infect. Dis. 213, 579–583. doi: 10.1093/infdis/jiv476

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z., Lu, L., Du, J., Yang, L., Ren, X., Liu, B., et al. (2018). Comparative analysis of rodent and small mammal viromes to better understand the wildlife origin of emerging infectious diseases. Microbiome 6:178. doi: 10.1186/s40168-018-0554-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Z., Ren, X., Yang, L., Hu, Y., Yang, J., He, G., et al. (2012). Virome analysis for identification of novel mammalian viruses in bat species from Chinese provinces. J. Virol. 86, 10999–11012. doi: 10.1128/jvi.01394-12

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Wu, Z., Ren, X., Yang, F., He, G., Zhang, J., et al. (2013). Novel SARS-like betacoronaviruses in bats, China, 2011. Emerg. Infect. Dis. 19, 989–991. doi: 10.3201/eid1906.121648

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Wu, Z., Ren, X., Yang, F., Zhang, J., He, G., et al. (2014). MERS-related betacoronavirus in Vespertilio superans bats, China. Emerg. Infect. Dis. 20, 1260–1262.

Google Scholar

Yang, X. L., Hu, B., Wang, B., Wang, M. N., Zhang, Q., Zhang, W., et al. (2015). Isolation and characterization of a novel bat coronavirus closely related to the direct progenitor of severe acute respiratory syndrome coronavirus. J. Virol. 90, 3253–3256. doi: 10.1128/JVI.02582-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X. L., Tan, C. W., Anderson, D. E., Jiang, R. D., Li, B., Zhang, W., et al. (2019). Characterization of a filovirus (Mengla virus) from Rousettus bats in China. Nat. Microbiol. 4, 390–395. doi: 10.1038/s41564-018-0328-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, G., Cowled, C., Shi, Z., Huang, Z., Bishop-Lilly, K. A., Fang, X., et al. (2013). Comparative analysis of bat genomes provides insight into the evolution of flight and immunity. Science 339, 456–460. doi: 10.1126/science.1230835

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., Fan, H., Lan, T., Yang, X. L., Shi, W. F., Zhang, W., et al. (2018). Fatal swine acute diarrhoea syndrome caused by an HKU2-related coronavirus of bat origin. Nature 556, 255–258. doi: 10.1038/s41586-018-0010-9

PubMed Abstract | CrossRef Full Text | Google Scholar