Molecular Characterization of Leishmania RNA virus 2 in Leishmania major from Uzbekistan.

Here we report sequence and phylogenetic analysis of two new isolates of Leishmania RNA virus 2 (LRV2) found in Leishmania major isolated from human patients with cutaneous leishmaniasis in south Uzbekistan. These new virus-infected flagellates were isolated in the same region of Uzbekistan and the viral sequences differed by only nineteen SNPs, all except one being silent mutations. Therefore, we concluded that they belong to a single LRV2 species. New viruses are closely related to the LRV2-Lmj-ASKH documented in Turkmenistan in 1995, which is congruent with their shared host (L. major) and common geographical origin.

1. Introduction

Leishmaniavirus is a genus of the family Totiviridae. In addition to this, the family includes similar viruses found in other protists (Giardia, Trichomonas, Eimeria) and various fungi [1]. Their virions are icosahedral (T = pseudo 2), non-enveloped, and approximately 40 nm in diameter [2,3,4]. The double-stranded RNA (dsRNA) genome of Leishmaniavirus is not segmented, is approximately 5.2 kb in length, and contains two open reading frames (ORFs) for the capsid protein and the RNA-dependent RNA polymerase (RDRP) [5,6]. The capsid ORF has a dedicated start-codon, whereas RDRP is translated as a C-terminal extension of the capsid. In different leishmania viruses, the mechanism of fusion Gag-Pol protein production varies, with RDRP located in +1 or -1 frameshift relative to or in-frame with the capsid [7,8].

Based on their phylogeny, leishmania viruses were divided into two major groups: Leishmania RNA virus 1 (LRV1), infecting New World Leishmania (representatives of the subgenus Viannia: L. guyanensis, L. braziliensis) and LRV2, documented in the Old World leishmanias (representatives of the subgenus Leishmania: L. major, L. aethiopica, L. infantum). The first LRV2 was isolated from L. major in Turkmenistan [9]. The phylogenies of LRV1/2 and their respective hosts are congruent, suggesting the coevolution between the virus and its Leishmania host [8,10]. It was assumed that the entire lifecycle of totiviruses takes place in the cell cytoplasm and they never produce extracellular virions [1,11]. This view was supported by the coevolution pattern observed in LRVs. However, our recent study provided phylogenetic evidence for the LRV transfer between Leishmania and Blechomonas, two separate genera of the family Trypanosomatidae [12]. In addition, it has been recently demonstrated that LRV1 can exit the cell via extracellular vesicles, explaining the assumed cell-to-cell virus transfer [13].

LRV1 contributes to the pathogenicity of the New World L. guyanensis by interfering with the immune response of a vertebrate host [14]. Viral dsRNA interacts with the Toll-like receptor 3 (TLR3) inside parasitophorous vacuole, leading to the overexpression of pro-inflammatory cytokines (TNF-α, IL-6, INF-γ) and, subsequently, chronic inflammation in the primary lesion. This, in turn, facilitates the parasites’ dissemination to the secondary site (nasal mucosa), manifesting itself as a muco-cutaneous leishmaniasis (MCL) [15,16]. These crucial findings spurred the research on prevalence and diversity of Leishmania viruses [17], as well as viruses infecting other trypanosomatids [18,19]. The primary focus of this research was LRV1 in Leishmania (Viannia) spp. [20,21]. It was reported that LRV1 is more common in Leishmania isolates originating from Brazilian Amazonia, whereas isolates from the southern states of Brazil are mostly virus-free [21]. This implies that other factors, such as co-infection with other RNA viruses [22] or glycoconjugates [23], may also contribute to the Viannia pathogenicity.

The vast majority of the Old World leishmanias causes localized, slowly healing sores (cutaneous leishmaniasis (CL)). The notable exception is Leishmania donovani, which causes disseminated visceral leishmaniasis affecting many internal organs, such as the spleen, liver, and bone marrow [24]. Despite high pathogenicity, no LRVs were detected in 22 clinical isolates of L. donovani. However, 15 of these isolates contained another virus-infected trypanosomatid, Leptomonas seymouri, alongside L. donovani [25]. Leptomonas seymouri bears a capsid-less Narna-like virus and lacks RNA interference pathways, leading to accumulation of high amounts of dsRNA in the cytoplasm [26,27]. Thus, the influence of viral RNA on the development of symptoms, caused by L. donovani, is plausible. In Leishmania aethiopica, LRV2 was associated with elevated levels of TNF-α and IL-6 in vitro, suggesting that LRV1 and 2 may have a similar impact on the vertebrate immune system [8]. Besides originally described LRV2 from L. major in Turkmenistan and L. aethiopica in Ethiopia, a handful of virus positive isolates have been recently documented in Iran (L. infantum and L. major) [28] and Turkey (L. major) [29]. For these samples, only partial sequences were reported, precluding their robust phylogenetic analysis. In general, LRV2s appear to be less abundant compared to LRV1s.

