The samples used in this study were previously collected from three geographically distinct locations within the Loreto region of the Peruvian Amazon (Supplementary Table S1). The first location, Nueva Jerusalen (NJ), is an indigenous community located northeast of Loreto, approximately 100 km from Iquitos, accessible primarily by river transport for 2 days. The P. vivax population in NJ showed modest genetic diversity, while P. falciparum had a lower diversity. There were also significant levels of genetic differentiation between parasites from NJ and other remote areas with some temporal clustering observed in the P. falciparum population, linked to an outbreak in February 2020 (Cabrera-Sosa et al., 2024). Samples from NJ were collected through active case detection (ACD) in November 2019 (P. vivax, n = 48; P. falciparum, n = 17) and passive case detection (PCD) from December 2019 to May 2020 (P. vivax, n = 20; P. falciparum, n = 49) (Cabrera-Sosa et al., 2024).

The second location, Mazan (MZ), is a riverine district situated around 50 km southeast of Iquitos. Due to its proximity to the urban center, MZ experiences higher human mobility, which may influence malaria transmission. Referring to the previous studies, the P. vivax population showed little differentiation between MZ and other communities, suggesting a high level of gene flow and connectivity (Manrique et al., 2019). In contrast, the P. falciparum population exhibited great genetic differentiation, indicating that transmission between different areas is limited (Kattenberg et al., 2023a). Samples from MZ were collected by two population-based cross-sectional surveys during 2018 (P. vivax, n = 13; P. falciparum, n = 10) (Villasis et al., 2021; Rosado et al., 2022).

Lastly, SE is a remote riverine community located 30 km south of Mazan, accessible only by boat. Its isolation makes it an important site for studying malaria in remote communities. In this community, the P. vivax population displayed moderate genetic diversity. In contrast, P. falciparum population showed slightly lower diversity with low to moderate levels of population differentiation (Manrique et al., 2015; Cabrera-Sosa et al., 2024). Samples were collected from SE by ACD and PCD in 2016 (P. falciparum, n = 14).

Blood samples (dried filter papers for NJ and SE, or packed red blood cells for MZ) were collected for molecular diagnosis. Also, thick and thin blood smears were taken for light microscopy diagnosis in the health post of each community. Regardless of the symptoms, treatment was provided for any positive microscopy result, following the national guidelines (Ministerio de Salud, 2015).

All sample collections were approved by the Institutional Ethics Committee at Universidad Peruana Cayetano Heredia (UPCH) (SIDISI codes: 64024, 101518 and 102725). Written informed consent was obtained from all participants in MZ and SE. Samples from NJ were collected as part of the Ministry of Health (MINSA) activities and then transferred to our team for research purposes. This study was also approved by Institutional Ethics Committee at UPCH (SIDISI 207543). All methods were performed following the MINSA guidelines and regulations.

DNA was extracted using EZNA^® Blood DNA (Omega Bio-tek, United States), following the manufacturer’s protocol with an elution volume of 50 µL. Molecular diagnosis was performed by real-time PCR (qPCR). Samples from MZ were diagnosed using a SYBR Green-based protocol (Mangold et al., 2005). A TaqMan-based assay (Rougemont et al., 2004) was used for samples from NJ and SE.

Eighty-one P. vivax and 86 P. falciparum samples with parasitemia greater than 5 parasites/μL were genotyped with the Pv and Pf AmpliSeq Peru assays (Supplementary Figure S1). Library preparation was performed as previously described (Kattenberg et al., 2023a; Kattenberg et al., 2023b; Kattenberg et al., 2024). In brief, 7.5 µL of DNA was amplified with two primer pools. Then, digestion, indexes ligation, library amplification and cleaning steps were performed. Libraries were quantified using Qubit High sensitivity DNA kit (Invitrogen), pooled, diluted to 2 nM and denatured with NaOH to 7pM. Finally, the denatured library pool was loaded on a MiSeq system (Illumina) for 2 × 300 paired-end sequencing (Miseq Reagent Kit v3, Illumina) with 5% PhiX spike-in (Illumina).

