Ethics statement

The Leishmania isolates used in this study were obtained from a previously published study [21]. Briefly, these isolates had been successfully cultured from biopsies collected between 2011 and 2014 from 20 adult patients (16 males and four females; 19–66 years old) diagnosed with ACL and enrolled in the study before treatment with pentamidine isethionate (Pentacarinat; Rhone Poulenc) (Table 1). All patients, in consultation for a suspicion of leishmaniasis at the dermatology unit of the Cayenne hospital (Centre Hospitalier Andrée Rosemon, CHAR), Cayenne, French Guiana, were informed during consultation by the clinician that case records and biological data might be further used in research and that they had the right to refuse. Biological samples were taken for diagnostic purposes and oral informed consent was documented in the case records by the clinician during consultation. Patient data was anonymized at the Laboratory of Parasitology and Mycology of the Cayenne hospital, which carried out diagnostic analyzes and initial parasite culture. The project did not raise any concerns and was approved by the Ethical Committee at Cayenne Hospital. Ethical approval was granted based on the human experimentation guidelines of the “Comité consultatif sur le traitement de l'information en matière de recherche” (CCTIRS: 2012–42). The monocentric audit of retrospective anonymized case record data was permitted by the “Commission nationale de l’informatique et des libertés” (CNIL: DR.2014-091) as well as the regulations of Cayenne Hospital (http://www.ch-cayenne.net/Droits-et-Devoirs.html).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Epidemiological features of the 20 patients with ACL, Leishmania species identification and characteristics of the 24 LRV1 strains sequenced.

https://doi.org/10.1371/journal.pntd.0005764.t001

LRV1 status and parasite species confirmation

All 20 isolates of Leishmania spp. had formerly been tested for their LRV1 status by RT-PCR and Leishmania spp. were identified using PCR-RFLP, as previously described (Table 1) [21, 36].

Parasite cultures

Parasites were cultured in Schneider’s medium (Sigma-Aldrich) complemented with 10% heat-inactivated fetal calf serum (FCS) (Gibco), 0.6 mg/L Biopterin (Sigma-Aldrich) and 5 mg/L folate hemin (Sigma-Aldrich) without antibiotics. Dry pellets were constituted by centrifugation 5 min at 2500 rpm and removal of supernatant. The pellets were stored at −80°C until RNA extraction.

Biological cloning

For parasite cloning, 0.5 mL of culture medium containing 2–4.10⁷ parasites was spread onto freshly prepared Schneider plates (Schneider’s medium, 10% FCS, Biopterin 0.6 mg/L and folate hemin 5 mg/L, 2% low-melting-point agarose). After 7 days, the colonies were harvested and cultured separately in tubes containing complemented Schneider’s medium for an additional 5–7 days. Separate cultures were then processed for RNA isolation and conventional RT-PCR using primer pairs specific to each virus (see below).

RNA extraction and sequence-independent amplification

Total RNA was extracted from dry pellets of parasites (10⁸ parasites) using the RNeasy mini kit (QIAGEN) as recommended by the manufacturer, treated with Turbo DNase (Life Technologies) for 1 h at 37°C and then measured on a Biophotometer (Eppendorf). Then the generation of double-stranded cDNA and SPIA amplification was performed from 50 ng of total RNA using the Ovation RNA-Seq system V2 (NuGEN Technologies, Inc.) as specified by the manufacturer’s protocol.

Generation of libraries and high-throughput sequencing

The libraries were generated using the NextFlex PCR-free DNA Seq Kit (Bioo Scientific, Austin, TX, USA) with (15 strains, with a two-letter and two-number code) or without (five strains, with a 20XX code) a 10-cycle PCR enrichment before quantification and validation. They were then sequenced on an Illumina MiSeq instrument in 250-base paired-end reads (Illumina, San Diego, CA, USA). Sequence files were generated using Illumina Analysis Pipeline version 1.8 (CASAVA).

Bioinformatics analysis

Raw sequences were filtered and trimmed using fastq-mcf with the following parameters: duplicates were removed when 50 bases were identical between reads (-D parameter) while minimum read length was set to 50 [37]. Adapter sequences were removed using -t (% of occurrence) set to zero. The other parameters were left to their default settings. Then Deconseq was used to eliminate phiX174 phage sequences [38]. Fastq clean files were assembled using SPAdes (version 3.0.0) with the following range of k-mers: 21, 55, 77, 99 and 127 [39]. We used blastn for sequence comparison between assembled contigs and reference LRV1 genomes (LRV1-1 and LRV1-4, accession numbers M92355 and U01899, respectively) [6]. The LRV1-positive contigs were then extracted and used as references for mapping using the bwa mem algorithm (version 0.7.10-r789) with a kmer set to 55 and a minimum alignment score of 40, using a minimum fastq quality of 30 [40]. Some assembly abnormalities were manually corrected. The contigs and known LRV genomes were then multiple-aligned with MAFFT alignment software (version 7.037b) [41].

