Encephalomyocarditis virus is potentially derived from eastern bent-wing bats living in East Asian countries.

Bats are reservoir hosts of many zoonotic viruses and identification of viruses that they carry is important. This study aimed to use high throughput screening to identify the viruses in fecal guano of Taiwanese insectivorous bats caves in order to obtain more information on bat-derived pathogenic viruses in East Asia. Guano samples were collected from two caves in Taiwan, pooled, and then subjected to Multiplex PCR-based next generation sequencing for viral identification. Subsequently, encephalomyocarditis virus (EMCV) sequence was detected and confirmed by reverse transcription PCR. EMCV is considered as rodent virus and thus, animal species identification through cytochrome oxidase I (COI) barcoding was further done to identify the viral source. Finally, determination of distribution and verification of the presence of EMCV in guano obtained from Japanese and South Korean caves was also done. We concluded that the guano collected was not contaminated with the excrement of rodents which were reported and presumed to live in Taiwan. Also, EMCV genome fragments were found in guanos of Japanese and South Korean caves. It is possible that the eastern bent-wing bat (Miniopterus fuliginosus) is one of the natural hosts of EMCV in East Asia.

Introduction

Bats are considered natural reservoirs for several zoonotic viruses including lyssaviruses (Davis et al., 2013a, b; Albas et al., 2011; Almeida et al., 2011, 2005), henipaviruses (Roche et al., 2015; Breed et al., 2011; Rahman et al., 2010), coronaviruses (Anthony et al., 2017; Chen et al., 2016; Anindita et al., 2015; Annan et al., 2013; Anthony et al., 2013), and/or picornaviruses (Wu et al., 2012; Lau et al., 2011). Most viruses, which originate from bats and are pathogenic in humans and other mammals, appear to cause no clinical signs of disease in bats under natural or experimental infection (Calisher et al., 2006). Bats are also considered the origin of the Middle East respiratory syndrome coronavirus (MERS-CoV) that was first reported in Saudi Arabia (Wang et al., 2014), and severe acute respiratory syndrome coronavirus (SARS-CoV) that caused outbreaks in China in 2003 (Ge et al., 2016, 2013; Balboni et al., 2012; Wang et al., 2006). A novel coronavirus, which is closely related to SARS-CoV and uses Angiotensin Converting Enzyme 2 (ACE-2) receptor from humans, civets, and Chinese horseshoe bats for cell entry, was isolated from Chinese horseshoe bats suggesting these bats might be natural reservoirs of SARS-CoV (Ge et al., 2013).

In 2012, an outbreak of rabies occurred among ferret badgers in Taiwan (Chiou et al., 2016; Tsai et al., 2016; Chiou et al., 2014) after the country had been considered rabies-free since the last reported human case in 1959 and non-human animal in 1961 (Huang et al., 2015; Chiou et al., 2014). During 2016–2017, a putative new lyssavirus was found in the Japanese pipistrelle bat in the same country (Hu et al., 2018). SARS-related coronavirus and Alphacoronavirus were also detected in the country’s bat population (Chen et al., 2016). To obtain more information on bat-derived pathogenic viruses in East Asia, this study aimed to determine the viruses present in fecal guano of Taiwanese insectivorous bats caves.

Materials and methods

Sampling

Sampling was done in the two caves of Taiwan located in Taoyuan City, Dasi district: 24°50′35″N, 121°17′58″E (Taiwan 1) and New Taipei City, Ruifang district: 25°07′34″N, 121°49′59″E (Taiwan 2) wherein fecal guanos were collected randomly at several points in the cave and then pooled. During the process of viral identification, the encephalomyocarditis virus (EMCV) genome was found. To determine the distribution of the virus and confirm the presence of EMCV in fecal guanos, sampling was also performed in other countries including the main island of South Korea (36°39′33″N, 128°11′34″E); Jeju Island, South Korea (33°26′08″N, 126°50′15″E); Nishi-Yoshino, Nara, Japan (34°18′58″N, 135°46′28″E); Kawakami, Nara, Japan (34°17′02″N, 136°00′09″E); Sabae-shi, Fukui, Japan (35°56′54″N, 136°08′08″E); and Wakasa, Fukui, Japan (35°31′52″N, 135°50′22″E).

Primer design (Multiplex PCR)