FASTQ files were processed using an analysis algorithm as previously described (Kattenberg et al., 2023b). After trimming, reads were aligned to the reference genome (PvP01 v46 or Pf3D7 v44 from PlasmoDB, https://plasmodb.org/plasmo/app). Then, variants were called, generating individual gVCF files per each sample, which were combined to call genotypes jointly. After, variants were hard filtered and annotated. Finally, Variant Call Format (VCFs) containing only the Peru SNP barcode in each AmpliSeq assay (28 for P. falciparum, 41 for P. vivax) were created (Kattenberg et al., 2023a; Kattenberg et al., 2024).

The diploid genotype table (GT) was extracted from the SNP barcode VCF and manually reformatted into a haploid GenAlEx file (Peakall and Smouse, 2012). We constructed multi-locus genotypes (MLG) by combining the alleles across all loci from the barcode and assumed 2 MLG per sample. Then, when a sample had at least one heterozygous position, two MLG were created for GenAlEx file. The reference allele of each heterozygous position was kept in the first “dominant” MLG, and the alternative allele in the other “secondary” MLG, maintaining the homozygous positions in both.

In total, 62 P. vivax and 90 P. falciparum samples were genotyped by using Plasmodium species-specific MS panels (Supplementary Figure S1). For P. vivax, a 16 MS panel (11.162, Ch2.121, 14.297, Ch2.152, Ch14.3021, Ch14.2986, Ch14.3010, Ch2.122, MS9, 13.239, Ch14.2981, MS6, MS4, MS15, 3.502, MS20), previously developed by our team, was used (Manrique et al., 2019). Samples were classified according to their parasitemia and followed a different processing before MS genotyping. Samples with >230 parasites/μL (par/μL) were diluted to 75–100 par/μL for the MS amplification. On the other hand, samples between 40 and 230 par/μL were amplified by selective whole genome amplification (sWGA) using the enzyme phi29 DNA Polymerase (NEB) before the genotyping. Finally, DNA from samples between 15 and 40 par/μL was reextracted, mixed with the first extraction and used for posterior PCRs (Manrique et al., 2019). Each MS was amplified by conventional PCR. In each reaction, one of the primers was labeled with a fluorophore for later identification. The master mix consisted of 1U of the AccuStart II Taq DNA polymerase enzyme (QuantaBio), 2.5 μL of PCR buffer, 1.5 mM MgCl₂, 200–400 nM of each primer (depending on the panel) and 2.5 μL of DNA. The amplification protocol consisted of an initial denaturation at 95°C for 5 min, followed by 40 cycles of 30s at 94°C, 30s at 55°C–67°C (Tm depending on the MS) and 30s at 72°C. The last step was a final extension at 70°C for 15 min (Manrique et al., 2019).

In the case of P. falciparum, a panel of seven MS (TA1, Poly-α, PFPK2, TA109, 2490, C2M34, C3M69), previously applied to Peruvian samples (Akinyi et al., 2013; Bendezu et al., 2022) was used. Three amplification protocols were used, including two direct PCR (for TA109, Poly-α, PFPK2, C2M34 and C3M69) and one semi-nested PCR (for TA1 and 2490) (Nair et al., 2003; Roper et al., 2003). The master mix consisted of 1X Promega Master-mix (Promega), 400 nM of each primer and 2 µL of DNA or 1 µL of the product of the first PCR. For all cases, extracted DNA was used directly, without dilution or sWGA.

In both P. vivax and P. falciparum, positive controls and no-template controls were included in the reactions. After MS amplification, the PCR products were diluted and mixed with formamide and the GeneScan 500 LIZ marker (Applied Biosystems). Finally, the mixes were analyzed by capillary electrophoresis on an ABI Prism 3130 (Applied Biosystems). Chromatograms were analyzed with the Microsatellite Analysis software (Applied Biosystems). A quality control of the size marker was performed and the list of peaks per sample was downloaded.

With this data, a semi-automated analysis was performed to select true alleles. First, the peak with the highest relative fluorescence units (RFU) value (height) in the no-template controls was used as a threshold and all peaks with lower height in samples were discarded. After, the peak with the highest height was considered the main allele for each MS and sample. Secondary alleles were considered if their peak height was at least 1/3 of the main allele, and their size was not within ± 2bp of the main allele. As a final validation step, we manually compared the selected alleles in each sample to the positive control, discarding any allele when a discrepancy in the shape or pattern of peaks was detected. If a previously selected main allele was discarded in this process, the remained allele with the highest height was reclassified as the new main allele.

The final allele table was reformatted manually into a haploid GenAlEx file (Peakall and Smouse, 2012). For subset analysis, we separated the data of samples with secondary alleles into 2 MLG in the GenAlEx formatted file, with the main alleles kept in the first “dominant” MLG, and the secondary alleles in the other MLG, maintaining the rest of MS alleles in both.

We applied some inclusion criteria for the population genetics analyses (Supplementary Figure S1). First, to ensure comparability, we selected samples with genotyping data from both AmpliSeq and MS (P. vivax, n = 62; P. falciparum, n = 86). Then, samples with missing data in more than 25% of panels (SNP or MS) were excluded. 80/86 P. falciparum and 51/62 P. vivax samples passed this filter. Next, we kept markers with less than 25% missing data. Using this criteria, 2/28 positions from the P. falciparum SNP barcorde and 1/16 MS from the P. vivax panel were excluded. At the end, the final genotyped data included 80 P. falciparum samples with 7 MS (105 MLG) and 26 SNP (117 MLG), and 51 P. vivax samples with 15 MS (86 MLG) and 40 SNP (68 MLG) (Supplementary File S1).

We assessed the inclusion of the secondary MLG in the analysis by comparing the genetic differentiation of MLG sets containing: Set 1) samples with only one MLG, Set 2) dominant MLG of all samples, Set 3) dominant MLG of all samples and secondary MLG of samples with in only one secondary allele, and Set 4) all dominant and secondary MLG of all samples. Fixation index (Fst) was calculated using hierfstat package in R (Goudet et al., 2022). Genetic differentiation (Fst) between all these sets were close to 0 (Supplementary Figure S2), indicating no genetic differences between them. Therefore, all MLG for all samples were included.

