The First Non-LRV RNA Virus in Leishmania

In this work, we describe the first Leishmania-infecting leishbunyavirus—the first virus other than Leishmania RNA virus (LRV) found in trypanosomatid parasites. Its host is Leishmania martiniquensis, a human pathogen causing infections with a wide range of manifestations from asymptomatic to severe visceral disease. This virus (LmarLBV1) possesses many characteristic features of leishbunyaviruses, such as tripartite organization of its RNA genome, with ORFs encoding RNA-dependent RNA polymerase, surface glycoprotein, and nucleoprotein on L, M, and S segments, respectively. Our phylogenetic analyses suggest that LmarLBV1 originated from leishbunyaviruses of monoxenous trypanosomatids and, probably, is a result of genomic re-assortment. The LmarLBV1 facilitates parasites’ infectivity in vitro in primary murine macrophages model. The discovery of a virus in L. martiniquensis poses the question of whether it influences pathogenicity of this parasite in vivo, similarly to the LRV in other Leishmania species.


Introduction
Bunyavirales is an order of negative-sense single-stranded RNA (-ssRNA) viruses [1]. They typically have three genomic segments (large, L; medium, M; small, S) encoding a viral RNA-dependent RNA polymerase L (RDRP L), a surface glycoprotein precursor, and a nucleoprotein, respectively [2]. Additional ORFs, usually involved in counteracting the host antiviral response, may be present in S or M segments [3,4]. Each viral segment has terminal complementary sequences governing its interaction with the polymerase. Furthermore, multiple molecules of a nucleoprotein wrap around genomic RNA following helical symmetry [4]. Together, an RNA molecule, a polymerase, and the nucleoproteins form a functional viral ribonucleoprotein (vRNP) capable of transcription and replication [5]. Virions are usually 90-100 nm in diameter and consist of vRNPs of each genomic segment enclosed by a lipid membrane with incorporated viral glycoproteins [3]. Many bunyaviruses (a generic term for Bunyavirales) are causative agents of arthropod-borne diseases of vertebrates and plants [6]. Total genomic DNA was isolated from 10 mL of log-phase trypanosomatid cultures with the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Small subunit rRNA gene was amplified using primers S762 and S763 [40], following the previously described protocol [41]. The obtained PCR fragments were sequenced directly at Macrogen Europe (Amsterdam, The Netherlands) using the primers 883F, 907R, S757, and A757 [42]. The identity of the strains was confirmed by BLAST analysis [43].

DsRNA Isolation and Next-Generation Sequencing
Total RNA was extracted from 10 8 cells using TRIzol (Thermo Fisher Scientific), following the manufacturer's guidelines. The dsRNA fraction was isolated from 200 µg of total RNA using the previously described DNase-S1 nuclease method [38] and visualized in 0.8% agarose gels. The abundance of fragments was analyzed using GeneTools v. 4.3.9 (Syngene, Cambridge, UK). RiboMinus libraries, prepared from the dsRNA sample, were sequenced on the Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) at Macrogen Inc. (Seoul, South Korea).

