Phylogenetic study of treatment failure clinical isolates of Leishmania major

Background: Molecular characteristics are necessary for designing programs for the control, prevention, and treatment against cutaneous leishmaniasis. In the present study, treatment failure (TF) clinical isolates of Leishmania major were phylogenetically analyzed using the barcode genes of cytochrome oxidase II (COXII) and 13A/B. Methods: Samples were collected from 126 patients referred to the Diagnosis Laboratory Center from October 2017 to December 2019. All TF cases were assessed using COXII and 13A/B, and phylogenetic analysis was created using BLASTn, T-COFFEE, and MEGA 7.0.21. Results: All 126 isolates were L. major, in which 5 isolates were TF. All isolates had successfully amplied by the specic primer for COXII region. The alignment analysis with all the 5 TF isolates and the standard Friedlin strain showed more than 98% similarity. Phylogenetic analysis with 13A/B showed that all 5 TF isolates are clustered in one group, although phylogenetic analysis showed different clustering after comparison of 5 TF isolates with 16 TR isolates and the same regions in other isolates in GenBank, NCBI. Conclusions: All 21 isolates studied in this study, comprising TF and TR isolates, were clustered in near groups. However, 13A/B could differentiate the 5 TF isolates from selected 6 TR isolates, completely.


Background
Leishmaniasis is a common disease in the tropical and subtropical regions of the world that has a wide range of clinical symptoms caused by a protozoan from Leishmania spp. and transmitted by Phlebotomus spp. The symptoms vary from spontaneous cutaneous lesion to lethal visceral form [1,2]. According to reports from the World Health Organization (WHO), more than 350 million people in over 98 countries are exposed to leishmaniasis. A total of 2 million new infections occur annually, with an estimated 0.7 to 1.5 million cases of cutaneous leishmaniasis [3]. CL is endemic in various parts of Europe, Africa, and Asian countries, particularly in Central Asia and the Middle East. Leishmaniasis is prevalent in many parts of Iran such as Tehran, Mashhad, Tabriz, Kerman, Isfahan, Neyshabur, Bam, and Yazd. CL mainly occurs in urban and rural forms, mainly caused by L. tropica and L. major, respectively.
Initially, Leishmania detection and identi cation were based on the size, the shape of the wound, geographic location, and injection to the laboratory animals [4]. These methods were unable to accurately identify the actual parasite species, so molecular techniques were used as a suitable alternative method.
Molecular techniques also help to nd out the genetic characterization of the parasite, which can be considered to result in precise and appropriate planning and control, as well as e cient vaccine and drug development [5]. Various target genes used in this eld include a variety of mitochondrial and nucleus barcodes [6,7].
Treatment of CL is considered one of the best strategies for control of the disease [8]. Recently, some cases of treatment failure (TF) have been reported from different foci of the world [9][10][11], which reported different genetic analysis. It seems that molecular characterization and phylogenetic studies of TF isolates can help to nd suitable strategies for vaccine development, prevention, and drug design. In this study, TF clinical isolates obtained from patients with CL were characterized molecularly and assessed phylogenetic analysis using one of the important mitochondrial barcodes named Cytochrome oxidase II (COXII) and the other important nucleous barcode gene named 13A/B.

Ethical statement
The informed consent was completed by each patient participating in this study before recording the information and sampling based on the Helsinki declaration. This study was approved by the Ethics Committee of Shahid Sadoughi University of Medical Sciences, Yazd, Iran.

Population study
This study included a population of patients (126 cases) with zoonotic CL caused by L. major that ve of them had no response to Glucantime, which referred to the laboratory of Isfahan Diagnosis Center, Iran, from October 2017 to December 2019. The patients with L. major were followed for three months. The cases with no responses to Glucantime treatment were considered as TF isolates. The cases with the response to Glucantime were considered as treatment response (TR).

Sampling
Samples were collected by scrubbing the lesion edge after disinfecting by using 70% alcohol. Two slides were prepared for either direct microscopic examination or molecular analysis for each patient.