Conventional RT-PCR and Sanger sequencing

Conventional methods were used to extend sequences when needed. Three sets of degenerate consensus primers were designed according to the NGS sequences obtained. Primers were degenerated so as to amplify products from all samples regardless of the sequence. Full information on the primers is available in Table 2. These sets of primers were used in RT-PCR reactions carried out using the Transcriptor one-step RT-PCR kit (Roche), as recommended by the manufacturer, using 500 ng of input total RNA. RT-PCR products at the expected size were sent to Beckman Coulter Genomics (http://www.beckmangenomics.com/) for direct sequencing. Furthermore, to specifically amplify the different viruses in case of multiple-infected Leishmania isolates, one set of nondegenerate specific primers per virus was designed and used for RT-PCR (Table 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Sequences of oligonucleotide primers used for consensus and specific RT-PCRs.

https://doi.org/10.1371/journal.pntd.0005764.t002

Phylogenetic analyses

Contigs from HTS as well as raw sequences downloaded from the Beckman Coulter Genomics website were analyzed and edited in MEGA 5.05 [42]. Multiple sequence alignments were constructed using ClustalW with all published complete and partial LRV1 sequences. Alignments were checked manually. Sequences were translated into amino acids and both nucleotide and amino acid sequences were checked for irregularities. For phylogenetic trees inferred from the aligned nucleotide sequences, the MrModeltest2.3 program was used to determine the optimal model of nucleotide evolution [43]. The GTR model, with a gamma distribution shape parameter (G) and invariable sites (I), was identified and used for the Bayesian approach, which was performed with Mr. Bayes 3.2.2 to infer phylogenetic relationships [44]. Markov Chain Monte Carlo (MCMC) simulations were run for 10,000,000 generations, with four simultaneous chains, using a sample frequency of 100 and a burn-in of 25,000. Majority rule consensus trees were obtained from the output. Validation of the inference was assessed based on the standard deviation of split frequencies, which was less than the expected threshold value of 0.01. For phylogenetic trees inferred from the aligned amino acid sequences, the ProtTest3 program was used to determine the optimal model of amino acid evolution. The JTT model, with gamma (G) distribution, was identified and used for the Bayesian approach [45], which was performed with Mr. Bayes 3.2.2 [44]. MCMC simulations were run for 10,000,000 generations, with four simultaneous chains, using a sample frequency of 100 and a burn-in of 25,000. Majority rule consensus trees were obtained from the output. Validation of the inference was assessed based on the standard deviation of split frequencies, which was less than the expected threshold value of 0.01.

Naming assignation

LRV1 strains were named according to a recent proposal to the International Committee on Taxonomy of Viruses, i.e., “LRV1” followed by two letters corresponding to the parasite species to which the strain designation was then assigned [46]. For example: LRV1-Lg-LF94 for LRV1 from L. (V.) guyanensis strain LF94 (Table 1).

Accession numbers

Sequences reported in this paper were deposited in the GenBank nucleotide database under accession numbers KY750607 to KY750630. The raw sequencing data are available in the Sequence Read Archive under accession number PRJNA371487 (www.ncbi.nlm.nih.gov/bioproject/371487).

Patients

In the context of epidemiological screening, we tested 129 isolates of Leishmania sp. Collected in French Guiana between 2011 and 2014 for the presence of LRV1 [21]. Here, we analyzed a subset of previously identified LRV1-positive isolates obtained from patients diagnosed with ACL. A total of 20 parasite cultures, obtained from human biopsies isolated from 16 men and four women, aged 19–66 years with a median age of 35 years, were selected for high-throughput sequencing. The data (gender and age) collected for each patient are listed in Table 1. In all, we sequenced 19 LRV1-positive L. (V.) guyanensis isolates and one LRV1-positive L. (V.) braziliensis isolate.

Global analysis of high-throughput sequencing data

Using Illumina MiSeq sequencing technology, 24 full-length or almost full-length genomic sequences were de novo assembled into individual contigs with a mean 49.8 depth of coverage. Sequences were 4.8–5.3 kb long. Alignment of the full-length genome sequences obtained with those available in the database made it possible to design three sets of consensus-degenerate primers (Table 2). These primer pairs were used for conventional RT-PCR and Sanger sequencing to extend the sequences to the 5’ and 3’ ends that were missing for some strains. High-throughput sequencing combined with conventional RT-PCR applied to the 20 LRV1-positive isolates generated almost complete genomic sequences (from nucleotide 60 to 5248 of LRV1-1, accession number NC_002063) from each viral strain. This covers the full-length coding sequence of the virus. In addition, Sanger sequencing of the RT-PCR products also allowed verifying the correctness of the sequences obtained. A total of 24 different viral sequences were obtained from the 20 LRV1-positive parasite isolates tested. All sequences were 5.2 kb long with an average G+C content ranging from 44.7% to 46.6%. Detailed genome statistics are outlined in Table 1. Unexpectedly, co-infections by two or three LRV1 viruses were identified in three parasite isolates, 2028, WF69 and XJ93. LRV1 full-length coding sequences obtained from the same isolate were named 2028_G1, 2028_G2 and 2028_G3, WF69_G1 and WF69_G2, XJ93_G1 and XJ93_G2.

Sequence features

The 24 sequences obtained had the same structural organization as partially overlapping open reading frames surrounded by untranslated regions similar to the sequence already described [10, 12]. In addition, they showed an identical genome size. All but one sequence (strain LF94) showed a double-nucleotide insertion in the 5’UTR region at positions 201–202 relative to LRV1-1 (numbers given here and below refer to the nucleotide positions in the published sequence of LRV1-1, acc. # NC_002063). At this exception of this shared insertion, only a few other nucleotide insertions and deletions were seen in the 5’UTR region of four viral sequences relative to LRV1-1. Thus, strain 2014 exhibited one nucleotide deletion at position 363, strain 2001 possessed two nucleotide insertions at positions 142 and 198 (these two insertions were also present in LRV1-2, LRV1-8 and LRV1-9 5’UTR sequences), strain XJ93_G2 showed a double-nucleotide insertion at position 198–199 and one nucleotide deletion at position 446, while strain YA70 presented one nucleotide insertion at position 198 and one deletion at position 393. In addition to these indels in the 5’UTR, of the 24 sequences, only the coding sequence of the RNA-dependent RNA polymerase gene of the YA70 strain presented one codon insertion at positions 73–75 (D25), one codon deletion at positions 169–171 and a double codon deletion at positions 493–498. These codon indels did not alter the reading frame. YA70 shared the first amino acid insertion at position 25 (aspartic acid or lysine, respectively) and the first codon deletion with strain LRV1-Lb2169 (international code MHOM/BO/2011/2169, acc. # KC862308). These two strains were identified from L. (V.) braziliensis parasites. The YA70 strain also possessed a double stop codon at the end of the coat and RNA-dependent RNA polymerase genes. Three other strains also possessed a double stop codon at the end of the RNA-dependent RNA polymerase gene: LRV1-1, LF94 and 2001. Except for YA70, no deletions or insertions were detected in the coding sequences of the other strains. According to these results, orf2 was 2229 bp in length and orf3 was 2631–2637 bp, encoding a capsid protein of 742 aa and a RNA-dependent RNA polymerase of 876–878 aa, respectively.

Sequence identities

The almost complete genome sequences obtained for each strain as well as the nucleotide and amino acid sequences of the capsid and RNA-dependent RNA polymerase were compared to one another and to those published in the databases. All sequences were unique. On the basis of the common ≈ 5260 bp sequences obtained for each strain, pairwise comparisons showed that nucleotide sequence identities ranged from 76.5% (2001 vs. LRV1-4) to 99.5% (PD46 vs. VL91). CP coding sequences exhibited among themselves from 76.3% (Lg1398 vs. LRV1-Lb2169) to 99.7% (2008 vs. LF98) nucleotide identity and from 89.5% (LRV1-4 vs. LRV1-Lb2169) to 100% (PD46 vs. VL91 and 2008 vs. LF98 and 2028_G1 vs. LL28/LV11) amino acid identity. For RdRp coding sequences, the nucleotide identities ranged from 72.4% (2001 vs. LRV1-Lb2169) to 99.4% (PD46 vs. VL91) and the amino acid identities ranged from 79.9% (2001 vs. LRV1-4 vs. LRV1-Lb2169) to 99.5% (2008 vs. LF98). Based on the genetic distance (Table 3) and phylogenetic relationship results (see below), six clades were clearly identifiable. We tentatively called them clades A–F. With the exception of clade A sequences for which the nucleotide divergence ranged from 0.3 to 10.4%, the intra-clade variability of all clades was below 10% at both the nucleotide and amino acid levels. Moreover, except for the A-B inter-clade variability, which ranged from 95.5 to 98.2% on CP and from 92.3 to 94.9% on RdRp on amino acid sequences, the inter-clade variability was over 10% at the nucleotide and amino acid levels. It is noteworthy that the two LRV1 sequences obtained from isolate WF69, WF69_G1 and WF69_G2, had a 91.8% nucleotide identity. Strains 2028_G2 vs. 2028_G3 gave the same percentage of nucleotide identity (91.8%), while 2028_G1 vs. 2028_G2 and 2028_G1 vs. 2028_G3 showed 86.9 and 87.4% nucleotide identity, respectively. Finally, strains XJ93_G1 vs. XJ93_G2 exhibited 83.4% nucleotide identity. Nucleotide substitutions occurred throughout the sequences.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Nucleotide and amino acid identities between the novel LRV1 sequences and those available in the databases, as reported by clades (A-F).

https://doi.org/10.1371/journal.pntd.0005764.t003

5’UTR sequence comparison

5’UTR sequences (from nucleotide 60 to 450) exhibited 86.2 (YA70 vs. 2001) to 100% (VL91 vs. YZ58 and LF98 vs. 2008) nucleotide identities between them and with the LRV1-1 and LRV1-4 5’UTR sequences. Of interest, it should be noted that YA70, which corresponded to the only sequence obtained from a L. (V.) braziliensis isolate, was the most divergent sequence. If we exclude it from the data set, the nucleotide sequence conservation between LRV1 5’UTR sequences from L. (V.) guyanensis isolates ranged from 88.8 to 100%. The same analysis restricted to the LRV1 strains from L. (V.) guyanensis from French Guiana (excluding strain 2001, for which the geographic location of contamination was reported as the region of Manaus, Amazonas state, Brazil) demonstrated a nucleotide sequence conservation ranging from 91.1% to 100%. In addition, taking into account all the sequences analyzed, 74.7% of the 5’UTR nucleotide positions were identical. By excluding YA70, this percentage increased to 79.1%. Again, by restricting our analysis to LRV1 5’UTR sequences obtained from L. (V.) guyanensis strains from French Guiana, 82.1% of the nucleotide positions were identical. These 100% conserved nucleotide positions were distributed into blocks of different sizes. The longest one, a 50-bp motif, covering nucleotides 313–362, encompassed the previously identified pyrimidine-rich region (333-UUUCUUGUUACUAUU-347), whose first nine nucleotides are complementary to a purine-rich region of the Leishmania 18S rRNA [47], as well as the nucleocapsid endoribonuclease cleavage site 314-GAUCˇCGAA-321 [13, 14]. In addition, analysis of the stem-loop IV sequences located upstream, from nucleotide 266 to 281, showed that all but two sequences were 100% conserved and identical to the published one [15]. The only two divergent sequences corresponded to strains 2001 and XJ93_G2. Strain XJ93_G2 presented two complementary substitutions at positions U269C and A278G of the stem, while strain 2001 presented one substitution at position A278G able to form a wobble base pair with the uracil at position 269 (all positions are given relative to LRV1-1 sequence, accession number NC_002063). This shows that LRV1 5’UTR sequences are strongly conserved in both nucleotide sequences and predicted RNA secondary structures.

Phylogenies

Phylogenetic analyses were done on the almost-complete genome sequences (nucleotides 60–5248 relative to the LRV1-1 sequence) (S1 Fig), as well as on either nucleotide or amino acid sequences of the capsid (Fig 1A) and RNA-dependent RNA polymerase (Fig 1B). They all gave strictly identical tree topologies. As mentioned above, we identified six monophyletic “clades”, A–F, all supported by high posterior probabilities. Most of our strains clustered in clades A and B. Indeed, nine of the strains, with LRV1-4 and LRV1-LgM5313, belonged to the monophyletic clade A while clade B was composed of 11 of our strains. All LRV1 strains belonging to clades A and B were from L. (V.) guyanensis parasites. Clades A and B belonged to a monophyletic lineage supported by a posterior probability value of 1. In addition, in clade A, LRV1-LgM5313 and LRV1-4 sequences formed a monophyletic well-supported branch in clade A (Fig 1A). These two sequences correspond to Brazilian L. (V.) guyanensis isolates, whereas the others were from French Guiana. Clade C contained only a new divergent strain, XJ93_G2, which is well supported phylogenetically in all analyses, while clade D was composed of strains 2001 and LRV1-Lg1398. In addition, a phylogenetic analysis restricted to a common 219-bp fragment of the 5’UTR region (S2 Fig), which allowed including sequences from most of the LRV1-1 to LRV1-13 strains (except LRV1-3, -5, -6 and -12 for which no sequence is available), demonstrated that strain 2001 was close to LRV1-8 and LRV1-9 strains. These four strains were identified from L. (V.) guyanensis isolates from Brazil. Clade E comprised two strains, one of our sequences, LF94, and LRV1-1. Finally, YA70 and LRV1-Lb2169, identified from L. (V.) braziliensis strains, belonged to a separate clade, F, distinct from the other clades composed of LRV1 sequences of L. (V.) guyanensis (Fig 1). Interestingly enough, the two viral sequences identified from WF69 belonged to two distinct groups of sequences within clade A, two of the three strains from 2028 (2028_G2 and 2028_G3) were also part of two distinct groups of sequences within clade A, while the third sequence, 2028_G1, belonged to clade B. Finally, XJ93_G1 belonged to clade B and XJ93_G2 was the only representative of clade C.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Phylogenetic analysis of Leishmania RNA virus 1 isolates.

The phylogenetic trees were inferred from A-) the complete amino acid sequences of the capsid gene (757 amino acids) and B-) the complete amino acid sequences of the RdRp gene (879 amino acids) using the Bayesian method with the JTT + G model of amino acid evolution. New LRV1 sequences generated in this study are in boldface. Sequence identifiers include the strain ID and the NCBI accession number. The major clades representing the different LRV1 genotypes identified (A–F) as well as the clades of LRV1 isolates from L. (V.) guyanensis (Lg) or L. (V.) braziliensis (Lb) parasites are labeled. Posterior probabilities of the Bayesian analysis (>90%) are shown next to each node. The scale bars indicate amino acid substitutions per site.

