Design and implementation of multiplexed amplicon sequencing panels to serve genomic epidemiology of infectious disease: a malaria case study

Panel designs

We developed a small multiplex of four highly polymorphic antigenic loci, dubbed ‘4CAST’: CSP, AMA1, SERA2, and TRAP (Fig. 2). All four amplicons use previously published primer sequences (Miller et al., 2017; Neafsey et al., 2015), as no modification was required for successful multiplexing.

• Download figure
• Open in new tab

Figure 2. Global characterization of loci in the 4CAST and AMPLseq panels.

Estimates of diversity of each locus in the 4CAST and AMPLseq panels, with one locus per row. We estimated haplotypic diversity from monoclonal P. falciparum WGS data from each country. The top 4 loci shown represent the 4CAST loci, which are also included in the AMPLseq panel. All 128 P. falciparum loci in the AMPLseq panel are shown; the single P. vivax locus is not shown.

In designing the larger ‘AMPLseq’ multiplexed amplicon panel, we first built a large pool of candidate loci, anticipating significant attrition of candidates due to primer incompatibility. We prioritized four classes of loci: loci within antigens of interest (Helb et al., 2015), loci with high population diversity for relatedness inference (Taylor et al., 2019), loci included in the SpotMalaria v1 panel (Chang et al., 2019; Jacob et al., 2021), and known drug resistance markers. We contracted the services of GTseek LLC (https://gtseek.com) to design multiplexed oligo panels according to the criteria previously described for the Genotyping-in-Thousands by sequencing (GT-seq) protocol (Campbell, Harmon, & Narum, 2015) (S1 Supporting information). We optimized the final primer set and and reaction conditions through several sequencing runs and determined that the primers for the four 4CAST loci (CSP, AMA1, SERA2, TRAP) could be added to the panel without impacting amplification of the other loci. We also successfully added primers amplifying known markers of antimalarial drug resistance within the genes dhfr, dhps, mdr1, and kelch13 (S2 Table). Furthermore, we successfully added previously described primers for PvDHFR (Lefterova, Budvytiene, Sandlund, Färnert, & Banaei, 2015) in order to identify P. falciparum / P. vivax co-infections that have gone undetected in preliminary screening by microscopy or rapid diagnostic test (RDT). The final panel, dubbed ‘AMPLseq’ (short for Assorted Mix of Plasmodium Loci) contains this single P. vivax locus and 128 P. falciparum loci (Fig. 2), with a median length across all amplicons of 276 bp (S1 Fig.).

Panel protocols

To create the primer pool used in 4CAST PCR1, we combined 100 µM of each 4CAST primer (S3 Table) and diluted the combined primer mix to 6.25 µM per primer in nuclease-free water (NF dH₂O). Each 10.5 µl PCR1 reaction incorporated 1.5 µl combined primer mix, 5 µl KAPA HiFi HotStart ReadyMix (2x), and 4 µl sample template. PCR1 amplification consisted of an initial incubation step at 95 °C (3 min); 25 amplification cycles at 95 °C (20 s), 57 °C (15 s) and 62 °C (30 s); and a final extension step at 72 °C (1 min). Each 12.2 µl PCR2 reaction (which adds sample indices and sequencing adapters) incorporated 2.2 µl unique dual index (10 µM Illumina Nextera DNA UD Indexes), 5 µl KAPA HiFi HotStart ReadyMix (2x), 2 µl NF dH₂O and 3 µl PCR1 product. PCR2 consisted of an initial incubation step at 95 °C (1 min); 10 amplification cycles at 95 °C (15 s), 55 °C (15 s) and 72 °C (30 s); and a final extension step at 72 °C (1 min). We combined PCR2 products in equal volumes and performed double-sided size selection using Agencourt AMPure XP beads (Beckman Coulter): we incubated 100 µl library with 55 µl beads, immobilized beads via magnet rack, and transferred the supernatant to a new tube. We incubated the transferred supernatant with 20 µl beads and washed immobilized beads twice with 80% ethanol. We eluted the library in 25 µl EB buffer (10 mM Tris-Cl, pH 8.5), subsequently adding 2.5 µl EB buffer with 1% Tween-20. We verified size selection via Agilent BioAnalyzer 2100 and sequenced the selected library at 6 pM with >10% PhiX in paired-end, 500-cycle format using MiSeq Reagent Kit v2 (S1 Protocol).

We followed a similar nested PCR and pooled clean-up procedure for AMPLseq library construction. Primer sequences, input volumes and concentrations are listed in S3 Table and PCR conditions and size selection steps are shown in S2 Protocol. As detailed therein, library construction for AMPLseq library construction differs in a few minor aspects. For example, primer input quantities vary slightly (800 pmol +/- 33%) to account for amplification rate differences among loci. PCR1 products are diluted 1:12 in NF dH₂O prior to PCR2 and only single-sided (left-tailed) bead-based size selection is used to enhance yield. Sequencing also occurs via paired-end, 500-cycle MiSeq but with a higher final library loading concentration (12 pM) and a lower fraction of PhiX (8%).

Mock samples

We generated mock samples from parasite lines 3D7 and Dd2, cultured at 3% hematocrit in commercially obtained red blood cells as previously described (Trager & Jensen, 1976). We extracted genomic DNA using the Qiagen Blood and Tissue Kit on cells previously lysed with 0.15% saponin. We generated positive control template representing DNA extractions from whole human blood infected with 10,000 monoclonal 3D7 parasites/µl by diluting genomic DNA from 3D7 to 0.25 ng/µl in 10 ng/µl human genomic DNA, and storing in 10 mM Tris-HCl (pH 8.0), 1 mM EDTA (Promega, Madison, WI). We generated further control templates representing 1000, 100, and 10 3D7 parasites/µl by serial 1:10 dilution of the 10,000 3D7 parasites/µl control, likewise using 10 ng/µl human genomic DNA as diluent. We also generated a 10,000 parasites/µl positive control as described above but using Dd2 instead of 3D7 strain genomic DNA. We generated mixed-strain control templates by combining the 10,000 3D7 parasites/µl control with this 10,000 Dd2 parasites/µl control at 1:1, 3:1, and 10:1 ratios (respectively). We serially diluted the 1:1 ratio to 1000, 100, and 10 parasites/µl concentrations and diluted the 3:1 and 10:1 ratios to 1000 and 100 parasites/µl concentrations using 10 ng/µl human genomic DNA diluent as before. We also applied selective whole genome amplification (sWGA) to all above control templates representing ≤ 1000 parasites/µl. The 50 µl sWGA reaction followed Oyola et al. 2016 (Oyola et al., 2016) with the exception of fixing template input volume to 10 µL. We purified sWGA products with Agencourt AMPure XP beads (Beckman Coulter) on the KingFisher Flex (S3 Protocol) and verified amplification success via NanoDrop (ThermoFisher Scientific).

Clinical samples

We tested the panels on clinical dried blood spot (DBS) samples from Mali and Guyana. Tran et al. collected samples in Kalifabougou, Mali in 2011 – 2013 as previously described (Tran et al., 2013). The Kalifabougou cohort study was approved by the Ethics Committee of the Faculty of Medicine, Pharmacy and Dentistry at the University of Sciences, Technique and Technology of Bamako, and the Institutional Review Board of the National Institute of Allergy and Infectious Diseases, National Institutes of Health (NIH IRB protocol number: 11IN126; https://clinicaltrials.gov/; trial number NCT01322581). Written informed consent was obtained from participants or parents or guardians of participating children before inclusion in the study. The Guyana Ministry of Health collected samples from Port Kaituma and Georgetown, Guyana between May – August 2020, by spotting participants’ whole blood onto Whatman FTA cards and storing the samples with individual desiccant packets at room temperature. Informed consent (or parental assent for minors) was obtained for all subjects according to protocols approved by ethical committees.

We punched DBS samples 3 – 5 times into a 96-well deep well plate using the DBS pneumatic card puncher (Analytical Sales and Services, Inc.) equipped with a 3 mm cutter. We then extracted gDNA following the DNA purification from buccal swab section of the KingFisher Ready DNA Ultra 2.0 Prefilled Plates for KingFisher Flex instruments protocol (ThermoFisher Scientific) with minor modifications (S4 Protocol). We used the same sWGA procedure as above on the extracted gDNA.

Whole genome sequencing and variant calling

We performed whole genome sequencing on clinical samples collected in Guyana to validate 4CAST and AMPLseq outcomes. We performed sWGA on DNA samples as described above to enrich parasite DNA. We used the enriched DNA to construct Illumina sequencing libraries from the amplified material using the NEBNext Ultra II FS DNA prep kit (NEB #E6177) prior to sequencing on an Illumina HiSeqX instrument at the Broad Institute, using 150 bp paired-end reads, targeting a sequencing depth of at least 50X coverage. We aligned reads to the P. falciparum v3 reference genome assembly using BWA-MEM (Li, 2013) and called SNPs and INDELs using the GATK HaplotypeCaller (DePristo et al., 2011; McKenna et al., 2010; Van der Auwera et al., 2013) according to the best practices for P. falciparum as determined by the Pf3k consortium (https://www.malariagen.net/resource/34). Analyses were limited to the callable segments of the genome (Miles et al., 2016) and excluded sites where over 20% of samples were multiallelic. Data from these samples were submitted to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under accession PRJNA758191.

Amplicon data analysis

We developed an application named AmpSeQC (S2 Supporting information) to assess sequence quality and amplicon/sequence run success (S2 Fig.). We also used AmpSeQC for P. vivax detection by applying a concatenation of the P. falciparum 3D7 and P. vivax PvP01 reference genomes during the BWA-MEM alignment step. For in-depth assessment of P. falciparum sequence variation, we processed paired-end Illumina sequencing data in the form of FASTQ files using a custom analysis pipeline (S2 Supporting information) that leverages the Divisive Amplicon Denoising Algorithm (DADA2) tool designed by Callahan et al. 2016 (Callahan et al., 2016) to obtain microhaplotypes (S2 Fig.). We mapped microhaplotypes obtained from DADA2 against a custom-built database of 3D7 and Dd2 reference sequences for each amplicon locus and filtered microhaplotypes based on edit distance, length, and chimeric identification, using a custom R script. (S2 Supporting information). We summarized observed sequence polymorphism into a concise format by converting individual microhaplotypes into pseudo-CIGAR strings using a custom python script. Microhaplotypes were discarded if supported by fewer than 10 read-pairs or by less than 1% total read-pairs within the locus, or if they exhibited other error features (S3 Supporting information).

We analyzed native and pre-amplified mock samples to determine precision and sensitivity of the DADA2 pipeline and filters. We defined a true positive (TP) as a microhaplotype with a pseudo-CIGAR string identical to the reference strain (either 3D7 or Dd2). We defined a false positive (FP) as a microhaplotype with a pseudo-CIGAR string not matching 3D7 (in the case of samples containing only 3D7) or not matching 3D7 or Dd2 (in the case of the mixes), and we defined a false negative (FN) as a locus without any correct microhaplotype representation. We defined precision as TP/(TP+FP), and sensitivity (recall) as TP/(TP+FN). Forty-two of 128 P. falciparum loci in AMPLseq exhibit identical 3D7 and Dd2 reference sequences; we only included these in precision and sensitivity calculations for pure 3D7 controls (i.e., TP+FN = 128); precision and sensitivity calculations for strain mixtures considered only the 86 loci that differ between 3D7 and Dd2 reference sequences (i.e., TP+FN = 86).

All amplicon sequencing data were submitted to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under accession PRJNA758191.

Comparator Panels

We compared AMPLseq and 4CAST to two previously published AmpSeq panels for malaria molecular surveillance, Paragon HeOME v1 (Tessema et al., 2020) and SpotMalaria v2 (Jacob et al., 2021).

Paragon HeOME v1, designed via CleanPlex algorithm (Paragon Genomics Inc, USA), contains 100 primer pairs in a single pool. These target 89 P. falciparum loci selected based on high diversity and differentiation (Jost D ≥ 0.21) among clinical isolates from Africa as well as eleven drug resistance-associated loci. A distinctive feature of HeOME library construction involves its requirement for bead-based clean-up and CleanPlex digestion of each sample between PCR1 and PCR2. The protocol therefore does not require sWGA prior to PCR1.

SpotMalaria v2, designed via Agena BioScience and MPrimer design software, contains 136 primer pairs divided into three different pools. These target loci were considered best able to recapitulate pairwise genetic distance, population differentiation, and sample heterozygosity inferences from global WGS data (MalariaGEN Plasmodium falciparum Community Project, 2016). Primers also target various drug resistance-associated loci (some known to exhibit copy number variation) and mitochondrial loci with conserved primer binding sites among Plasmodium spp. Library construction requires sWGA prior to PCR1 but no special processing between PCR1 and PCR2.

We also compared our amplicon panels to a 24-SNP molecular barcode assay (Daniels et al., 2008). The SNPs targeted by this Taqman qPCR-based assay were chosen principally for their high minor allele frequency (average MAF > 0.35) in parasite sample collections from Thailand and Senegal (Daniels et al., 2008).

paneljudge and in silico data simulations

We used WGS data to simulate genotypic panel data for simulations. This publication uses data from the MalariaGEN Plasmodium falciparum Community Project as described online pending publication and public release of dataset Pf7 (https://www.malariagen.net/resource/34). Specifically, we used genomic data from monoclonal samples collected in Mali, Malawi, Senegal, and Thailand (Zhu et al., 2019), and from Colombia and Venezuela (ENA accession numbers in S4 Table). We also used previously published monoclonal genomic data from Guyana (SRA BioProject PRJNA543530) (Mathieu et al., 2020) and French Guiana (SRA BioProject PRJNA242182) (Pelleau et al., 2015). We used the scikit-allel library (Miles et al., 2020) to process the data and then estimate microhaplotype frequency and diversity. Specifically, we used the read_vcf, is_het, and haploidify_samples functions as described (S1 Supporting information), and we estimated haplotype frequencies with the distinct_frequencies function.

We assessed the performance of different panels for relatedness inference using simulated data. We generated data on pairs of haploid genotypes (equivalent to pairs of monoclonal malaria samples) using paneljudge (Taylor & Jacob, 2020), an R package that we built to simulate data under a hidden Markov model (HMM) (Taylor et al., 2019) (S2 Supporting information). For each panel, we calculated inter-locus distances from the median nucleotide position of each locus and set distances as infinite between chromosomes. For each panel and population of interest, we calculated haplotype frequency estimates using scikit-allel, as described above. Given these distances and frequency estimates, we simulated data using relatedness parameter values of 0.01 (unrelated), 0.50 (siblings), and 0.99 (clonal), and switch rate parameter values of 1, 5, 10, and 50. For each combination of panel, population, relatedness parameter, and switch rate parameter, we simulated data on 100 haploid genotype pairs. For each haploid genotype pair, we then generated estimates of the relatedness parameter and the switch rate parameter using paneljudge, with 95% confidence intervals (CIs), under the same model used to simulate the data. We next performed relationship classification from these estimates and CIs. For estimates of unrelated pairs (relatedness parameter of 0.01), we generously classified estimates as correct if the lower limit of the 95% confidence interval (LCI) was below or equal to 0.01 and the upper limit of the 95% confidence interval (UCI) was below 0.99. We classified estimates of sibling-level relatedness (0.50) as correct if the LCI was above 0.01 and the UCI was below 0.99. We classified estimates of clonal pairs (0.99) as correct if the LCI was above 0.01 and the UCI was above or equal to 0.99. In all relatedness levels, if the 95% confidence interval spanned both 0.01 and 0.99 (i.e., LCI < 0.01 and UCI > 0.99), then we denoted the estimate as unclassified.

To evaluate panel performance in COI estimation, we combined monoclonal WGS data to engineer in silico polyclonal samples using vcftools (Danecek et al., 2011). We then estimated microhaplotype frequencies for each locus of a given panel, using scikit-allel as described above, and counted the number of distinct microhaplotypes observed at each locus per sample. We estimated COI as the maximum number of distinct microhaplotypes observed at any locus within a sample.

To evaluate panel performance for geographic attribution, we identified microhaplotypes at loci as described above. We used the microhaplotype sequences themselves and visualized these data using the Rtsne package (Krijthe, 2015), with 5000 iterations, Θ of 0.0, and perplexity parameters of 10 (for 4CAST and the 24 SNP barcode) or 30 (for the remaining panels).

4CAST and AMPLseq validation

We validated assay precision (defined as TP/(TP+FP)), sensitivity (defined as TP/(TP+FN)), and depth of coverage using 3D7 mock clinical samples representing parasitemia levels between 10 and 10000 parasites/µl in 10 ng/µl human DNA. Both 4CAST and AMPLseq generated 3D7 microhaplotype calls with 100% precision for all parasitemia levels assessed, both with and without pre-amplification by sWGA. 4CAST achieved high sensitivity and depth without preliminary sWGA, generating a median of 43 read-pairs per locus from native templates representing 10 parasites/µl (Fig. 3A). Median depth increased to 443 and 1312.5 read-pairs per locus for native templates representing 100 and 1000 parasites/µl, respectively. Read-pair counts were also evenly distributed among 4CAST loci using native DNA (Fig. 3A).

• Download figure
• Open in new tab

Figure 3. 4CAST and AMPLseq panel validation with mock and clinical samples.

A) Boxplots of read-pairs per locus in the 4CAST panel. The first facet shows read-depth per locus across mock samples ranging from 10 - 10000 p/µl, both native DNA and sWGA DNA (n=4 per condition). The second and third facets show read-depth per locus across two sets of clinical samples, <1 year old and ∼10 year old dried blood spots, respectively (n=32 per sample set). B) The same 3D7 mock sample sets were used to assess AMPLseq sensitivity (top left panel of B) and sequencing depth (bottom left panel of B). Each point in the bottom left panel of B represents read-pair support for one AMPLseq locus. Positions on the y-axis indicate median read-pair support across replicate samples. Low sensitivity observed using native templates (grey) representing ≤100 3D7 parasites/µl suggests that clinical samples should be pre-amplified with sWGA (results at right). C) Ratio of read-pairs from microhaplotypes assigned to 3D7 (x-axis) or Dd2 (y-axis) from mock mixtures of these DNAs in ratios of 1:1 (tan), 3:1 (pink), and 10:1 (dark red). All samples contained 1000 parasites per µl total, across both DNA sources. Dashed lines represent the expected ratio, and each point represents a 4CAST locus per sample (n=4 per condition). D) AMPLseq read-pair ratios observed in native 3D7+Dd2 mock mixtures (1000 parasites/µl) are plotted as above (C) for 4CAST.