Detection and identi cation
A microscopic examination was carried out to nd the Leishman body. The DNA was extracted from stained smears using the DNA extraction kit (GeneAll, South Korea) based on manufacturer instructions. The quantity of the isolated DNA was evaluated using a spectrophotometer by nanodrop (ABI, USA). With the aim of veri cation of Leishmania spp., ITS1-PCR was used using LITSr-F 5'-CTGGATCATTTTCCgATg-3' and L5.8s-R 5'-TGATACCACTTATCGCACTT-3' [12]. Each ampli cation reaction mixture included 1x PCR buffer, 0.2 mM dNTP, 1.5 mM MgCl 2 , 1.5 U Taq DNA polymerase, 0.5 µM each primer, and 100 ng DNA.
The ampli cation temperature included the rst denaturation 94 °C for 5 min that followed by 35 cycles of 94 °C for 45 s, 50°C for 45 s, and 72 °C for 45 s. The nal elongation was done by 72 °C for 5 min. Positive and negative controls were considered in each run using L. major (MRHO/IR/75/ER) and ddH 2 O, respectively. The ampli cation analysis was done using agarose gel electrophoresis (1%) alongside with 50 bp DNA ladder. The expected amplicon was about 300-350 bp for veri cation of Leishmania spp. detection. Molecular identi cation was performed on positive amplicons digested using Hae III (Bsu RI) for 3 h at 37 °C. The digestion banding pattern analysis was assessed using 3% agarose gel electrophoresis alongside with 50 bp DNA ladder. L major (MRHO/IR/75/ER) was considered as positive control. A banding pattern revealing two main bands of 220 and 140 bp was considered as L. major. All tests were done in triplicate.

Molecular analysis of barcode gene of COXII
The primer pair of COXII was COXII-F 5'-GCATAAATCCATGTAAAACACCACA-3' and COXII-R 5'-TGGCTTTTATWTTATCATTTTGAATG-3' [13]. Each ampli cation reaction included 1x PCR buffer, 0.2 mM dNTP, 1.5 mM MgCl 2 , 1.5 U Taq DNA polymerase, 0.5 µM each primer, and 100 ng DNA. The ampli cation protocol was the rst denaturation at 94 °C in 1 cycle, followed by 35 cycles of 94 °C for 60 s, 57 °C for 20 s, and 72 °C for 30 s. The nal elongation was done by 72 °C for 5 min. Positive and negative controls were considered in each run with L. major (MRHO/IR/75/ER) and ddH 2 O, respectively. The ampli cation analysis was done using agarose gel electrophoresis (1.5%) alongside with 50 bp DNA ladder. The expected amplicon was 602 bp.

Molecular analysis of barcode gene of 13A/B
The primers of 13A/B and the reaction conditions were obtained from the study by Aghai-Maybodi et al. [14]. Positive and negative controls were considered in each run with L. major (MRHO/IR/75/ER) and ddH 2 O, respectively. The ampli cation analysis was done using agarose gel electrophoresis (3%) alongside with 50 bp DNA ladder. The expected amplicon was 120 bp.

Sequencing and molecular analysis
The amplicons of the target barcode gene of COXII in 21 cases, including 5 TF and 16 TR cases were puri ed and sequenced. The amplicons of the target gene of 13A/B regarding 5 TF isolates and 6 TR isolates were puri ed and sequenced. The sequences were repaired and BLAST. All the studied sequences were submitted in the GenBank database using BankIt (NCBI). The related sequences were searched from GenBank (NCBI) using BLASTn, https://blast.ncbi.nlm.nih.gov/Blast.cgi [15]. Multiple alignments were done using T-COFFE. The phylogenetic analysis was done using MEGA 7.0.21 [16]. The evolutionary history was inferred using the Neighbor-Joining method [17]. The evolutionary distances were computed using the Maximum Composite Likelihood method [18] and were in the units of the number of base substitutions per site. All cases were eliminated if they had gaps and missing data. Bootstrap of 1000 was performed to estimate the node reliability for COXII and 500 for 13A/B. Sequences of COXII were compared to seven entries retrieved from GenBank, including KU680818.1, KU680819.1, KF815210.1, KF815211.1, EU140338.1, AF287688.1, and EF633106.1. The sequences of 13A/B regarding to 5 TF and 6 TR isolates were compared each other and phylogenetic tree was drawn.

Molecular analysis of COXII
All 126 cases were successfully ampli ed using the primer pair of COXII (Supplementary le 1). The 21 isolates, including 5 TF and 16 TR, which were selected for sequencing, were sequenced successfully. All of them were submitted at GenBank, NCBI, using BankIt with the accession numbers of MH443402. Phylogenetic tree analysis Sequences of the COXII region from all 21 sequences were analyzed using nBLAST, T-COFFE, multiple alignments, and Mega 7.0.21 software. The phylogenetic trees of each isolate were drawn by the UPGMA method. The multiple alignments of ve TF isolates' sequences were done with the Friedlin strain of L. major (Supplementary le 2). Accordingly, all ve isolates have conserved areas in most sites compared with the European standard L. major (Friedlin), being signi cantly different from the standard European strain in several nucleotides except in terminal regions. Isolate 2 with the dispersed point in two nucleotides differs from the European standard strain and all isolates are dissimilar with the European standard strain in four nucleotides. Using the similarity matrix indicated that isolates 1-5 had the highest genetic similarity of nearly 100% (99.996%) with the European standard strain. The phylogenetic tree drawn based on the COXII locus shows this high genetic similarity by placing all the isolates in a single cluster (Supplementary le 2).  (Figure 1).