Viral Sequence Assembly
Transcriptome assembly was carried out essentially as described earlier [32]. In brief, reads were trimmed with Trimmomatic v. 0.36 [44] and assembled de novo using Trinity v. 2.4.0 [45]. Minimal k-mer was set to 5, and other parameters were not changed. Read mapping was performed in Bowtie2 v. 2.3.4.1 [46] and SAMtools v. 1.8 [47], and the coverage was calculated using BEDTools v. 2.25 software [48]. Viral segments were identified by BLAST searches of the 100 most abundant transcripts. Borders of viral segments were determined based on coverage value (with 10 reads per base as the threshold) and presence of specific terminal sequences. To obtain the terminal complementary sequences, original reads were trimmed with BBduk and mapped with BBmap (https://jgi.doe.gov/data-and-tools/bbtools/) to viral contigs assembled previously. GenBank accession numbers for the L, M, and S segment sequences are MK356554, MK356555, and MK356556, respectively.

Prediction of Functional Elements
The search for ORFs in the viral contigs was performed using NCBI ORFfinder [49] with the minimal ORF length set to 150 nt. The identification of the RDRP domain was done using the NCBI Conserved Domain Search [50]. Predictions of the transmembrane domains and membrane-targeting signal peptides were made using the TMHMM v. 2.0 (www.cbs.dtu.dk/services/TMHMM/), TMPred [51], Phobius [52], MEMSTAT3 on PSIPRED server [53], and SignalP v. 4.1 [54] software packages. N-glycosylation sites were identified with NetNGlyc 1.0 Server (www.cbs.dtu.dk/services/NetNGlyc/).

Phylogenetic Analyses
Full-length amino acid sequences of Leishbunyaviridae and Phenuiviridae RDRPs were aligned using MAFFT v. 7.313 E-INS-i algorithm [55]. The alignment was trimmed in TrimAl v. 1.4 with "automated1" algorithm [56], producing a matrix with 1772 amino acid positions that was used for Viruses 2020, 12, 168 4 of 17 phylogenetic reconstructions. Maximum likelihood analysis was performed in IQ-TREE v. 1.6.1 [57]. The best amino acid substitution model, LG with rate heterogeneity across sites approximated using proportion of invariant sites and 4 categories of discrete Γ distribution (+ I + G4), as well as the empirical amino acid frequencies (+ F), was selected by both corrected Akaike information criterion and Bayesian information criterion in the built-in ModelFinder [58]. Statistical supports for the branches were generated by running 1000 thorough bootstrap replicates. Bayesian inference was accomplished in MrBayes v. 3.2.6 [59] with the same substitution model and estimated during the run using "mixed" prior (resulting in 1.0 posterior probability of LG) and other model parameters specified above. The analysis was run for 1,000,000 Monte-Carlo Markov chain generations with default settings. For the comparison of nucleoproteins and glycoproteins, the respective alignments were prepared and trimmed in the same way as described above resulting in 163 and 190 aa data matrices. Maximum likelihood analysis for the nucleoproteins was performed similarly to the RDRPs. For the glycoproteins, pairwise p-distances were estimated in MEGA X [60].

Negative-Stain Transmission Electron Microscopy
In brief, gradient-purified virus samples were applied to a carbon-coated copper grid, stained with molybdenum acetate, and examined under a Philips 201C transmission electron microscope as described previously [38].

Treatment with Ribavirin
Virus-positive L. (M.) martiniquensis culture was treated with 2 mM of ribavirin (Sigma-Aldrich) for 4 weeks. The cultures were passaged weekly and the viral loads were measured by RT-qPCR in the LightCycler480 (Roche Life Science, Penzberg, Germany) as described previously [61,62] using the SYBR Green Master mix (Roche Life Science) and the following primer pairs: LBV_RDRP_for 5'-ggatcagcaaacaggagtcag-3', LBV_RDRP_rev 5'-acatccaaaggctggcataca-3'; and 18S_for 5'-ttatggagctgtgcgacaag-3', 18S_rev 5'-agtacgttctcccccgaact-3'. The cDNA was synthesized with random hexamer primers using the Super Script III-First strand synthesis kit (Thermo Fisher Scientific) following the manufacturer's instructions. Then, 18S rRNA expression was used for normalization. The anti-viral treatment was stopped after 4 weeks, but the viral load was followed for 2 more weeks to ensure stable depletion.

Macrophage Infection
Mouse bone-marrow derived macrophages were infected as described previously [63] with modifications [64]. In brief, differentiated macrophages were cultured in complete RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS), 50 units/mL of penicillin, 50 µg/mL of streptomycin, 2 mM of L-glutamine, and 0.05 mM of 2-mercapto-ethanol (all from Sigma-Aldrich) at 37 • C with 5% CO 2 . These cells were plated into CellStar 24-wells (Greiner Bio-One GmbH, Kremsmünster, Austria) at 4 × 10 5 cells/mL. The stationary-phase Leishmania cells were added at a parasite to macrophage ratio of 6 promastigotes to 1 macrophage. After 2 h, cells were left either in complete RPMI-1640 or in the media combined with 50 U/mL IFN-γ (Bio-Rad) and 0.5 µg/mL LPS (Sigma-Aldrich) (classically stimulated macrophages) or with 25 ng/mL IL-4 (eBioscience/Thermo Fisher Scientific) (alternatively stimulated macrophages). Then, 72 h post infection, macrophages were lysed and amastigotes were counted by a hemocytometer after resuspension in the complete RPMI medium. All experiments were performed in two independent biological replicates and samples were analyzed in triplicate. Statistical analysis was done with a generalized linear model of the negative binominal distribution.
Ethics statement: Animals were maintained and handled in the animal facility of Charles University in Prague in accordance with institutional guidelines and Czech legislation (Act No. 246/1992 and 359/2012 coll. on protection of animals against cruelty in present statutes at large), which complies with all relevant European Union guidelines. All the experiments were approved by the Committee on the Ethics of Laboratory Experiments of the Charles University and were performed under permission No. MSMT-31114/2015-13 of the Czech Ministry of the Environment. All efforts were made to minimize the number and suffering of experimental animals during the study.