We calculated population genetic parameters using each SNP and MS panel for P. vivax and P. falciparum samples. Genetic diversity was measured as expected heterozygosity (He) and calculated using the formula: He = [n/(n − 1)][1 − Σp²], where n is the number of genotyped samples and p is the frequency of each allele at a given locus (Nei and Roychoudhury, 1974; Serrote et al., 2020). Genetic differentiation was expressed as pairwise Fst values and 95% confidence interval (95% CI), calculated based on Weir and Cockerham’s Method using the diveRsity package in R (Weir and Cockerham, 1984; Keenan et al., 2013). Population structure was explored using principal component analysis (PCA) (Jolliffe and Cadima, 2016) and the software STRUCTURE v2.3 (Pritchard et al., 2000). PCA was performed using the “prcomp” function in stats R-package. STRUCTURE determined the most likely number of clusters (K). Runs were performed exploring K from 1 to 10 (20 iterations each), assuming a mixture model, correlated allele frequencies, and using sampling location as priors (LOCPRIOR option). Each run had a burn-in period of 50,000 iterations followed by 150,000 Markov Chain Monte Carlo (MCMC) iterations. The most likely K was defined by calculating the rate of change of K (ΔK) (Evanno et al., 2005), using the R package pophelper (Francis, 2017). The complexity of infection (COI) was used to assess the proportion of polyclonal infections (multiple clones in one sample). Samples with 2 MLG (i.e., at least one heterozygous SNP genotype for the AmpliSeq data or at least one secondary allele in any MS) were considered polyclonal. For comparison, we also determined COI with the within-sample F (Fws) statistic using all biallelic SNPs detected by the AmpliSeq assays (P. vivax: 2923 SNPs, P. falciparum: 963 SNPs) using the moimix package in R. Samples were considered a monoclonal infection when the Fws was ≥ 0.95 (Auburn et al., 2012).

He values across different locations were compared performing the Kruskal-Wallis rank sum test. Dunn’s post hoc tests were used to perform pairwise comparisons between all possible pairs of groups to determine which specific groups are significantly different from each other. Comparison of He results obtained with MS and SNP were calculated using the correlation coefficient using Spearman’s rank correlation. Pairwise Fst from SNP and MS data were compared using Wilcoxon Signed-Rank Test. Proportions of mono/polyclonal infections according to SNP barcode, MS and Fws were compared for each species (P. falciparum and P. vivax) using the chi-squared test. Cohen’s Kappa coefficient was used to assess the agreement between SNP and MS in determining the complexity of infection. All statistical analyses were performed in R Studio (version 2022.12.0) using R version 4.2.2. P-values < 0.05 were considered significant. The Bonferroni correction was applied to adjust the p-values when performing multiple pairwise comparisons.