In Uzbekistan, over 6000 cases of CL, caused by Leishmania major, were reported throughout the mid-20th century. Some of these cases were associated with a highly aggressive strain of L. major [30]. In the years following the collapse of the Soviet Union, the absence of appropriate preventive and monitoring activities led to the reemergence of leishmaniasis in Uzbekistan and other countries in Central Asia, making it a major public health problem in this region [31,32]. CL in Central Asia is mostly zoonotic, with gerbils serving as a natural reservoir of parasites [33,34]. In this study, we conducted a survey and sequence analysis of LRVs found in different Leishmania spp. isolated from both wild rodents and human patients.

2. Materials and Methods

2.1. Parasite Culture and RNA Isolation

Cryopreserved stocks of 10 Leishmania spp. isolates (3 L. major, 4 L. turanica, 2 L. gerbilli, and 1 L. infantum), collected from wild gerbils and lesions of infected patients in Central Asia, were initially grown on a biphasic blood agar overlaid with the M199 medium (Sigma–Aldrich, St. Louis, MO, USA) for 1 week. Leishmania promastigotes were then transferred to the M199 medium, supplemented with 10 mg/mL of hemin (Jena Bioschience Gmbh, Jena, Germany), 10% fetal bovine serum, 500 units/mL of penicillin, and 50 µg/mL of streptomycin (all from Thermo Fisher Scientific, Waltham, MA, USA) at 25 °C.

2.2. dsRNA Isolation and Next-Generation Sequencing

Total RNA extraction from 108 promastigotes was performed using TRIzol isolation reagent (Thermo Fisher Scientific) as described previously [19]. Fifty mg of total RNA from each sample were treated with RNase-free DNase I (New England Biolabs, Ipswich, USA) and S1 nuclease from Aspergillus oryzae (Sigma–Aldrich) according to the manufacturer’s instructions. Viral dsRNA bands were visualized on 0.8% agarose gel and stained with ethidium bromide. Individual dsRNA bands were gel purified using a Zymoclean Gel RNA Recovery Kit (Zymo Research, Irvine, CA, USA). RiboMinus libraries were generated and sequenced using the Illumina NovaSeq platform (Illumina, San Diego, CA, USA) at Macrogen Inc. (Seoul, South Korea).

2.3. Bioinformatics Analysis

The raw sequence reads were trimmed with Trimmomatic v. 0.36 [35]. Read mapping was done in Bowtie2 v. 2.3.4.1 [36] and SAMtools v. 1.8 [37], and assembled de novo with Trinity [38]. Coverage was calculated using BEDTools v. 2.25 software [39]. ORFs of the assembled viral contigs were analyzed with NCBI ORFfinder [40]. The mutational analysis was accomplished in the HIVE-hexagon Population Analysis Tool program [41]. For the phylogenetic analyses, amino acid sequences of the capsid and RDRP were separately aligned by MAFFT v. 7.313 E-INS-i algorithm [42], trimmed with TrimAl v. 1.7 using “automated1” settings [43] and concatenated in FaBox v. 1.5 [44]. For both partitions, LG + F + I + G4 was selected as a best fit model by ModelFinder [45]. Maximum likelihood analysis with standard bootstrap method (1000 replicas) was performed in IQ-TREE v. 1.6.1 [46].

3. Results

Ten Leishmania strains (L. major, L. turanica, L. gerbilli) isolated mostly from the CL-infected gerbils and human patients in different endemic areas of Turkmenistan, Kazakhstan, and Uzbekistan (Table 1) were screened for the presence dsRNA (indicating viral infection [19]). Two L. major isolates, both isolated form the soars of humans in Muborak, Qashqadaryo region in southern Uzbekistan, were found to be positive. Other analyzed samples from other districts of Uzbekistan (Karakalpak, Jambyl region, and Termez, Surxondaryo region), Turkmenistan (Serdar, Balkan region, and Tejen, Ahal region,) and Kazakhstan (Embi, Aktobe region) appeared negative on dsRNA gel electrophoresis. Both positive samples were characterized by a distinct 5.2 kbp band, which corresponded in size to the genomic dsRNA of the Leishmania RNA virus (Figure 1).

Table 1

Leishmania spp. analyzed in this work.

