Divergent bufavirus harboured in megabats represents a new lineage of parvoviruses

Bufavirus is a recently recognized member of the genus Protoparvovirus in the subfamily Parvovirinae. It has been reported that human bufavirus was detected predominantly in patients with diarrhoea in several countries. However, little is known about bufavirus or its close relatives in nonhuman mammals. In this study, we performed nested-PCR screening and identified bufavirus from 12 megabats of Pteropus spp. in Indonesia. Furthermore, we determined nearly the full genome sequence of a novel megabat-borne bufavirus, tentatively named megabat bufavirus 1. Phylogenetic analyses showed that megabat bufavirus 1 clustered with known protoparvoviruses, including human bufavirus but represented a distinct lineage of bufavirus. Our analyses also inferred phylogenetic relationships among animal-borne bufaviruses recently reported by other studies. Recombination analyses suggested that the most common recent ancestor of megabat bufavirus 1 might have arisen from multiple genetic recombination events. These results characterized megabat bufavirus 1 as the first protoparvovirus discovered from megabats and indicates the high genetic divergence of bufavirus.

Bats are traditionally divided into two suborders: Megachiroptera (megabats) and Microchiroptera (microbats). Megabats and microbats mainly consist of frugivorous bats and insectivorous bats, respectively 15 . In addition to their biological differences, they are phylogenetically distinguishable 16,17 . Megabats are a potential reservoir of emerging zoonotic viruses, including Nipah virus, Hendra virus, SARS coronavirus and Ebola virus 18,19 . In addition, a metagenomic study showed that megabats harbour various mammalian viruses, some of which are genetically related to human viruses 20 . However, BuVs or related parvoviruses were detected in microbats but not in megabats. In the present study, we investigated BuVs in flying foxes (Pteropus spp.), which belong to the megabat family. Based on the obtained sequence data, we examined the phylogenetic relationships between the identified-BuVs and previously identified BuVs. One of the BuVs identified in this study was further characterized and we obtained a nearly full genome sequence.

