Genetic characterization of a novel picornavirus in Algerian bats: co-evolution analysis of bat-related picornaviruses

Bats are reservoirs of numerous zoonotic viruses. The Picornaviridae family comprises important pathogens which may infect both humans and animals. In this study, a bat-related picornavirus was detected from Algerian Minioptreus schreibersii bats for the first time in the country. Molecular analyses revealed the new virus originates to the Mischivirus genus. In the operational use of the acquired sequence and all available data regarding bat picornaviruses, we performed a co-evolutionary analysis of mischiviruses and their hosts, to authentically reveal evolutionary patterns within this genus. Based on this analysis, we enlarged the dataset, and examined the co-evolutionary history of all bat-related picornaviruses including their hosts, to effectively compile all possible species jumping events during their evolution. Furthermore, we explored the phylogeny association with geographical location, host-genus and host-species in both data sets.

Sequence data selection. All known bat picornavirus sequences, representing either complete or partial coding sequences were retrieved from GenBank, including the novel mischivirus sequence presented in this study. A total of 70 sequences were analysed (3 kobuvirus, 9 mischivirus, 1 crohivirus, 1 kunsagivirus, 1 sapelovirus, 4 hepatovirus, 1 shanbavirus, 50 unassigned viruses), plus 1 amphibian ampivirus which was represented as an outgroup. The sampling location, collection date, and host genus were listed as indicated in GenBank sequence annotation, and/or the literature 19 (Supplementary Table S1). The BtPVs sequences were sampled from 26 bat species belonging to 9 genera (Table 1). Sequences provenance hail nearly entirely from Europe n = 35, Asia n = 25, Africa n = 9 and America n = 1. Sequences with unknown hosts were discarded. In order to represent mammalian host evolution, we downloaded both complete and partial mitochondrial cytochrome b gene (CYTB) 20 . In reference to the M. schreibersii bat, we obtained additional CYTB sequences from the Hungarian Natural History Museum in Budapest.
Sequence data editing and phylogenetic analysis. Both BtPVs and mischiviruses, and their respective hosts' sequences were aligned using the MAFFT alignment tool 21 . Sequence length adjustment was acquired in the use of GeneDoc 22 . The size ranged between 1419-1838 bp regarding the P1 region, 1724-2082 bp representative of the P2 region and 343-1404 bp in reference to the RdRp. Host sequences were not modified. Prior to applying the datasets for phylogenetic reconstruction, we implemented the finest substitution model selection using Mega v6 23 . The GTR + G substitution model was applied for phylogenetic construction based on 3D pol gene from Mischiviruses and their hosts. Additionally, the same substitution model was used to create a tree from P1 region of all BtPVs. Likewise, the GTR + G + I substation model was used to implement phylogeny based on P2 region and 3D pol gene of all BtPVs and their hosts. An amphibian picornavirus Ampivirus A sequence was used to root the viral phylogenies 9 , while Furipterus horrens CYTB sequence was used as a representative of the outgroup in the host tree. Non-clock Bayesian phylogenetic trees were constructed using MrBayes v3.2.4 software 24 . Each analysis operated for 10 million generations (25% were discarded as burn-in) and sampled every 1000 generations, and the resultant trees were then edited using iTOL 25 .
Pairwise genetic distances were calculated between RdRp nucleotide and amino acid sequences, using the MegAlign pro program (DNASTAR v15.2.0) with uncorrected pairwise distances as a metric, and P distance in MEGA v6 23 .
To assess the temporal signal in the viral data above, a regression method of root-to-tip distances against dates of sampling was implemented regarding the RdRp Bayesian trees, of which, TempEst was used 26 .
Phylogeny-trait association analysis. Phylogeny-trait statistics were performed, using the association index (AI), parsimony score (PS), and maximum monophyletic clade (MC) index statistics available in the BaTS package 27 . Mischiviruses and all bat picornaviruses (based on the 3D pol gene) were inspected using BaTS software. Both of these exhibited significant bunching by the following character states of interest: bat host genus, species, or sampling location. The values obtained were interpreted according to Parker and colleagues 27 . This analysis compared the posterior distribution of trees regarding our data formerly mentioned, to a null distribution of 1000 trait-randomized trees. The results were interpreted in accordance with Parker and colleagues 27 . Prior, the trace files generated by MrBayes were analyzed in Tracer v1.6 26 , with the aim of discarding the burn-in trees.
Co-evolution analysis. In order to estimate the virus-host co-divergence scope, we simultaneously analysed picornaviruses (RdRp) and their hosts' phylogenies along with mischiviruses and their hosts' phylogenies, all in the operational use of Jane v4.0 28 . It deduces the nature and the frequency of different evolutionary events, by determining the congruence with the least costly reconstructions of the host-parasite connection, using the tree topologies. Thus, the parameters for the entire event costs (co-speciation, duplication, duplication and host switch, loss and failure to diverge) were set to 0, then 1, and after co-speciation, equal to 0 and other events equal to 1 with a population size equal to 100 and 100 generations for both datasets mentioned above.
Tanglegrams for all bat PVs and their hosts, in addition to mishiviruses and their hosts were created using Dendroscope v3-9-5 29 .
Recombination analysis. To effectively detect recombination, phylogenetic trees generated from various regions of BtPVs genomes (P1, P2 and 3D pol ) were examined regarding tree structure incongruities. Subsequently, all BtPVs aligned nucleotide sequences were imported in the Recombination Detection Program (RDP 4). Recombination events, parental and recombinant sequences as well as putative breakpoints, all underwent analysis using, GENECONV, BOOTSCAN, GENCONV, SISCAN and MAXCHI methods aligned to default settings 30 , and RDP with internal references only as a parameter.
Statistical analysis. The strength of the correlation among picornavirus diversity (the number of detected PV clusters) and the number of PV species for each host genus was estimated using the Spearman coefficient (r). The value interpretation was as follows: 0.00-0.39 "weak" correlation, 0.40-0.59 "moderate", 0.60-0.79 "strong" and 0.80-1.0 "robust".