Unlike 4CAST, AMPLseq required sWGA for 3D7 mock samples representing ≤ 100 parasites/µl (Fig. 3B). Following sWGA on mock samples representing 10 parasites/µl, the assay generated ≥ 10 read-pairs at a median of 126 loci, with a median of 465 read-pairs after excluding loci with fewer than 10 reads. Values were statistically similar for pre-amplified samples representing 100 parasites/µl and increased to 692 read-pairs for pre-amplified samples representing 1000 parasites/µl (S3 Fig.).

We also validated the sensitivity of 4CAST and AMPLseq for genotyping polyclonal infections by using mock samples containing both 3D7 and Dd2 templates (likewise in 10 ng/µl human DNA). These mixtures featured Dd2 at 50% (i.e., 1:1 3D7:Dd2 ratio), 25% (3:1), and 9% (10:1) relative abundance. Total parasitemia levels ranged between 10 and 10000 parasites/µl. Both 4CAST and AMPLseq generated microhaplotype calls with 100% precision at the 86 loci that are dimorphic between the 3D7 and Dd2 references (including all four 4CAST loci and an additional 82 loci in AMPLseq). This perfect precision was observed at all parasitemia levels in both native and pre-amplified mock mixtures of the two strains.

4CAST showed high sensitivity for Dd2 without the need for sWGA. At 1000 parasites/µl, the assay detected Dd2-specific microhaplotypes at each of its four loci in all 1:1, 3:1, and 10:1 mixture replicates (Fig. 3C). At 100 parasites/µl, median Dd2 sensitivity remained 100% at 1:1 and 3:1 ratios but lowered slightly to 94% at 10:1. Sensitivity for each strain decreased at 10 parasites/µl, with 1:1 ratios yielding a median of 3 target loci for 3D7 and a median of 2 targets for Dd2; median sensitivity in these samples rose to 3.5 and 3 loci (respectively) following pre-amplification with sWGA, but this led to unbalanced read-pair support between the two strains (S4A Fig.), possibly due to differential sWGA success on low-quality Dd2 vs. high-quality 3D7 templates. 4CAST read-pair ratios generated from native templates, by contrast, showed a very high correlation with input ratios at 100 p/µl (S4A Fig.) and 1000 p/µl (Fig. 3C). Ratios became less informative at 10 parasites/µl (S4A Fig.).

