Genomic and epidemiological evidence for the emergence of a L. infantum/L. donovani hybrid with unusual epidemiology in northern Italy

| Genetics and Molecular Biology | Research Article Genomic and epidemiological evidence for the emergence of a L. infantum/L. donovani hybrid with unusual epidemiology in northern Italy F. Bruno,1 G. Castelli,1 B. Li,2 S. Reale,1 E. Carra,3 F. Vitale,1 S. Scibetta,1 M. Calzolari,3 S. Varani,4 M. Ortalli,4,5 E. Franceschini,6 W. Gennari,7 G. Rugna,3 G. F. Späth8 AUTHOR AFFILIATIONS See affiliation list on p. 16. ABSTRACT Leishmania (L.) infantum is one of the main causative agents of animal and human leishmaniasis across many endemic areas in South America, Europe, North Africa, and Asia. Despite its clinical significance, little is known about the genetic diversity of L. infantum circulating in a given endemic area. Here, we investigate this important open question by applying a comparative genomics approach to seven L. infantum isolates from different hosts and Italian regions, including the northern part of the country (Emilia-Romagna, RER), Sicily, and Sardinia, as an initial attempt to explore the breadth of parasite genetic heterogeneity in Italy. Additionally, microsatellite analysis was carried out to compare the isolates from RER with other 70 L. infantum strains from the same region as well as 65 strains belonging to the L. donovani complex from other countries. We revealed important karyotypic instability and identified strain-specific changes in gene dosage, which affected important virulence factors such as amastins and surface antigen-like proteins. Single nucleotide polymorphism-based clustering analysis of these genomes together with over 80 publicly available L. infantum and L. donovani genomes placed the Italian isolates into three geographically distinct clusters within the Mediterra nean basin and uncovered three isolates clustering with putative L. infantum/L. donovani hybrids isolated in Cyprus. As judged by microsatellite profiling, these hybrid isolates are representative of a sub-population of parasites circulating in northern Italy that preferentially infect humans but not dogs. Our results place Italy at the crossroads of L. infantum infection in the Mediterranean and call attention to the public health risk represented by the introduction of non-European Leishmania species. IMPORTANCE This study closes important knowledge gaps with respect to Leishmania (L.) infantum genetic heterogeneity in a given endemic country, as exemplified here for Italy, and reveals genetic hybridization as a main cause for re-emerging human leishmaniasis in northern Italy. The observed high diversity of Leishmania parasites on the Italian peninsula suggests different geographical origins, with genomic adaptation to various ecologies affecting both pathogenicity and transmission potential. This is documented by the discovery of a putative L. infantum/L. donovani hybrid strain, which has been shown to preferentially infect humans but not dogs. Our results provide important information to health authorities, which need to consider the public health risk represented by the introduction of new Leishmania species into EU countries due to population displacement or travel from countries where exotic/allochthonous parasite species are endemic. KEYWORDS Leishmania, comparative genomics, genetic hybrid, visceral leishmaniasis July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 1 Invited Editor David L. Sacks, NIH, NIAID, Bethesda, Maryland, USA Editor L. David Sibley, Washington University in St. Louis School of Medicine, St. Louis, Missouri, USA Address correspondence to G. F. Späth, [email protected], or G. Rugna, [email protected]. F. Bruno and G. Castelli contributed equally to this article. Author order was determined based on alphabetical order. The authors declare no conflict of interest. See the funding table on p. 16. Received 4 April 2024 Accepted 30 April 2024 Published 4 June 2024 Copyright © 2024 Bruno et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license. Leishmaniases are a group of sand fly-transmitted parasitic diseases that can affect mammals, including humans. In humans, these diseases are characterized by cutaneous, muco-cutaneous, or visceral tissue damage, the latter often being fatal if left untreated. The pathogenesis and clinical outcome of Leishmania infection rely on a series of parasite-specific traits that ensure proliferation and fitness gain in different host systems. First, these parasites have evolved different life cycle stages, including motile promastigote forms that are adapted for survival inside the midgut of phlebotomine sand flies and intracellular amastigotes that infect immune cells of the reticuloendothe lial system inside the mammalian host. Second, aside stage differentiation, Leishmania can further adapt to environmental variations encountered inside both insect and vertebrate hosts; this phenomenon has been linked to the intrinsic plasticity of the Leishmania genome and the genetic recombination that can follow the formation of parasite hybrids between strains and species inside the insect vector. In the wake of next-generation sequencing and the development of powerful computational pipelines for sequencing data analysis (1), comparative genomics approaches have provided important insight into the epidemiology of Leishmania infection in the field or parasite fitness gain in experimental settings. For example, comparative genomics has informed on (i) the mechanisms underlying drug resistance and parasite evolution in treated patients (2), (ii) the origin, population structure, and phylogenetic relationship of clinical isolates, (iii) hybridization events and their effect on tissue tropism or clinical outcome of infection (3, 4), (iv) the role of dynamic gene dosage changes in parasite fitness gain (5), or (v) the genomics of geographic adaptation, which is the scope of the current study. Leishmania infantum causes leishmaniasis in Mediterranean Europe, with significant emergence of human cases observed in localized areas of Spain (6) and Italy (7). Italy is hypo-endemic for human leishmaniasis, with less than 100 cases of cutaneous and visceral leishmaniasis (VL) reported to the WHO in 2020 (8), while dogs are consid ered the main reservoir hosts. Historically, leishmaniasis in Italy was restricted to the Tyrrhenian littoral, the southern peninsula, and the islands. However, in the last 20 years, sand fly vectors as well as human and canine Leishmania infections have been detected in northern Italy, traditionally classified as a cold area unsuitable for sand fly survival (9). In the last decade, a surge of human leishmaniasis cases has been observed in the Emilia-Romagna region (RER), northern Italy (7). However, there was no concurrent increase in the prevalence of canine leishmaniasis in the same period. In addition, molecular typing studies in RER showed that the Leishmania strains circulating in dogs belonged to a different population compared to strains isolated from human VL cases and sand flies (10). Furthermore, the biting preference of Phlebotomus (Ph) perfiliewi, the suspected vector responsible for the increase of VL cases in RER, suggested the presence of a peri-urban or sylvatic reservoir other than dogs, which is capable of spreading the infection in such areas (11). Consistent with this observation, a high frequency of Leishmania infection was recently identified in other potential reservoir hosts, including roe deer, hares, wolves, or red foxes (12). While some studies have evaluated the epidemiology of VL in Italy, genetic characteri zation of L. infantum across different regions is lacking. Here, we combined comparative genomics and microsatellite profiling approaches on L. infantum field isolates that were obtained in Italy to shed light on parasite genetic heterogeneity. We provide evidence for a remarkable genetic diversity of L. infantum across the Italian peninsula and identify the strains from RER as putative L. infantum/L. donovani hybrids, which are closely related to Cypriot strains, thus correlating their unusual epidemiology and transmission pattern to their unique genetic constitution. Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 2 RESULTS Cultured L. infantum isolates show strain-specific and convergent aneuploi dies We applied a comparative genomics approach to gain insight into the genetic diversity of L. infantum field isolates in Italy. We deliberately chose initially five isolates from different hosts (human, dog, cat, and marten) and geographic regions (RER in northern Italy, Sicily, and Sardinia) to maximize the assessment of parasite genetic heterogeneity (Table 1; Fig. 1A). Following culture expansion, DNA extraction, and Illumina short read sequencing, we applied our genome instability pipeline (GIP) (1) to map the reads on the JPCM5 L. infantum reference genome (13) and assess read depth variations for each individual sample. We then compared chromosome read depth variations between samples using per chromosome boxplots of normalized coverage (somy score) in 300 bp genomic bins (see Materials and Methods). As judged by the changes in the somy score, all samples showed important karyotypic variations, thus confirming the intrinsic instability of the L. infantum genome as previously reported (2) (Fig. S1). As expected from previous studies, all samples were tetrasomic for chromosome (chr) 31 (14) while displaying unique karyotypic profiles considering the other chromosomes. This diversity may be largely attributed to the strain-specific selection of different aneuploidies during culture adaptation, which we previously linked to in vitro fitness gain (5). The heat map shown in Fig. 1B reveals leish4, a strain isolated from a Sicilian dog from Palermo (Table 1), as the karyotypically most stable sample, with most chromosomes showing a median somy score between 2 and 2.5. All other samples showed more important somy variations, including full trisomies or higher, mosaic aneuplodies. Our data revealed pentasomies for chr 8 and 23 in the leish16 sample. Given that this sample represents the first L. infantum genome sequenced from a feline isolate, we investigated if this unusual amplification pattern may be selected in this particular host. This possibility was ruled out by the quantification of chr 8 copy number by quantitative PCR (qPCR) analysis of three independent feline isolates, which showed a disomic amplification score comparable to disomic chr 28 included as a control (Fig. 1C). Thus, the pentasomy for chr 8 may either have been uniquely selected in leish16 during feline infection or underwent important amplification during the 16 culture passages of this sample (see Table 1). Of note, amplification of chr 8 and 23 is quite common in cultured promastigotes of other species. They have been described as among the most frequent polysomal chromosomes in L. donovani isolates (15) and have been shown to arise during L. major self-hybridization (16). Thus, the Italian L. infantum field isolates show important somy variations in culture compared to the haploid reference genome. Convergent amplification of a small number of chromosomes in most or all strains suggests that these aneuploidies are a potential TABLE 1 Overview of the L. infantum isolates that were analyzed in this study by whole-genome sequencing Lab ID International code Host/location (province)/date/ culture passage (p) Information Sequencing reads: number/% mapped leish3 MCAN /IT/2021/51327 Dog/Sardinia (Cagliari)/2021/p8 Glucantime resistant 27,739,245/0.98 leish4 MCAN/IT/2020/1265 Dog/Sicily (Palermo)/2020/p8 –a25,856,894/0.91 leish5 MMST/IT/2006/V2921 Marten/Sicily (Palermo)/2006/p15 – 21,183,149/0.98 leish7 MHOM/IT/2014/IZSLER-MO22 Human/Emilia-Romagna (Modena)/ 2014/p8 VL case; comorbidity: solid neoplasia 26,421,239/0.98 leish16 MFEL/IT/2018/10816 Cat/Sicily (Palermo)/2 018/p16 – 25,578,672/0.81 leishMO23 MHOM/IT/2014/IZSLER-MO23 Human/Emilia-Romagna (Modena)/ 2014/p8 VL case; comorbidity: hematologic disorder 22,041,069/0.98 leishMO38 MHOM/IT/2016/IZSLER-MO38 Human/Emilia-Romagna (Bologna)/ 2016/p8 VL case; comorbidity: HIV positive 20,902,208/0.92 a"–" indicates that information is not available. Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 3 driver of L. infantum fitness gain in culture. Surprisingly, similar chromosomes were linked to culture adaptation of L. donovani isolates from the Indian sub-continent and Sudan (e.g., chr 5, 20, and 26) (5), suggesting conserved fitness mechanisms between both parasite species when selected for growth in culture. Aside from these convergent aneuploidies, we further observed strain-specific, karyotypic changes that may fine-tune adaptation in response to other genetic differences between the strains, which we further analyzed in the following section. FIG 1 Karyotypic analyses of the L. infantum isolates (n = 5). (A) Map indicating the geographic location and origin of the L. infantum isolates investigated in this study. (B) Heatmap showing the median somy score for the different strains across the 36 chromosomes. Darker tonalities reflect higher somy values. (C) Quantitative PCR analysis with SYBR Green assay to quantify chr 8 copy number in feline Leishmania strains. CT values are shown for the histone deacetylase gene on chr 8, and the glycerophosphoryl diester phosphodiesterase gene on disomic chr 28 used as a control. Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 4 Gene copy number variations in the L. infantum isolates affect known Leishmania virulence genes We next examined the gene-level coverage information generated by GIP to assess strain-specific and convergent gene dosage changes with respect to the JPCM5 reference genome (supplemental material; Fig. S2; and Data Set S1A). Unlike karyotypic changes that are rapidly selected during short-term culture adaptation (5–20 passages), gene copy number variations (CNVs) are much less dynamic and thus can reveal gene dosage changes that were selected in the field. Plotting the normalized read depth coverage per gene, we observed a wavy pattern for all isolates caused by minor read depth variations across the chromosomes (Fig. 2A; supplemental material; Fig. S2). A series of gene CNVs were observed, which correspond to the amplification or deple tion of gene copies with respect to the reference genome, but also in between the five samples. Significantly, these variations were not scattered across the genome but localized in defined genomic regions that were similar in all samples, thus revealing possible hot spots of gene dosage changes (Fig. 2A). Closer inspection of the 50 most significant gene CNVs identified both single gene and sub-chromosomal amplifications, many of which have been previously linked to infectivity, including amastins (17), surface antigen-like proteins (18), or genes implicated in the synthesis of phosphogly can virulence factors (19), such as phosphoglycan beta 1–3 galactosyltransferase genes (LINF_020006900 and LINF_020007000) or ppg3 (LINF_350010100 and LINF_350010200) encoding filamentous proteophosphoglycan (20) (supplemental material; Fig. S3A; Data Set S1B). Significantly, most gene amplifications are observed across several isolates. Such convergence is not only observed at the gene level but also at the level of gene function, with amplified genes often sharing the same annotation even though they are on different chromosomes and encoding different proteins, as shown, for example, for the genes encoding surface antigen proteins and amastins (supplemental material; Fig. S3B). We also identified 116 gene depletions (normalized mean coverage < 0.5) and deletions (normalized mean coverage = 0) that were shared or strain-specific in the five L. infantum isolates (Fig. 2B). Most of these changes occurred in multi-copy gene arrays (Data Set S1C), where recombination events between identical sequences allow for dynamic gene copy number changes. In particular, such dynamic variations between the strains were observed for (i) five gene copies on chr 9 encoding for the autophagy gene ATG8, which has been implicated in Leishmania stress response and infectivity (21), (ii) six gene copies on chr 22 encoding for an uncharacterized, hypothetical protein, (iii) a cluster of nine amastin genes encoded on chr 34 corresponding to two multi-copy gene arrays (Fig. S3B) and that were previously classified as δ-amastins by Jackson (22) known to control amastigote infectivity (17), (iv) 10 HSP70 gene copies encoded on chr 28 that we previously revealed as a potential hot spot of environment-genotype interaction in L. donovani field isolates (23), and (v) 10 genes on chr 34 encoding for paraflagellar rod protein (Fig. 2C). The reduced read depth that we observed for individual gene copies could indicate mosaic gene CNVs, i.e., mixed populations of wild-type and deleted genotypes, or represent a populating-wide, heterozygous state. Only a few genes show true deletions as judged by the absence of reads, including LINF_270010900 likely encoding for a calpain-like protease as suggested by a BLAST search (data not shown), which is deleted in leish16, leish3, leish5, and leish7, or the ATG8 gene LINF_190013600 that is deleted in leish5 (Fig. 2C). The presence of these deletions within non-clonal, heterogeneous parasite populations suggests selection against these genes, even though loss due to genetic drift cannot be ruled out. Discovery of leish7 as a highly divergent L. infantum strain We next assessed the evolutionary relationship between our five L. infantum isolates by comparing single nucleotide polymorphism (SNP) numbers, positions, and frequencies from the GIP output. Based on the number of alternative (alt) alleles, leish16 is the most Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 5 FIG 2 Gene copy number analyses of the L. infantum isolates (n = 5). (A) Scatterplot showing the normalized mean coverage (y-axis) for each gene across the 36 chromosomes (x-axis). Different strains are shown in different colors (see legend in the figure). (B) Heatmap of normalized mean coverage. The plot shows gene depletions defined by a normalized mean coverage < 0.5 (green, dark green corresponding to deletions with read depth = 0) and gene amplifications defined by a normalized mean coverage > 1.5 (red). The column on the left indicates the MAPQ score (red, MAPQ < 40; black, MAPQ > 40). (C) Radar plots of the normalized mean coverage values for the genes indicated. Different strains are shown in different colors (see legend in the figure). The expected normalized coverage with respect to the reference genome is 1. Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 6 FIG 3 Genome-wide SNP analyses of the L. infantum isolates (n = 5). (A) Upset plot of SNPs identified in the different strains. The bars visualize the number of SNPs common to a given combination of samples. The numbers in the bar plot indicate shared SNPs for the given combination. The numbers on the left side of the strain ID indicate the total number of SNPs in each isolate. (B) Violin plot representing density estimates of the distribution of alt-allele frequencies for the indicated samples. Strip plots are superimposed on the violin plots in order to show individual loci at intermediate frequencies. Only those loci where at least one sample had an alt-allele frequency between 0.2 and 0.8 are represented. (C) Cluster analysis based on pairwise Euclidean distances between alt-allele frequencies across SNPs. The tree is constructed from these distances using the UPGMA method. (D) The scatterplots display individual SNPs as dots. The x-axis indicates the relative position of each SNP along the genome, while the y-axis indicates the variant allele frequency (see Fig. S4 for individual panels). The different chromosomes are displayed one after the other, and their boundaries are visualized as vertical lines. The individual SNPs are colored according to the samples as indicated by the legend in the graph. The dense superposition of dots corresponding to SNPs at frequency 0 or 1 results in blurry lines at y = 0 and y = 1. There are no SNP in the (0, 0.1) frequency interval due to the filtering applied during the SNP calling step of the GIP pipeline (the corresponding frequency is set to 0 when data pertaining to various samples are merged together). Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 7 similar to the JPCM5 reference genome with only 451 SNPs identified, while leish7 is the most divergent with 49,053 SNPs (Fig. 3A and B; Data Set S1D). Of note, including Illumina reads of the JPCM5 reference strain itself (24) revealed a total of 387 SNPs, of which 267 were unique, while 120 were shared with one or more of the isolates. Cluster analysis based on these variants revealed the highest similarities between JPCM5, leish4, and leish16 on one hand, and leish3 and leish5 on the other hand, while leish7 represen ted its own branch, with similar distance from either of the two other clusters (Fig. 3C). Plotting frequency against location revealed defined patches of heterozygous SNPs at a frequency of about 50% in leish7, suggesting the presence of distinct haplotypes that may be remnants of former hybridization events (Fig. 3D, see also individual plots shown in Fig. S4). Investigating the possible origin of the leish7 putative hybrid We next established the genome sequence for two additional isolates (leishMO23 and leishMO38) obtained from human VL cases from RER, which have been previously shown to cluster with leish7 (originally named MO22) (10). Cluster analyses revealed their close relationship with leish7, confirming the transmission of a distinct parasite sub-popula tion—likely a hybrid—in RER (Fig. 4A, top right inset, and Data Set S1D). We further expanded the cluster analysis including 11 genomes that we previously generated as part of the LeiSHield project (www.leishield.org) (2) and 63 genomes analyzed by Franssen et al. (15). While four of our Italian L. infantum strains clustered with isolates from Tunisia (leish3 and leish5) and Spain (leish4 and leish16), the three RER hybridlike strains (leish7, leishMO23, and leishMO38) clustered with the known L. infantum/L. donovani hybrids first described in Cyprus (Linf_CH33, 35, 36 and Ldo_CH33). Compari son of SNPs between these hybrids revealed over 34,000 shared SNPs (Fig. 4B), indicating that these geographically distinct hybrids are likely the result of similar hybridization events between L. infantum and L. donovani parental strains. This common hybrid nature is further supported using the Kraken software, a sequence classification tool developed for metagenomic analyses, which assigns taxonomic labels to DNA sequences using the exact alignment of k-mers (25). This approach allowed us to classify all strains (including the Italian and Cypriot hybrids) as belonging to the L. donovani complex, with only the hybrids revealing a significant match at the species level to L. donovani genomes (Fig. 4C). However, whether the RER and Cypriot hybrids are of common origin remains to be established. Population structure and phylogenetic analysis of putative hybrid strains from northern Italy In order to better define the prevalence of the leish7 hybrid-like genotype in northern Italy, microsatellite analysis was performed on 73 L. infantum strains, including 22 strains obtained from human VL cases, 8 strains obtained from sand flies, and 40 strains obtained from canine cases from various areas of RER (Fig. 5). Bayesian clustering as implemented in STRUCTURE assigned the 73 multilocus microsatellite typing (MLMT) profiles from these L. infantum strains to two main populations (the highest value of ΔK was for K = 2; Fig. S5), namely, population A (PopA) that consisted of canine L. infantum strains only and population B (PopB) that consisted of all the L. infantum strains from VL (including leish7, leishMO23, and leishMO38) and sand flies (Fig. 5). A minor ΔK peak at K = 5 (Fig. S5) indicated a substructure in the whole data set, which was first explored by increasing the number of populations from K = 3 to K = 5: PopA exhibited subdivision into two groups at K = 3 and three groups at K = 4, while PopB split into two groups at K = 5 (Fig. S6). Subsequently, we analyzed each of the two main populations (PopA and PopB) separately; the calculation of ΔK values (Fig. 5) showed the presence of four subpopulations: PopA1, PopA2, PopB1 (including leishMO38), and PopB2 (including leish7 and leishMO23), while a fifth sub-cluster that was indicated by the minor peak at K = 5 in the first STRUCTURE analysis (all strains) was not confirmed. Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 8 The spatial distribution of the genetic groups showed that the two main popula tions (PopA and PopB) overlapped (Fig. 5), while partial overlap was observed for the four subpopulations: PopA1 was widespread all over RER, PopB1 was prevalent in the central-western areas, PopB2 was distributed in central-eastern areas, and PopA2 was limited to an eastern area of RER (Fig. S7). The measures of genetic diversity for each main population defined by STRUCTURE are summarized in Table 2. PopA and PopB showed a low expected heterozygosity (mean He = 0.258), which was higher than the observed heterozygosity (mean Ho = 0.056). Furthermore, the inbreeding coefficients (Fis) were very high (mean Fis = 0.807). The global phylogenetic analysis showing the position of human and sand fly strains from RER within the L. donovani complex is represented in Fig. 6. Canine L. infantum strains from RER belonging to the original PopA were grouped in a monophyletic clade together with VL strains from other Italian regions and Mediterranean human and canine strains belonging to L. infantum zymodemes MON-1, MON-72, and MON-98. On the contrary, L. infantum strains from VL patients and sand flies from RER belonging to PopB grouped in a monophyletic clade that was intermediate between L. infantum and L. donovani clusters. In addition, the PopB clade was close to the cluster, which included two strains from Cyprus, namely, CH35 and CD44cl.1, both classified as L. donovani MON-37 (26). DISCUSSION The present study combines comparative genomics with molecular epidemiology approaches to deliver a comprehensive, high-resolution investigation of intraspecific genetic diversity and heterozygosity within the L. infantum species in Italy, one of the most important Leishmania species worldwide from both a veterinary and public health perspective. We revealed a remarkable heterogeneity in Italian L. infantum isolates that are related to various parasite genotypes of different geographic origins, including a hybrid-like genotype first described in isolates from Cyprus. Our results place Italy at the crossroads of L. infantum infection in the Mediterra nean and inform on the possible origin of regional sub-populations and their local transmission cycles. We identified a series of L. infantum strains from Sardinia and Sicily that are genetically close to those from Spain (southwestern[SW] Europe) and Tunisia (North Africa), highlighting Italy’s central location in the Mediterranean Basin and its historical and current role as one of the main routes of human migration. Based on genetic similarity, the strains leish3 and leish5 likely belong to the North African MON-1 population endemic in southern and eastern Mediterranean regions. This population is genetically different from the European MON-1 population found in SW Europe and South America, which corresponds to the strains leish4 and leish16 and includes the JPCM5 strain used as a reference in this study. The complex Leishmania epidemiology is the result of the parasite’s transmission dynamics, the distribution and geographic abundance of competent vectors, past and current exposure of human populations to the parasite, environmental parameters that impact the seasonal activity and distribution of sand flies (e.g., air temperature, relative humidity, ground cover, altitude, and vegetation), and the presence of reservoir hosts. In Italy, as in the whole Mediterranean basin, dogs (Canis lupus familiaris) are incrimi nated as the primary reservoir hosts of zoonotic human leishmaniasis, even if other domestic and wild mammals infected with L. infantum could play a possible role in the biological cycle of the parasite (27). When compared to L. donovani species, L. infantum MON-1, the most common zymodeme in dogs and humans, has been considered a genomically rather homogeneous population regardless of geographic distribution (15). The parasite’s adaptation to dogs may have been an evolutionary advantage, and the frequent mobility of the canine reservoir may have contributed to the spread of the L. infantum MON-1 group (28). However, microsatellite analyses and WGS showed that the MON-1 “core” L. infantum clade displays a certain degree of isolation by distance. Three distinct populations have been shown to circulate within the Mediterranean Region, Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-24 9 AUTHOR AFFILIATIONS 1WOAH Leishmania Reference Laboratory, Istituto Zooprofilattico Sperimentale della Sicilia, Centro di Referenza Nazionale per le Leishmaniosi (C.Re.Na.L.), Palermo, Italy 2Bioinformatics and Biostatistics Hub, Institut Pasteur, Université Paris Cité, Paris, France 3Istituto Zooprofilattico Sperimentale della Lombardia e dell'Emilia Romagna "B. Ubertini", Brescia, Italy 4Department of Medical and Surgical Sciences, University of Bologna, Bologna, Italy 5IRCCS Azienda Ospedaliero-Universitaria di Bologna, Bologna, Italy 6Infectious Disease Unit, Azienda Ospedaliera Universitaria di Modena, Modena, Italy 7Virology and Molecular Microbiology Unit, University Hospital of Modena, Modena, Italy 8Unité de Parasitologie moléculaire et Signalisation, INSERM U1201, Institut Pasteur, Université Paris Cité, Paris, France AUTHOR ORCIDs M. Calzolari http://orcid.org/0000-0001-7489-2556 G. Rugna http://orcid.org/0000-0001-7437-1929 G. F. Späth http://orcid.org/0000-0002-0256-2029 FUNDING Funder Grant(s) Author(s) EC | Horizon 2020 Framework Programme (H2020) LeiSHield-MATI-REP-778298-1 G. F. Späth Ministero della Salute (Italy Ministry of Health) E54I19002870001 G. Rugna Ministero della Salute (Italy Ministry of Health) IZS SI RC SI 04/21 F. Bruno Ministero della Salute (Italy Ministry of Health) GR‐ 2019‐ 12369134 G. Castelli European Commission (EC) Project no. PE00000007 INF-ACT S. Varani DATA AVAILABILITY All relevant data are within the paper and its supplemental files. Genome sequence data have been deposited in BioProject at PRJNA1008390. ETHICS APPROVAL Samples from human cases were collected for diagnostic purposes. Storing and molecular evaluation of retrospectively collected samples received approvals from both the Ethics Committee of the St. Orsola-Malpighi University Hospital, Bologna, Italy (protocol No. 3729/2017, last amendment approved: EM414-2023 97/2017/O/Tess/ AOUBO) and the Modena Ethics Committee, Modena, Italy (protocol No. 239/13). Samples were coded and anonymized. All animal samples were obtained for diagnostic purposes with no unnecessary invasive procedures. Evaluation of these samples did not require ethical approval according to the European Directive 2010/63/EU. ADDITIONAL FILES The following material is available online. Supplemental Material Data Set S1 (mBio00995-24-s0001.xlsx). WGS supporting information. Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-2416 Data Set S2 (mBio00995-24-s0002.xlsx). MLMT supporting information. Supplemental material (mBio00995-24-s0003.docx). Supporting methods and supplemental figures. REFERENCES 1. Späth GF, Bussotti G. 2022. GIP: an open-source computational pipeline for mapping genomic instability from protists to cancer cells. Nucleic Acids Res 50:e36. https://doi.org/10.1093/nar/gkab1237 2. Bussotti G, Benkahla A, Jeddi F, Souiaï O, Aoun K, Späth GF, Bouratbine A. 2020. Nuclear and mitochondrial genome sequencing of NorthAfrican Leishmania infantum isolates from cured and relapsed visceral leishmaniasis patients reveals variations correlating with geography and phenotype. Microb Genom 6:mgen000444. https://doi.org/10.1099/ mgen.0.000444 3. Lypaczewski P, Matlashewski G. 2021. Leishmania donovani hybridisation and introgression in nature: a comparative genomic investigation. Lancet Microbe 2:e250–e258. https://doi.org/10.1016/S26665247(21)00028-8 4. Van den Broeck F, Savill NJ, Imamura H, Sanders M, Maes I, Cooper S, Mateus D, Jara M, Adaui V, Arevalo J, Llanos-Cuentas A, Garcia L, Cupolillo E, Miles M, Berriman M, Schnaufer A, Cotton JA, Dujardin J-C. 2020. Ecological divergence and hybridization of neotropical Leishmania parasites. Proc Natl Acad Sci U S A 117:25159–25168. https://doi.org/10. 1073/pnas.1920136117 5. Prieto Barja P, Pescher P, Bussotti G, Dumetz F, Imamura H, Kedra D, Domagalska M, Chaumeau V, Himmelbauer H, Pages M, Sterkers Y, Dujardin J-C, Notredame C, Späth GF. 2017. Haplotype selection as an adaptive mechanism in the protozoan pathogen Leishmania donovani. Nat Ecol Evol 1:1961–1969. https://doi.org/10.1038/s41559-017-0361-x 6. Horrillo L, Castro A, Matía B, Molina L, García-Martínez J, Jaqueti J, García-Arata I, Carrillo E, Moreno J, Ruiz-Giardin JM, San Martín J. 2019. Clinical aspects of visceral leishmaniasis caused by L. infantum in adults. Ten years of experience of the largest outbreak in Europe: what have we learned? Parasit Vectors 12:359. https://doi.org/10.1186/s13071-0193628-z 7. Todeschini R, Musti MA, Pandolfi P, Troncatti M, Baldini M, Resi D, Natalini S, Bergamini F, Galletti G, Santi A, Rossi A, Rugna G, Granozzi B, Attard L, Gaspari V, Liguori G, Ortalli M, Varani S. 2024. Re-emergence of human leishmaniasis in northern Italy, 2004 to 2022: a retrospective analysis. Euro Surveill 29:2300190. https://doi.org/10.2807/1560-7917.ES. 2024.29.4.2300190 8. Ruiz-Postigo JA, Jain S, Mikhailov A, Maia-Elkhoury AN, Valadas S, Warusavithana S, Osman M, Lin Z, Beshah A, Yajima A. 2021. Global leishmaniasis surveillance: 2019-2020, a baseline for the 2030 roadmap/ surveillance mondiale de la leishmaniose: 2019-2020, une periode de reference pour La Feuille de route a L'horizon 2030. WER 96:401–420. 9. Gradoni L, Ferroglio E, Zanet S, Mignone W, Venco L, Bongiorno G, Fiorentino E, Cassini R, Grillini M, Simonato G, Michelutti A, Montarsi F, Natale A, Gizzarelli M, Foglia Manzillo V, Solari Basano F, Nazzari R, Melideo O, Gatti D, Oliva G. 2022. Monitoring and detection of new endemic foci of canine leishmaniosis in northern continental Italy: an update from a study involving five regions (2018-2019). Vet Parasitol Reg Stud Rep 27:100676. https://doi.org/10.1016/j.vprsr.2021.100676 10. Rugna G, Carra E, Bergamini F, Calzolari M, Salvatore D, Corpus F, Gennari W, Baldelli R, Fabbi M, Natalini S, Vitale F, Varani S, Merialdi G. 2018. Multilocus microsatellite typing (MLMT) reveals host-related population structure in Leishmania infantum from northeastern Italy. PLoS Negl Trop Dis 12:e0006595. https://doi.org/10.1371/journal.pntd. 0006595 11. Calzolari M, Romeo G, Bergamini F, Dottori M, Rugna G, Carra E. 2022. Host preference and Leishmania infantum natural infection of the sand fly Phlebotomus perfiliewi in northern Italy. Acta Trop 226:106246. https:/ /doi.org/10.1016/j.actatropica.2021.106246 12. Taddei R, Bregoli A, Galletti G, Carra E, Fiorentini L, Fontana MC, Frasnelli M, Musto C, Pupillo G, Reggiani A, Santi A, Rossi A, Tamba M, Calzolari M, Rugna G. 2022. Wildlife hosts of Leishmania infantum in a re-emerging focus of human leishmaniasis, in Emilia-Romagna, Northeast Italy. Pathogens 11:1308. https://doi.org/10.3390/pathogens11111308 13. Peacock CS, Seeger K, Harris D, Murphy L, Ruiz JC, Quail MA, Peters N, Adlem E, Tivey A, Aslett M, et al. 2007. Comparative genomic analysis of three Leishmania species that cause diverse human disease. Nat Genet 39:839–847. https://doi.org/10.1038/ng2053 14. Rogers MB, Hilley JD, Dickens NJ, Wilkes J, Bates PA, Depledge DP, Harris D, Her Y, Herzyk P, Imamura H, Otto TD, Sanders M, Seeger K, Dujardin JC, Berriman M, Smith DF, Hertz-Fowler C, Mottram JC. 2011. Chromo some and gene copy number variation allow major structural change between species and strains of Leishmania. Genome Res 21:2129–2142. https://doi.org/10.1101/gr.122945.111 15. Franssen SU, Durrant C, Stark O, Moser B, Downing T, Imamura H, Dujardin J-C, Sanders MJ, Mauricio I, Miles MA, Schnur LF, Jaffe CL, Nasereddin A, Schallig H, Yeo M, Bhattacharyya T, Alam MZ, Berriman M, Wirth T, Schönian G, Cotton JA. 2020. Global genome diversity of the Leishmania donovani complex. Elife 9:e51243. https://doi.org/10.7554/ eLife.51243 16. Ferreira TR, Inbar E, Shaik J, Jeffrey BM, Ghosh K, Dobson DE, Beverley SM, Sacks D. 2022. Self-hybridization in Leishmania major. mBio 13:e0285822. https://doi.org/10.1128/mbio.02858-22 17. Rochette A, McNicoll F, Girard J, Breton M, Leblanc E, Bergeron MG, Papadopoulou B. 2005. Characterization and developmental gene regulation of a large gene family encoding amastin surface proteins in Leishmania spp. Mol Biochem Parasitol 140:205–220. https://doi.org/10. 1016/j.molbiopara.2005.01.006 18. Bhandari V, Kumar D, Verma S, Srividya G, Negi NS, Singh R, Salotra P. 2013. Increased parasite surface antigen-2 expression in clinical isolates of Leishmania donovani augments antimony resistance. Biochem Biophys Res Commun 440:646–651. https://doi.org/10.1016/j.bbrc.2013. 09.113 19. Rogers ME. 2012. The role of Leishmania proteophosphoglycans in sand fly transmission and infection of the mammalian host. Front Microbiol 3:223. https://doi.org/10.3389/fmicb.2012.00223 20. Depledge DP, Evans KJ, Ivens AC, Aziz N, Maroof A, Kaye PM, Smith DF. 2009. Comparative expression profiling of Leishmania: modulation in gene expression between species and in different host genetic backgrounds. PLoS Negl Trop Dis 3:e476. https://doi.org/10.1371/ journal.pntd.0000476 21. Giri S, Shaha C. 2019. Leishmania donovani parasite requires Atg8 protein for infectivity and survival under stress. Cell Death Dis 10:808. https:// doi.org/10.1038/s41419-019-2038-7 22. Jackson AP. 2010. The evolution of amastin surface glycoproteins in trypanosomatid parasites. Mol Biol Evol 27:33–45. https://doi.org/10. 1093/molbev/msp214 23. Drini S, Criscuolo A, Lechat P, Imamura H, Skalický T, Rachidi N, Lukeš J, Dujardin J-C, Späth GF. 2016. Speciesand strain-specific adaptation of the HSP70 super family in pathogenic trypanosomatids. Genome Biol Evol 8:1980–1995. https://doi.org/10.1093/gbe/evw140 24. González-de la Fuente S, Peiró-Pastor R, Rastrojo A, Moreno J, CarrascoRamiro F, Requena JM, Aguado B. 2017. Resequencing of the Leishmania infantum (strain JPCM5) genome and de novo assembly into 36 contigs. Sci Rep 7:18050. https://doi.org/10.1038/s41598-017-18374-y 25. Wood DE, Lu J, Langmead B. 2019. Improved metagenomic analysis with Kraken 2. Genome Biol 20:257. https://doi.org/10.1186/s13059-0191891-0 26. Gouzelou E, Haralambous C, Amro A, Mentis A, Pratlong F, Dedet J-P, Votypka J, Volf P, Toz SO, Kuhls K, Schönian G, Soteriadou K. 2012. Multilocus microsatellite typing (MLMT) of strains from Turkey and Cyprus reveals a novel monophyletic L. donovani sensu lato group. PLoS Negl Trop Dis 6:e1507. https://doi.org/10.1371/journal.pntd.0001507 27. Azami-Conesa I, Gómez-Muñoz MT, Martínez-Díaz RA. 2021. A systematic review (1990-2021) of wild animals infected with zoonotic Leishmania. Microorganisms 9:1101. https://doi.org/10.3390/ microorganisms9051101 28. Maia C, Cardoso L. 2015. Spread of Leishmania infantum in Europe with dog travelling. Vet Parasitol 213:2–11. https://doi.org/10.1016/j.vetpar. 2015.05.003 Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-2417 29. Antoniou M, Haralambous C, Mazeris A, Pratlong F, Dedet J-P, Soteriadou K. 2009. Leishmania donovani leishmaniasis in Cyprus. Lancet Infect Dis 9:76–77. https://doi.org/10.1016/S1473-3099(09)70004-0 30. Calzolari M, Carra E, Rugna G, Bonilauri P, Bergamini F, Bellini R, Varani S, Dottori M. 2019. Isolation and molecular typing of Leishmania infantum from Phlebotomus perfiliewi in a re-emerging focus of leishmaniasis, northeastern Italy. Microorganisms 7:644. https://doi.org/10.3390/ microorganisms7120644 31. Dujardin J-C, Campino L, Cañavate C, Dedet J-P, Gradoni L, Soteriadou K, Mazeris A, Ozbel Y, Boelaert M. 2008. Spread of vector-borne diseases and neglect of leishmaniasis, Europe. Emerg Infect Dis 14:1013–1018. https://doi.org/10.3201/eid1407.071589 32. Lex A, Gehlenborg N, Strobelt H, Vuillemot R, Pfister H. 2014. UpSet: visualization of intersecting SETS. IEEE Trans Vis Comput Graph 20:1983– 1992. https://doi.org/10.1109/TVCG.2014.2346248 33. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, Burovski E, Peterson P, Weckesser W, Bright J, et al. 2020. SciPy 1.0: fundamental algorithms for scientific computing in python. Nat Methods 17:261–272. https://doi.org/10.1038/s41592-019-0686-2 34. McKinney W. 2010. “Data structures for statistical computing in python” Python in Science Conference; Austin, Texas: p 51–56. https://doi.org/10. 25080/Majora-92bf1922-00a 35. Reback J, McKinney W, Den Bossche J, Augspurger T, Cloud P, Klein A, Hawkins S, Roeschke M, Tratner J, She C. 2020. Pandas-Dev/Pandas: Pandas 1.0. 5. Zenodo 36. Hunter JD. 2007. Matplotlib: A 2D graphics environment. Comput Sci Eng 9:90–95. https://doi.org/10.1109/MCSE.2007.55 37. Waskom ML. 2021. Seaborn: statistical data visualization. JOSS 6:3021. https://doi.org/10.21105/joss.03021 Research Article mBio July 2024 Volume 15 Issue 7 10.1128/mbio.00995-2418