Leishmania spp.	WHO Code	LRV2	Source 1	Origin
L. turanica	MRHO/KZ/87/MARZBK7	-	R. opimus	Embi, Aktobe region, Kazakhstan
L. major	MHOM/UZ/1998/Isv M-01h	+	H. sapiens	Muborak, Qashqadaryo region, Uzbekistan
L. turanica	MRHO/UZ/87MarzKK-52R	-	R. opimus	Karakalpak, Jambyl region, Uzbekistan
L. turanica	MRHO/SU/Marz KD-51	-	R. opimus	Karakalpak, Jambyl region, Uzbekistan
L. gerbilli	MRHO/TM/95/9520	-	R. opimus	Serdar, Balkan region, Turkmenistan
L. gerbilli	MRHO/TM/95/9518	-	R. opimus	Serdar, Balkan region, Turkmenistan
L. major	MHOM/UZ/1998/Isv M-27h	+	H. sapiens	Muborak, Qashqadaryo region, Uzbekistan
L. infantum	MHOM/KZ/75/MarzDzha	-	H. sapiens	Karakalpak, Jambyl region, Uzbekistan
L. turanica	I/TM/95/Ph-82	-	P. papatasi	Tejen, Ahal region, Turkmenistan
L. major	MRHO/UZ/2003/Isv T-38g	-	R. opimus	Termez, Surxondaryo region, Uzbekistan

1 Species are abbreviated as follows: R. opimus is Rhombomys opimus (great gerbil), H. sapiens is Homo sapiens, P. papatasi is Phlebotomus papatasi (sandfly), LRV2 is Leishmania RNA virus 2. “+” and “-” denote presence and absence, respectively.

Figure 1

Detection of viral double-stranded RNA (dsRNA) on agarose gel. M; Gene Ruler 1 kb Plus ladder. Sizes on the left are in kb.

RNA-sequencing was performed on positive samples and sequences were deposited to GenBank as LRV2-Lmj-Uzb1 and 2 (GenBank Acc. No. MN418974 and MN418975, respectively). Both sequences were 97% similar to the previously described LRV2-Lmj-ASKH isolated from L. major in Turkmenistan and had the same in-frame arrangement of capsid and RDRP ORFs with a single stop codon between the frames. Nucleotide sequences around the stop codon, where RNA pseudoknot and ribosome shunting sites were previously predicted [9], were invariably conserved in all three viral isolates. In-depth mutational analysis showed 170 nt/40 aa and 174 nt/41 aa differences between Uzb1/ASKH and Uzb2/ASKH sequences, respectively. In turn, Uzb1 and Uzb2 isolates varied by 19 SNPs, with only one resulting in amino acid substitution (T197I) in the RDRP (Table 2). Out of 19, there were three mutations with less than 100% frequencies in RNA-seq data, indicating that several viral subpopulations may be present in one isolate. Clonal analysis revealed heterogeneity in the Uzb1 isolate (Figure 2A), whereas the Uzb2 isolate was uniform (Figure 2B). The Uzb1 isolate showed two minor alternatively assembled regions, around 450 and 320 bp long (Figure 2A), which differed by 2 and 1 SNPs, respectively. Quantitatively, these alternative sequences accounted for 4.3% and 2.7% of the main contig abundance, respectively.

Table 2

SNP comparison between LRV2-Lmj-Uzb1 and LRV2-Lmj-Uzb2. Coverage indicates a total number of sequence reads covering each position on the genome.

Position, nt	LRV2-Lmj-Uzb1	LRV2-Lmj-Uzb1	Frequency, %	Coverage
389	C	T	100	384
539	T	A	40	372
1181	C	A	100	1525
1415	T	C	52	430
1517	C	T	100	732
1772	A	G	100	3616
1835	A	G	100	4188
1892	G	A	100	3890
1922	T	G	100	1434
2603	T	A	100	836
2648	G	A	100	593
3064	T	C	100	345
3125	G	A	100	493
3434	A	G	31	224
3464	A	G	100	289
3806	A	T	100	698
3809	G	A	100	686
3881	G	A	100	704
4892	C	T	100	1154

Figure 2

Clonal analysis revealed heterogeneous viral populations in LRV2-Lmj-Uzb1. Nucleotide positions are on the X axis, the Y axis serves as a ruler in the Sankey graph and represents depth of coverage. (A). LRV2-Lmj-Uzb1 has additionally two minor clones (green colour inserts). (B). LRV2-Lmj-Uzb2 is homogeneous.