AMPLseq (with sWGA) was also successful in detecting Dd2-specific microhaplotypes, but only at a maximum of 77 of 86 dimorphic loci (in the 1:1 ratio at 10000 parasites/µl). Dd2-specific sequences were detected at a minimum of two dimorphic loci for all three input ratios (1:1, 3:1, 10:1) and parasitemia levels (≥ 10 p/µl) assessed. Like with 4CAST, however, the use of sWGA decorrelated read-pair ratios from input ratios (S4B Fig.). Dd2 sensitivity was also reduced relative to 3D7 sensitivity (median Dd2 sensitivity / 3D7 sensitivity = 58%) with sWGA. These discrepancies were not observed with native templates (Fig. 3D) at ≥1000 parasites/µl for which AMPLseq achieves high read-pair support without the use of sWGA.

We also tested both panels on genomic DNA extracted from dried blood spots collected by the Guyana Ministry of Health in 2020 from individuals diagnosed as P. falciparum-positive via a rapid diagnostic test (RDT). Ten Guyanese samples were tested with both panels, and an additional six were tested with AMPLseq. Using 4CAST, we observed coverage across all loci in all samples, with a median read-pair depth per locus of 1162 read-pairs without sWGA (Fig. 3A). Using AMPLseq (with sWGA), we observed a median of 122 loci with ≥10 read-pairs and a median read-pair depth of 298 read-pairs per covered locus (Fig. 3B).

