A preliminary study of viral metagenomics of French bat species in contact with humans: identification of new mammalian viruses

Abstract

The prediction of viral zoonosis epidemics has become a major public health issue. A profound understanding of the viral population in key animal species acting as reservoirs represents an important step towards this goal. Bats harbor diverse viruses, some of which are of particular interest because they cause severe human diseases. However, little is known about the diversity of the global population of viruses found in bats (virome). We determined the viral diversity of five different French insectivorous bat species (nine specimens in total) in close contact with humans. Sequence-independent amplification, high-throughput sequencing with Illumina technology and a dedicated bioinformatics analysis pipeline were used on pooled tissues (brain, liver and lungs). Comparisons of the sequences of contigs and unassembled reads provided a global taxonomic distribution of virus-related sequences for each sample, highlighting differences both within and between bat species. Many viral families were present in these viromes, including viruses known to infect bacteria, plants/fungi, insects or vertebrates, the most relevant being those infecting mammals (Retroviridae, Herpesviridae, Bunyaviridae, Poxviridae, Flaviviridae, Reoviridae, Bornaviridae, Picobirnaviridae). In particular, we detected several new mammalian viruses, including rotaviruses, gammaretroviruses, bornaviruses and bunyaviruses with the identification of the first bat nairovirus. These observations demonstrate that bats naturally harbor viruses from many different families, most of which infect mammals. They may therefore constitute a major reservoir of viral diversity that should be analyzed carefully, to determine the role played by bats in the spread of zoonotic viral infections.

Figure 1. Distribution of unassembled read sequences after BLASTn analysis.

(A) Percentage of sequences related to the main categories of existing viruses: vertebrate (blue), plant/fungal (green), invertebrate (brown), protozoan (yellow) viruses and bacteriophages (gray), and unassigned viral sequences (no taxonomic data concerning the family available, indicated in red). The total number of viral read sequences is indicated below each pie chart. (B) The percentage of sequences related to the most abundant viral families, indicated in the same colors for each main viral category as in (A): blue = vertebrate, brown = invertebrate, gray = phage. Viral families accounting for less than 2% of total sequences were pooled and represented as the “other” category (in purple), and read sequences with no available data concerning the taxonomic family were considered to be unassigned (in red).

Figure 2. Distribution of contig sequences after BLASTx analysis.

(A) Percentage of sequences related to the main categories of existing viruses: vertebrate (blue), plant/fungal (green), invertebrate (brown), protozoan (yellow) viruses and bacteriophages (gray), and unassigned viral sequences (no data available concerning the taxonomic family, indicated in red). The total number of viral contigs is indicated below each pie chart. (B) Percentage of sequences related to the most abundant viral families, indicated by the same color range for each of the main viral categories as in (A) : blue = vertebrate, green = plant/fungal, brown = invertebrate, yellow = protozoan, gray = phage. Viral families accounting for less than 2% of all sequences are pooled in the “other” category (in purple), and read sequences with no available data regarding the taxonomic family are considered to be unassigned (in red).

Figure 3. Phylogenetic analysis of the bat nairovirus-related sequences.

(A) Schematic representation of the large (L) segment (almost 12,260 nt encoding the RNA-dependent RNA polymerase of almost 4,040 aa) from the genome of Dugbe virus (GenBank number U15018), with black and gray bars corresponding to the longest contig sequences (>300 nt) of the bat nairovirus (named Ahun nairovirus) identified in specimens b8 and b9, respectively. The genomic region amplified by PCR is represented by a dashed bar, and the sequence used for phylogenetic analysis is indicated with an asterisk. (B) Phylogenetic tree produced from the amino-acid alignment of the partial polymerase fragment (396-nt sequence, translated into a 132-aa sequence, aa positions 2317 to 2448 of the L protein of Dugbe virus, and aligned in accordance with previous studies [86], [87]). The name of the bat nairovirus is indicated in bold. The various serogroups of nairoviruses and prototype species of the other genera belonging to the family Bunyaviridae are indicated on the right of the tree. The scale bar indicates branch length, and bootstrap values ≥70% are shown next to the relevant nodes. The tree is midpoint-rooted for purposes of clarity only.

