Whole-Genome Sequencing and Comparative Genomic Analysis of Leishmania (Viannia) naiffi and L. (Viannia) shawi Reveal Species-Specific Genes and Novel Potential Drug Targets

TL;DR

This study sequenced and compared the genomes of two Leishmania species, identifying unique genes and potential drug targets for treating leishmaniasis.

Contribution

The study provides novel insights into species-specific genes and potential drug targets in Leishmania (Viannia) naiffi and L. (Viannia) shawi.

Findings

L. naiffi and L. shawi have 46 and 25 unique genes, respectively, identified through pangenome analysis.

Functional screening identified calpains, ABC transporters, and GSK-3 as potential drug targets.

L. naiffi showed greater genomic diversity with 34,480 SNPs compared to 26,562 in L. shawi.

Abstract

This study presents the complete sequencing and comparative genomic analysis of Leimania (Viannia) naiffi and Leishmania (Viannia) shawi, species of epidemiological relevance in the Brazilian Amazon. Genome assemblies yielded sizes of 32.13 Mb and 32.51 Mb, with 8170 and 7767 annotated genes, respectively. Predicted gene functions were primarily related to catalytic, binding, and ATP-dependent activities. Pangenome analysis revealed a core genome of 6256 genes alongside notable species-specific differences, including 46 and 25 unique genes in L. naiffi and L. shawi. Functional screening identified pharmacologically promising proteins such as calpains, ABC transporters, and notably, GSK-3. Ploidy analysis indicated tetraploidy on chromosome 8 in L. naiffi and chromosome 2 in L. shawi. Genetic variability assessment detected 34,480 SNPs in L. naiffi and 26,562 in L. shawi, indicating greater genomic diversity in the former. Phylogenetic inference based on the polA1 gene confirmed the placement of both species within the Leishmania (Viannia) subgenus. These findings advance Leishmania genomics knowledge by highlighting unique genetic signatures, regions of high variability, and potential therapeutic targets. This work establishes a foundation for future research on evolution, pathogenicity, and drug development for leishmaniasis.

Figures

Full text

1. Introduction

Leishmaniasis is a neglected disease caused by flagellated protozoa of the genus Leishmania, responsible for infecting thousands of people worldwide each year [1,2]. Species of this genus belong to the family Trypanosomatidae and are traditionally classified into the subgenera Leishmania (Leishmania) and Leishmania (Viannia), which exhibit relevant epidemiological, clinical, and genetic differences [3].

Transmission of leishmaniasis occurs through the bite of sand flies belonging primarily to the subtribes Brumptomyiina, Lutzomyiina, and Psychodopygina [4,5]. It is a complex zoonosis whose main reservoirs include dogs, rodents, and other wild mammals. Clinical manifestations in humans vary according to the parasite species involved and may present as localized cutaneous, diffuse cutaneous, mucocutaneous, and visceral forms—the latter being considered the most severe [6].

It is estimated that approximately 30,000 new cases of visceral leishmaniasis and more than 1 million new cases of cutaneous leishmaniasis occur annually worldwide [7]. In Brazil, recent data indicate a substantial impact of the disease, particularly in the Northern Region, where the diversity of species of the subgenus Leishmania (Viannia) is high, including L. (V.) braziliensis, L. (V.) guyanensis, L. (V.) shawi, L. (V.) lainsoni, L. (V.) naiffi, and L. (V.) lindenbergi [8,9].

Leishmania (V.) naiffi, initially described by Lainson and Shaw in 1989, was first identified in a vertebrate host, the nine-banded armadillo (Dasypus novemcinctus), which is widely distributed throughout the Amazon region. Although human cases were subsequently confirmed, infections caused by this species have likely been historically underdiagnosed due to their generic classification as Leishmania (Viannia) spp. [10,11]. Clinical manifestations associated with L. naiffi tend to be localized cutaneous forms, usually presenting as a single lesion, and display particular features, such as low pathogenicity in classical animal models [8,12].

