Coexistence of multiple coronaviruses in several bat colonies in an abandoned mineshaft

Coronaviruses are divided into the genera Alphacoronavirus, Betacoronavirus, Gammacoronavirus, and a proposed Deltacoronavirus genus(Woo et al., 2012a), and belong to the subfamily Coronavirinae, family Coronaviridae, of the order Nidovirales. Alphacoronaviruses infect pigs(Kusanagi et al., 1992), humans(Macnaughton and Madge, 1978; van der Hoek et al., 2004), cats(Pedersen et al., 1984), bats(Poon et al., 2005), and other mammals. Betacoronaviruses infect humans(Bruckova et al., 1970; Peiris et al., 2003; Woo et al., 2005; Zaki et al., 2012), mice(Mizutani, 2001), bats(Li et al., 2005), and other mammals. Gammacoronaviruses mostly infect birds though some infect sea mammals(Mihindukulasuriya et al., 2008; Woo et al., 2012a). The proposed deltacoronaviruses primarily infect birds but are also reported to infect mammals(Woo et al., 2012a).

Human coronaviruses used to have low pathogenicity and cause mild respiratory symptoms. However, two coronaviruses, severe acute respiratory syndrome corona-virus(SARS-CoV)(Peiris et al., 2003) and Middle East respiratory syndrome coronavirus(MERS-CoV), can cause severe respiratory diseases(Zaki et al., 2012). Phylogenetic evidence indicates that SARS-CoV and MERS-CoV originated from bats(van Boheemen et al., 2012; Ge et al., 2013; Lau et al., 2013). Since the 2002-2003 outbreak of SARS in China, coronaviruses have been discovered in 71 bat species throughout Asia, America, Africa, and Europe(Tang et al., 2006; Lau et al., 2012; Woo et al., 2012b; Anthony et al., 2013; De Benedictis et al., 2013; Geldenhuys et al., 2013; Goes et al., 2013; Graham et al., 2013; Lelli et al., 2013; Shi, 2013). Among the 15 established alpha-and betacoronav-irus species, eight were identified in bats in China, including Miniopterus bat coronavirus 1(BtCoV 1), Miniopterus bat coronavirus HKU8(BtCoV HKU8), Rhinolophus bat coronavirus HKU2(BtCoV HKU2), Scotophilus bat coronavirus 512(BtCoV 512), SARS-related or SARS-like CoVs(SL-CoV; WIV1, Rp3, HKU3, etc.), Pipistrellus bat coronavirus HKU5(BtCoV HKU5), Rousettus bat coronavirus HKU9(BtCoV HKU9), and Tylonycteris bat coronavirus HKU4(BtCoV HKU4)(Lau et al., 2005; Li et al., 2005; Woo et al., 2006; Lau et al., 2007; Chu et al., 2008; Woo et al., 2009; Cavanagh and Britton, 2012; Woo et al., 2012b). Viruses similar to the above bat coronavirus species or strains(BtCoV 1, HKU2, HKU9, 512, HKU8, and SL-CoV)have been found in southeast Asia and Europe(Watanabe et al., 2010; Gouilh et al., 2011; Wacharapluesadee et al., 2015; Drexler et al., 2010; Lelli et al., 2013).

Considering the species diversity and geographical distribution of bats, it is likely that novel bat coronaviruses remain to be discovered. Furthermore, continuous surveillance of bat coronaviruses will provide insights into the evolution of coronaviruses as well as help in evaluating the potential for interspecies transmission of these viruses from bats to humans and other animals(Drexler et al., 2011). In this study, we conducted a two-year RT-PCR surveillance of bat coronaviruses in an ab and oned mineshaft in Mojiang County, Yunnan Province, China. We found a high infection rate with a genetically diverse array of coronaviruses in different bat species in the mineshaft.