https://doi.org/10.1371/journal.pntd.0005764.g001

The last phylogenetic analysis, restricted to partial capsid sequences (299bp/99aa), allowed including sequences recently published by Adaui et al. obtained from L. (V.) braziliensis isolates from Peru and Bolivia [29]. All of the LRV1 sequences from L. (V.) braziliensis isolates formed a distinct monophyletic clade, highly supported by a posterior probability of 0.91 (Fig 2). LRV1 sequences from L. (V.) braziliensis isolates were subdivided into distinct clades as observed for those from L. (V.) guyanensis. The two sequences previously identified as clade F (YA70 and LRV1-Lb2169) grouped with the other sequences from L. (V.) braziliensis isolates. Interestingly, YA70, which corresponded to the unique strain obtained from a L. (V.) braziliensis isolate from French Guiana, possessed a basal position to the clade. Nevertheless, phylogenetic relationships of LRV1 isolates from L. (V.) guyanensis strains were less robust given that clades D and E clustered with clade F of L (V.) braziliensis viruses.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Phylogenetic analysis of partial capsid gene sequences (alignment of 299 nucleotides).

The tree was inferred using the Bayesian method with the GTR + G + I model. Novel sequences generated in this study are shown in bold. Sequence identifiers include the strain ID and the NCBI accession number. The major clades representing the different LRV1 genotypes identified (A–F) as well as the clades of LRV1 isolates from L. (V.) guyanensis (Lg) or L. (V.) braziliensis (Lb) parasites are labeled. Posterior probabilities of the Bayesian analysis (>90%) are shown next to each node. The scale bar indicates nucleotide sequence divergence among sequences.

https://doi.org/10.1371/journal.pntd.0005764.g002

Biological cloning

To determine if the multiple LRV1 strains identified in three Leishmania isolates were real viral co-infections or due to the simultaneous presence of distinct parasite populations in our cultures, we proceeded to the biological cloning of the parasites. Two of the three parasites, XJ93 and WF69, were tested. The initial parasite cultures were divided for either RNA isolation or biological cloning. Seven Leishmania clones and three subclones of one of the clones were obtained from XJ93. Five clones were obtained from WF69. RNA extracted from the initial cultures and from the different clones and subclones was then submitted to RT-PCR using different combinations of primers, each one specific to a virus (Table 2). RT-PCR proved the presence of co-occurring viruses in the initial cultures as well as in all the different clones and subclones.