By contrast, Leishmania (V.) shawi was first isolated in 1984 from Cebus apella (now Sapajus apella), an arboreal New World monkey. It is mainly associated with cutaneous leishmaniasis and may present with either single or multiple lesions, possibly related to lymphatic dissemination. This species has been isolated from several wild reservoirs, including primates, coatis, and sloths, reinforcing its ecological and epidemiological importance in the Amazon region [11].

Despite the clinical and epidemiological relevance of leishmaniasis, the available therapeutic options remain limited and are often associated with high toxicity, elevated costs, prolonged treatment regimens, and the emergence of parasite resistance [13]. In this context, the identification of new therapeutic targets represents an urgent need for the development of more effective and safer strategies.

Advances in genomic approaches have made a decisive contribution to the understanding of the biology of parasites of the genus Leishmania, revealing high levels of genomic plasticity, variations in gene content, chromosomal aneuploidy, and intraspecific diversity, in addition to enabling the rational identification of potential pharmacological targets [14]. Comparative genomics and pangenomics studies are therefore fundamental tools for elucidating evolutionary mechanisms, host adaptation, and functional differences among closely related species.

In this context, the present study aimed to present the complete genome sequencing and comparative genomic analysis of the strains L. naiffi (MDAS/BR/1979/M5533) and L. shawi (MCEB/BR/1984/M8408), exploring similarities and differences related to gene annotation, pangenome composition, chromosomal ploidy, genetic variability, genomic coverage, variant calling, phylogenetic inference, and the identification of potential therapeutic targets based on protein homology.

2. Materials and Methods

2.1. DNA Extraction and Sequencing

The reference strains Leishmania naiffi and Leishmania shawi were provided by the Leishmania Collection (CLIOC), which is part of the Network of Biological Resource Centers for Conformity Assessment of Biological Material at the Instituto Oswaldo Cruz, Rio de Janeiro, Brazil. The parasites were maintained in Schneider’s medium supplemented with 20% fetal bovine serum at 25 °C. Information regarding the number of in vitro passages was not available. Genomic DNA from the L. naiffi (MDAS/BR/1979/M5533) and L. shawi (MCEB/BR/1984/M8408) strains was extracted using the Wizard^®^ Genomic DNA Purification Kit (Promega, Madison, WI, USA). The DNA samples were extracted and rehydrated to a final volume of 50 µL for each sample. Genomic libraries were prepared using the Nextera XT DNA Library Preparation Kit (Illumina, Inc., San Diego, CA, USA) and sequenced on the NextSeq 500 platform (Illumina, Inc., San Diego, CA, USA) using the NextSeq 500/550 High Output v2.5 Kit (300 cycles) in a paired-end sequencing format for both samples. All laboratory procedures were performed according to the manufacturer’s instructions, unless otherwise specified.

2.2. Read Quality Assessment

Read quality was assessed using FastQC version v0.12.1 [15], and adapter removal as well as filtering of low-quality reads were performed using Fastp (version v0.12.1) [16].

2.3. Genome Assembly

Genome assembly was performed using a de novo strategy with MEGAHIT v1.2.9 [17], applying default parameters and automatic k-mer selection (21, 29, 39, 59, 79, 99, 119, and 141). MEGAHIT employs succinct de Bruijn graphs (SdBGs), which reduce memory usage and enable efficient processing of large datasets. Its iterative multi–k-mer approach favors the recovery of low-coverage regions and the resolution of repetitive sequences, while the use of mercy k-mers ensures the inclusion of low-abundance regions, which are common in metagenomic data.

2.4. Taxonomic Classification

The contigs resulting from the assembly of L. naiffi and L. shawi were subjected to a local alignment against a custom database created from complete genomes of the genus Leishmania available on NCBI. For this analysis, the BLASTn (v.2.5.0) algorithm [18] was used.

Taxonomic classification of the contigs was performed using Kraken2 software (v.2.1.2) [19], which allowed the identification and separation of contigs related to the Leishmania genus and the assignment of taxonomic labels to each sequence. This classification enabled the extraction of all contigs related to L. naiffi and L. shawi.

To visualize and analyze the taxonomic classification data of the contigs, the Pavian application [20] was used, which is a tool for exploring metagenomic classification results for pathogen detection, in this case, Leishmania.

