Genetic diversity of coronaviruses in bats in Lao PDR and Cambodia.

South-East Asia is a hot spot for emerging zoonotic diseases, and bats have been recognized as hosts for a large number of zoonotic viruses such as Severe Acute Respiratory Syndrome (SARS), responsible for acute respiratory syndrome outbreaks. Thus, it is important to expand our knowledge of the presence of viruses in bats which could represent a risk to humans. Coronaviruses (CoVs) have been reported in bat species from Thailand, China, Indonesia, Taiwan and the Philippines. However no such work was conducted in Cambodia or Lao PDR. Between 2010 and 2013, 1965 bats were therefore sampled at interfaces with human populations in these two countries. They were tested for the presence of coronavirus by consensus reverse transcription-PCR assay. A total of 93 samples (4.7%) from 17 genera of bats tested positive. Sequence analysis revealed the presence of potentially 37 and 56 coronavirus belonging to alpha-coronavirus (αCoV) and beta-CoV (βCoV), respectively. The βCoVs group is known to include some coronaviruses highly pathogenic to human, such as SARS-CoV and MERS-CoV. All coronavirus sequences generated from frugivorous bats (family Pteropodidae) (n=55) clustered with other bat βCoVs of lineage D, whereas one coronavirus from Pipistrellus coromandra fell in the lineage C of βCoVs which also includes the MERS-CoV. αCoVs were all detected in various genera of insectivorous bats and clustered with diverse bat αCoV sequences previously published. A closely related strain of PEDV, responsible for severe diarrhea in pigs (PEDV-CoV), was detected in 2 Myotis bats. We highlighted the presence and the high diversity of coronaviruses circulating in bats from Cambodia and Lao PDR. Three new bat genera and species were newly identified as host of coronaviruses, namely Macroglossus sp., Megaerops niphanae and Myotis horsfieldii.

Introduction

Coronaviruses (CoVs) are enveloped, single stranded positive-sense RNA viruses displaying a large genome of 26 kb to 32 kb (Masters, 2006). They are classified within 4 groups based on genotype and serological characteristics: Alphacoronaviruses (αCoVs), Betacoronaviruses (βCoVs), Gammacoronaviruses, and the recently created group of the Deltacoronaviruses (Woo et al., 2009, Woo et al., 2012). Their host range is very wide and includes both mammalian and avian species. Coronaviruses can cause acute and chronic respiratory, enteric, neurological and hepatic diseases in their hosts (Weiss and Navas-Martin, 2005). Coronaviruses of animal origin were responsible for the Severe Acute Respiratory Syndrome (SARS) outbreak in 2003–2004, which was associated with deaths in Hong Kong, China, South East Asia and North America (Centers for Disease Control and Prevention (CDC), 2003; Peiris et al., 2003) and the current epidemics of MERS in the Arabian Peninsula (Alsahafi and Cheng, 2016) and Korea (Choi, 2015).

Bats have been identified as natural reservoirs for several zoonotic viruses, such as henipaviruses (Young et al., 1996, Chua et al., 2000, Halpin et al., 2000) lyssaviruses variants (Banyard et al., 2014) and Ebola virus (Leroy et al., 2005, Leroy et al., 2009). Bats have been identified as the natural host of the SARS-CoV, (Wang et al., 2006), and recently, the NeoCoV from the clade the Middle Eastern Respiratory Syndrome (MERS-CoV) was detected in a sub-Saharan bat (Neoromicia capensis) besides its camel host (Corman et al., 2014). A growing number of coronaviruses have been detected in bats since the SARS-CoV outbreak (Chu et al., 2006, Chu et al., 2008, Lau et al., 2007, Lau et al., 2010, Lau et al., 2012, Watanabe et al., 2010, Gouilh et al., 2011, Tsuda et al., 2012, Wacharapluesadee et al., 2013, Wacharapluesadee et al., 2015, Anindita et al., 2015, Xu et al., 2016, Chen et al., 2016, Kim et al., 2016) including a high diversity of coronaviruses, recently detected in five Thai provinces neighboring Cambodia (Wacharapluesadee et al., 2013, Wacharapluesadee et al., 2015).

The order Chiroptera represents approximately 20% of all living mammal species (Teeling et al., 2005). Over 25% of the worlds bat diversity is found in South-East Asia, established in many natural, urban and suburban environments (Kingston, 2013). Seventy species of bats have been described so far in Cambodia whereas ninety species are known in Lao PDR, including Yangochiroptera and Yinpterochiroptera (Matveev, 2005, Sarak et al., 2013). Apart from Singapore, bats are hunted for food or preparation of traditional medicines and are found in food markets throughout South-East Asia, despite bats are protected by the law in these countries (Lee et al., 2014, Mildenstein et al., 2016). Bat farms, where artificial roosts are erected to facilitate bat guano harvest to serve as agricultural fertilizer, are becoming common in South-East Asia, including in Cambodia (Thi et al., 2014) In Thailand, coronaviruses belonging to the lineages B and C of betacoronaviruses were detected in bat guano (Gouilh et al., 2011, Wacharapluesadee et al., 2013). Moreover, several studies evidenced the presence of coronavirus belonging to the betacoronavirus group in close areas (He et al., 2014, Anindita et al., 2015, Wacharapluesadee et al., 2015)