Results
Detection and genome organization of the novel mischivirus from bats. Out of the 35 sequenced metagenomic libraries regarding viral discovery, picornavirus (PV) was found in one library with 1179 reads. Bat DNA barcoding revealed the virus was detected from a M. schreibersii bat. Based on genome sequence identity level, the novel sequence (MG888045) described in this study is grouped within the Mischivirus genus. Nearly the entire genome (6961 nt) was obtained (some 1400 bp are missing, 5′ UTR and the beginning of the L protein). This demonstrates the typical PV characteristic genome organization of UTR [L-P1(VP0, VP3, VP1)-P2(2 A, 2B, 2Chel)-P3(3 A, 3BVPg, 3Cpro, 3Dpol)] UTR-poly(A); and the conservative motifs were very similar to the Hungarian Mischivirus B described by Kemenesi and colleagues 11 . Genome organization pattern, hypothetical cleavage sites and conserved motifs according to the first start codon in the obtained sequence are indicated in Fig. 1.

Phylogeny and PVs bunching by host and sampling location.
According to the Blast results, the novel Algerian sequence shared 85% of nucleotide identity with the Hungarian virus and 73% identity with the Chinese strain. Moreover, it shared between 91-94% identity with shorter sequences available from the Bulgarian tentative mischiviruses.
The phylogenetic analysis predicated on RNA-dependent RNA polymerase gene (RdRp) of mischivirus sequences (Fig. 2) revealed how the novel Algerian BtPV formed a monophyletic clade together with the Hungarian Mischivirus B sequences 11 ; in addition, the Bulgarian and Romanian tentative Mischivirus B 12 , including the Chinese Mischivirus A 10 . The phylogenetic relationship is supported with high posterior probability values (>90%).
An extended phylogenetic analysis to all bat picornaviruses (Fig. 3, Supplementary Fig. S1) exhibited a grouping, related in some areas of the tree with both host genus and large-scale sampling location (continent), while in other areas it is interspersed. Furthermore, bat mischiviruses clustered according to host genus, and sampling location for each virus species (Fig. 2).
Regression analyses ( Supplementary Fig. S2) exhibited no association between sampling times and root-to-tip genetic distances, neither for bat mischivirus dataset R 2 = 0.0829 ( Supplementary Fig. 2a) nor for all BtPVs R 2 = 0.0218 ( Supplementary Fig. 2b), and due to this, any molecular clock dating was excluded.
Phylogeny-trait association analysis tests (AI = 0.028, PS = 2) offered statistical support (p < 0.05) regarding the clustering of bat mischivirus (n = 16) when considering host genus. The null hypothesis of no association between phylogeny and host species character trait was accepted based on the association index test (AI = 0.38, p = 0.065), while it was rejected based on the parsimony score test (PS = 3, p = 0). The MC statistic supported the association for both Miniopterus (MC = 12, p = 0.001) and Hipposideros (MC = 3, p = 0.001) genera, in addition to the M. schreibersii species (MC = 12, p = 0.001). Regarding the geographical macro area of sampling, the association index (AI = 0.029, p = 0.003) suggested a strong phylogeny-trait association, once the parsimony score (PS = 3, p = 0.23) exhibited no significant relationship. Furthermore, the MC permitted an inspection in each geographic region alone and provided connecting proof for character-trait Europe (MC = 10, p = 0.009). Note how individual traits (single countries, species) consistently provided non-significant results (see Supplementary Table S2).   Table S3). Based on the phylogenetic analyses (Fig. 3) executed using the 69 RdRp gene sequences, 28 different bat genus-specific clusters were identified: three in Kubovirus, three in Mischivirus, one in Crohivirus, one in Shanbavirus, one in Kunsagivirus, one in Sapelovirus, four in Hepatovirus, and fourteen clusters in different unassigned viruses. Likewise, we tallied six PV phylogenetic clusters for Miniopterus genus, five for Myotis, five for Rhinolophus, five for Eidolon, two for Hipposideros, two for Nyctalus, one for Coleura, and one for Ia. Since these BtPVs are related to different host species and sampled from varying locations, the number of RdRp sequences and/or host species within each cluster are different. Strikingly, several of the clusters are composed of a single BtPV RdRp sequence. A very strong correlation was observed among the number of BtPV specific clusters and both its species richness (r = 0.94; p = 0.0001), or geographical sampling area (r = 0.84; p = 0.002). The accurate classification of BtPVs is related to the length of the viral sequences, nonetheless, even short RdRp sequences permitted us to acquire a primary classification of the unassigned sequences.
BtPV sequences are very diverse and pairwise distances calculation resulted in relevant data. The mean nucleotide divergence between sequences from cluster C5 (Myotis) and C4 (Hipposideros) suggest they are closely related, thus C5 (unassigned) may belong to the genus Mischivirus (Table 2 and Fig. 3  www.nature.com/scientificreports www.nature.com/scientificreports/ Coevolutionary analyses. In a dataset comprising ICTV classified mischiviruses and our sequence, topological similarities were observed mischiviruses and their hosts' phylogenies ( Fig. 4) suggesting first, a co-divergence evolutionary scenario, which interestingly, was not supported by reconciliation analysis (using Jane), when considering the least costly event. When taking into account both all putative and classified mischiviruses, incongruence topology was observed, in addition to reconciliation analysis which revealed how host-jumping events were involved in the evolution of mischiviruses (Fig. 5).
When extending the analysis to all BtPVs and their hosts' phylogenies, the history of their evolution was explained by more host-jumping events than co-speciation events, generating incongruent tree topologies ( Supplementary Fig. S3). Jane analysis displayed more duplication and host switch events than co-divergence events, independently of co-divergence costs ( Table 3). The least costly events represented the best evolutionary scenario, Jane results helped us to determine hypothetic donors and receptors in cross-species transmission events. Recombination analysis. The phylogenetic trees built from P1, P2 and 3Dpol regions of BtPVs genomes displayed discordances in their structures indicating potential recombination events (Fig. 6a-c). Out of the 69 recombination events detected with the RDP 4 program, only 11 unique putative recombination events identified with two or more methods were retained (Table 4).

Discussion
Numerous factors (such as geographical, demographical or ecological) likely influence the occurrence of spill-over events and act as driver factors in viral evolution 31,32 . Destruction of the natural habitat of bats worldwide facilitated the urbanization of several dedicated species or simply established more opportunities for human-bat contact events 33 . This constitutes novel possible factors regarding viral emergence as already seen with coronaviruses 34 and rabies 35 . Understanding the evolutionary mechanisms of BtPVs is a prominent direction of research.  www.nature.com/scientificreports www.nature.com/scientificreports/ In addition, the zoonotic potential of all documented BtPVs is clearly unidentified yet and fairly misunderstood 36 . Although animal originated PVs were previously exemplified as potential zoonotic agents, as in the case of the encephalomyocarditis virus, which was clearly revealed through experimental infections on human tissues and primary cell cultures 37 .
Understanding viruses and their hosts' coevolution is crucial to show up the evolutionary character and understand potential disease emergence factors 38,39 . Thus far, several phylogenetic and systematic evolutionary studies were performed with regards to the Picornaviridae family, as the virus members of this family are increasingly discovered 40,41 . However, little is still relatively unknown in reference to the virus-host coevolution patterns for this group. In this study, we describe the first picornavirus from Algerian bats, and we include this novel sequence in detailed coevolutionary analysis, focusing on BtPVs.
Picornaviruses exhibit a high genetic diversity (quasi-species) similar to other RNA viruses 42,43 , and a broad geographical spectrum in Chiroptera order. Moreover, a positive relationship was previously described between the viral richness and geographical distribution of bats 44,45 . Similarly, in the present study, a positive correlation between PV diversity, species richness, and geographical distribution was revealed. Furthermore, more than one virus cluster and larger viral clusters were obtained for genera in which more species have been sampled in a large geographical range.
BtPVs originating from the same host genus were very diverse and not closely related, hence, both BtPVs and bat mischiviruses did not demonstrate any specificity, whether on the species level nor on the genus level. Likewise, for tortoise picornavirus 40 and contrariwise for coronaviruses 46 .
Several BtPVs clusters possess a single viral sequence whereas several sequences were considerably short. For instance, mischiviruses displayed three host genus related clusters. Miniopterus genus virus related cluster comprising the Asian mischivirus sequence (8468 bp), the new Algerian sequence (6,961 pb), the six Hungarian sequences (6,855 bp), and four partial 3Dpol Bulgarian sequences (three sequences 343 bp, one 983 bp). Myotis genus cluster composed of three Romanian partial 3D pol sequences (993 bp), while the Hipposideros genus cluster consisting of a single Mischivirus C sequence from the Congo (8096 bp). Both the number and the length of sequences often limit the analyses. By way of illustration through the use of BaTS software if and when the www.nature.com/scientificreports www.nature.com/scientificreports/ number of sequences regarding a character trait is lower than three, the result is declared insignificant (case of America). Therefore, short sequences may likely be misclassified.
According to Lewis-Rogers and Crandall 43 , PVs and their hosts evolve through host-jumping events and not via co-speciation. Chiropteran order was not included in the previous study. The main finding of the present study was the frequency regarding host-jumping events occurring in the evolutionary history of BtPVs, and reflected as incongruence between the virus and the host phylogeny. Furthermore, we could observe nearly all bat genera host phylogenetically divergent viruses and an absence of species specificity. According to studies performed on coronaviruses 14 it may likely be due to multiple introductions of PVs, this supposition is coherent with the detection of highly related PVs in humans and different animal species, such as the aichivirus 1 in humans 47  Cross-species transmission event occurrence is multifactorial, nevertheless, we concluded in some cases of BtPVs, sympatry may increase host-jumping events. Co-roosting of Myotis and Miniopterus bats may likely explain the detection of closely related kobuvirus associated among these bats 14 . The migratory ability of Miniopterus bats with the longest distance of 883 km recorded in Europe, which possibly eases the spread of viruses among bat populations spread out in geographically distant areas. For example, M. schreibersii, originating from Europe and Algeria, which possess different geographical distribution and share the same virus species Mischivirus B 54 .
Despite the fact in which no interaction is known regarding M. oxygnathus and H. gigas in accordance with their ecologies, hypothetical cross-species transmission events were detected among these two species. The length of the viral sequence obtained from M. oxygnathus (993 bp) is a limitation regarding accurate classification.  www.nature.com/scientificreports www.nature.com/scientificreports/ Based on the phylogenetic analyses (Fig. 2), the putative mischivirus sequence obtained from M. oxygnathus and M. myotis were closely related to Mischivirus C, identified from H. gigas. Although, they revealed an identical branching pattern upon the structure of the phylogenetic tree as Mischivirus A from M. schreibersii sampled in China and Mischvirus B from Algeria and both Bulgaria and Hungary, indicating BtPVs from Romania may not belong to Mischivirus C. Furthermore, it has been demonstrated in which M. schreibersii bats from Hainan were actually M. fuliginosus 10,55 . This is due to the cryptic nature of several Miniopterus bat species. Distinctly, the identification based on their morphology is not sufficient, and typically requires the use of DNA barcoding or echolocation studies 56 .
Our results emphasized the evolution among PVs within bats horizontally, through host jumping mechanism, rather than co-speciation. Additionally, host specificity was not observed, suggesting it is not involved in the evolutionary history regarding BtPVs. Moreover, we found the occurrence of cross-species transmission events Bayesian phylogenetic reconstruction of RdRp region for all BtPVs genomes included in this study. Ampivirus A sequence (NC027214) used as an outgroup for the three viral phylogenetic trees. Genus-specific clusters are colored, based on bat genus. Recombination may be reflected in tree structure incongruities between phylogenetic trees (a-c).