Costs associated with AmpliSeq and MS genotyping for 96 samples (and controls) were estimated, including materials (reagents and kits) and time-personnel costs. All material prices were calculated based on quotations and invoices from Lima, Peru, from 2022. The cost does not consider DNA extraction, molecular diagnosis and plastic labware because those methods are similar for both procedures. For P. vivax MS genotyping, we assumed that 60% of samples required sWGA (Manrique et al., 2019). Total costs included portion of materials used exactly for 96 samples. Separately, we calculated the cost of sequencing (MiSeq) and capillary electrophoresis (ABI Prism 3130), including equipment and annual maintenance.

Personnel cost per day (60 USD/working day) was calculated considering monthly salary of 1,200 USD and 20 working days in one month. Number of working days in each activity was based on practical experience. Costs in Peruvian Nuevos Soles (PEN) were converted into United States Dollars (USD) using the exchange rate of 1 PEN = 0.26 USD.

For P. vivax, only the reference allele was detected for the SNP PvP01_13_v1_32509. The median minor allele frequency (MAF) for the rest of 40 SNPs was 0.27 [range: 0.02–0.46, mean ± standard deviation (sd): 0.25 ± 0.14]. All P. vivax samples had ≤25% missing genotypes in the barcode. In the case of P. falciparum, all alternative alleles at 28 SNPs were genotyped successfully. The median MAF was 0.19 (range: 0–0.48, mean ± sd: 0.20 ± 0.15). Two SNPs (Pf3D7_02_v3_694307 and Pf3D7_09_v3_231065) had >25% of missing data (57% and 66%, respectively) and were excluded for downstream analyses. A low proportion of samples (4.65%, 4/86) had >25% of missingness in the SNP barcode and were subsequently excluded.

The amplification efficiency of both MS panels is shown in Supplementary Figure S3. For P. vivax, MS had low proportion of missing data (3.2%–22.6%), except for Ch14.3010 (43.6%), which was excluded for later analyses. 52 (83.4%) samples had ≤25% missing genotypes. On the other hand, amplification efficiency in the P. falciparum MS panel was ≥90% (% missing data: 0%–10%). Only 5 samples (5.6%) with >25% of missingness were excluded.

We found similar trends of population diversity (expressed as He) with both panels (MS or SNP) for both species (Figure 1; Supplementary Figure S4). P. vivax populations in MZ and NJ had similar genetic diversity when using MS (Figure 1A, median He: 0.68–0.73, p = 0.23) or SNP (Figure 1B, median He: 0.37–0.38, p = 0.80). We also noted that 2/40 (5%) of the SNP positions were fixed in all population (Supplementary Table S2).

On the other hand, P. falciparum isolates from MZ (median He: 0 with both markers) had lower diversity than parasites from NJ using MS (Figure 1C, median He: 0.36–0.37, p = 0.03–0.04) or SNP (Figure 1D, median He: 0.04–0.09, p = 0–0.03). Also, He values in P. falciparum samples from MZ and NJ were similar to SE isolates, determined by MS (median He: 0.28, p = 0.78–1) and by SNP (median He: 0, p = 0.36–1). Finally, around 40% of SNP positions (10/26) were fixed in these populations (Supplementary Table S2).

Despite the similar trends in the genetic diversity, we found no significant correlation between He values obtained by MS and SNP in the P. vivax (Figure 1E, R = 0.5, p = 1) and P. falciparum (Figure 1F, R = 0.74, p = 0.26) populations, although the correlation lacked statistical power as the degrees of freedom were low (1 and 2, respectively).

We observed no significant differences between the pairwise Fst values obtained using MS and SNP across the P. vivax (p = 0.34–0.85) and P. falciparum (p = 0.24–0.83) populations (Figure 2). In case of P. vivax, low differentiation of parasites within NJ was noted when using MS (Fst = 0.04, 95% CI: 0.01–0.07) or SNP (Fst = 0.03, 95% CI: 0–0.09) (Figure 2A). Moderate levels of differentiation were observed between MZ and NJ with MS (Fst = 0.13–0.14, 95% CI: 0.08–0.22) and SNP (Fst = 0.09–0.12, 95% CI: 0.02–0.22) (Figure 2A). In the P. falciparum population, the lowest Fst values were found between the two groups in NJ when using MS (Fst = 0.14, 95% CI: 0.07–0.24) and SNP (Fst = 0.19, 95% CI: 0.09–0.32). On the other hand, the highest Fst values was observed between MZ and SE isolates (Fst = 0.65, 95% CI: 0.41–0.84 with MS, Fst = 0.61, 95% CI: 0.31–0.83 with SNP) (Figure 2B). The consistency between the two markers suggests that both MS and SNP markers are capable of reliably capturing genetic differentiation in the Peruvian setting.