Due to evolving land-use such as deforestation, infrastructure development, urban development, and agricultural expansion, bat populations are settling in areas closer to human dwellings (Jung and Threlfall, 2016), increasing the likelihood of contact between bats and humans. Socio-economic-driven changes of the environment are also impacting the bats and thus may affect virus biodiversity (Looi and Chua, 2007, Turmelle and Olival, 2009, Brierley et al., 2016). These factors may also contribute to an increased occurrence of contacts between bats and humans, which in turn, may increase the potential transmission of zoonotic pathogens, including a coronaviruses. Despite the prolificacy of such interfaces in the region, so far no study has been conducted on bat coronaviruses in Cambodia and Lao PDR. To our knowledge, this is the first study to investigate the presence and the diversity of coronaviruses in bats in Cambodia and Lao PDR.

Materials and methods

Ethics

The study was approved by the National Veterinary Research Institute and Forest Administration department of the Ministry of Agriculture Forestry and Fisheries in Cambodia, as well as by the National Animal Health Laboratory of the Ministry of Agriculture Forestry and Fisheries, in Lao PDR, and under the Institutional Animal Care and Use Committee at the University of California, Davis (protocol number: 16048). Since Cambodia and Lao PDR have no ethic committee overseeing animal experimentation, animals were treated in compliance with the guidelines of the American Society of Mammalogists and within the framework of the European Union legislation guidelines (Directive 86/609/EEC).

Collection of bat samples

Bat samples from various species were collected from different locations in Lao PDR and Cambodia, over a 3-year period (from November 2010 to December 2013). Bat species were identified to genus or species by trained biologists and veterinarians. The sampling was performed in two phases: Phase 1 was carried out in 2010 by the Institut Pasteur in Cambodia (IPC) and the Muséum national d'Histoire naturelle (MNHN; Paris, France) in cooperation with Cambodian government partners. Bats were captured and humanely euthanized in full compliance with local ethical and legal guidelines. The capture sites were located in Ratanakiri, Stung Treng and Preah Vihear provinces in Cambodia. They were chosen because of their easier accessibility to the collection team. Rectal swabs were stored in viral transport medium solution (VTM; containing tryptose phosphate Broth 2.95%, 145 mM of NaCl, 5% gelatin, 54 mM Amphotericin B, 106 U of penicillin-streptomycine per liter, 80 mg of gentamycine per liter [Sigma-Aldrich, Irvine, UK]). Tissue specimens (i.e. lung, liver, spleen, kidney, heart) were placed in separate cryotubes. All specimens were immediately transferred into liquid nitrogen containers before being transported to the Institut Pasteur laboratory where they were stored at − 80 °C prior to testing. The second sampling phase was performed by the Wildlife Conservation Society (WCS) from 2011 to 2013. During this phase 2, no animals were killed and samples were collected at four key interfaces for increased contact (and potential disease transmission) between bats and humans: fresh-food markets; freshly hunted in rural communities for subsistence or sale; wild meat restaurants; and bat guano farms. In Lao PDR the majority of samples were obtained from freshly dead bats found in food markets, with a small number collected from bats captured by subsistence hunters; in Cambodia the majority of samples were collected at wild meat restaurants where bats were butchered, prepared and served, with additional samples obtained from freshly trapped bats (live and dead) held by hunters and middle-men in rural communities, and at a bat guano farm. In this study, sterile swabs were used to collect freshly voided fecal samples from tarpaulins placed under the guano farm roosts. Rectal and oral swabs and tissue samples were also collected from individual animals that had died of natural causes and were found fresh beneath bat guano farm roosts. Oral and rectal swabs were collected opportunistically from fresh dead and live bats at the market, hunting and restaurant locations. Swab and tissue samples were placed in separate cryovials in VTM and immediately stored on liquid nitrogen in dewars for transport to the laboratory where they were stored in − 80 °C freezers until testing.

RNA extraction and nested-RT-PCR