Bat fecal swabs were collected in August and September of 2012 and in April and July of 2013 as described previously(Li et al., 2005). The collected samples were kept in viral transport medium(1 × Hank's balanced salt solution, 1% bovine albumin, pH 7.4, 15 μL/mL amphotericin B, 100 units/mL penicillin G and 50 ug/mL streptomycin) and stored at -80 ℃ before processing. The bat species were first identified based on morphology. PCR amplification and sequencing of cytochrome b(Cytb)or NADH dehydrogenase subunit 1(ND1)(Irwin et al., 1991; Mayer and von Helversen, 2001)from DNA extracted from the fecal swabs was used to confirm the bat species. Viral RNA was extracted from 140 μL of the fecal swab samples with Viral RNA Mini Kit(Qiagen, Westburg, The Netherl and s)per the manufacturer's instructions. RNA was eluted in 60 μL of AVE(RNase-free water with 0.04% sodium azide)buffer, aliquoted, and stored at -80 ℃.

One-step RT-PCR(Invitrogen, San Diego, USA)was used for the amplification of a 440-bp fragment targeting the RNA-dependent RNA polymerase gene(RdRp)of all known alpha-and betacoronaviruses(de Souza Luna et al., 2007). Amplification of an 816-bp fragment extending the 440 bp was performed using published methods(Drexler et al., 2010). Spike(S)genes were amplified using degenerate primers designed based on the alignment of known coronaviruses(sequences provided upon request). St and ard precautions were taken to avoid PCR contamination, and no false positives were observed in negative controls. PCR products were gel-purified and sequenced on the ABI PRISM^® 3100 Genetic Analyzer(Applied Biosystems, Foster City, USA). The sequencing chromatograms were carefully inspected for overlapping multicolor peaks, which are an indicator of sequence heterogeneity in the amplicons. The PCR products of corresponding samples were cloned and 10 clones of each sample were sequenced to confirm sequence heterogeneity. To avoid PCR contamination, all positive samples were checked through two independent PCRs by different experimenters. The positive samples detected in this study were named using the abbreviated bat species name plus the bat sample number abbreviation. For example, a virus detected from Rhinolophus sinicus in sample number 4017 was named RsBtCoV/4017. If the bat was co-infected by two different coronaviruses, numbers were appended to the sample names, such as RsBtCoV/4017-1 and RsBtCoV/4017-2.

Preliminary sequence management and analysis were carried out using Geneious(version 5.5.9, Biomatters Limited, New Zeal and), sequence alignment and editing were performed using ClustalW(version 2.0) and BioEdit(version 7.1.9)(Hall, 1999; Larkin et al., 2007). Potential in vitro recombinant sequences(i.e., PCR artifacts)were screened for and discarded using Recombination Detection Program v2.0(Martin et al., 2005). Obtained consensus sequences were compared with known sequences of the coronavirus RdRp genes in GenBank. Phylogenetic trees were constructed using MrBayes(version 3.2) under the GTR model and the Maximum Likelihood algorithm in Geneious(Huelsenbeck and Ronquist, 2001).

In total, 276 bats(83 in August 2012; 97 in September 2012; 52 in April 2013, and 44 in July 2013) were sampled in a mineshaft in Mojiang(Figure 1). Six bat species were identified based on morphology and confirmed by sequencing of the Cytb or ND1 gene(Table 1). Bat species or most closely related species were R. sinicus, Rhinolophus affinis, Hipposideros pomona, Miniopterus schreibersii, Miniopterus fuliginosus, and Miniopterus fuscus. Full-length Cytb and ND1 gene sequences were deposited in GenBank under accession numbers KP876547 to KP876557.

Figure 1.  Sample collection site in Mojiang, Yunnan Province. The mineshaft where bat samples were collected is indicated with a red dot, three county towns (Mojiang, Zhenyuan, and Ning'er) near the mineshaft are labeled with green dots.

Among the collected 276 fecal samples, 138(50%)were positive for coronavirus: 45.7%(37/81) in August and 74.7%(74/99) in September 2012, and 46.2%(24/52) in April and 6.8%(3/44) in July 2013. All six bat species showed high infection rates, from 35% to 73%(Table 1). R. sinicus, R. affinis, and M. schreibersii were observed to have the highest infection rates in September 2012. Due to the small sample size, it was not possible to accurately assess changes in the infection rate for H. pomona, M. fuliginosus, and M. fuscus.

From the 138 positive samples, 152 RdRp partial coronavirus sequences(approximately 400 bp)were obtained, indicating co-infections of two viruses. Two sequences(HiBtCoV/3740-2 and RaBtCoV/4991) were homologous to betacoronaviruses, all other 150 sequences were homologous to alphacoronaviruses, including the established bat coronavirus species BtCoV 1, BtCoV HKU2, and BtCoV HKU8, unassigned BtCoV HKU7 and BtCoV HKU10, and unclassified alphacorona-virus(Table 2, Figure 2). To better underst and the phylogenetic relationships between these bat coronaviruses, we selected 12 samples representing different branches of the tree and extended the sequence of the partial RdRp screening fragment to 816 bp(Figure 3). The partial RdRp sequences obtained in this study were submitted to GenBank under accession numbers KP876505 to KP876546 and KU343189 to KU343200.

Figure 2.  Phylogenetic analysis of bat coronaviruses detected in this study based on the partial RdRp gene sequences. The partial RdRp gene sequences (~370 bp) obtained in this study were aligned with those of known representative coronaviruses and used to construct the phylogenetic tree by MrBayes V3.2 under assumption of GTR model, using 1, 000, 000 trees sampled every 100 steps, annotated with a burn-in of 90% using TreeAnnotator V1.0.6. Detected coronaviruses were classified into seven alphacoronavirus groups and two betacoronavirus lineages separated with horizontal lines. Coronaviruses co-infecting an individual are indicated by sequential number after the name of each strain. All the detected coronaviruses here are in bold and co-infection coronaviruses are labeled with a black triangle to the left of the name.

Among the 150 RdRp partial sequences that were homologous to alphacoronaviruses, fifty-one were related to BtCoV 1, thirty-seven to BtCoV HKU2, one to BtCoV HKU7, fifty-one to BtCoV HKU8, four to BtCoV HKU10, and four to unclassified alphacoronaviruses(Table 2). Most sequences shared 87%-98% nucleotide(nt)identity and 95%-100% amino acid(aa)identity with the established bat coronavirus species, suggesting that some of these coronaviruses represent different strains of known viral species. Only a few sequences were related to unclassified bat coronaviruses. One of two sequences detected in a R. affinis bat(RaBtCoV/4307-1) was related to an unclassified alphacoronavirus Hipposideros BtCoV Ratcha-67/THA/2007 detected in Thail and, with 83% nt and 92% aa identities, based on the 816-bp sequence(Figure 3)(Gouilh et al., 2011). Three sequences(one in M. fuliginosus and two in R. sinicus)were related to another unclassified alphacoronavirus Cardioderma BtCoV Kenya/KY43/2006 detected in Kenya, with 80% nt and 93% aa identities(Tong et al., 2009).

Figure 3.  Phylogenetic analysis of coronaviruses based on the amino acid sequences of the partial (816-bp) RdRp gene. Sequences obtained in this study were aligned with those of known representative coronaviruses and used to construct the phylogenetic tree using MrBayes. Coronaviruses co-infecting an individual are indicated with a sequential number after the name of each strain. The aa identities between coronavirus sequences detected in this study and the closest reference sequence are labeled following the virus names. The coronaviruses detected in this study are in bold and co-infecting viruses are labeled with a black triangle to the left of the sample name.

Only two sequences detected in this study were homologous to betacoronaviruses. One of them(RaBtCoV/4991) was detected in a R. affinis sample and was related to SL-CoV. The conserved 440-bp RdRp fragment of RaBtCoV/4991 had 89% nt identity and 95% aa identity with SL-CoV Rs672(Yuan et al., 2010). In the phylogenetic tree, RaBtCoV/4991 showed more divergence from human SARS-CoV than other bat SL-CoVs and could be considered as a new strain of this virus lineage(Figure 2). Another sequence(HpBtCoV/3740-2) was detected in a H. pomona sample and shared the closest similarity(81% and 89.1% identities at nt and aa levels, respectively, based on the 816-bp sequence)with a Hipposideros bat coronavirus Hipposideros/GhanaBoo/348/2008 discovered in Ghana(Quan et al., 2010). These two strains, together with Zaria bat ZBCoV(host Hipposideros commersoni) and Thail and bat BtCoV/B56054(host Hipposideros larvatus)(Tong et al., 2009; Wacharaplue-sadee et al., 2015), formed an independent lineage distantly related to SL-CoVs(Figure 2, 3). The partial RdRp sequences of this lineage shared less than 90% aa identity compared to the closest CoVs in lineage B and thus, this lineage may represent a novel betacoronavirus species.

Co-infection with different coronavirus species was found in 14 out of 138 positive specimens. Co-infection events were detected in all six bat species(Table 2, Figure 2). The most frequent co-infections were with the BtCoV 1B and BtCoV HKU8 lineages, as found previously by Chu et al. (2008). The other co-infections involved BtCoV HKU8 or BtCoV HKU2 and the unclassified bat alphacoronaviruses. Most interestingly, we also found a co-infection with coronavirus species from two different genera, BtCoV HKU10 and an unclassified betacoronavirus in one H. pomona bat(no. 3740).

We detected high interspecies host diversity, with the same coronavirus species infecting different bat species within and /or across families. BtCoV 1 was detected in three out of six bat species of two families(Vespertilionidae and Rhinolophidae), including M. schreibersii, M. fuliginosus, and R. sinicus. BtCoV HKU2 was detected in two closely related species of the same family, R. sinicus, and R. affinis. BtCoV HKU8 was detected in five out of six bat species of two families, including M. schreibersii, M. fuliginosus, M. fuscus, R. sinicus, and R. affinis.

To further address interspecies infections of coronaviruses in bats, we amplified 8 full-length S genes from BtCoV 1-positive samples including M. schreibersii, M. fuliginosus, and R. sinicus. The full-length S genes ranged from 4, 125 to 4, 137 bp(Supplementary Table S1). In the phylogenetic tree, MfulBtCoV/3759-1, MsBtCoV/4068-1, MfulBtCoV/3736-1, MsBtCoV/4001-1, MfulBtCoV/3709, and RsBtCoV/3716 clustered with BtCoV 1B, and MsBtCoV/4056-1 and MsBtCoV/3710 clustered with BtCoV 1A(Figure 4). Sequence comparison showed that the S gene of RsBtCoV/3716, derived from R. sinicus, showed higher divergence from other strains among the BtCoV 1B group(Supplementary Table S1). Interestingly, the S genes of coronaviruses detected from the same bat family displayed high divergence. The S gene sequences obtained in this study were deposited in GenBank under accession numbers KU343201 to KU343208.

Figure 4.  Phylogenetic analysis of bat coronaviruses of BtCoV 1 based on full-length S genes. Eight full-length S gene sequences obtained in this study were aligned with 4 reference S genes and used to construct the phylogenetic tree using MrBayes. Sequences obtained in this study are in bold.

Bats and rodents are both considered important natural reservoirs of zoonotic viruses(Luis et al., 2013). A previous study reported the discovery of paramyxovirus(MojV)in rats(Rattus flavipectus)in the same mineshaft in Mojiang(Wu et al., 2014). Here, focusing on another host, we performed a surveillance of coronaviruses in bats. In contrast to the sole strain MojV detected in rats, we found highly divergent alphacoronavirus strains and two betacoronavirus strains in bats from the same mineshaft. Among the eight known bat coronavirus species, four were detected in the cave. In addition, a new putative coronavirus species was detected in a Hipposideros bat.

Most coronaviruses first replicate in epithelial cells of the respiratory and /or enteric tract. Unapparent enteric infection of animals maintains the virus in the population(Lai et al., 2007). Previous studies have indicated that bats host a large number of viruses, but seldom display clinical symptoms under natural or experimental infection conditions(Baker et al., 2013). In this study, we conducted an 18-month surveillance of bat coronaviruses in an ab and oned mineshaft cohabited by six bat species across three families. High infection rates of bat coronaviruses, particularly alphacoronaviruses, were detected in the four sampling periods. More interestingly, we observed a high rate of co-infection with two corona-virus species and interspecies infection with the same coronavirus species within or across bat families. These phenomena may be owing to the diversity and high density of bat populations in the same cave, facilitating coronavirus intra-and interspecies transmissions, which may result in recombination and acceleration of coronav-irus evolution. Alphacoronaviruses were found in all six collected bat species. The most frequently detected sequences were related to BtCoV 1, BtCoV HKU2, and BtCoV HKU8, BtCoV HKU10, BtCoV HKU7, and two unclassified alphacoronaviruses were less frequently detected. As previously reported, BtCoV 1 and BtCoV HKU8 are highly prevalent in Miniopterus species(Chu et al., 2006; Woo et al., 2006; Chu et al., 2008), while BtCoV HKU2 mainly infects Rhinolophus species(Woo et al., 2006; Lau et al., 2007). Due to the overlapping host range and cohabitation, high rates of co-infection with BtCoV 1 and BtCoV HKU8 in individual bats were detected in Miniopterus species. As BtCoV 1B were detected mainly in Miniopterus bats, it could be speculated that RsBtCoV/3716 may have been a spillover from Miniopterus bats and adapted to R. sinicus through S gene mutation.

Most viral sequences detected had high nt/aa identities with the established coronavirus species, indicating that these viruses had relatively stable genetic evolution. Several identified sequences were homologous to unclassified alphacoronaviruses previously detected in Thail and and Kenya, suggesting that these viruses represent novel groups of alphacoronaviruses(Tong et al., 2009; Gouilh et al., 2011; Wacharapluesadee et al., 2015). Further sequencing of full-length genomes will be helpful to determine their taxonomic positions.

Members of the genus Betacoronavirus are divided into four lineages: A, B, C, and D(previously termed groups 2A, 2B, 2C, and 2D). Lineage A includes two human coronaviruses, hCoVs OC43 and HKU1, which are typically associated with common colds. Lineage B includes the highly pathogenic human SARS-CoV and the diverse bat SL-CoVs(WIV1, SHC014, Rp3, HKU1, and others). Lineage C includes the highly pathogenic human MERS-CoV and several bat MERS-CoV-related coronaviruses(BtCoVs HKU4 and HKU 5). Lineage D only contains coronaviruses discovered in bats(BtCoV HKU9). In this study, we detected a SL-CoV-related sequence in R. affinis. This strain is distantly related to the previously discovered bat SL-CoVs in other Rhinolophus species and represents a new strain of SL-CoVs. Additionally, we detected a new betacoronavirus(RaBtCoV/3740-2) in a H. pomona sample that formed a distinct lineage in Betacoronavirus along with several other coronaviruses found in Hipposideros bats(Tong et al., 2009; Quan et al., 2010; Gouilh et al., 2011; Wacharaplu-esadee et al., 2015). Previous studies have indicated that these coronaviruses in Hipposideros bats are SARS-related and should be considered members of betacoronav-irus lineage B(Tong et al., 2009; Quan et al., 2010; Gouilh et al., 2011). Based on the sequence comparison and phylogenetic analysis in this study, we suggest that these coronaviruses found in Hipposideros bats of different countries may represent a novel betacoronavirus lineage(tentatively, lineage E). Considering that the two highly pathogenic human coronaviruses(SARS-CoV and MERS-CoV)in this genus originated from bats(Ge et al., 2013; Lau et al., 2013; Corman et al., 2014), attention should be particularly paid to these lineages of bat coronaviruses.

This work was jointly funded by the National Natural Science Foundation of China(81290341) and Scientific and Technological Basis Special Project(2013FY113500), and China Mega-Project for Infectious Disease(2014ZX10004001-003) from the Minister of Science and Technology of the People's Republic of China, and USNIAID(R01AI110964). We are very grateful to Adam Nitido for help in language editing.

The authors declare that they have no conflict of interest. This study was approved by the Ethics Committee of the Yunnan Institute of Endemic Disease Control and Prevention. All institutional and national guidelines for the care and use of animals were followed.

ZLS designed and coordinated the study. XYG, YZZ, JHZ, CML, XLY, LJW, BW, YZ, and ZXL collected samples. NW, XYG, WZ, BH, and BL performed molecular studies. XYG and NW analyzed the data. XYG and ZLS drafted the manuscript. All authors read and approved the final manuscript.

Supplementary Table S1 is available on the website of Virologica Sinica: www.virosin.org; link.springer.com/journal/12250.

Comparison of 8 S genes obtained here with 4 reference sequences of BtCoV 1.