Molecular analysis of 13A/B
Because COXII region could not differentiate the TF and TR isolates in this study, we used 13A/B analysis as one of the important nuclear barcode gene. For this reason, all 5 TE isolates were assessed with six TR isolates ( Figure 2). None of the sequences related to 13A/B were deposited in GenBank, NCBI.

Discussion
This study evaluated molecular characterization of 21 isolates of L. major obtained from patients with CL using PCR assays by mitochondrion gene of COXII. All 126 isolates in this study were L. major and ampli ed COXII region successfully using the speci c primer. Sequence analysis of the mentioned region in 21 isolates was clustered in a similar group in comparison with the same region in some other isolates from GenBank, NCBI.
Cytochromes are involved in the electron transport process of the mitochondrial respiratory chain. They are considered one of the most useful genes for molecular and phylogenetic studies [19,20]. The COXII has been sown a higher evolutionary rate in some other studies [6,21]. Because of its multicopy nature, it can be used for clinical material and environmental samples, such as sandflies. Nevertheless, only a few phylogenetic studies of the genus Leishmania have included this marker and its discriminative power has not been established in detail [6,19,21]. Pilot studies have shown it to be useful for the discrimination of subgenera or species complexes, but its applicability for species delimitation has to be further studied. In our study, the region of COXII was used for ampli cation, and the results showed that this region had a suitable sequence for assessing the phylogenetic differences, in addition to being capable of detecting Leishmania with high sensitivity. It exhibited a good replication quality. This region is located on mitochondrial DNA maxicircles with different levels of conservation on DNA, which was also used to determine the heterogeneity of Leishmania and the phylogenetic relationship of Leishmania species [6].
After sequencing the ampli ed region with COXII in our study, it was found that all isolates had the highest similarity with the European standard strain, along with the identi cation of genus Leishmania. Isolates had many conserved sites with the European standard L. major and are ultimately conserved in most sites. Except in the terminal regions, all ve TF isolates are signi cantly different from L. major (Friedlin) European standard strain in several nucleotides. The phylogenetic tree drawn based on the COXII locus shows this high genetic similarity by placing all the isolates in a single cluster. Although all isolates are placed in a single cluster suggesting more genetic similarity, T-COFFEE alignment analysis disclosed that these isolates contained many central conserved areas and little changes occurred in these areas during evolution. This nding is similar to that of Aghai-Maybodi et al. [14], who emphasized slight changes in the central conserved areas of Iranian isolates, as con rmed by its comparison with the standard strain. Although previous studies have further highlighted the ability of this primer pair to distinguish between species and subspecies and are used to con rm species differentiation, this was not achieved in this study. A possible reason could be the small size of the ampli ed sequence because species differentiation cannot be made in this area. Although the primer pairs revealed a high similarity between the ampli ed Leishmania isolates, it is used in studies on the differentiation between strains in different geographical regions [22]. But this barcode gene could not differentiate the TF and TR isolates in this study. However, the 13A/B region could differentiate them successfully (Figure 2). It seems that the barcode gene of 13A/B could be considered as one the important region for differentiation the isolates with different phenotypes. We recommended using this area for differentiation of TF isolates in some other different geography area.

Conclusions
The results showed that the region of COXII was suitable for the identi cation and phylogenetic analysis in Leishmania spp. All the 21 isolates studied in this study, comprising TF and TR isolates, were clustered in near groups

Consent for publication
Consent form was lled out by all patients.

Availability of data and materials
All data are available.

Competing interests
The authors declare that they have no competing interests Funding We declared that the funding body was used for collection, analysis, and interpretation of data, and writing the manuscript.
Authors' contributions SH collected the samples, applied the molecular experiments, and wrote the draft of the manuscript, VR analyzed the patient data, MT interpreted the bioinformatics analysis, GE designed the study and contributed in major manuscript writing, SSH supervised the laboratory tests, MJB helped in manuscript writing. All authors read and approved the nal manuscript. Phylogenetic tree of all 21 isolates mentioned in this study using the mitochondrion barcode gene of COXII. Five isolates with the accession number of MK972457-MK972461 are treatment failure isolates of Leishmania major obtained from patients with cutaneous leishmaniasis. Other 16 treatment response isolates in this study have the accession number of MH443389-MH443404. The other seven isolates were selected from the GenBank, NCBI.