Additionally, we tested both panels on gDNA extracted from 16 dried blood spot samples collected in Mali in 2011 (Tran et al., 2013) and subsequently stored at room temperature for ten years. Using 4CAST (without sWGA), we observed a median read-depth of 407 read-pairs per locus (Fig. 3A). Using AMPLseq (with sWGA), we observed a median of 75 loci with ≥10 read-pairs and a median read-depth of 112 read-pairs per covered locus (Fig. 3B).

Evaluation of panel performance for relatedness

We used the R package paneljudge to assess in silico the impact of choosing a specific genotyping panel for relatedness inference. Considering the choice of panel, we evaluated relatedness estimation from data simulated on our 4CAST and AMPLseq panels, the SpotMalaria v2 (Jacob et al., 2021) and Paragon HeOME v1 (Tessema et al., 2020) amplicon panels, and a barcode of 24 SNPs. When data were simulated using microhaplotype frequency estimates of Senegalese parasites, we found that almost all estimates of unrelated or clonal pairs were correctly classified, regardless of the panel (Fig. 4A). All three large panels also performed similarly well in accurately identifying partially-related parasite pairs, despite being the product of three distinct design processes. Neither 4CAST nor the 24 SNP barcode estimated relatedness for partially-related samples as well as the larger panels. We also evaluated panel performance in less diverse parasite populations (Colombia and Thailand), including a population not used in the panel designs (Colombia). We repeated the simulations using microhaplotype frequencies estimated with these data. Again, we found that all panels performed well for estimating relatedness of clonal pairs, and that the 24 SNP barcode and 4CAST were less likely to have correctly classified estimates of non-clonal pairs. With the data simulated using Colombian microhaplotype frequencies from the Pacific Coast region, all three large panels performed well for all three relatedness values, despite the Colombian data not having informed the design of any of the panels.