Viral dsRNA in Leishmania (M.) martiniquensis
Four isolates of four different species of the leishmanial subgenus Mundinia were screened for the presence of dsRNA viruses. In one of these isolates, L. (M.) martiniquensis MHOM/MQ/92/MAR1, we documented the presence of three major dsRNA bands designated as L, M, and S for large, medium, and short, respectively ( Figure 1A). This sample was sequenced using the Illumina HiSeq platform, yielding 5.4 Gbp of sequence data. The three viral contigs (6.1, 1.2, and 0.7 kb long) were highly abundant (60.1 to 354.2 fold above the average RPKM (Reads Per Kilobase per Million mapped reads) value), which facilitated their quick and reliable identification. Each contained a single ORF; 2012, 334, and 165 aa long in the L, M, and S fragments, respectively. As previously reported for other leishbunyaviruses [32], the proportions of particular viral segments were not even. As compared to the L RNA, the S segment was about six-fold more abundant (Table 1 and Table S1). This is in agreement with the higher demand for the S RNA-encoded nucleoprotein in vRNP formation. The sequences of all three viral segments were complete and included both 5 and 3 terminal "panhandle" inverted repeats (5 -acacaaaga tctttgtgt-3 , Figure 2) necessary for replication, transcription, and translation in bunyaviruses [3]. The sequences of the identified terminal repeats were identical to those of other known LBVs [32,38].

Viral dsRNA in Leishmania (M.) martiniquensis
Four isolates of four different species of the leishmanial subgenus Mundinia were screened for the presence of dsRNA viruses. In one of these isolates, L. (M.) martiniquensis MHOM/MQ/92/MAR1, we documented the presence of three major dsRNA bands designated as L, M, and S for large, medium, and short, respectively ( Figure 1A). This sample was sequenced using the Illumina HiSeq platform, yielding 5.4 Gbp of sequence data. The three viral contigs (6.1, 1.2, and 0.7 kb long) were highly abundant (60.1 to 354.2 fold above the average RPKM (Reads Per Kilobase per Million mapped reads) value), which facilitated their quick and reliable identification. Each contained a single ORF; 2012, 334, and 165 aa long in the L, M, and S fragments, respectively. As previously reported for other leishbunyaviruses [32], the proportions of particular viral segments were not even. As compared to the L RNA, the S segment was about six-fold more abundant (Table 1 and Table S1). This is in agreement with the higher demand for the S RNA-encoded nucleoprotein in vRNP formation. The sequences of all three viral segments were complete and included both 5' and 3' terminal "panhandle" inverted repeats (5'-acacaaaga tctttgtgt-3', Figure 2) necessary for replication, transcription, and translation in bunyaviruses [3]. The sequences of the identified terminal repeats were identical to those of other known LBVs [32,38].  BLASTp searches demonstrate that the ORF sequences within the L and S segments are very similar to the RDRPs (up to 43% identity with 96% coverage) and nucleocapsid proteins (up to 51% identity with 96% coverage) of leishbunyaviruses. The Conserved Domain search identified a bunyaviral RDRP domain (pfam04196) between aa 588 and 1306 in the L segment ORF with an Evalue = 2.24e -22 . The region between the aa 86 and 151 of the same ORF displayed organization typical for the endonuclease domain of leishbunyaviruses ( Figure 2).    BLASTp searches demonstrate that the ORF sequences within the L and S segments are very similar to the RDRPs (up to 43% identity with 96% coverage) and nucleocapsid proteins (up to 51% identity with 96% coverage) of leishbunyaviruses. The Conserved Domain search identified a bunyaviral RDRP domain (pfam04196) between aa 588 and 1306 in the L segment ORF with an E-value = 2.24e −22 . The region between the aa 86 and 151 of the same ORF displayed organization typical for the endonuclease domain of leishbunyaviruses ( Figure 2). Consistent with the previously published data on LBVs [38], the search for the homologs encoded in the viral M segment did not return any hits with BLASTp, Conserved Domain search, PHYRE2, and HHpred software. The analysis of the M segment-encoded glycoprotein with TMHMM, TMPred, and Phobius did not identify any transmembrane domains (TMDs) in the viruses under study. However, like in other LBVs and consistent with the glycoprotein annotation, SignalP detected the N-terminal membrane insertion peptide (with cleavage site between aa 21 and 22) and NetNGlyc predicted two N-glycosylation sites (at amino acid positions 34 and 237) in M segment sequences. Previously, similar results were obtained for LepmorLBV1, whereas in CabsLBV1, CotoLBV1, and the LBVs of Blechomonas spp., two TMDs were predicted in this segment [32,38]. We posit that such discrepancy could be explained by extreme sequence divergence preventing unambiguous identification of these elements. Application of a more sensitive algorithm, MEMSAT3, predicted one TMD in the virus investigated here, LepmorLBV1, as well as in Duke bunyavirus, which was not analyzed before. Similar to typical bunyaviruses, other LBVs have three TMDs in their glycoprotein ORFs. Analyses presented above suggest that the virus under investigation, as other bunyaviruses, can utilize host machinery for glycoprotein synthesis and virion assembly [65,66]. Indeed, the negatively stained transmission electron microscopy on purified virions from L. (M.) martiniquensis demonstrate the typical envelope with evenly spaced surface projections ( Figure 1B).
In summary, we demonstrate that the new virus possesses many characteristic features of leishbunyaviruses and, therefore, we named it Leishmania martiniquensis leishbunyavirus 1 (LmarLBV1).