Subsequently, the SSPACE v3.0 program was used to order and join contigs and generate scaffolds [21].

2.5. Genome Annotation

Genome annotation was performed using AUGUSTUS v3.5 [22], employing the genomic structure of Leishmania (Sauroleishmania) tarentolae as a reference. A local genomic database comprising 69 Leishmania genomes (Supplementary Table S1) was constructed from GenBank to improve genome annotations and to analyze and refine the positions and structures of open reading frames (ORFs). Geneious v8.1.4 [23] was used for genome visualization and for the manual sequence edition of the identified ORFs after assessing the quality and integrity of each ORF.

To identify and classify genes according to their molecular and biological functions, the PANTHER database [24] was used. Through its web service, genes were functionally annotated based on their associated biological processes and molecular functions.

2.6. Orthologous Gene

Genomic data from L. naiffi and L. shawi were used to analyze the pangenome and identify gene orthologs by comparison with assembled genomes available in the NCBI database for L. (V.) guyanensis (MHOM/BR/75/M4147) and Leishmania (L.) major (MHOM/IL/80/Friedlin). These analyses were performed using BLASTp v2.5.0 [25] and the R statistical software v4.1.2 with the venn package [26].

2.7. Search for Pharmacological Targets

To identify genes with potential to serve as pharmacological targets in L. naiffi and L. shawi, a sequence homology-based approach using amino acid sequences was adopted. Initially, a local database was constructed containing protein sequences derived from three-dimensional structures previously associated with the rational development of drugs against leishmaniasis, including calpains (PDB IDs: 1TLO, 1MDW, 1ZCN), DYRK1A (3ANQ), GSK-3 (7S6V), and an ABC family transporter (7OJ8).

Predicted protein sequences derived from the genomes of L. naiffi and L. shawi were compared against this database using BLASTp, from the BLAST+ suite (version 2.15.0) [27].

BLASTp results were subjected to a set of minimum filtering criteria to ensure consistency in the identification of potential functional homologs. Only alignments simultaneously meeting the following parameters were retained: (i) amino acid identity ≥ 35%, (ii) alignment length ≥ 60 residues, (iii) e-value ≤ 1 × 10^−3^, and (iv) bitscore ≥ 30. These thresholds are widely accepted for the detection of moderately conserved functional homologies, allowing the exclusion of short or statistically weak alignments [28].

Approved hits were subsequently ranked based on a combined score integrating information on percentage identity, alignment coverage, statistical significance (logarithm of the e-value), and bitscore. For each analyzed genome, the top 200 hits were selected according to this score. Finally, all results were consolidated into a single file, enabling a global comparative interpretation among the analyzed species.

2.8. Reference Mapping, Genomic Coverage, and Variant Calling

Reference mapping of the L. naiffi and L. shawi genomes was performed using Bowtie2 [29], enabling the assessment of chromosomal coverage. Read mapping also allowed the inference of chromosomal ploidy, calculated as the ratio between the coverage of each individual chromosome and half of the global mean coverage across the 35 chromosomes—characteristic of the Leishmania (Viannia) subgenus—as expected for a diploid organism.

Single-nucleotide polymorphisms (SNPs) and insertions/deletions (indels) were identified using BCFtools v1.2 [30] to detect nucleotide differences between the sequenced strains L. naiffi and L. shawi and the corresponding reference genomes available in GenBank (L. naiffi: GCA_962239345; L. shawi: GCA_962240455.1).

2.9. Phylogenetic Inference Methods

The alpha catalytic subunit of DNA polymerase (polA1) was selected as the molecular target for phylogenetic inference of Leishmania, based on its phylogenetic signal evaluated using the TREE-PUZZLE v5.3 algorithm [31] with the likelihood mapping method. A total of 65 Leishmania genomes, together with the genomes of L. naiffi (MDAS/BR/1979/M5533) and L. shawi (MCEB/BR/1984/M8408), were included in this analysis.