The viral genera used for primer design (Supplementary Table 1) were selected in reference to the ICTV Master Species List (King et al., 2018). Thirty-seven genera including viruses possibly infectious to mammals and birds were selected. The data sets for primer design of each genus were compiled using viral strains specified by RefSeq on GenBank. The CoCoMo primer algorithms (Jimba et al., 2010) and CLC Genomics Workbench 6.5.1 (CLC bio) were used for the design. If there was only one RefSeq sequence per genus, another sequence of the same genus was added to the data set. The product size of the PCR was designed at a maximum of 1000 bp. All primers used for the Multiplex PCR are presented in Supplementary Table 1. Previously reported universal primers for Herpesvirus, Adenovirus, Orthoreovirus, and Bornavirus (Weissenbock et al., 2009; Wellehan et al., 2004; Spinner and Di Giovanni, 2001; Mizutani et al., 1999; VanDevanter et al., 1996) were also used.

Animal species identification

Identification of animal species was done from the guano through DNA barcoding based on cytochrome c oxidase I (COI) sequence using the specific primers (5′-GGTCAACAAATCATAAAGATATTGG, 5′-TAAACTTCAG GGTGACCAAAAAATCA).

Data sets for mapping

For viral identification, the sequences used for primer design were used as the data set for mapping reads from next generation sequencing (NGS). For animal species identification, sequences from a total of 75 species of mammals (Supplementary Table 2) that had been recorded and/or presumed to be present in Taiwan, and/or distributed worldwide were used as reference (Vincenot, 2017; Heaney and Molur, 2016; Ruedas et al., 2016; Challender et al., 2014; Lin and Motokawa, 2014; Bates et al., 2008; Hai Yin and Richardson, 2008; Hutson et al., 2008; Smith and Johnston, 2008; Wilson and Reeder, 2005).

Total nucleic acid extraction and cDNA synthesis

Total viral nucleic acid from pooled guano from each cave was extracted using the High Pure Viral Nucleic Acid Kit (Roche, Basel, Switzerland); half of the samples were used for RNA extraction by purification with the TURBO DNA-free™ Kit (Life Technologies, CA, USA) according to the manufacturer’s instruction. Complementary DNA (cDNA) was generated using random hexamers and the SuperScript™IV First-Strand Synthesis System (Thermo Fisher Scientific, California, USA) by incubating for 10 min at 23℃, 10 min at 50℃, and 10 min at 80℃.

Multiplex PCR

A primer mix was prepared for each well so that multiplex PCR using approximately 700 primers could be performed on one 96-well plate. Primers were mixed for each viral genus making the maximum number of primers in one well as 10, occupying 94 wells. The cDNA library was amplified using KOD-Multi & Epi-® (TOYOBO, Osaka, Japan) with initial denaturation at 94℃ for 10 min., followed by 40 cycles with denaturation at 98℃ for 10 s, annealing at 45℃ for 30 ss, and extension at 68℃ for 30 s, and a final extension at 72℃ for 10 min. PCR products were then collected from all wells (94 wells), mixed, and purified with Agencourt AMPure XP (Beckman Coulter, Brea, Calif., USA, USA) using 80% ethanol to remove primer dimers of less than 100 bp.

Next generation sequencing (NGS) and data analysis

A NGS library was first constructed according to the manufacturer’s protocol using the NEBNext®Ultra™ RNA Library Prep Kit for Illumina® (New England Bio Labs, Massachusets, USA). Sequencing of paired ends was performed with the MiSeq Reagent Kit v2 (300 cycles) (Illumina®, San Diego, CA, USA) wherein each read was output as a FASTQ format file from the MiSeq reporter and analyzed using CLC genomics Workbench 6.5.1 (Filgen, Nagoya, Japan). Commands for trimming sequences with low quality and commands for removing 20 bases corresponding to primers from the 5’-end of each read were executed. Trimmed reads were mapped to each sequence of the prepared data sets for viral identification (Supplementary Table 1) with a 50% length fraction and 80% similarity fraction. Reads with a minimum consensus length of 200 and a minimum mapping read number of 100 were selected, and the reference sequence that was mapped at the forward and reverse end of the predicted fragment between the primers’ binding site was selected and judged as positive. For animal species identification, NGS data that showed less than 50% similarity to the COI reference in consecutive full-length sequence of 5’ and 3′ ends of NGS data (less than 151 bp each sequence) were eliminated. Since only one COI genomic sequence per species was included in the reference despite the intraspecies variation, the maximum variation in COI barcoding was computed with the formula: similarity + [(100- similarity) × 0.11]. The maximum similarity was used as the threshold to determine the species.

Amplification of EMCV genomic fragments