Phylogeny
The amino acid sequence of the RDRP was used in the phylogenetic inference of L. martiniquensis leishbunyavirus, using sequences of related Phenuiviridae as an outgroup (Figure 3). LmarLBV1 was nested within the clade Leishbunyaviridae with its closest relative being the Duke bunyavirus [67], which presumably infects a trypanosomatid from bees [38]. These two species proved to be sister to a big cluster of viruses from various monoxenous trypanosomatids. Judging by its phylogenetic position, we propose that L. martiniquensis acquired leishbunyavirus from a monoxenous trypanosomatid.
Although the glycoprotein sequences of LBVs are quite divergent, we perceived that the C-terminal part in some of them displayed conserved residues ( Figure S1). Of note, all these viruses were those with one predicted TMD. Moreover, the sequence of this region in LmarLBV1 is more similar to that in LepmorLBV1s, than in DuBV, its closest relative according to the RDRP tree ( Figure 3). Indels, rather than amino acid substitutions, distinguished these glycoprotein sequence fragments. LmarLBV1 had 28 and 77 indels compared to LepmorLBV1s and DuBV, respectively. The analysis of p-distances in trimmed alignments of the full glycoprotein sequences of all available species also showed markedly higher similarity between LmarLBV1 and LepmorLBV1s than between LmarLBV1 and DuBV (Table S3). Of interest, although nucleoprotein sequences were too short for reliable phylogenetic analysis, they grouped LmarLBV1 with DuBV, similarly to the RDRP-based tree ( Figure S2). The amino acid sequence of the RDRP was used in the phylogenetic inference of L. martiniquensis leishbunyavirus, using sequences of related Phenuiviridae as an outgroup (Figure 3). LmarLBV1 was nested within the clade Leishbunyaviridae with its closest relative being the Duke bunyavirus [67], which presumably infects a trypanosomatid from bees [38]. These two species proved to be sister to a big cluster of viruses from various monoxenous trypanosomatids. Judging by its phylogenetic position, we propose that L. martiniquensis acquired leishbunyavirus from a monoxenous trypanosomatid.  Table S2 [32].
Although the glycoprotein sequences of LBVs are quite divergent, we perceived that the Cterminal part in some of them displayed conserved residues ( Figure S1). Of note, all these viruses were those with one predicted TMD. Moreover, the sequence of this region in LmarLBV1 is more

