Selective whole-genome amplification reveals population genetics of Leishmania braziliensis directly from patient skin biopsies

Abstract

In Brazil, Leishmania braziliensis is the main causative agent of the neglected tropical disease, cutaneous leishmaniasis (CL). CL presents on a spectrum of disease severity with a high rate of treatment failure. Yet the parasite factors that contribute to disease presentation and treatment outcome are not well understood, in part because successfully isolating and culturing parasites from patient lesions remains a major technical challenge. Here we describe the development of selective whole genome amplification (SWGA) for Leishmania and show that this method enables culture-independent analysis of parasite genomes obtained directly from primary patient skin samples, allowing us to circumvent artifacts associated with adaptation to culture. We show that SWGA can be applied to multiple Leishmania species residing in different host species, suggesting that this method is broadly useful in both experimental infection models and clinical studies. SWGA carried out directly on skin biopsies collected from patients in Corte de Pedra, Bahia, Brazil, showed extensive genomic diversity. Finally, as a proof-of-concept, we demonstrated that SWGA data can be integrated with published whole genome data from cultured parasite isolates to identify variants unique to specific geographic regions in Brazil where treatment failure rates are known to be high. SWGA provides a relatively simple method to generate Leishmania genomes directly from patient samples, unlocking the potential to link parasite genetics with host clinical phenotypes.

Fig 1. SWGA primer design and evaluation.

(A) The number of exact match ‘hits’ per megabase (Mbp) for each of the 23 identified SWGA primers against Leishmania and host reference genomes, and (B) the fold difference in exact matches against L. braziliensis compared to human, mouse, or canine genomes. (C) Heatmap showing percent reads aligning to L. braziliensis for each of the four SWGA primer sets used to carry out SWGA on known ratios of L. braziliensis DNA spiked into human genomic DNA (0.1 and 1% final parasite DNA). (D) The number of reads is shown in relation to the percentage of the parasite genome covered at ≥1x (blue line), 5x (red line) and 10x (green line). Vertical dashed line indicates a sequencing effort of 100 x 10⁶ 150bp paired-end reads.

Fig 2. In vivo validation of SWGA for Leishmania.

Percentage of reads mapping to L. braziliensis genome in DNA from (A) infected mouse ears (n = 3 animals infected with L. major, triangles; n = 5 animals infected with L. braziliensis, circles) or (B) patient lesion biopsies, sequenced before (pre) and after (post) SWGA. Data shown are from the SWGA primer set the yielded the best amplification for each sample. (C) Genome coverage for SWGA data from a single patient sample (patient #7, blue point from panel B). (D) Coverage of four selected L. braziliensis chromosomes in SWGA data from a single patient (#7; blue lines) compared to whole genome sequencing (WGS) of pure, cultured L. braziliensis (orange lines). Data shown in panel C and D are merged from all SWGA primer sets to maximize coverage.

Fig 3. Allele frequency determined by SWGA.

Alternate Allele Read Depth Proportion (AARDP) histograms for L. braziliensis chromosomes 10, 23, 28, and 31, for (A) whole genome sequencing (WGS) of pure cultured parasites, (B) SWGA of pure cultures, (C-D) SWGA of synthetic controls consisting of 1% (C) or 0.1% (D) parasite DNA, (E-F) SWGA on two patient samples from Fig 2B. Peaks centered on 0.5 indicate disomic chromosomes, while peaks at approximately 0.25, 0.5 and 0.75 indicate tetrasomic chromosomes. Green, red, and blue dashed lines denote an AARDP of 0.25, 0.5, and 0.75, respectively.

Fig 4. Scalable SWGA profiling of patient samples.

(A) QPCR is used to prioritize samples that have the highest parasite burden and, therefore, the greatest likelihood of success for SWGA. (B) SWGA is carried out in 96-well plates using multiple primer sets and primer set combinations (plate rows) for each patient (plate columns). (C) Shallow sequencing is used to determine which samples showed the best amplification by SWGA. (D) All successful SWGA reactions are pooled for each patient and (E) subjected to deep sequencing. (F) Results of selective whole genome amplification of L. braziliensis from 18 primary patient samples.

Fig 5. Integrating SWGA and WGS genomes for population genomics.

(A) Map showing all 59 samples, from this study and four previously published reports, included in the analysis [7, 46, 8, 47]. (B) Zoomed in view of Bahia, Brazil showing region covered by samples from this study. White point indicates position of field hospital where patients were seen. (C-D) Principal component analysis of SNP data from 59 genomes, colored by country of origin. (E) Maximum likelihood tree constructed using 877713 variants from 59 L. braziliensis genomes and the L. guyanensis outgroup, compared to the L. braziliensis reference. Branch length of outgroup was shortened for figure preparation. Tree is rooted using the L. guyanensis outgroup. The same cultured laboratory clone of L. braziliensis from Brazil was sequenced either by traditional WGS (black circle) or SWGA (white triangle). Map data from Maps Mapbox (www.mapbox.com/about/maps) and OpenStreetMap (www.openstreetmap.org/about).

Fig 6. Identification of variants unique to Northeastern Brazil.

(A) Table showing variants identified by integrated analysis of WGS and SWGA genomes (top), and studies included (+) or excluded (-) from the analysis (bottom). Venn diagrams indicate how each of the five studies (labeled a-d; [7, 46, 8, 47]) were used in the integrated analysis to generate the variants shown in table column above. (B) Bubble chart showing results of Gene Ontology (GO) enrichment for Molecular Function terms associated with 149 genes containing frame-shift variants (left) or 152 genes identified with high-frequency missense mutations in Northeast (NE) Brazil (right). All terms shown were associated with ≥ 5 genes. FC = fold change; FDR = false discovery rate (Benjamini-Hochberg correction). (C) Four representative parasite genes that were enriched for high-frequency missense mutations in genomes from Northeast Brazil.