• Download figure
• Open in new tab

Figure 4. In silico relatedness estimation comparisons among panels and empirical AMPLseq validation against WGS.

A) Evaluation of relatedness estimation from data simulated on genotyping panels using the paneljudge R package. Pairs of haploid genotypes were simulated at each locus of a panel, using microhaplotype frequencies estimated from a given parasite population (Colombia, Thailand, or Senegal, as shown in columns from left to right). Genotype pairs were simulated at three levels of relatedness: unrelated (relatedness = 0.01), related (relatedness = 0.50), and clonal (relatedness = 0.99), as shown in the rows from top to bottom. Relatedness estimates of these pairs were classified using their 95% confidence intervals (LCI = lower limit of the 95% confidence interval, UCI = upper limit of the 95% confidence interval). Estimates could be correctly classified, misclassified, or unclassified, as described in the grey box. Each bar represents the proportion of simulations per condition (n = 400) classified in each category. Bars that are filled with a color represent correctly classified simulations, bars that are hashed represent misclassified simulations, and bars that are filled with white represent simulations that were unable to be classified. The colors of the bars represent the panel used in that set of simulations. B) Empirical AMPLseq results recapitulate WGS-based relatedness inference. Points represent relatedness estimates (hmmIBD (Schaffner et al., 2018) ‘fract_sites_IBD’ computed under default settings) for pairs of Guyanese samples using WGS (n= 9408 variants) vs. AMPLseq (n = 220 variants, from within 128 AMPLseq P. falciparum loci).