Viral RNA was extracted using the QIAamp Viral RNA Mini Kit (Qiagen, Hilden, Germany) according to the supplier's instructions. Reverse transcription (RT) was performed using SuperScript III (Invitrogen, San Diego, CA). The PCR mixture (final volume: 25 μl) contained 2 μl of cDNA, PCR buffer (50 mM Tris-HCl (pH 9.0); 50 mM NaCl; 5 mM MgCl2), 200 μM (each) deoxynucleoside triphosphates (dNTPs), 20 pmol of each primer targeting the RdRp gene (adapted from Watanabe et al., 2010) (Table 1 ), and 1 U of HOT FIREPol® DNA Polymerase (Solis BioDyne, Tartu, Estonia). The PCR mixture was incubated at 95 °C for 12 min, followed by 35 cycles at 95 °C for 30 s, 50 °C for 30 s, and 72 °C for 1 min, and by a final extension at 72 °C for 5 min. A nested PCR amplification using forward primer CoV-Fwd2 and the same reverse primer, was performed on 1 μl of the primary PCR products, using the same amplification conditions. The amplification of sequences specific to CoVs was attested by the visualization of a 440 bp and a 434 bp fragment after the first and second PCR round, respectively. To limit the risk of contamination, RNA extraction, reverse transcription-PCR (RT-PCR), nested-PCR and gel electrophoresis were carried out in separate rooms. In addition, negative controls (water) were included in each run of the nested-RT-PCR assay and results were validated only if these controls tested negative while the positive controls (plasmids prepared by cloning the gene of interest) had to test positive. Amplification of longer fragments of the RdRp gene (i.e. 1370 bp) was performed on cDNA of 38 positive samples, by a nested RT-PCR using gene-specific primers designed by multiple alignments of known CoV sequences from the same clusters. Amplified product were sequenced in both directions by direct Sanger sequencing in commercial facilities (Macrogen, Inc., Seoul, Korea). Sequence data were deposited in GenBank and accession numbers ranging from KX284902 to KX520662 and KY010629 to KY010666 are provided in Supplementary Table 1. Since many coronaviruses from different bat genera were detected, and for the sake of clarity, a short name of the strain, the sample code, and the host classification were abbreviated to be used in the sequence nomenclature (Supplementary Table 2). For example, the sequence “Bat coronavirus 512/2005/PREDICT-KHP13-PTR1-0109”, corresponding to the Bat coronavirus 512/2005 strain (BatCoV-512), detected in the sample KHP13-PTR1-0109 (P109), from a Scotophilus kuhlii (Sku) was abbreviated as: BatCoV-512_P109_Sku. The abbreviations used to code the bat classification, are listed in the Supplementary Table 2.

Table 1

Primers used for the detection of CoVs.

Gene	Round	Primer name	Sequence 5′-3′
RdRp	1st	CoV-Fwd1	GGTTGGGAYTAYCCHAARTGTGA
RdRp	1st and 2nd	CoV-Rvs2	CCATCATCASWYRAATCATCATA
RdRp	2nd	CoV-Fwd2	GAYTAYCCHAARTGTGAYAGAGC

Sequence analysis

Sequences were assembled and analyzed using CLC Genomics Workbench version 3.6.1 and BioEdit, version 7.0.9.1. Sequences were aligned with a representative set of CoV sequences retrieved from GenBank using Seaview, version 4.5.4 (Gouy et al., 2010). Phylogenetic trees based on RNA sequences were constructed using the Maximum Likelihood method with the GTR + G + I model and bootstrap values (BP) were calculated after 1000 replicates. The best evolutionary model was determined using MEGA version 6.06. Trees displaying protein sequences were constructed using the neighbor joining method and bootstrap values were calculated after 1000 replicates. Phylogenetic and molecular evolutionary analyses were conducted using Seaview version 4.5.4.

Geographic data

Land cover data was obtained from GlobeLand30 service operated by the National Geomatics Center of China (NGCC, 2014). Initial data was produced in 2010 with an update in 2014. Images used for GlobeLand30 (GLC30) classification were multispectral images with a 30-meter resolution. Six classes of land cover were considered: cropland, forest, grassland, wetland, water bodies and human settlement area. Land cover structure for each sampling location of sampling is described in Supplementary Table 3. Data mapping was conducted with Quantum GIS, version 2.8.2.

Results

Sampling location, land cover use and hosts

A total of 1965 bats were sampled in 44 locations from nine provinces in Cambodia and eight provinces in Lao PDR (Supplementary Fig. 1). The characteristics of the environments around sampling site are described in Supplementary Table 3. Bats originating from 5 districts (i.e. Kasi and Vang Vieng districts in Lao PDR and Choam Khsant, Kean Svay, and Moung Russei districts in Cambodia) represented 60% of all the bats collected. Sampled bats were wild, live captured during phase 1 (n  = 322), where as in phase 2 (n  = 1643) they were: hunted by villagers for local consumption (n  = 455); for sale in markets (n  = 791); being butchered and prepared for sale at wild meat restaurants (n  = 392) or on a bat guano farm (n  = 5) (Supplementary Table 3). Apart from bats trapped during phase 1 or animals sampled on the guano farm (site reference C6) during phase 2, the exact location where the animals were captured remained imprecise. However, the chiroptera were always captured by hunters in areas close to the location where the animals were sampled, i.e. mostly in areas at the border of deep forests, in mixed agricultural zones with sparse forests, in suburban zones close to sparse forest, in naturals protected forest areas, in places close to water surfaces or in limestone karst areas with mountain forests (Supplementary Fig. 1, Supplementary Table 3). Landscape analysis around each sampling location led to a typology comprising four main groups of land cover composition (Supplementary Table 3, Supplementary Fig. 1): 1) Deep forest area, well isolated from human settlements (sites: C3, C25, C29, C30, C32, C33, C34, C41, C42 and C43), 2) Isolated pockets of croplands surrounded by forest, very often forest is under fragmentation (sites: C2, C12, C16, C18, C19, C21, C22, C23, C24, C26, C27, C28, C30, C35, C36, C37, C38, C40 and C44), 3) Forest edge, mixed agricultural zones with sparse forests (sites: C5, C7, C13 and C14) and 4) typically agricultural zones, villages and suburban zones, often close to highly fragmented forests (sites: C1, C4, C6, C8, C9, C10, C11 and C20). For sites C1, C4, C6, C8 and C9 a key parameters was the presence of wetlands in the immediate vicinity.