The polA1 gene was identified and extracted using Geneious v8.1.4 [23], and phylogenetic analysis was performed using the maximum likelihood (ML) method implemented in IQ-TREE2 [32], following selection of the most appropriate nucleotide substitution model. Parameters were optimized to identify the most likely tree topology, with robustness assessed by bootstrap analysis (1000 replicates). The best phylogenetic tree was selected based on genetic identity thresholds among Leishmania species, inferred from a sequence identity matrix constructed using the software. All programs were run with default parameters unless otherwise specified.

3. Results

3.1. Quality Control, Genome Assembly, and Taxonomic Classification

The results of read trimming, the total number of assembled contigs, and the taxonomic classification of contigs for L. naiffi and L. shawi are shown in Table 1.

The genomes of L. naiffi and L. shawi were sequenced, yielding 16,008,124 and 10,855,510 reads, respectively. After the quality control (trimming) step, the final number of reads was 15,103,553 for L. naiffi and 9,741,698 for L. shawi.

During genome assembly, 20,542 contigs were obtained for L. naiffi and 13,849 for L. shawi. Following taxonomic classification, these numbers were refined to 20,446 and 13,816 contigs, respectively. The final assembled genome sizes were 32.13 Mb for L. naiffi and 32.51 Mb for L. shawi.

For both species, scaffolds had a minimum length of 200 bp, whereas the maximum scaffold length was 17,834 bp for L. naiffi and 49,040 bp for L. shawi. The N50 values, which indicate the median scaffold length, were 2583 bp for L. naiffi and 5882 bp for L. shawi.

3.2. Genome Annotation Analysis

Based on genome annotation performed using AUGUSTUS, open reading frames (ORFs) were established for L. naiffi, in which 8170 genes were identified, compared with 7767 genes identified in L. shawi. Of these, only 2935 genes in L. naiffi were assigned known functions, whereas this number was slightly higher in L. shawi, reaching 3033 genes.

In the prediction of molecular functions (Table 2), most proteins from L. shawi (62.0%) and L. naiffi (61.1%) could not be assigned to functional categories in the PANTHER database. Among the classified proteins, the most representative functions corresponded to catalytic activity (16.6% in L. shawi and 16.8% in L. naiffi) and binding activity (12.1% and 12.3%, respectively). Other categories, such as ATP-dependent activity, transporter activity, and structural molecule activity, were detected at lower proportions but showed similar distributions between the two genomes (Table 2).

Regarding biological processes (Table 3), a high proportion of proteins could not be classified (52.7% in L. shawi and 52.0% in L. naiffi). Among the identified categories, cellular processes (20.0% and 20.4%) and metabolic processes (14.6% and 14.7%) were the most prominent, respectively. Processes related to localization, biological regulation, and response to stimuli were observed at intermediate proportions, whereas categories such as reproduction, homeostasis, and development were poorly represented, with nearly identical frequencies between the analyzed species (Table 3).

3.3. Orthologous Gene Evaluation

Pangenome analysis of Leishmania species, represented by the Venn diagram, showed that the core genome—defined as the set of genes shared among all analyzed species—comprised 6256 genes. In addition, species-specific gene sets and genes partially shared among subsets of species were identified (Figure 1).

L. naiffi and L. shawi share a set of 6681 conserved genes. Despite this similarity, both species also exhibit particular features in this analysis: L. naiffi possesses 46 exclusive genes, whereas L. shawi presents 25 species-specific genes.

3.4. Identification of Pharmacological Targets

Protein homology analysis enabled the identification of 21 genes with therapeutic potential, of which 11 were identified in L. naiffi and 10 in L. shawi (Table 4). Among the candidate genes, members of the GSK-3, calpain, and ABC transporter families were particularly prominent—protein classes that are widely explored as pharmacological targets in trypanosomatids.

The distribution of quality scores among the candidate genes revealed quantitative differences between the two species (Figure 2). In L. naiffi, the highest score was observed for gene g6653.t1_1, annotated as GSK-3, whereas in L. shawi the highest score corresponded to gene g2262.t1_1, also belonging to the GSK-3 family. These findings reinforce the relevance of this protein family as a conserved and potentially exploitable target in both analyzed species.