Pairwise relatedness estimates (Schaffner et al., 2018) from AMPLseq correlated highly with those from WGS data available for the Guyanese sample set (Pearson’s r = 0.86, slope = 1.01, p < 0.001) (Fig. 4B). Despite patient travel history metadata suggesting infections to have occurred in various geographic regions of Guyana (S1 Table), AMPLseq relatedness estimates for the Guyanese sample set are significantly higher than those for the Malian sample set (Mann Whitney U, p < 0.001), consistent with anticipated lower transmission levels in Guyana. Nevertheless, the wide range of relatedness estimates (0.007 – 1) observed among Guyanese sample comparisons suggests AMPLseq capacity to indicate epidemiologically relevant microstructure even in relatively unstructured parasite populations. For example, the first (lowest) quartile of pairwise relatedness estimates from Guyana was enriched for comparisons involving A2-GUY and C5-GUY, two highly related samples that in whole genome analysis show 50% relatedness with a sample from Venezuela (SPT26229, see S4 Table).

Geographic attribution

We again engineered amplicon data in silico to evaluate the relative signal in genotyping panels for geographic attribution of samples (Fig. 5). By sub-sampling WGS variant calls, calling microhaplotypes, and visualizing these data using t-SNE plots, we found that results from both the 24 SNP barcode (Daniels et al., 2008) and 4CAST distinguished samples by continent of origin, though not by country. Results from all three larger panels additionally distinguished non-African samples by country, and these panels separated East African (Malawi) from West African samples (Mali/Senegal) to varying degrees; no panel was able to distinguish between Malian and Senegalese samples in this visualization. We also added empirical AMPLseq data from 5 Guyanese samples (C3-GUY, C4-GUY, C5-GUY, C7-GUY, and C8-GUY) and WGS data sub-sampled to AMPLseq coordinates for Venezuelan sample SPT26229 (S5 Fig.). The AMPLseq samples formed a small cluster beside the WGS-based Guyanese and French Guianese samples. The Venezuelan sample SPT26229 also placed on the perimeter of the Guyana/French Guiana sample cluster, sharing the same axis-2 position as the empirical AMPLseq points. Results show that empirical AMPLseq data can distinguish autochthonous samples from the Guiana shield, and we expect geographic attribution in the region to improve as more data are collected from infections originating in Venezuela and other undersampled localities.