When the EMCV genome was detected using NGS, reverse-transcription and PCR amplification of EMCV genes was done using SuperScript™IV First-Strand Synthesis System (Thermo Fisher Scientific, California, USA) first by incubating for 10 min at 23℃, 10 min at 50℃, and 10 min at 80℃. Amplification was carried out using the specific primers for the EMCV genome (Table 1 ), which were designed based on the highly conserved regions of EMCV (Liu et al., 2016; Duke et al., 1992), and KOD-Plus-Neo (Toyobo, Japan) at the following conditions: 2 min at 94℃; 35 cycles of 10 ss at 98℃, 30 ss at 58℃, and 30 ss at 68℃, with a final extension at 72℃ for 5 min.

Table 1

Primers for EMCV genome.

EMCV genes	Forward Primer (5`-3`)	Reverse Primer (5`-3`)	Amplicon size (bp)
Promoter (1st set)	TGAATGTCGTGAAGGAAGCAGT	ACCTCGACTAAACACATGT	315
Promoter (2nd set)	TGGAAGCTTCTTGAAGACA	AGATCAGATCCCATACAAT	240
1A	ATTCCACCTCCTCAGACAAGA	AGCTAGCAATGGAAGCATAT	206
2B (1st set)	TTCATGTTTAGACCAAGGAAACA	AGAGATTGTGGGAAACCGT	432
2B (2nd set)	AGACCCAAGGAGCGGCAGTGT	TGTCTTGAACTTAGCTGCTAT	361
3 A/B (1st set)	TCCAGTAGACGAGGTCAGT	TGAATGTCCAACAACTGCA	321
3 A/B (2nd set)	TCCAGCAGCTTAAAGCAAGACA	AACTCTCGCCGTCTCATTGTA	256

Results

Guanos from two caves in Taiwan were subjected to NGS using the designed primers and the reads were mapped to the viral sequences used for primer design. In two caves (Taiwan 1 and Taiwan 2), sequences matched to the reference were homologous to EMCV, enterovirus F strain, rabbit vesivirus, and Spanish goat encephalitis virus (Supplementary Table 3). To confirm the presence of these viral genomes, conventional RT-PCR followed by Sanger sequencing was performed using specific primers for EMCV (5’CACCATATTGCCGTCTTT, 5’CCACAACTATCCAACTCAC), enterovirus F strain (5’CCACGTGGCGCTAGTACC, 5’ACGGGGTACCGAAAGTAG), rabbit vesivirus (5’TCAAACCCAATTGGCATCTTC, 5’GCTTTGATGTAGTAGAATTGCCCA), and Spanish goat encephalitis virus (5’TGCGTGTACAACATGATGGGA, 5’CACGCAGTCATCACCACTGA). However, only the EMCV genomic sequences were amplified. To confirm the presence of EMCV, RT-PCR was performed using the specific primers for the EMCV genome (Table 1) wherein only the promoter region was amplified (Fig. 1 A). To eliminate the possibility that the EMCV genome is incorporated in the genome of rodents, Rhinolophinae, and Miniopterinae, the obtained sequence by RT-PCR was then subjected to BLAST wherein the result showed the absence of EMCV in the genome of these animals.

Fig. 1

Electrophoresis result of RT-PCR of EMCV genes. A, RT-PCR result of 4 EMCV genomic fragments using the primers for EMCV genome (Table 1): 1, promoter; 2, 1 A; 3, 2B; 4, 3 A/B; in samples collected from Taiwan 1 (Taoyuan city, Dasi district) and Taiwan 2 (New Taipei city) showing positive PCR-fragment for EMCV promoter after sequencing (arrows). B, RT-PCR result of EMCV promoter for samples collected in Korea (Main Island, South Korea), Jeju Is. (Jeju Island, South Korea), Nara 1 (Nishi-Yoshino, Nara, Japan), Nara 2 (Kawakami, Nara, Japan), Fukui 1 (Sabae-shi, Fukui, Japan), and Fukui 2 (Wakasa, Fukui, Japan); showing positive PCR-fragment for EMCV promoter after sequencing (arrow) in all the samples.