LmarLBV1 Has Minor Effect on Leishmania Infectivity In Vitro
To assess the role of LmarLBV1 in Leishmania biology, we first established an isogenic line of L. (M.) martiniquensis MHOM/MQ/92/MAR1 depleted of leishbunyavirus using ribavirin ( Figure 4A). After four weeks of treatment, the viral load (as judged by RT-qPCR) was significantly diminished in the treated, compared to the untreated cells. Importantly, it stayed low even after the treatment was stopped ( Figure 4A, asterisk), indicating that depletion was not transient.
Wild type and LmarLBV1-depleted L. (M.) martiniquensis were used to infect non-stimulated, classically (LPS/IFN-γ), or alternatively (IL4) stimulated primary murine macrophages to assess early stages of infection. As expected, the infection level in the classically stimulated macrophages was significantly lower compared to either non-stimulated or IL-4-treated cells. Importantly, parasites, which were depleted of virus, were less infective, compared to their wild type kin ( Figure 4B). The effect of viral presence is minor, yet it is statistically significant and may point out the potential role of classically (LPS/IFN-γ), or alternatively (IL4) stimulated primary murine macrophages to assess early stages of infection. As expected, the infection level in the classically stimulated macrophages was significantly lower compared to either non-stimulated or IL-4-treated cells. Importantly, parasites, which were depleted of virus, were less infective, compared to their wild type kin ( Figure 4B). The effect of viral presence is minor, yet it is statistically significant and may point out the potential role of LmarLBV1 in L. (M.) martiniquensis biology.

Discussion
In this study, we describe the first leishbunyavirus of Leishmania. Previously, these viruses were discovered in monoxenous trypanosomatids of the subfamily Leishmaniinae (mainly in Crithidia spp.), genus Blechomonas (subfamily Blechomonadinae), as well as in one plant-infecting dixenous Phytomonas sp. ( Figure 5) [32,38,39]. Leishbunyaviral sequences are also found in metatranscriptomes of insects infected by flagellates of other trypanosomatid genera, such as Strigomonas, Herpetomonas, and Trypanosoma ( Figure 5, light grey) [38]. This is the most widespread and species-rich group of RNA viruses in trypanosomatids known to date. This fact, along with the discordance of viral and trypanosomatid phylogenies documented in the previous studies [32,38], strongly suggests that host-to-host transition is significantly facilitated in this group of viruses. It is explained when taking into account two facts: (i) LBVs are able to form membrane-bound viral particles [68] and (ii) the flagellar pocket of trypanosomatids is an organelle-governing intensive exchange with the milieu by endo-and exocytosis [69][70][71][72]. The documented particles of LBVs measure about 100 nm [38] corresponding to the typical size of clathrin-coated endocytic vesicles in trypanosomatids [73]. Interestingly, clathrin-mediated endocytosis is the general route for an uptake of bunyaviruses [74]. Bunyaviruses evolved to utilize the eukaryotic endomembrane system for virus assembly and spreading. Apparently, LBVs use the same strategy in trypanosomatids. account two facts: (i) LBVs are able to form membrane-bound viral particles [68] and (ii) the flagellar pocket of trypanosomatids is an organelle-governing intensive exchange with the milieu by endoand exocytosis [69][70][71][72]. The documented particles of LBVs measure about 100 nm [38] corresponding to the typical size of clathrin-coated endocytic vesicles in trypanosomatids [73]. Interestingly, clathrin-mediated endocytosis is the general route for an uptake of bunyaviruses [74]. Bunyaviruses evolved to utilize the eukaryotic endomembrane system for virus assembly and spreading. Apparently, LBVs use the same strategy in trypanosomatids.