• Download figure
• Open in new tab

Figure 5. In silico geographic attribution comparison among panels.

Visualization of WGS data subset to coordinates of genotyping panels. Microhaplotypes called at each locus were visualized using tSNE through the Rtsne package, with parameter Θ of 0.0, 5000 iterations, and a perplexity parameter of 30 (A, C, D) or 10 (B, E). Each dot represents a single sample, colored by its country of origin (which was not included in the tSNE algorithm). One genotyping panel is visualized in each plot: (A) AMPLseq, (B) 4CAST, (C) SpotMalaria v2, (D) Paragon v1, (E) 24 SNP barcode.

COI estimation

We evaluated COI estimation based on 4CAST as opposed to the single locus AMA1, which is commonly used for this purpose, alone or with another locus (Lerch et al., 2017; Miller et al., 2017; Nelson et al., 2019). We engineered in silico samples with COI ranging from 2 – 10 (100 simulations per COI level) and used the maximum number of unique microhaplotypes present at any locus as a simple objective method to estimate COI. 4CAST provided more accurate estimates of COI than AMA1 alone in these simulated data, especially at simulated COI levels between 5 and 7. (Fig. 6A). S6 Fig. indicates that estimation improves at simulated COI=8 using AMPLseq, but to reap the full benefit of the larger panel in practice will require a more elaborate probabilistic algorithm that accounts for variable coverage/sensitivity among loci and incorporates a multinomial approach for polyallelic loci.

• Download figure
• Open in new tab

Figure 6. In silico and empirical complexity of infection (COI) inference.

A) Scatter plots of estimated COI for samples simulated from combinations of monoclonal WGS data, subsetted to the loci of interest (AMA1 locus or 4CAST loci). The x-axis represents the number of monoclonal genomes combined into each simulation, and the y-axis represents the COI estimated using the simulated data. COI was naively estimated as the maximum number of unique microhaplotypes present at any locus per sample (n=100 samples per condition). Each dot represents a sample, jittered for visibility. The black bars represent the median and light grey boxes represent the 25th – 75th quantiles. B) Estimated COI for clinical samples sequenced using 4CAST. COI was again estimated as the maximum number of unique microhaplotypes present at any locus in the sample.

In the absence of such an algorithm, we proceeded with the naive estimation method above to classify COI in the Mali and Guyana clinical samples. Only a single polyclonal infection (C6-GUY) occurred among Guyanese samples assayed by 4CAST. The repeated detection of two CSP and SERA2 alternate alleles at depths ranging from 32 to 168 read-pairs enabled unambiguous COI=2 classification for the sample. WGS sequencing coverage, by contrast, detected only moderately elevated SNP heterozygosity (1.9%) in C6-GUY and this elevation was not sufficient to classify COI>1 via The Real McCoil (Chang et al., 2017) (S7 Fig.). AMPLseq also identified COI=2 for C6-GUY but without consistent support between replicates (2 vs. 6 biallelic loci). Six additional Guyanese samples were assayed by AMPLseq and one was classified as COI=2. This sample (A5-GUY) gave a stronger minor variant signal in both AMPLseq (15 biallelic loci in both replicates) and WGS data (10.9% SNP heterozygosity) (S7 Fig.).

For the Malian sample set, 4CAST and AMPLseq both classified samples E5-PST030 and C6-PST063 as monoclonal and all other samples as polyclonal based on presence/absence of multiallelic loci. While 4CAST detected as many as 10 alternate alleles (median = 3) per sample locus (Fig. 6B), AMPLseq detected at most 4 (median = 2). These results reaffirm 4CAST as a tool of choice for resolving higher COI levels and when parasitemia levels are low. AMPLseq may reach simulated performance levels (S6 Fig.) by reducing sample multiplexing or sequencing on higher output platforms (e.g., NovaSeq) for sample sets with low parasitemia.

Longitudinal sampling: distinguishing recrudescence vs. reinfection

We used 4CAST to examine longitudinal samples that were likely to be diverse and polyclonal. We sequenced samples from the same asymptomatic individual in the longitudinal Mali cohort over three consecutive visits (Fig. 7) (Tran et al., 2013). We detected a single microhaplotype at each locus that was present in the first two time points, suggesting a continued infection during the two weeks between samples. At the third time point, we detected a single, distinct microhaplotype at each locus, suggesting that a new infection had occurred and the original infection had disappeared or decreased below our limit of detection (<10 p/µl). In this particular case, the individual was asymptomatic and did not receive anti-malarial treatment between any samples; however, this simple example demonstrates the clarity that 4CAST can bring to tracking infection turnover in longitudinal studies, and suggests its utility in distinguishing recrudescence vs. reinfection in therapeutic efficacy studies.