Maximum likelihood phylogenetic analysis positioned newly identified isolates sister to the LRV2-Lmj-ASKH with maximal statistical support, which correlates with their common geographical origin in Central Asia (Figure 3). Currently, three host-specific LRV lineages can be defined in Leishmania: LRV1 (New World clade), infecting L. gyanensis and L. braziliensis, LRV2 (Ethiopian clade) found in L. aethiopica, and LRV2 (Central Asian clade) from L. major. When comparing nucleotide sequences of Uzb1 and Uzb2 to the reference strain LRV2-Lmj-ASKH, we noticed that both isolates had atypical frame-shift regions in the capsid and RDRP ORFs. There were four such regions with the length of 51 bp (capsid), 23 bp, 28 bp, and 21 bp (RDRP) long (Figure S1). Importantly, closely related LRV sequences, belonging to the same clade, lack indels and differ only by SNPs. Large indels are starting to be apparent only when LRVs from different clades are compared, at which point sequence similarity drops below 68%. Visual examination of the frame-shift regions revealed the presence of the singleton indels in the LRV2-Lmj-ASKH sequence, which can be explained by the sequencing errors (it was analyzed in 1995 [9]).

Figure 3

Maximum likelihood tree of LRVs built on amino acid sequence alignment using LG + F + I + G4 model. Newly identified isolates are highlighted in black. Standard bootstrap support (1000 replicas) values above 85 are shown. Victorivirus and Eimeriavirus (Totiviridae) were used as an outgroup. Scale bar represents the number of substitutions per site.

4. Discussion

In this paper, we present sequence and phylogenetic analysis of two newly identified leishmania viruses (LRV2) of Leishmania major isolated from human CL patients in southern Uzbekistan. This is the second report of LRV2 in L. major in Central Asia, following the original description of LRV2-Lmj-ASKH isolate from Turkmenistan in 1995 [9]. Two viruses were documented in the same host species (L. major) in the same geographical region, arguing that they may represent just one viral species.

Phylogenetic and mutational analyses suggest that LRV2-Lmj-Uzb1 and LRV2-Lmj-Uzb2 were most closely related to each other and the previously described LRV2-Lmj-ASKH. The phylogenetic position of these viruses correlates with their common geographical origin in Central Asia. Also of note, the leishmanian LRV1 and LRV2, infecting the same or closely related species of trypanosomatids, are monophyletic. Together, these observations support the currently prevailing view of predominantly vertical transmission and tight coevolution of the virus with its Leishmania spp. host [8,10]. Assembly of the RNA-seq data of the Uzb1 isolate yielded alternative contigs with the abundance 20–30 times lower than that of the major contig. This indicates the presence of the viral subpopulation in the Uzb1 sample. Given the low proportion and small number of mutations, it is parsimonious to suggest that the observed minor sequences are the result of clonal evolution of the virus. To date, it appears that the “tight coevolution” rule is applicable only to LRVs infecting Leishmania. The recent finding of LRV3 and 4 in Blechomonas spp. [12] suggests that there were at least two events of horizontal viral transfers between Blechomonas and Leishmania spp. In general, horizontal viral transfers are quite common within monoxenous trypanosomatids, with arthropods serving as mediators of viral exchange [19]. It is possible that the transition of Leishmania to a dixenous lifecycle, which necessitated adaptation to a specific arthropod vector [47,48], isolated these flagellates from their monoxenous kin. This, in turn, has cut down the viral flow between Leishmania and other trypanosomatids, resulting in the strict coevolutionary patterns observed nowadays.

Our analysis revealed only two positive isolates out of 10 samples analyzed, which is consistent with the overall low prevalence of Leishmaniavirus in the Old World leishmanias [16]. In other studies, only 2 LRV-positive isolates (one L. infantum and one L. major) were documented out of 50 tested in Iran [28] and no LRV-positive L. donovani isolates were found among 22 tested in India [25]. Old World Leishmania spp. do not cause MCL, in which LRV1 was shown to be involved [14,16]. On the other hand, LRV2 from L. aethiopica was found to display similar immunological effects as LRV1 in vitro [8]. Thus, LRV2 might influence development of the visceral disease, which is a predominant form of disseminated leishmaniasis in the Old World. Currently, there is only one report of LRV2 in L. infantum isolated from a visceral leishmaniasis patient, but the potential role of this virus in etiology and progression of the disease was not investigated [28]. Studies of this kind are further complicated by the absence of isogenic virus-free isolates for Old World leishmanias [49]. Furthermore, it might be difficult to obtain such isolates as these parasites do not have a functional RNA-interference pathway [50,51] that was harnessed for curing LRV1 in L. guyaynensis [52]. Consequently, broader sampling and more in-depth molecular studies are needed to elucidate the diversity of LRV viruses, as well as their interaction with Leishmania and vertebrate hosts.