Figure 5.
A schematic phylogenetic tree of the family Trypanosomatidae (modified from [13]), demonstrating the distribution of leishbunyaviruses (triangles) and leishmaniaviruses (circles) over the genera of these flagellates. Viruses identified in metatranscriptomes of trypanosomatid-infected insects [38] are shown in grey.
Infectivity and formation of viral particles in bunyaviruses depend on glycoproteins, type I transmembrane proteins that are proteolytically processed and glycosylated in the ER [3,8]. Their Cterminal cytoplasmic domains are thought to bind viral ribonucleoproteins and play a crucial role in genome packaging [75][76][77], whereas the N-terminal ectodomains are involved in receptor recognition and membrane fusion [65,78,79]. In leishbunyaviruses, the M segments and putative glycoproteins encoded within them are significantly reduced in size, extremely divergent, and sometimes contain Figure 5. A schematic phylogenetic tree of the family Trypanosomatidae (modified from [13]), demonstrating the distribution of leishbunyaviruses (triangles) and leishmaniaviruses (circles) over the genera of these flagellates. Viruses identified in metatranscriptomes of trypanosomatid-infected insects [38] are shown in grey.
Infectivity and formation of viral particles in bunyaviruses depend on glycoproteins, type I transmembrane proteins that are proteolytically processed and glycosylated in the ER [3,8]. Their C-terminal cytoplasmic domains are thought to bind viral ribonucleoproteins and play a crucial role in genome packaging [75][76][77], whereas the N-terminal ectodomains are involved in receptor recognition and membrane fusion [65,78,79]. In leishbunyaviruses, the M segments and putative glycoproteins encoded within them are significantly reduced in size, extremely divergent, and sometimes contain a reduced number of transmembrane domains. We hypothesize that such a layout reflects a reduced functionality of these proteins and potential broad specificity of viral infection, which explains their facilitated host-to-host transition. It was demonstrated that extended deletions in the bunyaviral glycoprotein N-terminus (ectodomain) do not prevent cell fusion and transport to the Golgi, but lead to attenuation of viruses [80]. This illustrates the propensity of these proteins to undergo reduction. The opportunity to be inherited vertically probably removes the need for efficient proliferation and infection and may even make such properties undesirable.
Transfer of viruses between different species of trypanosomatids is possible because of coinfections, which are quite common in these parasites [81][82][83][84][85]. Coinfections were previously reported for Leptomonas moramango [39]. Here we did not observe viral coinfection but revealed a putative consequence of such an event-re-assortment of genomic segments. This assumption arises from discordance of phylogenies of proteins from the L and S segments on one hand, and the M segment on the other. The RDRP and nucleocapsid of LmarLBV1 are closely related to their counterparts in DukeBV, whereas its glycoprotein is more similar to the corresponding proteins of LepmorLBV1a and LepmorLBV1b.
LmarLBV1 is the first non-LRV virus discovered in Leishmania. It was found in one of the members of the most enigmatic subgenus of these flagellates-Mundinia. Although the first species was characterized over 70 years ago, the subgenus itself was established only recently [21]. For the moment, this taxon contains four described species: L.  [86][87][88][89][90][91]. The first two infect guinea pigs and kangaroos, respectively, while the remaining two are isolated from humans. In contrast to other human-infecting Leishmania, which use sand flies as their vectors, these flagellates may be transmitted by biting midges [92,93]. Host switching may have shaped the genome evolution in these flagellates [94]. While the parasitofauna of biting midges is understudied, several species of monoxenous trypanosomatids are documented in these insects [95][96][97]. This is in agreement with our proposal that LmarLBV1 originates from LBVs of monoxenous trypanosomatids.
Leishmania martiniquensis is frequently found in skin lesions of immunocompromised patients indicating that it may be an opportunistic pathogen [98][99][100][101]. However, recent analysis of multiple records in Thailand and Myanmar reveals that neither the presence nor the severity of the infection is necessarily associated with HIV [87]. Notably, the clinical manifestations range from asymptomatic infection and various types of lesions to visceral disease. Previously, it was demonstrated that LRV1 boosts virulence of Leishmania guyanensis in humans [33,34,102,103]. The discovery of a virus in L. martiniquensis poses an important question on whether it also influences the pathogenicity of this parasite. We demonstrate that the presence of LmarLBV1 is slightly beneficial for Leishmania. The molecular mechanism of such facilitation may be non-specific, since it was recently shown that simultaneous inoculation of virus-negative L. guyanensis and Toscana virus (Bunyavirales, Phenuiviridae) increases footpad swelling and parasite burden in mice, reminiscent of the reaction to the LRV1-positive L. guyanensis [35]. Although it was not shown experimentally, the presence of the membrane-bound viral particles in LBVs suggest that they can be shed by trypanosomatid cells. This way, LmarLBV1 can interact with the immune system of a vertebrate host, increasing the severity of leishmanial infection. Our results signify the need for a systematic exploration of trypanosomatid viromes.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4915/12/2/168/ s1, Figure S1: Alignment of C-terminal part of the putative glycoproteins, possessing only one predicted transmembrane domain. Columns with ≥3 functionally similar amino acids are shaded, those with four functionally similar amino acids are marked with asterisk. Figure S2: Maximum likelihood phylogenetic tree of leishbunyaviral nucleoproteins. Numbers at branches indicate bootstrap support, values below 50 are not shown. The tree is rooted in agreement with RDRP-based reconstruction. Table S1: Summary statistics for RNA-seq data. Table S2: RDRP sequences of Bunyavirales with working abbreviations of viral names used in phylogenetic inferences. Table S3: Pairwise p-distances between glycoprotein sequences in Leishbunyaviridae.