Results
A total of 183 megabats consisting of four different species were included in this study (Table 1). Among these, Pteropus sp. and Dobsonia sp. were not identified at species level, but other genetically close species Pteropus hypomelanus and Dobsonia moluccensis, were included 21 . In an earlier study, BuV was detected from the spleen tissues and faeces of shrews 9 . Therefore, we used nested-PCR targeting for the NS1 gene of BuV using 183 spleens and 96 faeces samples (faecal samples were unavailable for a subset of megabats). As shown in Table 1, BuV was identified in three spleen and nine faeces specimens from Pteropus vampyrus, and one spleen specimen from Pteropus sp. Among these, one Pteropus vampyrus in Magelang was positive for BuV in both the spleen and faeces. With a limited number of samples, the detection rate of BuV in spleen (2%) was lower than that in faeces (9%).
Currently, DDBJ/EMBL/GenBank databases contain complete or partial nucleotide sequences of BuV from humans, nonhuman primates, shrews, rats, pigs and microbats. A phylogenetic tree was generated on the basis of the identified sequences of NS1 gene in this study and a corresponding region of the available reference sequences of protoparvoviruses including BuVs (Fig. 1). Known BuVs formed subgroups by host were distinguishable from other protoparvoviruses (groups 1-6, Fig. 1). MgBuV1, MgBuV2 and MgBuV3 diverged from a common ancestor of known protoparvoviruses and were likely to form a distinct subgroup from other BuVs, including microbat-borne BuVs (group 7, Fig. 1).
To reveal the whole genome sequence and genetically characterise further the identified BuVs, we subjected DNA samples isolated from spleens and faeces to next-generation sequencing. A small number of sequence reads from MgBuV1, but not MgBuV2 and MgBuV3, enabled to amplify the overlapping large fragments by PCR and determine nearly the whole genome sequence, including the complete protein-coding region. The obtained genome sequence of MgBuV1 (accession number LC085675) was 4765 nt in size and encoded three complete ORFs: the non-structural protein (NS1) and the viral capsid proteins (VP1 and VP2). These proteins showed 46-61% amino acid sequence identities to those of other BuVs of which complete genome sequences were available ( ( Fig. 2). NS1 contained a putative endonuclease metal coordination motif "HIH" and the helicase motif Walker A "GPASTGKS". The VP1 of most parvovirus species possesses a phospholipase A2 (PLA2) motif that is required for viral entry 22 . Three potential sites responsible for the enzymatic activity of PLA2 were strictly conserved in the VP1 of MgBuV1: a calcium-binding loop site "YLGPG" and two catalytic sites "HDLEY" and "D". The N-terminus of the VP1 and VP2 of MgBuV1 had a glycine-rich motif (G-rich) required for the cellular entry of parvovirus 23 .
Recently, six partial genome sequences of human BuV-like parvovirus were identified from insectivorous microbats in Hungary and China 13,14 . These sequences covered a partial coding region of NS1 and almost complete to complete coding regions of VP1 and VP2 (Fig. 2). Therefore, using the coding region of VP2, we performed a phylogenetic analysis to infer the deeper phylogenetic relationships between these microbat-borne BuVs and megabat-borne MgBuV1 (Fig. 3). All microbat-borne BuVs formed a single cluster (group 5, Fig. 3). MgBuV1 branched from the common ancestor of human BuVs, non human primate-borne BuV (WUHARV parvovirus) and pig-borne BuV, and represented a distinct lineage of BuV (linage 4, Fig. 3). WUHARV parvovirus (linage 3, Fig. 3) fell into the cluster of human BuVs (group 2, Fig. 3), consistent with a previous study of past recombination events behind the origin of WUHARV parvovirus 13 .  We searched for putative recombination events on the genome of MgBuV1 using two different approaches based on sliding window analyses, one using the Bootscan approach and the second using Bayesian inference (Fig. 4A,B). Both analyses suggested a mosaic genome structure in MgBuV1 and possible recombinations among ancestors of MgBuV1, WUHARV parvovirus and porcine parvovirus Zsana, around nucleotide positions 778 to 2384 of the MgBuV1 genome (Fig. 4C). These observations were supported by independent phylogenetic trees built from five regions bordered by the inferred recombination breakpoints (Fig. 5A-E). Among the BuVs and related protoparvoviruses, MgBuV1 branched from the ancestor of porcine parvovirus Zsana in the phylogenetic  tree at position 778 to 1219 (Fig. 5B) and position 1856 to 2384 (Fig. 5D). However, MgBuV1 branched from its ancestor WUHARV parvovirus at position 1219 to 1856 in the phylogenetic tree (Fig. 5C). In contrast, MgBuV1 formed a distinct lineage in phylogenetic trees constructed from position 2384 to 4765 (Fig. 5E). Based on this evidence, we inferred the parsimonious origin of the mosaic genome in MgBuV1 from both recombination analyses (Fig. 4C).

Discussion
Since the discovery of human BuV in 2012, its genetic divergence and worldwide distribution have been updated by a series of reports [1][2][3][4][5][6][7] . Although recent studies described BuVs in nonhuman mammals, little is known about the epidemiology of animal-borne BuVs . [9][10][11][12][13][14] . In this study, we report for the first time, a BuV identified in megabats. MgBuV1 is the fourth animal borne-BuV where the whole genome sequence of the protein-coding region has been determined. In addition, the detection of partial NS1 gene sequences different from MgBuV1 (called MgBuV2 and MgBuV3) indicates the presence of MgBuV1 variants. The discovery of megabat-borne BuVs facilitates a better understanding of the diversity and host-range of BuV.
MgBuV1 has a similar genome organization and is phylogenetically related to other mammalian-borne BuVs. According to the species demarcation criteria of the International Committee on Taxonomy of Viruses, parvoviruses in the same species should share > 85% amino acid sequence identity in their whole NS1 8 . As shown in Table 2, the NS1 of MgBuV1 showed < 60% identity to those of other related parvoviruses. Taken together, these results suggest that MgBuV1 is a novel species of protoparvovirus harboured by megabats.
Recently, complete or partial genome sequences of BuVs related to human BuV were identified from nonhuman primates, shrews, rats, pigs and microbats [9][10][11][12][13][14] . However, the phylogenetic relationships among these viruses remains to be elucidated. Our phylogenetic analyses covered the available sequences of these BuVs and showed a large degree of genetic divergence of BuV in the genus Protoparvovirus. Both phylogenetic analyses based on NS1 and VP2 genes showed that a group of BuVs comprised host-specific lineages and megabat-borne BuVs were clearly distinguished from microbat-borne BuVs. While the phylogeny of BuV is likely to be associated with their host, further studies of BuVs in animals are needed to infer the virus-host coevolution.
Genome recombination is a compelling argument to explain the evolution of a variety of viruses 24,25 ; furthermore, evidence for recombination events have been also observed in the genome of some protoparvoviruses 13,[26][27][28] . In this study, we revealed putative recombined regions and phylogenetic incongruences across the MgBuV1 genome, suggesting genetic recombination between MgBuV1 and porcine parvovirus Zsana (event 1, Fig. 4C) and other recombination events between MgBuV1 and WUHARV parvovirus (event 2, Fig. 4C). These results imply that the most common recent ancestor of MgBuV1 might have arisen from multiple recombination events. Such events also imply the occurrence of dual infection with ancestors of MgBuV1, WUHARV parvovirus and porcine parvovirus Zsana in nature.
In conclusion, we demonstrated the presence of megabat-borne BuV, the first protoparvovirus discovered from megabats. Further genetic characterization revealed the phylogenetic distinctiveness of MgBuV1 and suggested a molecular evolution of MgBuV1 driven by multiple recombination events. These results provide new insights into the high genetic divergence of BuV and will help us to investigate further BuV in animals.  21,[29][30][31] . Species were identified based on the morphological features and nucleotide sequence analysis of both mitochondrial 16S rRNA and cytochrome b as described previously 21 . Sample information is summarized in Table 1.

Methods
Nested-PCR screening. Faeces suspended in RNAlater reagent (Ambion, Life Technologies, Carlsbad, CA) were subjected to DNA extraction using a High Pure Viral Nucleic Acid kit (Roche Diagnostics GmbH, Mannheim, Germany) or QIAamp Viral RNA Mini Kit (Qiagen, Valencia CA). DNA from spleen tissues was prepared using DNAzol reagent (Molecular Research Center, Cincinnati, OH) or QIAamp DNA Mini kit (Qiagen). For nested PCR screen targeting for NS1 gene of BuV, we followed the method of Sasaki et al. 9 .
Complete genome sequencing. To determine the complete genome sequence of megabat BuV, a DNA library was prepared from the faecal-DNA samples using the Nextera XT DNA sample preparation kit (Illmina, San Diego, CA) and were sequenced on an illumina MiSeq instrument using the MiSeq Reagent Kit v3 600-cycle (Illumina) according to the manufacturer's instructions. The sequence reads homologous to known BuVs were identified by read mapping using CLC Genomics Workbench (CLC bio, Aarhus, Denmark). The overlapping large fragments were amplified by conventional PCR with specific primer sets designed from each identified sequence read. The 3′ end of virus genome was amplified by the modified rapid amplification of cDNA ends as reported previously 32 . All the amplified fragments were assembled manually into a contiguous sequence.
Phylogenetic analyses. Sequences were aligned with MAFFT 33 using the FFT-NS-I method. Phylogenetic inferences were performed on multiple sequence alignments with MrBayes 3.2 34 using a General Time Reversible model with a rate variation between sites modelled with a gamma distribution (+ Γ ) with six categories and allowing for a percentage of invariable sites (+ I). Additionally, this study implemented a second analysis approach over the same sliding window that incorporated the power of Bayesian inference to support or reject hypotheses regarding the topology of the phylogenetic tree. For each window, a hypothesis of the clustering of two particular sequences was tested against the null hypothesis that the cluster was less likely than those sequences independently clustering with any other under the same model parameters. To maximize the inferential power of the approach, the marginal likelihoods of the models assuming opposite topological hypotheses were compared with the stepping stone 37 sampling approach implemented in MrBayes 3.2. Then, the difference of the calculated marginal likelihoods for both models, i.e., the model assuming a clustering and the model without the clustering, represented the Bayes factor, which usually is interpreted as supporting the first model if the difference exceeds three and strongly supportive when it exceeds five 38 . All hypotheses were tested under the same model assumed for the phylogenetic analyses.

Recombination analyses.
The dataset in both analyses was composed of sequences in the BuV cluster. Adopting a parsimonious approach, this study accepted only those events suggested by Bootscan that also had support from the Bayesian hypotheses testing.
Nucleotide sequences. All the nucleotide sequences determined in this study were deposited in the DDBJ/ EMBL/GenBank databases under accession numbers LC085665− LC085676.