Detection of coronaviruses

Coronavirus RNA was detected in 93 bats (4.7%) out of 1965 animals tested. Coronavirus RNA was detected in 21 of the 44 sites with detection rates varying depending upon sites (Table 2 ). However, it is not possible to run an in depth statistical analysis owing to the fact that sampling was not developed for that purpose and data are therefore biased. Detailed data are presented in Supplementary Table 3. Nevertheless, in highly transformed areas (mostly croplands and urbanized areas), bat biodiversity was lower while the detection rate of coronavirus was higher (Table 2, Supplementary Table 3). The phylogenetic analysis based on 320-bp nucleotides of the 93 sample sequences and 34 reference sequences from GenBank indicated that 37 sequences belonged to the αCoV genus while the 56 others fell into the βCoV genus (Fig. 1 ). Out of the 56 bat-βCoVs, 55 belonged to the lineage D. The remaining strain detected (PREDICT-CoV-34_GT1-3_Pisp) belonged to the lineage C which comprises MERS-CoV and MERS-CoV-related viruses. The 55 lineage D bat-βCoVs clustered into four subclusters (D1, D2, D3 and D4) out of the five subclusters identified in this group D (Fig. 1). The subcluster D1 (BP = 100) comprised 7 different βCoVs, 2 sequences detected in pteropodid bats from Lao PDR and 5 from Cambodia. The subcluster D2 (BP = 80), contained 17 βCoVs detected from the Cambodian provinces of Preah Vihear, Stung Treng and Battambang. The subcluster D3 (BP = 45) contained 24 bat-βCoVs, which were mostly obtained from animals in Lao PDR (n  = 22) while only 2 were detected in bats from Cambodia. The subcluster D4 comprised the remaining bat-βCoVs (BP = 78), detected in bats from Lao PDR only. Bat-αCoVs were unevenly distributed within three subclusters (Fig. 1). The αCoV subcluster 1 included 32 sequences from Cambodian bats. The αCoV subcluster 2 included two αCoVs from Cambodia, while the αCoV subcluster 4 cluster contained one sequence from Cambodia and two from Lao PDR. The same phylogenetic clustering was seen with analyses using the amino-acid sequences. A similar topology was obtained when analyzing the protein sequence (Supplementary Fig. 2). The phylogenetic analysis based on 38 longer sequences of 1370 bp (Fig. 2 ) and the corresponding 454 amino acid sequence (Supplementary Fig. 3) confirmed the tree topology described on shorter sequences in Fig. 1 and Supplementary Fig. 2, respectively. Coronavirus RNA was detected in bats belonged to 10 distinct genera. The proportion of CoVs-positive bats significantly varied depending upon the bat genus (Pearson's Chi2, p  < 0.0001). The highest detection rates were found for the long-tongue fruit bats (Macroglossus sp.: 20%), the nectar cave bat (Eonycteris spelaea: 11.2%) and the lesser Asian bat (Scotophilus sp.: 9.3%). The 37 αCoV-positive bats belonged to five genera of the Yangochiroptera suborder, i.e. Rhinolophus (n  = 1), Myotis (n  = 2), Scotophilus (n  = 31), Pipistrellus (n  = 1) and Hipposideros (n  = 2). A total of 56 animal tested positive for a βCoV, out of which 55 belonged to five genera of Yinpterochiroptera suborder, family Pteropodidae: i.e. Rousettus (n  = 25), Cynopterus (n  = 15), Eonycteris (n  = 10), Macroglossus (n  = 4), Megaerops (n  = 1) whereas the bat infected by the PREDICT-CoV-34 strain was a Yangochiroptera of the genus Pipistrellus (Fig. 2). All αCoV RNA were detected in insectivorous bats while most of the βCoV RNA were detected in fruit bats.

Table 2

Detection rates of coronavirus RNA in bats and geographical origin of the positive animals.

Bat	Total no. of samples (percentage (%) of positive)	Origin of positive samplesaNo. of samples (no. of positive)			CambodiaNo. of samples (no. of positive)	Site	Lao PDRNo. of samples (no. of positive)	Site
		H	W	G
F/family
Genus or species
Emballonuridae	148 (0)				148 (0)		0 (0)
Taphozous sp.	148 (0)	118 (0)	30 (0)-	–	148 (0)	C9, C1, C32, C33	0 (0)
Hipposideridae	62 (3.2)				4 (0)		58 (2)
Aselliscus sp.	7 (0)	7 (0)	–	–	0 (0)		7 (0)	C25
Hipposideros sp.	55 (3.6)	51 (2)	4 (0)	–	4 (0)	C31, C33	51 (2)	C19, C14(+), C21, C23, C24
Megadermatidae	21 (0)				21 (0)		0 (0)
Megaderma sp.	21 (0)	–	21 (0)	–	21 (0)	C32, C33, C41	0 (0)
Pteropodidae	1124 (4.9)				416 (25)		708 (30)
Cynopterus sp.	341 (4.4)	183 (8)	158 (7)	–	318 (15)	C5(+), C7, C10(+), C12(+), C11, C29, C30, C32(+), C34, C35, C36, C37(+), C38(+), C39(+), C40, C41, C42, C44	23 (0)	C21, C23, C24, C26, C27, C28
Eonycteris spelaea	89 (11.2)	61 (6)	28 (4)	–	28 (4)	C37(+), C38(+), C39(+), C40	61 (6)	16, C21(+), C23(+), C24(+)
Macroglossus sp.	28 (14.3)	28 (4)	–	–	21 (4)	C7, C10(+)	7 (0)	C21, C24
Megaerops niphanae	130 (0.8)	122 (1)	8 (0)	–	16 (1)	C5, C7, C29, C30, C34, C39	114 (0)	C13, C14, C15, C16, C20 C21 C23, C24, C26, C27
Megaerops sp.	12 (0)	–	12 (0)	–	12 (0)	C36, C37	0 (0)
Pteropus sp.	10 (0)	10 (0)	–	–	10 (0)	C11	0 (0)
Rousettus sp.	514 (4.9)	506 (24)	8 (1)	–	11(1)	C10, C11, C33(+), C34, C37	503 (24)	C14, C15(+), C16(+), C17, C20, C21(+), C22, C23(+), 24(+)
Rhinolophidae	154 (0.7)				52 (1)		102 (0)
Rhinolophus sp.	154(0.7)	102 (0)	52 (1)	–	52 (1)	C31, C33(+), C37, C41, C42, C43	102 (0)	C14, C21, C23, C25
Vespertilionidae	456 (7.7)				418 (35)		38 (0)
Harpiocephalus sp.	1 (0)	1 (0)	–	–	0 (0)		1 (0)	C20
Ia io	32 (0)	32 (0)	–	–	0 (0)		32 (0)	C18
Myotis horsfieldii	50 (4)	50 (2)	–	–	50 (2)	C4(+)	0 (0)
Myotis ricketti	5 (0)	5 (0)	–	–	0 (0)		5 (0)	C20
Pipistrellus coromandra	29 (6.9)	29 (2)	–	–	29 (2)	C2(+), C3(+)	0 (0)
Scotophilus sp.	338 (9.3)	333 (30)	–	5(1)	338 (31)	C1(+), C6(+), C8, C9(+)	0 (0)
Tylonycteris sp.	1 (0)	–	1 (0)	–	1 (0)	C35	0 (0)
Total	1965 (4.7)	1638 (79)	322 (13)	5(1)	1059 (61)		906 (32)

Bat families are shown in bold.

Sites where bats tested positive for coronavirus are in bold and marked with (+).

H: hunted bats, sold in markets of restaurants; W: wild bats, caught in their natural environment; G: bats collected in bat guano farms.

Fig. 1

Phylogenetic analysis of coronavirus sequences based on a 320 nucleotides fragment on the RdRp gene.

The tree was constructed using a Maximum Likelihood method with the GTR + G + I model. Bootstrap values (BP) were calculated after 1000 replicates. Tree is rooted to an avian coronavirus (FJ376622). The sequences detected in this study are shown with a bullet. Accession numbers related to the sequences are presented in Supplementary Table 1.

The tree was constructed using a Maximum Likelihood method with the GTR + G + I model. Bootstrap values (BP) were calculated after 1000 replicates. Tree is rooted to an avian coronavirus (FJ376622). The sequences detected in this study are shown with a bullet. Accession numbers related to the sequences are presented in Supplementary Table 1.

Supplementary Fig. 2

Distribution of coronavirus sequences based on a 115 amino acid fragment encoded by the RdRp gene.

The tree was constructed by neighbor joining and bootstrap values were determined by 1000 replicates. The tree is rooted to an avian coronavirus (FJ376622). The sequences detected in this study are shown with a bullet. Accession numbers related to the sequences are presented in Supplementary Table 1.

Fig. 2

Sequence analysis of 1370 nucleic acid fragments of the RdRp gene of coronavirus sequences.

The tree was constructed using a Maximum Likelihood method with the GTR + G + I model. Boostrap values (BP) were calculated after 1000 replicates. Tree is rooted to an avian coronavirus (FJ376622). The sequences detected in this study are shown with a bullet. Accession numbers related to the sequences are presented in Supplementary Table 1.

The tree was constructed using a Maximum Likelihood method with the GTR + G + I model. Boostrap values (BP) were calculated after 1000 replicates. Tree is rooted to an avian coronavirus (FJ376622). The sequences detected in this study are shown with a bullet. Accession numbers related to the sequences are presented in Supplementary Table 1.

Supplementary Fig. 3

Distribution of coronavirus sequences based on a 115 amino acid fragment encoded by the RdRp gene.

The tree was constructed by neighbor joining and bootstrap values were determined by 1000 replicates. The tree is rooted to an avian coronavirus (FJ376622). Bootstrap values (BP) were calculated after 1000 replicates. The sequences detected in this study are shown with a bullet. Accession numbers related to the sequences are presented in Supplementary Table 1.

Phylogenetic clustering of CoVs according to host and location

1965 bats were from 17 different genera and 5 families. Eighty percent of the samples belonged to two bat families only, Vespertilionidae (n  = 456) and Pteropodidae (n  = 1124). In Lao PDR, 95.7% of the samples belonged to the Yinpterochiroptera suborder, including families Pteropodidae (78.1%), Hipposideridae (6.3%), and Rhinolophidae (11.3%) while only 4.3% were Yangochiroptera suborder from the family Vespertilionidae. In Cambodia, 35% of the bats collected belonged to suborder Yinpterochiroptera, and suborder Yangochiroptera accounted for the remainder of the collected samples (75%), including families Emballuronidae and Vespertilionidae. Only sequences from fruit bats from Cambodia and Lao PDR from the genera Rousettus, Eonycteris, Macroglossus, Cynopterus and Megaerops were found in the lineage D of βCoVs. Subcluster D1 contained two sequences from Macroglossus sp. from Cambodia (PREDICT-CoV-22 strains) and five virus sequences from Eonycteris spelaea from Lao PDR and Cambodia (PREDICT-CoV-22, R91, R77, R74, R58). Seventeen sequences were found in Cynopterus (PREDICT-CoV-24 strains, R96, R75, R72, R65, R59, R71), Rousettus sp. and Megaerops niphanae (PREDICT-CoV-24) collected in Preah Vihear and Battambang in 2010 and 2013, respectively. They fell into subcluster βCoV-D2 and displayed 98% to 100% of amino acid similarity with CoV sequences previously detected in Cynopterus sphinx and Hipposideros lekaguli in South East Thailand and China, respectively (Wacharapluesadee et al., 2015, Xu et al., 2016). Similarly, the βCoV-D3 subcluster comprised 24 sequences found in Rousettus, Eonycteris and Macroglossus and CoVs previously identified in pteropodids from Hong Kong, Kenya, Thailand and Indonesia (genera Rousettus and Dobsonia) (Lau et al., 2010, Wacharapluesadee et al., 2015, Anindita et al., 2015). They displayed 96.9% to 100% of amino acid identity. Seven sequences detected in Rousettus bats from Lao PDR and displaying 100% of amino acid identity, formed their own branch in subcluster D4. Additionally, the PREDICT-CoV-34 strain, detected in Pipistrellus coromandra fell into lineage C of βCoV (BP = 100) which contained sequence from MERS-CoV as well as sequences detected in Pipistrellus from Hong Kong, Neoromicia from South Africa and Myotis from China (Annan et al., 2013, Woo et al., 2006, Corman et al., 2015, Corman et al., 2014, Xu et al., 2016). The PREDICT-CoV-34 strain showed the highest similarity to the bat coronavirus isolate JPDB144 recently described in Myotis daubentonii in China (89.5% nucleic acid identity and 95% amino acid identity) (Xu et al., 2016). Sequences from the subcluster 1 displayed similar traits and were only detected in Cambodian Scotophilus (n  = 31), with the exception of one sequence detected in Cambodian Pipistrellus (BatCoV512_SL2-9_Pisp). The αCoV subcluster 4 contained sequences related to the strain HKU10, detected in various bat genera in Thailand and Hong Kong (Gouilh et al., 2011, Lau et al., 2012, Wacharapluesadee et al., 2015). In this subcluster, two sequences found in Hipposideros larvatus (PREDICT-CoV-53 strains) from Lao PDR formed their own branch (subcluster αCoV_4c), supported by a bootstrap of 100. They were related to the αCoV BatCoV Ratcha-67 (96.9% amino acid identity) previously detected in the same bat genus in Thailand (Gouilh et al., 2011). Additionally, the sequence BatCoV-25 PA201 detected in Rhinolophus shameli, showed the highest similarity with a sequence detected in the same host species in Thailand (BatCoV B554) (Wacharapluesadee et al., 2015). Interestingly, two αCoV sequences detected from Myotis horsfieldii in Cambodia, PEDV_Kch1-04 and PEDV_Kch1-13, fell with the CoV responsible for porcine epidemic diarrhea episodes in pigs (PEDV), in the αCoV subcluster 2. They displayed 93.4% to 97.2% nucleotide identity and 96.13% to 100% amino acid identity with the porcine CoV strain PEDV. All sequences of viruses identified in Scotophilus fell in the subcluster 1 of αCoVs, those from Hipposideros belonged to the subcluster αCoV_4c, the sequences found in Myotis were all from the subcluster 2 of αCoV and the sequences of the viruses identified in Rousettus were all grouped in the lineage βCoV_D. In Pipistrellus bats from Mondulkiri province caught in 2013, the viruses detected belonged to both the αCoV and the βCoV groups.

Discussion

In this study, we report for the first time, the detection and description of coronaviruses in chiroptera from Cambodia and Lao PDR. This report is of importance for public health as these countries host a high biodiversity of bats and the interface at which humans and chiroptera are in close proximity. Bats are used for food and reared as a source of guano. Searching for bat-borne viruses known to be potentially pathogenic for humans such as coronaviruses is therefore of high importance. A significant diversity of CoVs was found in pteropodids both in Cambodia and Lao PDR. Overall, CoV detection rates (6.5% in Lao PDR and 4.85% in Cambodia) were in the same range as those found in Thailand (7%) (Wacharapluesadee et al., 2015), but lower than those reported in Hong Kong (12%) or Philippines (29.6%) (Tsuda et al., 2012, Woo et al., 2006). Interestingly, for the first time, coronaviruses were detected in Megaerops niphanae and Myotis horsfieldii.

In this study, αCoVs were mostly associated with Yangochiroptera bats whereas and βCoVs were found in Yinpterochiroptera bats, with the exception of three viruses, i.e., 2 CoVs from the subcluster 4 of found in hipposiderids, and one βCoV from lineage C, detected in Pipistrellus. In Lao PDR, βCoVs from lineage D were all associated with pteropodids whereas αCoVs were detected from hipposiderids. In Cambodia, an association between host and virus clade was observed between bats from the genus Myotis and the αCoVs genetically related to PEDV strains. These results are in line with conclusions of previous studies: all CoVs detected in bats belonging to the Myotis genus were always αCoVs (Tang et al., 2006, Woo et al., 2006, Dominguez et al., 2007, Gloza-Rausch et al., 2008, Osborne et al., 2011, August et al., 2012, Kemenesi et al., 2014, Fischer et al., 2016). βCoVs from lineage D have frequently been found in frugivorous bat species from Madagascar, Kenya, Thailand, and Hong Kong (Woo et al., 2007, Anindita et al., 2015, Razanajatovo et al., 2015, Wacharapluesadee et al., 2015, Xu et al., 2016). A similar trend was observed in this study, as βCoVs from the lineage D were only affiliated with frugivorous bats.

Although only part of the samples could yield a longer sequence, phylogeny based on these sequences was similar to that based on short sequences. This suggests that the phylogeny based on short sequences is reliable. Unfortunately sequencing of coronavirus can be challenging, due in part, to limited nucleic acid in field samples, as well as the high genetic diversity of the viruses (Drexler et al., 2010, King et al., 2012). Although some findings in this study suggest the possibility of host specificity, in particular for the genera Hipposideros and Myotis (associated to αCoVs from the subclusters 4c and 2 respectively), we also report shared hosts from different families for both αCoV and βCoV. For example, Pipistrellus bats from the same location were found to harbor αCoV from the subcluster 1 and lineage C βCoVs sequences. Coronaviruses from the latter lineage, which also includes the highly pathogenic MERS-CoV (Reusken et al., 2016), have been detected in humans, camels, insectivorous bats (Vespertilionidae) and frugivorous bats (Phyllostomidae) (Reusken et al., 2010, Annan et al., 2013, Wacharapluesadee et al., 2013, Corman et al., 2014, Wang et al., 2014, Munster et al., 2016). There have also been found differences in the percentage of seropositive between Rhinolophus ferrumequinum and Myotis myotis in studies of lyssavirus of European bat colonies (Serra-Cobo et al., 2013).

Pipistrellus coromandra is a synanthropic bat and the presence of CoVs in the same clade with pathogenic viruses raises the question of the potential risk for human health resulting from deforestation and urbanization that creates habitats for these bats. Pipistrellus bats, like Myotis, comprise species with differing habitats. However, in this work, Pipistrellus coromandra was found in deforested and agricultural regions, confirming thus its synanthropic behavior in the study area considered. Further risk-assessment studies should focus on the correlation between landscape change, land use and deforestation that may affect the distribution of human-dwelling bats and therefore the coronaviruses they harbor. This risk in Cambodia might be worsened by the fact that Pipistrellus bats are hunted for food. Another related risk might be the development of guano farms, also associated with agriculture (Stibig et al., 2007, Broadhead and Izquierdo, 2010). Bats reared for guano in Thailand have shown to harbor lineage C βCoV (Wacharapluesadee et al., 2013). This activity is also developing in Cambodia (Thi et al., 2014) and the perception of risk from workers is weak. Furthermore, hunters and restaurant workers at sampling sites described being bitten by bats and were exposed to urine and feces (and in some cases blood). Bats on Cambodian guano farms are wild and free-ranging, with farmers constructing artificial roosts adjacent to their homes in order to attract bats. Farmers, who wear no protective equipment or clothing, collect guano voided onto tarpaulins or nets laid beneath the roosts and are regularly urinated and defecated upon and sometimes bitten. The risk of contamination by direct contact through urine and feces, or aerosols must therefore be considered.

The hypothesis of a potential case of horizontal transmission between livestock and bats can be raised with the PEDV-like viruses, detected in Myotis horsfieldii. Another αCoV closely related to PEDV-CoV was recently detected in Brazil in Mexican free-tailed bats (Tadarida brasiliensis) (Simas et al., 2015). Not only is this first report of the presence of this coronavirus in Myotis horsfieldii, but the strain detected seems genetically closely related to PEDV strains that infect swine and cattle (Song and Park, 2012). Myotis bats also belong to the Vespertilionidae family and dwell in dark places including houses, farms and barns. Investigating the presence of pig farms in areas where coronavirus are detected, and the circulation of Myotis bats in the surroundings, would provide data to explore potential relations between bats and livestock. The hypothesis of a possible origin of PEDV from bats as well as a potential cross-species transmission has been raised by previous studies (Huang et al., 2013, Tang et al., 2006) and would benefit from further investigation. To date, only a limited number of studies have investigated the capacity of bat coronaviruses to be infectious to other mammals, including humans. This is notably due to a lack of data on the spike protein of these bat viruses. Indeed, the spike protein is the primary determinant for the cell tropism and pathogenesis (Belouzard et al., 2012).

Our study had some limitations related to sampling procedure that may have affected results. Many of the samples (n  = 1838) were collected from dead bats intended for human consumption. If rectal swab is the sample which represents the highest probability for detecting coronaviruses in bats (Watanabe et al., 2010), it is possible that virus survival in dead animals might be affected. The health status of the animals sampled in our study was unknown but previous studies suggest that bats might not develop disease during coronavirus infections. RNA viruses seem to have little pathogenic effect on bat's life cycle (Li et al., 2005) which may explain that bats are excellent reservoirs for zoonotic viruses including CoVs (Omatsu et al., 2007, Brook and Dobson, 2015, Han et al., 2015). Seasonality has been shown for some zoonotic virus infection rates in Chiroptera. Reproduction periods, female status and resource availabilities have been proven to affect the prevalence in bats infected by other RNA viruses (Middleton et al., 2007, Plowright et al., 2008, Wacharapluesadee et al., 2010, George et al., 2011, Amman et al., 2012, Hayman, 2015, Amengual et al., 2007). Greater shedding over such periods increase the probability of transmission of viruses to humans and the risk of emerging zoonoses (Serra-Cobo et al., 2013). Seasonality was not explored in our study due to the heterogeneity of the sampling session over time, but it would be interesting to include these parameters in further investigations. Another limitation in this work is the lack of molecular identification of bats. The collection of bat samples for barcoding or Cyt identification was not part of the field procedure and thus no bat sample has been collected during the field work. The identification was therefore limited to the genus.

Almost all the bats which tested positive for coronaviruses were meant to be consumed by local human populations. The examples of Ebola or SARS-CoV outbreaks already suggested that wildlife hunting and consumption provided opportunities for human contamination (Bengis et al., 2004, Xu et al., 2004, Leroy et al., 2009). In Cambodia and Lao PDR, bat consumption is widespread (Mickleburgh et al., 2009, Lee et al., 2014, Mildenstein et al., 2016). These practices may increase the risk of human exposure to viruses through hosts that may be reservoirs for pathogens. Practices such as hunting, selling or cooking bats might represent efficient interfaces for virus transfer from bats to humans, and therefore need to be further investigated. It demonstrates the importance to develop guidance for rural communities exposed to bats, on how to deal with them and the potential virus threat.

The diversity of bat coronaviruses found circulating in Cambodia and Lao PDR suggests a correlation may exist between coronaviruses and host diversity which is at least for part the consequence of anthropogenic natural environment change, mostly deforestation. The natural landscapes surrounding the collection points correlated with the known biology and ecology of the bats sampled (Supplementary Table 3). Pteropodids and vespertilionids, known to roost in caves (Eonycteris or Rousettus, and Ia respectively), were collected close to karst areas in Lao PDR. In Cambodia, the Yinpterochiroptera bats belonged to 5 genera of pteropodids known to mostly live in trees (Cynopterus, Macroglossus, Megaerops, Pteropus, Rousettus). Forest coverage throughout Continental Southeast Asia has drastically evolved between 1990 and 2015 (Supplementary Fig. 4). Cambodia is the country where forest loss has been the most intensive with 20% of the forest surface lost since 1990, while about 50% of the remaining surface is fragmented and 70% of the forest is in a “perforation” state, associated with agricultural development (WWF, 2013; World Bank, 2016). In Lao PDR, the forest coverage looksstable (73% in 1990, 81% in 2015), but the structure of the forest has changed. A typical feature is the fragmentation of the natural forest and its replacement by cultivated areas. However, compared to other countries in the region, the forest coverage remains the highest in Lao PDR and shows signs of resilience. This is mostly due to the geology of the country which is mainly karstic and thus less useful for agriculture. For both countries the demographic growth remains a key point for natural environment evolution, forest loss and dynamic of urban and suburban expansion. This environmental change is affecting the biodiversity of bats and therefore that of their coronaviruses. Continued deforestation, agricultural expansion and suburban growth might facilitate encounters between humans and coronaviruses from human-dwelling bats that could potentially become harmful for humans, specific attention and studies should therefore be devoted to this aspect.

The following are the supplementary data related to this article.

Supplementary Fig. 1

Geographic distribution of samples in Lao PDR and Cambodia. Sampling sites are shown by triangles and labeled from C1 to C44 (Supplementary Table 2). The 21 sites where bats tested positive for coronaviruses are shown with a red triangle.

Supplementary Fig. 4

Forest cover change in South East Asia from 1990 to 2015. a. Map of forest cover change. b. Dynamic of forest cover change.

Supplementary Table 1

Accession numbers of sequences deposited in Genbank, and corresponding named used for sequence analysis.

Supplementary Table 2

Mammals and provinces abbreviations used for coding coronavirus samples.

Supplementary Table 3

Origin of the samples collected and description of the environment around the sampling sites.

Conflict of interest

Philippe Buchy is currently an employee of GSK vaccines.