Genes classified with intermediate scores predominantly corresponded to partial fragments or conserved domains of GSK-3 and calpains, whereas the sole representative of the ABC transporter family showed a lower score, reflecting a shorter alignment length and moderate sequence identity, although still within the established minimum criteria.

3.5. Chromosomal Coverage, Ploidy, and Variant Calling

For L. naiffi, chromosomal coverage ranged from 22.55× (chromosome 28) to 68.53× (chromosome 2), with most chromosomes showing coverage values concentrated between approximately 24× (chromosomes 6, 14, 15, 21, 23, and 32) and 26× (chromosomes 1, 5, 7, 11, 18, and 22). In contrast, L. shawi exhibited higher chromosomal coverage, ranging from 50.76× (chromosome 35) to 165.46× (chromosome 8), with a predominance of values around 63× (chromosomes 14, 20, and 24) (Figure 3).

Regarding ploidy (Figure 4), L. naiffi exhibited variation ranging from 1.59 (chromosome 28) to approximately 4.85 (chromosome 2), whereas L. shawi showed values ranging from approximately 1.24 (chromosome 35) to 4.04 (chromosome 8). In both genomes, most chromosomes displayed a ploidy pattern close to 2, consistent with a predominantly diploid state (Figure 4).

The identification of the number of SNPs (Figure 5) and indels (Figure 6) across the chromosomes of L. naiffi and L. shawi revealed differences between the two species. In L. naiffi, SNP counts varied among chromosomes, ranging from 243 in chromosome 8 to 3041 in chromosome 31. In L. shawi, SNP counts ranged from 224 in chromosome 4 to 2313 in chromosome 35 (Figure 4). Across the entire genome, genetic variability analysis identified 34,480 SNPs in L. naiffi and 26,562 SNPs in L. shawi, indicating greater genomic diversity in the former species.

Similarly, the number of indels also varied across chromosomes in both species. In L. naiffi, indel counts ranged from 91 in chromosome 2 to 704 in chromosome 31. In L. shawi, the number of indels ranged from 75 in chromosome 2 to 791 in chromosome 35 (Figure 6).

In the genome of L. naiffi, chromosome 31 stood out for having the highest variant counts, with 3041 SNPs and 704 indels. Other chromosomes, such as 23 and 27, also showed elevated numbers of these genetic markers. In L. shawi, the distribution of SNPs and indels across the genome was similarly heterogeneous, with chromosome 31 highlighted for having 2074 SNPs and 462 indels.

The duplicated chromosomes (20.1 and 20.2) in L. naiffi also exhibited a high number of genetic variations, particularly chromosome 20.1, with 1721 SNPs and 452 indels. In L. shawi, the duplicated chromosomes showed a similar pattern, with chromosome 20.1 concentrating 1477 SNPs and 415 indels. Across the whole genome, 34,480 SNPs and 9104 indels were identified in L. naiffi, whereas L. shawi presented 26,562 SNPs and 8046 indels. Analyses of genomic coverage, ploidy, and variant calling for L. naiffi and L. shawi are provided in Supplementary Table S2.

3.6. Phylogenetic Tree Reconstruction and Topology

Phylogenetic analysis of the genomes of 67 Leishmania strains, using the polA1 gene as a molecular marker, produced a maximum likelihood tree (Figure 7). This analysis showed that L. naiffi clustered within the Leishmania (Viannia) subgenus with strong bootstrap support (bootstrap >90%), although positioned outside the L. guyanensis and L. braziliensis complexes. Conversely, L. shawi grouped within the clade corresponding to the L. guyanensis complex, also with strong bootstrap support (>90%).

4. Discussion

4.1. Genome Assembly, Taxonomic Classification, and Scaffold

The genome sizes were very similar (32.13 Mb and 32.50 Mb), yet the assemblies of L. naiffi and L. shawi revealed notable structural differences. The L. naiffi assembly was more fragmented, exhibiting a higher total number of contigs and scaffolds (20,446) compared to L. shawi (13,816). This pattern of increased genomic fragmentation in L. naiffi is consistent with prior observations from a comparison involving L. guyanensis [33] (Supplementary Table S3).