Figure 4. Phylogenetic analysis of VP1 bat rotavirus-related sequences.

(A) Schematic representation of the VP1 segment (almost 3,300 nt encoding the RNA-dependent RNA polymerase of almost 1,090 aa) of the genome of the lamb rotavirus strain Lamb-NT (GenBank number FJ031024), with black bars corresponding to the longest contig sequences (>300 nt) of the bat rotavirus (named Maule rotavirus) identified in specimen b8 (Myotis mystacinus). The genomic region amplified by PCR is represented by a dashed bar, and the sequence used for phylogenetic analysis is indicated with an asterisk. B) Phylogenetic tree produced from the amino-acid alignment based on the partial VP1 sequence (119 aa, positions 964 to 1082 of the VP1 protein of lamb rotavirus strain Lamb-NT) translated from one of the longest HSPs. The bat rotavirus-related sequence is indicated in bold within the various rotavirus groups. The scale bar indicates branch length, and bootstrap values ≥70% are shown next to the relevant nodes. The tree is midpoint-rooted for purposes of clarity only.

Figure 5. Phylogenetic analysis of the bat gammaretrovirus-related sequence.

(A) Schematic representation of the partial genome structure encompassing the pol (almost 3,580 nt encoding, the polymerase of almost 1,190 aa) gene of the porcine endogenous retrovirus (GenBank number Y17013), with black bars corresponding to the longest contig sequences (>900 nt) of bat gammaretrovirus (named Sers gammaretrovirus) identified in samples from b7 (Eptesicus serotinus). The genomic region amplified by PCR is represented by a dashed bar, and the sequence used for phylogenetic analysis is indicated with an asterisk. (B) Phylogenetic tree produced from the amino-acid alignment based on a selected region (155 aa) of the translated sequence obtained from the PCR product (almost 288 aa, approximate positions 173 to 423 of the pol protein of the porcine endogenous retrovirus). The bat gammaretrovirus-related sequence is indicated in bold, with circles in black indicating bat gammaretroviruses. The scale bar indicates branch length, and bootstrap values ≥70% are shown next to the relevant nodes. The tree is midpoint-rooted for purposes of clarity only. REV, reticuloendotheliosis virus; FeLV, feline leukemia virus; GALV, gibbon ape leukemia virus; F-MuLV, Friend MuLV; R-MuLV, Rauscher murine leukemia virus; M-MuLV, Moloney MuLV; M-CRV, murine type C retrovirus; PreXMRV-1/2, pre-xenotropic MuLV-related virus 1 and 2; PERV-A, porcine endogenous type C retrovirus class A; PERV-B, porcine endogenous retrovirus B; PERV-C, porcine endogenous retrovirus C; RD114, feline RD114 retrovirus; MDEV, Mus dunni endogenous virus; KoRV, koala retrovirus; OOEV, Orcinus orca endogenous retrovirus; BaEV, baboon endogenous virus; RpuRV, Rhinolophus pusillus; RmRV, Rhinolophus megaphyllus; RaRV, Rhinolophus affinis retrovirus; MrRV, Myotis ricketti retrovirus; PaRV, Pteropus alecto retrovirus and MlRV, Megaderma lyra retrovirus.

Figure 6. Phylogenetic analysis of the bat bornavirus-related sequences.

(A) Schematic representation of the L gene (almost 5,140 nt encoding the RNA-dependent RNA polymerase of almost 1,710 aa) of the genome of Borna disease virus isolate bil (GenBank number ACG59353), with black and dark gray bars corresponding to the three bat bornavirus HSPs identified in samples b1 and b6, respectively. The position of the genome region amplified by PCR is represented by a dashed bar, and the sequence used for phylogenetic analysis is indicated with an asterisk. (B) Phylogenetic tree produced from the nucleotide alignment based on a selected region (214–216 nt, approximate position 3233–3446 of the L gene sequence of Borna disease virus isolate bil) of the sequence obtained from the PCR product. Bat bornavirus-related sequences are indicated in bold. The scale bar indicates branch length, and bootstrap values ≥70% are shown next to the relevant nodes. The tree is midpoint-rooted for purposes of clarity only.