Rodents have been known to be the natural reservoir of EMCV. To determine if the guanos collected also contained the excreta of some rodents, identification of animal species was done through DNA barcoding based on the 5’-end and 3′-end regions of the mitochondrial gene, COI sequence. Using the next generation sequencer with CLC genomics software, the number of reads in Taiwan 1 and Taiwan 2 caves were 5’-end: 3256, 3′-end: 3526, and 5’-end: 3198, 3′-end: 3198, respectively. After elimination of reads whose similarity was less than 50% to the COI reference, the number of mapped reads was reduced to 5’-end: 2167, 3′-end: 1936 (Taiwan 1), and 5’-end: 880, 3`-end: 518 (Taiwan 2). Among these read sequences, a total of 189 kinds of consensus sequences were found. The 189 consensus sequences were mapped to 48 out of 76 reference species. Of the 48 reference species, 11 were from the order Chiroptera (Miniopterus fuliginosus, Myotis adversus, Myotis formosus, Nyctalus noctula, Murina puta, Tadarida teniotis, Rhinolophus Monoceros, Hipposideros armiger, Eptesicus serotinus, Rhinolophus Monoceros, Rhinolopus ferrumequinum), 19 from Rodentia (Rattus exulans, Rattus rattus, Rattus tanezumi, Rattus tiomanicus, Rhizomys pruinosus, Euchoreutes naso, Mesocricetus auratus, Eospalax baileyi, Petaurista alborufus, Callosciurus erythraeus, Tamiops maritimus, Apodemus agrarius, Mus caroli, Mus musculus, Eothenomys chinensis, Micromys minutus, Rattus losea, Demonys pernyi), three from Eulipotyphla (Condylura cristata, Mogera robusta, Hylomys suillus), five from Soricomorpha (Crocidura attenuata, Crocidura suaveolens, Suncus murinus, Chimarrogale himalayica, Soriculus fumidus), seven from Carnivora (Prionailurus bengalensis, Paguma larvata, Viverricula indica, Melogale moschata, Canis familiaris, Ursus thibetanus, Mustela nivalis), and three from Artiodactyla (Cervus nippon taiouanus, Cervus unicolor swinhoe, and Capricornis swinhoei) (Supplementary Table 3). Since the maximal intraspecific distances of COI is between 0.00–11.00% (Yassin et al., 2010), we considered 89.00% similarity, with the computed maximum variation, as the maximum threshold for identifying the species. Exactly 10/189 consensus sequences were above the threshold similarity to the reference species. Four of these, total of six reads in 5’ and 3′ ends, were 96.44%–100% matched to the reference animal species Miniopterus fuliginosus. These sequences were 97.60%–100% similar to the same species when subjected to BLAST (Supplementary Table 4). One (1 read) was 91.63% similar to Petaurista alborfus, two (1 read/consensus) were 96.71% similar to Prionailurus bengalensis, one (1 read) was 93.68% similar to Melogale moschata, and two (1 and 7 reads) were 91.28% similar to Mogera robusta. When subjected to BLAST for confirmation, all of them were 100% homologous to species of insects (Cochyllis hospes, Fosterina rustica, Alysiinae sp.) and a spider (Argyroides zonatus). The remaining consensus sequences (150/189) that were below 89.00% similarity, were also subjected to BLAST to determine the species, where 120 of 150 sequences showed ≥ 89.00% similarity to various species of insects and other species such as scorpions, fungi, gastropods, and even a nematode (Supplementary Table 5). Moreover, in 23 of the 150 sequences consensus was below 89.00%, with regard to both the reference species and to other species when BLASTed while 7 of the 150 species had no matching species when BLASTed and thus, these sequences could not be identified. All of the highest homologies of 23 sequences did not belong to vertebrates.

To further confirm the presence of the EMCV genome in the feces of insectivorous bats, guano samples were also collected in the bat caves in Japan and Korea. The samples from all of the caves examined contained EMCV promoter region (Fig. 1B), which was confirmed by sequencing. This region was reported to be 100% homologous to Chinese, Korean, and Belgian swine isolates; and other strains from chimpanzee and mice (Liu et al., 2016).

Discussion

This study showed that the EMCV genome was widely found in fecal guanos of East Asian countries, and the possibility of interfusion of rodents’ excretion to the guanos was very low. It is possible that Miniopterus fuliginosus is one of the natural reservoirs of EMCV.

EMCV, under the family Picornaviridae, genus Cardiovirus, has a single-stranded, positive sense RNA and the virus has worldwide distribution (Billinis et al., 2004; King et al., 2000). Although it is considered a rodent virus (Liu et al., 2017, 2016; Spyrou et al., 2004), it has wide-spectrum of hosts including monkeys and porcupines (Cardeti et al., 2016), elephants (Lamglait et al., 2015), wild boar (Liu et al., 2017, 2016), dogs (Luo et al., 2017), tapir (Vercammen et al., 2017), with pigs being the most sensitive (Liu et al., 2017; Maurice et al., 2016; Vansteenkiste et al., 2016). It causes myocarditis, encephalitis, neurological diseases, reproductive disorders, and diabetes (Luo et al., 2017; Carocci and Bakkali-Kassimi, 2012). As a potentially zoonotic virus, its discovery in these caves could pose a threat to animals and humans living nearby because this virus is resistant to a hostile environment and thus can remain infectious for days (Carocci and Bakkali-Kassimi, 2012).

Although the sample areas were mainly inhabited by Miniopterus fuliginosus and Rhinolophus ferrumequinum, detection of EMCV in fecal guanos could not conclude that the source of the viral genome was indeed these bats since the samples were not collected directly from their bodies. Guanos might be composed of excreta from rodents, the natural reservoir of EMCV, aside from bats. To determine the species that contributed most to the guano formation collected, animal species identification was done through NGS using the 1545-bp, at the 5’-end of the mitochondrial gene COI as the DNA barcode region. This COI gene locus, being easily amplified using universal primers, shows low intra-specific and high interspecific divergence (Dentinger et al., 2011), and possesses a higher range of phylogenetic signal than any other mitochondrial genes (Hebert et al., 2003a), and has been proven to be highly effective in identifying the species of various organisms (Andujar et al., 2018; Bingpeng et al., 2018; Changbunjong et al., 2018; Yassin et al., 2010; Hubert et al., 2008; Ward et al., 2005; Hebert et al., 2004, 2003b).

In both sample areas, a total of 189 consensus sequences were mapped to 48 out of 75 reference sequences wherein 10/189 were ≥ 89.00% similar to the reference species (Supplementary table 3). Four of these mapped sequences, containing two and four reads in Taiwan 1 and Taiwan 2, respectively, were 96.44%–100% matched to the species Miniopterus fuliginosus, and were 98.00%–100% similar to the same species when subjected to BLAST (Table 3). The majority of other sequences belonged to the order Diptera followed by Lepidoptera, Hymenoptera, Coleoptera, Trichoptera, Orthoptera, and Thysanoptera, which are the main diet of the insectivorous bats including Miniopterus fuliginosus (Hu et al., 2011; Lumsden and Bennett, 2005; Funakoshi and Takeda, 1998; Bernard and D’avidson, 1996). Also, it was unlikely that other unidentified sequences were derived from vertebrates. That is, the possibility that the guanos contained excreta of vertebrate other than Miniopterus fuliginosus was very low. The relative frequency of the reads from Miniopterus fuliginosus was very low compared to that from insects. Insectivorous bats can consume 3.60–8.00 g/day, 4.20–7.60 insects/minute in a foraging time of 20.30–462.70 min (Encarnação and Dietz, 2006). Thus, the proportion of the DNA of this bat, from its intestinal epithelium, in fecal guano could be low in comparison to the insects that had been ingested and excreted and/or insects that directly jumped into the guano and finally died. These results suggest that the main source of the guano is the Miniopterus fuliginosus and that the possibility of interfusion of rodents’ excreta to the guano was very low.

Apart from Miniopterus fuliginosus, a few numbers of Rhinolophus ferrumequinum were also seen in the caves; however, species identification though COI barcoding did not confirm this bat’s genome. This does not eliminate the possibility of this species being a host of EMCV. Since the population of this species was overwhelmed by Miniopterus fuliginosus, it is possible that the feces from Rhinolophus ferrumequinum were not sampled in the representative collection points in the cave. Even in the same cave, the group of Miniopterus fuliginosus and Rhinolophus ferrumequinum were segregated because bats usually inhabit separate rooms of the roosting site (Hengjan et al., 2017; Schmidly and Bradley, 2016). Bats are known to have a high viral sharing network connection, suggesting that viruses may pass easily among bat species (Luis et al., 2015). Thus, it is important to also consider Rhinolophus ferrumequinum as a possible host of EMCV.

One of the many traits of bats that makes them suited to hosting and transmitting viruses is their ability to fly long distances during seasonal migrations in search for food (Calisher et al., 2006). Microchiropterans can travel from 10 km up to 80 km per night (Mildenstein and De Jong, 2011; Wong et al., 2007). This ability is a major factor of viral transmission to various locations. For instance, rabies virus associated with some species of bats have been identified throughout the geographic range and migratory routes of these bats (Calisher et al., 2006). Aside from flight ability, accidental and intentional translocation of bats has also been reported both intra- and intercontinentally. In a review by Constantine (2003), translocation of bats happens accidentally through ships and shipping containers where they land and/or roost, through aircrafts such as the eastern pipistrelle bat found in a plane that travelled from Mexico to Texas, or intentionally when these bats are used for captivity or for release after being rescued. All of these translocation incidents are modes of viral expansion. This is supported in this study as the EMCV genome was also found in samples from Korea and Japan. Moreover, migration of the bats between these countries was previously reported according to genetic analysis (Iida et al., 2017).

This study suggests the presence of EMCV sequences in Miniopterus fuliginosus in Taiwan, Korea, and Japan. It is possible that Miniopterus fuliginosus is one of the natural hosts of EMCV in East Asia.