Despite higher coverage for L. naiffi (70.5×) compared to L. shawi (45×), this did not translate into greater assembly continuity for the former. This observation reinforces the previously documented structural constraints characteristic of the subgenus Viannia, which include significant genomic plasticity and structural variations [34]. Such features have also been confirmed in the phylogenetically related species L. braziliensis and L. panamensis [35].

The higher N25, N50, and N75 values for L. shawi compared to L. naiffi indicate a more contiguous assembly for the former. This pattern mirrors that observed for L. guyanensis when compared to the L. naiffi assembly in a prior genomic study [33].

The maximum contig sizes, 49,040 bp in L. shawi versus 17,834 bp in L. naiffi, further highlight this difference, clearly demonstrating the greater continuity of the former assembly compared to the latter.

The GC content was nearly identical between the two species, at approximately 57% of the genome. This value falls within the typical range reported for the subgenus Leishmania (Viannia) (57–57.5%) and aligns with data reported for L. peruviana [36].

4.2. Comparison of Genome Annotation Analysis

Functional analysis of genes via BLASTp and PANTHER revealed a highly conserved profile between L. shawi and L. naiffi, characterized by a predominance of unclassified genes (62.0% and 61.1%, respectively). This pattern, previously noted in other genomes of the subgenus Leishmania (Viannia) [33], reflects the large proportion of proteins that remain without detailed functional annotation. Catalytic functions represented the most abundant class (16.6% in L. shawi and 16.8% in L. naiffi), indicating conservation of essential metabolic processes. This finding is consistent with reports of amplified enzymatic activities in Leishmania (Viannia) subgenus, as also observed for L. panamensis [35].

Minor differences, such as the higher absolute number of structural proteins in L. naiffi (84 versus 76), may reflect genomic particularities previously associated with the structural plasticity of Leishmania, as also discussed elsewhere [33,35,36].

The analysis of biological processes corroborated this pattern of conservation: unclassified proteins were predominant (52.7% and 52%), as previously reported [33]. Cellular processes (20% and 20.5%) and metabolic processes (14.6% and 14.7%) showed similar proportions between the species, consistent with the flexible metabolism of Leishmania (Viannia) discussed elsewhere [35]. Categories associated with response to stimuli, regulation, and localization remained stable, reinforcing that functional differences likely arise mainly from mechanisms such as amplifications, aneuploidies, and copy-number variations, widely recognized in Leishmania (Viannia) subgenus [33,34]. The numerical proximity observed in these categories is compatible with the homogeneous life cycle of the genus and aligns with annotations described in the literature [37,38].

4.3. Orthologous Genes Comparison

Pangenome analysis revealed a well-conserved core genome of 6256 genes shared between L. naiffi, L. shawi, and two additional reference species (L. major and L. guyanensis). This conserved core underscores the essential functional repertoire maintained across the Leishmania genus. The size of this core is aligned with previous genomic studies focusing on the Viannia subgenus, which reported core sizes of 6635 [39] and 6784 [38] genes when analyzing different species combinations, including those from the Leishmania (Leishmania) subgenus. This consistency reinforces the concept of strong functional and syntenic conservation within the genus.

However, the core genome size identified here is smaller than the 7392 and 7157 genes reported in other studies [34,35]. This variation is likely attributable to differences in methodological pipelines, stringency of orthology assignment, or the distinct phylogenetic breadth of the genome datasets analyzed.

Notably, the analysis also delineated species-specific gene sets, with 46 unique genes in L. naiffi and 25 in L. shawi. These accessory genes may underpin adaptive traits particular to each species, potentially influencing mechanisms of host–parasite interaction, niche-specific environmental responses, or metabolic adaptations.

4.4. Potential Pharmacological Targets

The recurrent identification of genes belonging to the GSK-3 family among the candidates with the highest quality scores in L. naiffi and L. shawi reinforces the relevance of this kinase as a high-priority therapeutic target in Leishmania species. Recent studies have highlighted GSK-3 as a core regulatory protein for essential cellular processes in trypanosomatids, including cell cycle control, parasite survival, and host adaptation. These studies also provide robust genetic and pharmacological validation supporting its potential as a promising drug target [40,41].