Clustering patterns obtained from STRUCTURE were strikingly similar between MS and SNP in both species (Figure 3), although the most likely number of clusters (K) differed between panels. For P. vivax, the predicted best K were 2, 4, and 6 according to the SNP barcode, whereas the best K were 3 and 6 using MS (Supplementary Figure S5). For P. falciparum, the predicted best K using SNP were 3, 5 and 7, meanwhile the most likely K with MS were 2, 4, and 7 (Supplementary Figure S6). For comparison between MS and SNP, we plotted the assignation to 3, 4, and 6 clusters for P. vivax isolates (Figure 3A; Supplementary Figure S5) and to 4, 5 and 7 clusters for P. falciparum parasites (Figure 3B; Supplementary Figure S6).

For P. vivax populations, high admixture was observed in NJ, with all clusters present in both years. In contrast, all samples from MZ belonged to predominantly one cluster according to SNPs, while MS did not attribute to any single cluster but showed the same admixture pattern (Figure 3A; Supplementary Figure S5).

In the case of P. falciparum, strong clustering observed in all population, except in NJ 2019, which presented high admixture. Predominant clusters in NJ 2020 and MZ were different, and SE shared some clusters with NJ 2019 and MZ. SNPs also detected some admixed samples in NJ 2020, which were predicted to belong to a single cluster by MS (Figure 3B; Supplementary Figure S5).

Previously, we identified three genetic clusters (Pop 1 to Pop 3) in the P. falciparum population in NJ using all biallelic SNPs detected in the AmpliSeq assay based on results from PCA and supported by other analyzes such discriminant analysis of principal components (DAPC), identity by descent (IBD) network and phylogenetics (Cabrera-Sosa et al., 2024). Using STRUCTURE analysis, all those genetic clusters (Pop 1 to Pop 3) had distinct patterns with only the 28-SNP barcode or MS (Figure 3C; Supplementary Figure S6). In Pop 1, similar admixture clustering was frequent, except when K = 4 with SNPs were plotted. Different clusters were predominant in each Pop 2 and Pop 3 (Figure 3C; Supplementary Figure S6).

PCA analysis confirmed patterns observed in the STRUCTURE analysis, and overall similar clustering in both MS and SNP (Figure 4). For P. vivax, no clear geographical or temporal clustering was observed in either panel (Figures 4A, B). On the other hand, a stronger clustering of P. falciparum isolates of NJ 2020 was observed in SNP than with MS (Figures 4C, D). The variance explained by the first two principal components were similar for P. vivax (20% and 25% for MS and SNP, respectively) and P. falciparum (63% for MS, 62% for SNP).

We found differences in the proportion of polyclonal infections determined by the SNP barcode, MS and the Fws statistics (Table 1). For P. vivax, this proportion was higher with MS (69%) than with SNP and Fws (33% and 28%, respectively; p = 3.3 × 10⁻⁵). Five MS (MS9, MS15, Ch2.121, Ch2.122, Ch2.152) accounted for high proportion of samples with secondary alleles (26%–43%). For P. falciparum, the SNP barcode detected similar polyclonal infections as MS (46.2% and 31.2%, respectively, p = 0.21), but higher than Fws (29%, p = 0.043). Agreement between SNP barcode and MS was higher in P. falciparum (κ = 0.54, 95% CI: 0.36–0.72) than in P. vivax (κ = 0.23, 95% CI: 0.04–0.43).

We estimated higher costs associated with AmpliSeq compared to MS genotyping (Table 2). The estimated cost of the AmpliSeq was 183 USD/per sample, whereas using the MS panels was 49 USD/sample for P. vivax (16 MS) and 27 USD/sample for P. falciparum (7 MS). Notably, a large difference was observed for the personnel costs, with the AmpliSeq been 2 to 4 times cheaper than MS genotyping (5 vs. 11–19 USD/sample), as it was much faster (AmpliSeq: 6 working days, MS: 18 to 30 working days) (Supplementary File S2). Total cost for equipment and annual maintenance was higher with the MiSeq than ABI 3130 (190k vs. 160k USD), although when referring as external service, cost for AmpliSeq and P. vivax MS genotyping was similar (33 vs. 30 USD/sample) (Supplementary File S2).