• Download figure
• Open in new tab

Figure 7. Longitudinal tracking of infections using 4CAST.

Identification of distinct microhaplotypes present in samples from an individual at three consecutive time points. The x-axis represents the three time points, and the y-axis represents the individual microhaplotypes identified, grouped by locus. Dots represent the presence of that microhaplotype, connected when present in consecutive visits.

Drug resistance profiling

AMPLseq loci in dhfr, mdr1, dhps, kelch13, and mdr2 contain ten sequence regions that code for various amino acid (AA) polymorphisms that have previously been associated with resistance to antimalarial drugs (Ariey et al., 2014; Miotto et al., 2015; Mita et al., 2007; Veiga et al., 2016)Thirteen of these 18 positions of interest contained nonsynonymous mutations in Malian and Guyanese clinical samples of this study (Fig. 8). Positions of interest that lacked mutations across both sample sets were DHFR AA 59 and 164; DHPS AA 613; KELCH13 AA 580; and MDR2 AA 484. All Guyanese sequences shared the same mutant alleles at many loci. Malian samples, by contrast, did not show fixed mutant alleles at any amino acid position of interest. A mix of mutant and wildtype alleles occurred among Malian samples for DHFR AA 51 and 108; MDR1 AA 86, 184, and 1246; DHPS AA 436 and 437; and MDR2 AA 492. A previously reported synonymous polymorphism was observed in one Malian sample at KELCH13 AA 589 (Taylor et al., 2015).

• Download figure
• Open in new tab

Figure 8. Drug resistance-associated sequence profiling in Guyanese and Malian clinical samples.

Bars in the top plot indicate the occurrence of various drug resistance-associated amino acid changes within AMPLseq loci. Positions of interest assayed by AMPLseq but without mutant alleles (see x-marks) in the clinical samples profiled here are labeled in parentheses. Positions 484 and 492 in mdr2 have been suggested to be involved in artemisinin resistance despite lack of experimental data showing an association with a clinical phenotype (Chenet et al., 2017; Miotto et al., 2015). Bottom plot indicates chromosomal and amino acid (AA) positions of each drug resistance-associated AMPLseq locus. Asterisk indicates synonymous mutation within kelch13.

P. falciparum and P. vivax co-infection detection

To test the ability of AMPLseq to detect P. vivax co-infections via co-amplification of PvDHFR, two additional Guyanese blood spot samples that had been diagnosed as P. vivax-only (G4G430) and P. vivax + P. falciparum co-infection (G4G180) via RDT were included in the sample set. These samples did not undergo sWGA.

PvDHFR was detected at high depth in both samples (1068 – 1822 read-pairs for G4G430 and 234 – 560 read-pairs for G4G180) (Fig. 9). Only G4G180 also showed read-pair support at P. falciparum panel loci (>10 read-pairs at 100 – 115 loci). PvDHFR was not detected in any native or pre-amplified 3D7 or mixed-strain (3D7 + Dd2) templates. This demonstrates high specificity of both PvDHFR and P. falciparum AMPLseq primers to their intended target species without any apparent amplification inhibition by the presence of congeneric DNA.

• Download figure
• Open in new tab

Figure 9. Plasmodium vivax detection by AMPLseq.

The left panel demonstrates strong read support for PvDHFR (green circle) in native control samples previously suggested to contain P. vivax (Pv; G4G430) and P. vivax + P. falciparum (Pv + Pf; G4G180) via RDT. The right panel shows native and sWGA results for C7-GUY, a clinical sample that appears to have been misdiagnosed as Pf-only prior to AMPLseq. Blue circles represent read support for P. falciparum loci. Positions on the y-axis indicate median read-pair support across two sample replicates. Box and whiskers indicate quartiles.

PvDHFR was also detected at low levels (16 – 30 read-pairs) in both native template replicates of C7-GUY, one of the sixteen Guyanese samples previously diagnosed as P. falciparum-only via RDT. Surprisingly, two PvDHFR read-pairs were also detected in one of the two sWGA replicates from the sample, despite the expectation that sWGA would primarily amplify P. falciparum sequencines. Sensitivity of PvDHFR detection in pre-amplified samples could be enhanced by adding PvDHFR primers to the P. falciparum sWGA primer pool. PvDHFR detection did not occur in any Malian sample, consistent with low prevalence of P. vivax in West Africa relative to the Guiana shield.