The high degree of conservation observed for GSK-3 genes in the two analyzed species is consistent with their functional essentiality, a characteristic frequently associated with therapeutic targets of high biological value. It is important to note that this high conservation does not preclude its pharmacological exploitation. Structural, regulatory, and kinetic differences relative to the mammalian host homolog have already been described, enabling the achievement of molecular selectivity [42]. In this context, the prioritization of GSK-3 family genes observed in this study, based on objective criteria of sequence homology and quality score ranking, finds strong support in the specialized literature. These findings underpin the potential of GSK-3 as a conserved and exploitable therapeutic target, with applicability in drug intervention strategies aimed at multiple Leishmania species. This reinforces its relevance in the landscape of the rational development of new antileishmanial agents [40,41,42].

4.5. Coverage, Ploidy, and Variant Calling

Analysis of chromosomal coverage and ploidy reinforces the genomic plasticity of the subgenus. L. naiffi exhibited coverages between 22.55× and 68.53×, whereas L. shawi showed higher values (50.76× to 165.46×), suggesting greater overall coverage and potential for variant detection. Both species displayed extensive aneuploidies (ploidies from ~1.2 to ~4.8), a pattern consistent with the chromosomal mosaicism characteristic of the subgenus, as documented in prior studies [33,34].

The predominance of ploidy near 2 reflects a baseline diploid karyotype, also reported for L. panamensis and L. braziliensis [35]. SNP and indel analyses revealed marked heterogeneity. In L. naiffi, chromosome 31 stood out with 3.041 SNPs, while in L. shawi, chromosomes 31 and 35 showed a high density of variants, suggesting regions subject to distinct selective pressures. The duplicated chromosomes (20.1 and 20.2) exhibited a strong accumulation of variants in both species, reinforcing the role of duplications and structural divergence, as proposed in existing models of genomic plasticity [33]. Such patterns corroborate that chromosomal instability in Leishmania is an adaptive functional mechanism, as reported elsewhere in the literature [43].

4.6. Evolutionary Implications and Phylogenetic Relationships

Phylogenetic analysis based on the polA1 gene consistently placed L. naiffi and L. shawi within the subgenus Leishmania (Viannia) subgenus, reinforcing the previously described evolutionary patterns for this lineage. L. naiffi was robustly grouped within this subgenus but occupied an external phylogenetic position relative to the L. guyanensis and L. braziliensis complexes, indicating a distinct evolutionary trajectory. This placement is compatible with the high genomic plasticity and structural particularities already described for Leishmania (Viannia) species [33,34].

In contrast, L. shawi consistently grouped within the L. guyanensis complex, with strong bootstrap support (>90%), corroborating its close phylogenetic relationship with other species in this complex. The clear separation between the Leishmania (Leishmania) and Leishmania (Viannia) subgenera observed in the tree is congruent with previously published multi-locus analyses [35] and with the structural and evolutionary differences described for South American species [36].

Collectively, these results indicate that although L. naiffi and L. shawi share a conserved genetic core, they follow distinct evolutionary paths, which are reflected in both their pangenome organization and their patterns of chromosomal variation.

5. Conclusions

This study provides two new genomic assemblies for L. naiffi (MDAS/BR/1979/M5533) and L. shawi (MCEB/BR/1984/M8408), expanding the landscape for comparative analyses within the Leishmania (Viannia) subgenus. These resources establish the foundation for developing species-specific diagnostic tools based on unique molecular markers and enable future research on gene expression, novel drug discovery, and pathogenicity. Thus, the genomes presented herein not only contribute to the fundamental knowledge of these species but also pave the way for targeted clinical interventions and integrated molecular surveillance strategies for controlling cutaneous leishmaniasis in the Amazon region. The results of this study should be interpreted as exploratory and hypothesis-generating. Future studies incorporating multiple isolates per species and long-read sequencing technologies will be required to validate species-level genomic features.