Population genetics of Plasmodium vivax with transmission and rebound in two endemic areas of Papua New Guinea

Background: Global efforts to control and eventually eliminate malaria have been less effective for Plasmodium vivax relative to Plasmodium falciparum due to its unique biology, including dormant liver stages that cause later relapse, and earlier commitment to transmission stages. After the nationwide distribution of long-lasting insecticide treated nets (LLIN) in Papua New Guinea (PNG), P. vivax initially reduced to low prevalence, but again resurged to levels similar to those before LLIN distributions.

Method: To explore changes in P. vivax population structure and identify sources of resurgence over this period, we applied a previously validated genome-wide SNP barcode to genotype 336 P. vivax isolates obtained from serial cross-sectional surveys conducted over a decade in East Sepik (2005, 2012, 2016) and Madang Province (2006, 2010, 2014).

Results: Population genetic analyses of the resulting parasite genotypes revealed contrasting spatiotemporal patterns between the two provinces. In Madang, the complexity of infection, genetic diversity, and population structure varied with prevalence, with a possible population bottleneck and early clonal expansion at low transmission, and rapid recovery of the population with resurgence. In East Sepik, there was a less dramatic impact on the parasite population after prevalence decline, and ongoing transmission of multiple residual lineages throughout the study period. P. vivax decline was also accompanied by an increase in genetic differentiation between the two areas, which reduced with resurgence suggesting changes in parasite migration between areas associated with prevalence.

Conclusion: The earlier implementation of LLIN in East Sepik, smaller rebound, heterogeneity in transmission and relative isolation, compared to Madang may have contributed to these differing patterns. The results demonstrate that long term sustained control efforts are essential to make a lasting impact on the P. vivax population, and that SNP barcodes can provide valuable insights into parasite transmission dynamics as a result of control efforts.

Background

In the past two decades, there has been a dramatic reduction in malaria morbidity and mortality globally using currently available tools such as long-lasting insecticide treated bed nets (LLIN), rapid diagnostic testing (RDTs) and highly effective artemisinin combination therapies (ACTs) (WHO, 2024). This achievement has led many endemic countries to work towards malaria elimination by 2030 (Cibulskis et al., 2016; Rabinovich et al., 2017). However, the emergence of artemisinin resistance (Ippolito et al., 2021), high transmission heterogeneity at local scales (Omedo et al., 2017; Kattenberg et al., 2020a; Gul et al., 2021; Fadilah et al., 2022) and importation of cases from higher transmission regions through human mobility (Chang et al., 2019; Buyon et al., 2020; Neafsey et al., 2021) threatens these elimination plans. Furthermore, in areas where both major human malaria parasites, Plasmodium vivax and Plasmodium falciparum co-circulate, P. vivax is more resilient to control efforts, owing to its ability to relapse, cause lower density infections, and transmit prior to clinical symptoms (Adams and Mueller, 2017; Ferreira and de Oliveira, 2015).

For malaria endemic countries to drive further reductions in transmission and progress towards elimination, new and more efficient surveillance approaches are needed. As infections become rarer and more clustered in space and time, traditional estimates of malaria transmission such as the entomological inoculation rate and infection prevalence, become more difficult and expensive to measure (Tusting et al., 2014). Genomic surveillance applies novel and sensitive molecular tools to parasite isolates collected through community surveys or clinical cases to address a number of malaria control “use cases” (Dalmat et al., 2019), including monitoring changing transmission dynamics (Auburn and Barry, 2017; Neafsey et al., 2021), to track the emergence and spread of new strains such as drug and diagnostic resistant strains (Miotto et al., 2020; Fola et al., 2023), and to monitor population connectivity and the risk of imported infections (Taylor et al., 2017). Population genomic surveillance can reveal whether the transmission is endemic (heterogeneous), epidemic (clonal expansion) or a chain (no recombination) (Colijn and Gardy, 2014; Fola et al., 2018). These genomic analyses provide insights into the transmission dynamics underlying prevalence estimates, and help to determine whether control programs have interrupted endemic transmission, to identify the factors sustaining ongoing or resurgent transmission, and whether existing efforts are adequate to eliminate malaria (Volkman et al., 2012; Barry et al., 2015).

Papua New Guinea (PNG) has the highest P. vivax transmission in the world (Cleary et al., 2022). Intensified nationwide malaria control activities since 2004, with nationwide implementation of LLIN by 2008, significantly reduced malaria transmission in PNG between 2008 and 2010 (Hetzel et al., 2015). However, since 2014, malaria has rebounded substantially in the north coast (Koepfli et al., 2017; Kattenberg et al., 2020a) and throughout the country (Hetzel et al., 2017). Earlier studies on P. vivax population genetics before malaria control was intensified i.e. pre-LLIN (2005-6) revealed high genetic diversity and a lack of population structure (Koepfli et al., 2013; Jennison et al., 2015). We also investigated the spatiotemporal variation of these P. vivax populations from 2005 to 2016 using ten highly polymorphic microsatellite markers, finding that despite a significant decline in parasite prevalence, there was no change in population structure and parasite populations maintained high genetic diversity (Kattenberg et al., 2020b). Nevertheless, some signals of transmission disruption were observed, such as an increase in genetic differentiation (F_ST) between provinces, and the emergence of significant multilocus linkage disequilibrium (LD) and focal clusters of phylogenetically related parasites (Kattenberg et al., 2020b).

Microsatellites are fast evolving with high allelic diversity per locus (Li et al., 2002) and previous studies typed less than one marker per chromosome and therefore will have a limited ability to accurately measure relatedness and detect more distant relationships amongst parasite isolates. We have developed an informative single nucleotide polymorphism (SNP) barcode with 176 SNPs spanning all P. vivax chromosomes, and demonstrated its improved performance and higher resolution to detect variability in diversity and geographic population structure compared to microsatellites using the 2005-6 baseline samples, collected pre-LLIN (Fola et al., 2020). We hypothesized that the small number of microsatellite markers, their high polymorphism and possibility of technical artifacts (Sutton, 2013; Ruybal-Pesántez et al., 2024), may limit the accurate detection of parasite population structure and Identity by Descent (IBD) as a measure of parasite relatedness. Therefore, we aimed to assess P. vivax population genetic signals associated with changing transmission using the SNP barcode (Fola et al., 2020). Using a subset of these SNPs, we measured complexity of infection and proportion of polyclonal infections, population level measures of genetic diversity and structure, and parasite relatedness. The results provide unique and unexpected insights into the spatio-temporal dynamics of P. vivax populations associated with declining transmission and rebound, highlighting the importance of marker panel choice and population genetics in characterising transmission dynamics, to monitor P. vivax populations.

Results

Sample genotyping

One hundred (100) SNPs from a previously validated 176 SNP barcode (Fola et al., 2020), Supplementary Table S1) were used to genotype P. vivax isolates from East Sepik (ESP) and Madang Provinces of PNG (Figure 1). Samples were derived from cross-sectional surveys conducted between 2005 and 2016 in several villages of ESP surrounding Maprik town, and from the Mugil, Malala and Utu catchments of Madang Province (Mueller et al., 2009a; Senn et al., 2012; Koepfli et al., 2017; Kattenberg et al., 2020a). From a total of 376 P. vivax qPCR positive samples, 336 (89.3%) were successfully genotyped with at least 70% of the target SNPs called per genotype (see Materials and Methods, Table 1). Among the genotypes, 12 out of the 100 SNPs were excluded because at least 30% of samples failed. The remaining 88 SNPs were used to construct infection haplotypes for population genetic analysis (Supplementary Table 1). More samples were available from Madang due to the larger geographic area covered and three defined catchment areas.

Figure 1

Figure 1. Map of the study area and P. vivax prevalence and polyclonal infections trends from 2005 to 2016. The coloured dots in the map indicate villages included in serial cross-sectional surveys in each Province. Blue, East Sepik villages; Pink, Madang villages. The inset graphs associated with each Province shows P. vivax malaria prevalence trend measured by qPCR and proportion of isolates that were polyclonal i.e. COI greater than one according to the Real McCoil analysis of SNP genotypes (polyclonality %, Table 1).

Table 1

Table 1. Summary of genotyped P. vivax samples and complexity of infection.

Complexity of infection and genetic diversity

All samples were screened prior to SNP genotyping for multiplicity of infection using capillary electrophoresis of highly polymorphic microsatellites msp1f3 and ms16 (Koepfli et al., 2011; Arnott et al., 2013). Only samples with one or two clones were selected, however 238 of the 336 samples (70.83%) were found to have heterozygous SNPs, including 121 of the originally classified single infections, showing that SNP barcoding with deep amplicon sequencing is more sensitive to detect multiple clones. The remaining 98 samples were confirmed as single infections. Complexity of Infection (COI) was then measured using the SNP data (Table 1).

There was significant temporal variation in mean COI and the proportion of polyclonal infections (polyclonality) using SNPs in both study areas (ESP: p < 0.01; Madang: p < 0.001, Mann-Whitney U test). However, in Madang, there was a large decline in mean COI and polyclonal infections from 2006 to 2010 then an increase in 2014, whereas in ESP these values where more stable showing a subtle decline from 2005 to 2012 and continuing to decrease even despite the increase in prevalence in 2016 Table 1; Figure 1).

As a measure of within-infection genetic diversity, the mean number of heterozygous SNP calls per population was calculated. There was a significant change in this metric with declining transmission in ESP (2005: 21.27, 2012: 16.70, p < 0.01, Mann-Whitney U test) but a much greater decrease for Madang (2006: 22.20, 2010: 8.60), p < 0.001, Mann-Whitney U test) (Figure 2A). In ESP, heterozygous SNPs continued to decline throughout the study period even after rebound in 2016 (14.23). In Madang, the decline of heterozygous SNPs in 2010, was followed by a substantial increase with rebound in 2014 (16.20, Figure 2A). Spatial variation in these values, i.e. comparing ESP to Madang, was only observed at the lowest prevalence time points (ESP 2012, Madang 2010) due to the high proportions of heterozygous SNPs maintained in ESP.

Figure 2

Figure 2. Spatiotemporal changes in genetic diversity of P. vivax populations of Papua New Guinea (A) Mean number of heterozygous SNP loci, (B) Nucleotide diversity (SNP π) (C) Effective population size (N_e). Error bars represent 95% CIs. (D) Multilocus Linkage disequilibrium (LD) as measured by the Standardized index of association (I_A^S).

Minor allele frequency (MAF) distributions demonstrate that at baseline the 88 barcode SNPs represented a wide range of allele frequencies (Supplementary Figure S1). With transmission decline, some SNPs shifted from the low allele frequencies to the mid-high frequencies which is consistent with a population bottleneck. This pattern was more pronounced in Madang 2010 and remained through to the third time point, whereas for ESP, the changes were more subtle and SNPs with low allele frequencies increased in the third time point, consistent with population expansion.

The nucleotide diversity statistic, (SNP π), measures genetic diversity among infections. SNP π was high in both Madang (0.45) and ESP (0.41) at baseline (Figure 2B). There were minimal changes in this parameter in ESP as transmission declined from 2005 (0.41) to 2012 (0.41) and upon rebound in 2016 (0.38, Figure 2B). Whereas, in Madang there was a significant reduction in SNP π from 2006 (0.45) to 2010 (0.38, p = 0.024, Mann-Whitney U test) and an increase with transmission rebound in 2014 (0.43, p = 0.062). The change in effective population size (N_e) over time was also greater in Madang than ESP (ESP 2005: 4,410, 2012: 3,886, 2016: 4,022; Madang 2006: 4,262, 2010: 2,289, 2014: 3,980, Figure 2C).

Low levels of multilocus linkage disequilibrium (LD) in both provinces at baseline are consistent with high levels of transmission and outcrossing (Anderson et al., 2000). LD increased with declining transmission, with a significant change in Madang from 2006 to 2010 (Figure 2D, p < 0.001). Whereas, there was a significant decrease after transmission rebound in Madang 2014 relative to the 2010 timepoint (0.005, p-value <0.01).

Spatiotemporal changes in P. vivax population structure

Low levels of genetic differentiation (F_ST = 0.075, Supplementary Figure S2) and overlapping genotypes in the Principle Components Analysis (PCA) were observed between provinces at baseline (Supplementary Figure S3A) and are consistent with microsatellite data (Koepfli et al., 2013; Jennison et al., 2015). Comparing the low transmission time points ESP 2012 and Madang 2010, there was a very high degree of spatial genetic differentiation (F_{ST =} 0.252, Supplementary Figure S2) and clustering of genotypes (Supplementary Figure S3B). After transmission rebound, the genetic differentiation between the two provinces decreased (F_ST = 0.182, Supplementary Figure S2) and genotypes showed more overlap in the PCA (Supplementary Figure S3C) indicating an increase in gene flow between provinces. Within ESP there was little genetic differentiation between samples collected in consecutive surveys (2005/2012: F_ST = 0.07, 2012/2016: F_ST = 0.08) but high genetic differentiation comparing first and last time points (2005/2016: F_ST = 0.15, Supplementary Figure S2). However, there was high genetic differentiation among Madang populations between all paired time points (2006/2010: F_ST = 0.16, 2010/2014 F_ST = 0.17, 2006/2014 F_ST = 0.18, Supplementary Figure S2).

Phylogenetic analysis confirmed the spatiotemporal population structure and genotype clustering in Madang, with increasing structure in the tree for post-LLIN populations both within and between time points (2010,2014), whereas there was no substantial structure amongst pre-LLIN (2006) Madang genotypes and all ESP genotypes (Figure 3A; Supplementary Figure S4). Small clades of pre-LLIN Madang genotypes (2006) also cluster with ESP genotypes (2005) confirming significant connectivity between provinces at that time. Whereas there was limited clustering of genotypes from the two provinces at the low transmission mid-point (2010/12) consistent with a reduction in gene flow and higher LD. At low transmission, Madang 2010 genotypes comprised two major clades (Figure 3A; Supplementary Figure S4A). The less frequent Madang 2010 clade containing only genotypes from one catchment area (Utu), clustered with the majority of haplotypes from Madang 2014, indicating this lineage was a major source of the Madang resurgence.

Figure 3

Figure 3. Spatiotemporal changes in population structure of P. vivax in Papua New Guinea (A) Phylogenetic tree, An unrooted neighbour-joining tree illustrates the genetic relatedness between 336 genotypes of 88 SNPs based on pairwise distance matrices, using the number of differences method, with branch lengths scaled accordingly (scale indicated). Colours indicate geographic area and year, while shapes for Madang, indicate the three catchment areas: circle = Malala, square = Mugil, triangle = Utu. (B) Bayesian cluster analysis. Individual ancestry coefficients are shown for K = 3, 4 and 6, each vertical bar represents an individual haplotype and the membership coefficient (Q) within each of the genetic populations, as defined by the different colours.

Further analysis using STRUCTURE software (Kaeuffer et al., 2007) at K = 3 showed one cluster predominant in ESP 2005 and Madang 2006 and 2014. The other two clusters were predominant in ESP 2012 and 2016, and Madang 2010 populations (Figure 3B). At K = 4, Madang genotypes formed three distinct sub-populations according to the year of collection (2006, 2010 and 2014). Whereas, ESP isolates form two genetic clusters with high proportions of ancestry shared between baseline (2005) and post-control (2012, 2016) populations. Approximately half of the ESP 2005 samples share ancestry with Madang 2006 genotypes. Additional sub-structure was detected at K = 6 for ESP at the village level with genotypes from Sengo clustering in 2005, Sunuhu in 2012 (which was also observed using microsatellites, Kattenberg et al., 2020b), and for Madang in Malala-Amiten village. At K = 4 and 6, a minor Madang 2010 cluster (dark green) accounts for the majority of the Madang 2014 population, consistent with the phylogenetic analysis. The phylogenetic tree also shows genotypes from 2006 in this clade suggesting it may have originated from a persistent lineage circulating throughout the study period (Figure 3A).

Clonal expansion in Madang at low transmission and ongoing endemic transmission in East Sepik

To further explore the spatiotemporal distribution of parasite lineages, we measured the expected Identity by Descent (eIBD) parameter and investigated relatedness networks using the “isoRelate” R package (Henden et al., 2018; Taylor et al., 2019). We first evaluated the proportion of genotype pairs with eIBD greater than 0.55 that are likely to be closely related (siblings) or clones, and pairs with at least moderate eIBD using a cutoff of 0.30 (Figure 4A). There was a higher proportion of closely related genotypes observed in the low transmission midpoints for both populations. However, this was more pronounced in Madang where almost 10% of pairs shared high IBD and 30% with moderate levels (Figure 4A). Whereas, in ESP, there was only a small change in the proportions of related pairs. After resurgence, in Madang the proportion of moderate eIBD pairs reduced substantially though remained elevated relative to baseline levels, whereas in ESP eIBD pairs remained elevated relative to baseline.

Figure 4

Figure 4. Spatio-temporal changes in genotype relatedness within P. vivax populations of PNG. (A) Pairwise expected Identity by Descent (eIBD) was calculated using IsoRelate and stratified by geographic location and year (A) eIBD > 0.55 and eIBD > 0.30 proportions are shown for each of the populations and time points, with the three catchments for Madang shown separately. (B) Spatio-temporal changes in relatedness of P. vivax genotypes from Papua New Guinea. Pairwise expected Identity by Descent (eIBD) was calculated using IsoRelate and networks drawn on the basis of relatedness thresholds of 0.55 and 0.30. The colours indicate the study year and location, as shown in the key. Different symbols for Madang represent the different catchments.

In the network analysis, most populations show small clusters of related parasites, but Madang 2010 genotypes mostly clustered into two distinct lineages (Figure 4B). Setting the threshold to illustrate moderate relatedness (eIBD greater than 0.30), genotypes from all years connected to the smaller Madang 2010 lineage, aligning with the STRUCTURE and phylogenetic analysis. In ESP, very few closely related pairs were detected, but there was a large cluster of moderately related individuals connecting all three time points. When the two locations were analysed together, related pairs within each site are connected into one large cluster, the majority with moderate eIBD values. The smaller Madang 2010 lineage links with several baseline (2006) lineages via a dispersed cluster of 2014 genotypes suggesting the 2010 lineage recombined with rarer (i.e. undetected in the 2010 dataset) residual genotypes upon resurgence.

A minimum spanning network analysis further demonstrates P. vivax had a homogeneous diversity pattern at baseline consistent with endemic transmission (Supplementary Figure S5) while post-LLIN, the cluster of related Madang 2010 genotypes showed a star-like configuration indicative of a recent clonal expansion. ESP showed a more homogeneous pattern throughout all time points, suggestive of sustained low levels of recombination amongst distinct genotypes, and connections to the resurgent populations detected in both Madang 2014 and ESP 2016 (Supplementary Figure S5).

Discussion

P. vivax parasites are more diverse than P. falciparum therefore more substantial and sustained reductions in transmission may be needed to reduce parasite diversity from pre-control levels (Neafsey et al., 2012). Population genomic changes in P. falciparum strongly correlate with transmission intensity where increasing multilocus LD, higher allele sharing, fewer polyclonal infections and lower genetic diversity have been observed as transmission declines (Anderson et al., 2000; Daniels et al., 2013; Nkhoma et al., 2013). For P. vivax, signals such as increasing population differentiation and reductions in effective population size with significant inbreeding have been observed in P. vivax populations of Thailand (Kittichai et al., 2017) and Solomon Islands (Waltmann et al., 2018) after more than a decade of intensive control interventions. Our previous investigation of P. vivax population structure using ten highly polymorphic microsatellite markers, showed a minimal impact on parasite population diversity and structure, but increased multilocus LD with declining transmission, suggesting an increase in related individuals in the dataset (Kattenberg et al., 2020b). Here we confirm that a 88 SNP barcode and deep amplicon sequencing not only provides higher resolution insights, it identifies significant population structure and increasing proportions of related individuals with declining transmission and highly variable changes between the P. vivax populations from the two geographic areas.

In Madang, with declining transmission, SNP barcodes showed a decrease in COI, genetic diversity, and effective population size, and an increase in multilocus LD with large clusters of closely related genotype pairs. At low transmission, Madang 2010 genotypes formed at least two distinct lineages. The lower COI, genetic diversity and effective population size, and allele frequency distributions are suggestive of a prior bottleneck, while the LD, STRUCTURE and IBD results suggest an early clonal expansion. Following resurgence (13%–20% prevalence), the parasite genetic parameters returned to near baseline levels indicating population recovery with continuing transmission of at least one of the Madang 2010 lineages, and importation from East Sepik. In contrast, in East Sepik, where P. vivax prevalence declined to only 6% in 2010 with high heterogeneity across villages (Kattenberg et al., 2020a), resurging to only 9% by 2016, more subtle changes were observed. All genetic diversity parameters trended in the same direction at resurgence (2012 vs. 2016), as with transmission decline (2005 vs. 2012). While there were increased proportions of related genotype pairs there were no dominant or large clusters of closely related individuals, suggesting transmission of multiple residual lineages, rather than expansion of specific lineages. Between the two populations, the lower prevalence was accompanied by an increase in population structure and higher genetic differentiation (F_ST) between provinces, which reversed upon resurgence and likely contributed to population recovery. These variable patterns in population genetic changes observed in the two provinces at the low transmission time points may depend on the prevalence and population diversity at the commencement of interventions, the timing of cross-sectional surveys, as well as the local epidemiology, and other socio-demographic factors.

In regions with high transmission, multiple P. vivax genotypes within a single patient are common (Juliano et al., 2010; Fola et al., 2017). A reduction in these polyclonal infections as transmission declines reduces the chance of recombination of distinct Plasmodium genotypes (outcrossing) inside the mosquito, reducing the genetic diversity of parasite populations (Zhong et al., 2018). For P. vivax however, a decline in COI is expected to occur slowly due to the contribution of relapses to blood stage infections (Robinson et al., 2015). COI trended with prevalence in Madang but not in East Sepik showing that the relationship between P. vivax transmission intensity and COI is complex. The lack of association between COI and parasite prevalence in East Sepik, could be explained by the heterogeneity of malaria transmission at village level (Karl et al., 2016; Kattenberg et al., 2020a). This is consistent with other studies (Gunawardena et al., 2014; Noviyanti et al., 2015) and our previous observations across different endemic regions in PNG (Juliano et al., 2010; Fola et al., 2017). Other studies have suggested an association between P. vivax COI and transmission intensity but only in areas that have had sustained levels of transmission (Waltmann et al., 2018; Koepfli et al., 2018), suggesting the contribution of relapses to the pool of genotypes may remain high in populations with recent reductions in transmission (Mueller et al., 2009b). Thus, for P. vivax, COI can be used as a reliable marker of population reduction to measure longer term transmission intensity patterns, but may be less reliable following more recent changes in transmission.

The detection of reduced within-sample genetic diversity, as indicated by proportions of heterozygous SNPs, suggests an increase in relatedness between clones within infections and/or a reduction in the diversity of new infections and the hypnozoite reservoir. Changes in this parameter were more pronounced in Madang relative to East Sepik over the study period. In addition, the emergence of multilocus LD, together with high genetic differentiation and lower genetic diversity indicates local inbreeding (Anderson et al., 2000; Nkhoma et al., 2013). In Madang, the pattern following resurgence indicated increased random mating amongst genetically distinct individuals, while in East Sepik, LD values continued to trend upwards, albeit to a limited extent, consistent with highly focal transmission and inbreeding (Kattenberg et al., 2020a; Kattenberg et al., 2020b).

The loss of rare alleles in Madang after transmission decline, and detection of a large cluster of related genotypes in the Madang population after intensification of control, suggests a strong bottleneck had occurred which may have been accompanied by a reduction in the hypnozoite reservoir, as indicated by dramatically decreased COI and heterozygous SNPs. Accordingly, there was a significant change in population level diversity and effective population size in Madang after control, but limited changes in East Sepik. The lack of significant change in parasite population size after intervention in East Sepik is consistent with other studies where minimal or no change in parasite genetic diversity was observed with declining transmission based on microsatellites (Waltmann et al., 2018; Kittichai et al., 2017) and genome wide SNP markers (Nkhoma et al., 2013).

The spatio-temporal structure of P. vivax populations was also highly variable, with major population structure changes over time in Madang, and limited changes in East Sepik. We confirmed the relatively panmictic P. vivax population on the north coast of PNG pre-LLIN, as previously described (Jennison et al., 2015). An increase in population structure at the low transmission mid-point post-LLIN, indicates a significant impact of control on parasite migration and consequent gene flow between provinces. However, there was a 1 to 2 year gap between cross-sectional surveys which may have influenced this result. At low transmission, Madang also had a somewhat fragmented parasite population structure at the local level with Utu parasites clustering separately to Malala and Mugil isolates in the phylogeny, STRUCTURE and IBD analyses. One source of resurgence is indicated by the connectivity of the lineage from Utu in Madang 2010, and the broader Madang 2014 population, which also connects to certain Madang 2006 genotypes in the IBD network. Persistence of genotypes across timepoints, or found only in the higher transmission baseline and final timepoints indicates residual parasite lineages could be a major source of resurgence (Daniels et al., 2013). Interestingly, there was an increased prevalence of P. vivax dihydrofolate reductase triple and quadruple mutants, in the Madang parasite population in 2010 (Barnadas et al., 2015) suggesting that drug resistance may also play a role (Barnadas et al., 2015). Drug resistance data however is not available for the 2014 P. vivax populations so would need further investigation.

The East Sepik genetic diversity at baseline was similar to post-LLIN levels in Madang, whereas there was a slightly lower baseline COI than Madang. ESP P. vivax populations may have already been impacted by control pressure, because there were differences in intervention coverage before upscaling control nationally, with earlier distributions of LLIN in East Sepik (Hetzel et al., 2015). Additionally ESP may have had a highly heterogeneous transmission prior to any control intervention. The heterogeneous transmission (Kattenberg et al., 2020a) and small sample size at low transmission (only 20 isolates in the 2012 population from two villages) may have also prevented the detection of clusters and fluctuations in allele frequencies. Topography and connecting roads, uneven sampling, sample size variation and presence of isolated catchment areas within provinces may explain the variable transmission of malaria across PNG (Cleary et al., 2021), and could influence parasite population structure. Variability in parasite population connectivity via human movement for social and other business activities to urban areas (like Madang town) likely contribute to the patterns observed in this study, with the East Sepik study sites being more remote and less populated. Movement of parasites from other endemic areas in PNG, which would be more likely in Madang, and resultant connectivity between parasite populations regionally may also contribute to maintaining a diverse gene pool (Fola et al., 2018; Tessema et al., 2019; Ferreira and de Oliveira, 2015).

The markers used for detecting population genetic signals of transmission dynamics are critically important (Sutton, 2013; Ruybal-Pesántez et al., 2024). The results we have observed with the 88 SNPs are striking relative to our previous observations with ten microsatellite markers where no population structure was observed, with a single panmictic population found across all timepoints and locations (Kattenberg et al., 2020b). The SNP barcode therefore provided more sensitive and meaningful insights into transmission dynamics, as we previously observed using only the baseline samples (Fola et al., 2020). The greater resolution is likely due to the higher density of markers across the genome which facilitates the detection of moderately related genotypes, whereas microsatellites, with less than one marker per chromosome, likely only detects very closely related genotypes (Cabrera-Sosa et al., 2024). This warrants the exploration of sources of resurgence through SNP barcoding of additional isolates from these and other ongoing epidemiological surveys in PNG (Mueller et al., 2022).

Our study has some limitations that may be addressed in future studies. Firstly, the sample size was variable across timepoints and study sites, with some low sample sizes. This was the result of either lower density parasitaemias, the quality of the DNA or limited availability due to low prevalence. While every effort was made to produce similar sample sizes across populations, this may impact the accuracy of statistical outputs, particularly in populations where the sample size is small. While the 88-SNP barcode used for this study provides useful population-level insights, a barcode with more SNPs may provide greater resolution to detect related individuals or inbreeding, and whole genome sequences (WGS) would provide details on a range of features including regions under selection, such as drug resistance, or adaptations to lower transmission timepoints. Currently however, the cost of WGS (approximately USD150-200) compared to SNP barcodes (approximately USD20) is too high to allow large-scale surveillance. In previous studies, relatedness of parasites between geographic locations have been linked to human mobility data which, as mentioned above, may explain some of the differences between provinces (Chang et al., 2019).

Overall, the results demonstrate the variable impact of malaria control on parasite population genetic structure within the same high transmission region, despite similarities in parasite prevalence and reductions due to control, underscoring the critical role of genomic surveillance to detect underlying transmission dynamics that cannot be inferred from prevalence alone. Whilst the nationwide distribution of LLIN and other control measures in PNG made a major impact on the P. vivax transmission dynamics in Madang with the initial reduction in prevalence, it was not sustained. Resurgence and expansion of the residual P. vivax populations and increasing movement of parasites via the higher prevalence and human movement appears to have been an important contributor to the patterns observed. In East Sepik, prior control efforts and high heterogeneity of transmission may have contributed to the smaller resurgence and more subtle changes in parasite population genetic parameters. This research highlights the power of genomic surveillance to track the impact of malaria control and identify the possible sources of malaria resurgence. These observations can aid in the prioritisation of targeted control and reduce the re-emergence and spread of residual infections. For P. vivax, maintaining control pressure for long periods will be essential to exert sustainable impacts on transmission and prevent rapid population recovery.

Methods and materials

Samples and study sites

Human blood samples for the study were collected during serial cross-sectional studies conducted in the Mugil, Malala and Utu areas of Madang Province in 2006, 2010 and 2014 and Maprik area of East Sepik Province (ESP) in 2005, 2012–13 and 2016 (Mueller et al., 2009a; Senn et al., 2012; Koepfli et al., 2017; Kattenberg et al., 2020a). These geographic areas are lowland, inland and coastal areas that are hyperendemic for malaria with a high prevalence of both major human malaria parasites P. falciparum and P. vivax. In Madang, catchments are found approximately 30–50 km from the urban centre of Madang Town, with Utu inland, and Mugil and Malala located along the North Coast Highway, whilst in ESP, all villages were located in an inland area to the south, east and west of Maprik town (Figure 1). In both provinces, there were major reductions in P. vivax prevalence from the first time point (2005/6) to the second time point (Madang 2010 and ESP 2012–13) due to intensified control efforts in the intervening years (Figure 1; Table 1). Following this (Madang 2014 and ESP 2016), there was a documented resurgence of malaria overall (Cleary et al., 2022), and an increase in P. vivax prevalence in the study sites (Koepfli et al., 2017, unpublished data) (Figure 1; Table 1).

Ethical approval for the study was provided by the Institutional Review Board of the PNG Institute of Medical Research No. 12.21, the Medical Research Advisory Council of PNG No. 13.08 and the Walter and Eliza Hall Human Research Ethics Committee No. 13.14.

Sample processing and genotyping

Processing of samples included extraction of genomic DNA using Qiagen or Favorgen kits, followed by screening for Plasmodium species infection using both light microscopy and species-specific qPCR targeting the 18s rRNA genes (Kattenberg et al., 2020a; Koepfli et al., 2017). In previous studies, P. vivax positive samples were subject to genotyping for multiplicity of infection (MOI) using two highly polymorphic microsatellite markers (e.g. msp1f3, ms16) (Arnott et al., 2013). Samples with only one or two clones detected based on this genotyping were prioritised to avoid issues associated with reconstructing haplotypes from the SNP barcoding.

To streamline the protocol and reduce costs, a subset of 100 the 176 validated SNPs (Fola et al., 2020) were selected to obtain relatively even coverage of the genome (Supplementary Table S1). A total of 376 P. vivax positive isolates were genotyped using these 100 SNP markers, with all procedures conducted as previously described SNPs (Fola et al., 2020) with modifications in the number of pools of multiplexed SNP loci to accompany the smaller number of markers. Following data cleaning and filtering of SNPs and samples due to missing data, a total of 88 SNPs were used for further analysis.

Bioinformatic analysis

The raw FASTQ files were demultiplexed by binning based on the MID index, the read quality was checked using FastQC (Version 0.8.0) (Andrews and Bittencourt, 2010) and combined FastQC output for all samples were visualized using MultiQC (Ewels et al., 2016). Low-quality reads (<Q30), adaptors, primers, and reads shorter or longer than expected size of amplicon were trimmed using Trimmomatic (Bolger et al., 2014). Only reads that passed stringent quality filters progressed for alignment and variant calling. Unmapped BAM files were generated from quality filtered and trimmed FASTQ files using the FastqtoSam function (http://broadinstitute.github.io/picard/). The combined pipeline was then used to generate indexed, mapped BAM files. This pipeline consists of SamToFastq, bwa-mem and MergeBamAlignment to map reads, and generates a clean and indexed mapped BAM file. In brief, the sequenced reads were mapped to the P. vivax Salvador I strain reference genome using BWA MEM (Vasimuddin et al., 2019). Overall quality and genome coverage of mapped bam files were checked using QualiMap v.2.2.1 (García-Alcalde et al., 2012) and detailed bioinformatics described in our previous publication (Fola et al., 2020).

In sequencing large sample sets on multiple MiSeq runs, batch effects are inevitable, thus we used the same sequencing protocol and bioinformatics tools to analyze data from different sequencing runs. However, there were still some batch effects observed in the data set when SNPs were called for each run separately due to allele frequency variation between batches (Supplementary Figure S6A). Combining the data for variant calling was used to minimize the batch effect introduced by the different runs. This helped in reducing induced batch effects as seen in the multidimensional scale plots (Supplementary Figure S6B).

Complexity of infection and heterozygous SNPs

In malaria endemic countries like PNG, individuals often harbour complex genetically distinct clones of the same species, defined by the detection of multiple alleles in genotyping assays. The number of clones is known as the multiplicity of infection (MOI) and can be measured using the proportion of heterozygous SNP calls (Auburn et al., 2012). Expecting to observe less heterozygous SNP calls with declining transmission, we determined the proportion of heterozygous SNP calls for samples genotyped from different time points. The most likely number of clones present in each sample was calculated using the Real McCOIL (Chang et al., 2017) and only monoclonal or dominant genotypes constructed from the major alleles in polyclonal samples were used for population genetic analysis. Moreover, assuming recent intensive malaria control activities may lead to a change in the frequency of minor alleles with a loss of rare variants with declining transmission, minor allele frequencies (MAF) were determined for all SNPs for all genotyped samples from different time points. The MAF was computed as the proportion of genotyped samples carrying the genotype that was least common.

Multilocus linkage disequilibrium and effective population size

Multilocus linkage disequilibrium (LD) was determined in each parasite population collected from different time points expecting evidence of inbreeding or increased LD over time with declining transmission. Multilocus LD was measured using the statistic I_AS (standardized index of association), which compares the observed variance in the numbers of alleles shared between genotypes with that expected when genotypes share no alleles at different loci (Haubold and Hudson, 2000). I_A^S was calculated using LIAN version 3.6 software by applying a Monte Carlo test with 100,000 re-sampling steps. LIAN compares the observed association of markers to the values expected for random association based on the population diversity. To assess the effect of control on parasite effective population size (N_e), we examined population size using the method implemented in NeEstimator v. 2.0 software (Do et al., 2014), we considered that P. vivax has a similar generation time as P. falciparum (2 months) (Nkhoma et al., 2013) with a maximum N_e value of 20,000 assigned in the software since the P. vivax population has a higher population size than P. falciparum (Koepfli et al., 2013; Jennison et al., 2015).

Genetic diversity and population structure

The SNPs were converted to numeric data (i.e., target loci contain reference alleles = “0” and target loci contain no reference alleles = “1”) and used as an input file to determine genetic diversity and population structure analysis. Then the creation of input files for genetic analysis was performed using CONVERT version 1.3.1 (Lischer and Excoffier, 2012). We calculated nucleotide diversity (SNP π) and Wright’s F_ST using DnaSP Version 5.0 (Librado and Rozas, 2009). We measured parasite population structure and genetic differentiation between parasites sampled before and after an enhanced intervention. The Bayesian clustering software, STRUCTURE version 2.3.4 (Pritchard et al., 2000; Pritchard et al., 2000) was used to determine the number of population clusters (K) and whether haplotypes clustered according to geographical origin. We performed STRUCTURE runs with a burn-in period of 100,000 followed by 100,000 Markov chains (MCMC iterations). Evanno’s method of ΔK statistics implemented in the STRUCTURE HARVESTER (Earl and vonHoldt, 2012) to determine best genetic clusters. CLUMPAK-Cluster Markov Packager Across K web-based server was used for summation and graphical representation of STRUCTURE results (Kopelman et al., 2015).

Phylogenetic and relatedness analysis

We calculated genotype relatedness in each year based on the similarity of their SNP barcodes within the population. An increase in genotype relatedness (carrying similar SNP barcodes) is expected since there is a low chance of recombination of different clones inside mosquito with declining transmission. Pairwise comparisons of the number of alleles shared (PS) was computed based on the metric 1-PS using the “Ape” R package ‘dist.gene’ command. This was used to conduct multidimensional scaling (MDS) and principle components analysis (PCA). The phylogenetic relationships were inferred using the Neighbor-Joining method (and computed using the number of differences method in MEGA version 11 (Tamura et al., 2021). The goeBURST algorithm within PHYLOViZ 2.0 (Francisco et al., 2012) was used to generate a minimum spanning tree.

To measure pairwise relatedness values, we calculated identity by descent (IBD) using IsoRelate (Henden et al., 2018). IsoRelate uses a hidden Markov model (HMM) to detect genomic regions that are identical by descent (IBD) and aids simultaneous detection of parasite population clustering (https://github.com/bahlolab/isoRelate). We assume that all isolates are monoclonal and retained only isolates with missingness less than or equal to 30%, and SNP loci with missingness less than or equal to 30% and MAF less than or equal to 0.01, and 1 cM was considered equivalent to 17 kbp as for P. falciparum, (Henden et al., 2018). IBD parameters were then measured for each pair of genotypes using the ‘getIBDparameters’ command. Pairwise posterior probabilities of IBD sharing was estimated at each locus (i.e. the probability that isolate A and isolate B are IBD at a particular SNP locus, under the implemented HMM of relatedness and the parameters inferred in step 1) using the ‘getIBDposterior’ command. We then computed the expected fraction IBD (eIBD), defined to be the mean probability of IBD sharing across all loci (Taylor et al., 2017) for each pair of isolates. Note that the isoRelate output was corrected by a factor of 2. Sparse SNP barcode data cannot reliably infer low-level IBD sharing, with previous validation of a 155-SNP barcode for P. falciparum against whole genome sequence data indicating a lower resolution limit of 0.25 (i.e. ∼30 SNPs, data not shown). We selected eIBD threshold values of 0.30 and 0.55 for exploration of relatedness within and between P. vivax populations.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found in the European Nucleotide Archive (ENA) at EMBL-EBI with the accession number PRJEB105390.

Ethics statement

The studies involving humans were approved by PNG Institute of Medical Research Institutional Review Board No. 12.21, Medical Research Advisory Council of PNG No. 13.08, Walter and Eliza Hall Human Research Ethics Committee No. 13.14. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.

Author contributions

AF: Data curation, Methodology, Formal Analysis, Investigation, Writing – original draft. SM: Visualization, Formal Analysis, Writing – review and editing, Methodology, Investigation. ZR: Supervision, Methodology, Writing – review and editing, Investigation. DL-G: Resources, Writing – review and editing. EN: Resources, Writing – review and editing. SL: Writing – review and editing, Resources, Methodology. JoK: Resources, Writing – review and editing. CK: Writing – review and editing, Resources. JaK: Project administration, Writing – review and editing, Funding acquisition, Resources. MO-K: Resources, Project administration, Writing – review and editing. ML: Resources, Funding acquisition, Project administration, Writing – review and editing. LR: Writing – review and editing, Funding acquisition, Resources, Project administration. IM: Writing – review and editing, Resources, Supervision, Conceptualization, Project administration, Funding acquisition. AB: Resources, Data curation, Conceptualization, Methodology, Writing – original draft, Supervision, Writing – review and editing, Project administration, Formal Analysis, Funding acquisition.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This research was funded by the National Institutes of Allergy and Infectious Disease, USA U19 AI129382 and AI089686 and the National Health and Medical Research Council (NHMRC) of Australia Project Grant GNT1027108. The authors are grateful for support from the Victorian State Government Operational Infrastructure Support and Australian Government NHMRC Independent Research Institute Infrastructure Support Scheme. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

Acknowledgements

We are grateful to all participants, their families and the communities of Madang and East Sepik Province who supported the study. In addition, we would like to thank the staff of the PNG Institute of Medical Research and other investigators who contributed to the collection of samples and data management for the field studies.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author AB declared that they were an editorial board member of Frontiers at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2025.1621920/full#supplementary-material

Supplementary File 1 | Processed and cleaned dataset.

Supplementary Datasheet 1 | Figures S1 - S6, Table S1.

References

Adams, J. H., and Mueller, I. (2017). The biology of Plasmodium vivax. Cold Spring Harb. Perspect. Med. 7, a025585. doi:10.1101/cshperspect.a025585

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, T. J., Haubold, B., Williams, J. T., Estrada-Franco, J. G., Richardson, L., Mollinedo, R., et al. (2000). Microsatellite markers reveal a spectrum of population structures in the malaria parasite Plasmodium falciparum. Mol. Biol. Evol. 17, 1467–1482. doi:10.1093/oxfordjournals.molbev.a026247

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrews, S., and Bittencourt, S., (2010). FastQC: a quality control tool for high throughput sequence data. Berlin, Germany: ScienceOpen, Inc.

Google Scholar

Arnott, A., Barnadas, C., Senn, N., Siba, P., Mueller, I., Reeder, J. C., et al. (2013). High genetic diversity of Plasmodium vivax on the north coast of Papua New Guinea. Am. J. Trop. Med. Hyg. 89, 188–194. doi:10.4269/ajtmh.12-0774

PubMed Abstract | CrossRef Full Text | Google Scholar

Auburn, S., and Barry, A. E. (2017). Dissecting malaria biology and epidemiology using population genetics and genomics. Int. J. Parasitol. 47, 77–85. doi:10.1016/j.ijpara.2016.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Auburn, S., Campino, S., Miotto, O., Djimde, A. A., Zongo, I., Manske, M., et al. (2012). Characterization of within-host Plasmodium falciparum diversity using next-generation sequence data. PLoS One 7, e32891. doi:10.1371/journal.pone.0032891

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnadas, C., Timinao, L., Javati, S., Iga, J., Malau, E., Koepfli, C., et al. (2015). Significant geographical differences in prevalence of mutations associated with Plasmodium falciparum and Plasmodium vivax drug resistance in two regions from Papua New Guinea. Malar. J. 14, 399. doi:10.1186/s12936-015-0879-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Barry, A. E., Waltmann, A., Koepfli, C., Barnadas, C., and Mueller, I. (2015). Uncovering the transmission dynamics of Plasmodium vivax using population genetics. Pathog. Glob. Health 109, 142–152. doi:10.1179/2047773215Y.0000000012

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics 30, 2114–2120. doi:10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Buyon, L. E., Santamaria, A. M., Early, A. M., Quijada, M., Barahona, I., Lasso, J., et al. (2020). Population genomics of Plasmodium vivax in Panama to assess the risk of case importation on malaria elimination. PLoS Negl. Trop. Dis. 14, e0008962. doi:10.1371/journal.pntd.0008962

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabrera-Sosa, L., Safarpour, M., Kattenberg, J. H., Ramirez, R., Vinetz, J. M., Rosanas-Urgell, A., et al. (2024). Comparing newly developed SNP barcode panels with microsatellites to explore population genetics of malaria parasites in the Peruvian Amazon. Front. Genet. 15, 1488109. doi:10.3389/fgene.2024.1488109

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, H.-H., Worby, C. J., Yeka, A., Nankabirwa, J., Kamya, M. R., Staedke, S. G., et al. (2017). The real McCOIL: a method for the concurrent estimation of the complexity of infection and SNP allele frequency for malaria parasites. PLoS Comput. Biol. 13, e1005348. doi:10.1371/journal.pcbi.1005348

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, H.-H., Wesolowski, A., Sinha, I., Jacob, C. G., Mahmud, A., Uddin, D., et al. (2019). Mapping imported malaria in Bangladesh using parasite genetic and human mobility data. Elife 8, e43481. doi:10.7554/eLife.43481

PubMed Abstract | CrossRef Full Text | Google Scholar

Cibulskis, R. E., Alonso, P., Aponte, J., Aregawi, M., Barrette, A., Bergeron, L., et al. (2016). Malaria: global progress 2000 - 2015 and future challenges. Infect. Dis. Poverty 5, 61. doi:10.1186/s40249-016-0151-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Cleary, E., Hetzel, M. W., Siba, P. M., Lau, C. L., and Clements, A. C. A. (2021). Spatial prediction of malaria prevalence in Papua New Guinea: a comparison of Bayesian decision network and multivariate regression modelling approaches for improved accuracy in prevalence prediction. Malar. J. 20, 269. doi:10.1186/s12936-021-03804-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cleary, E., Hetzel, M. W., and Clements, A. C. A. (2022). A review of malaria epidemiology and control in Papua New Guinea 1900 to 2021: progress made and future directions. Front. Epidemiol. 2, 980795. doi:10.3389/fepid.2022.980795

PubMed Abstract | CrossRef Full Text | Google Scholar

Colijn, C., and Gardy, J. (2014). Phylogenetic tree shapes resolve disease transmission patterns. Evol. Med. Public Health 2014, 96–108. doi:10.1093/emph/eou018

PubMed Abstract | CrossRef Full Text | Google Scholar

Dalmat, R., Naughton, B., Kwan-Gett, T. S., Slyker, J., and Stuckey, E. M. (2019). Use cases for genetic epidemiology in malaria elimination. Malar. J. 18, 163. doi:10.1186/s12936-019-2784-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Daniels, R., Chang, H.-H., Séne, P. D., Park, D. C., Neafsey, D. E., Schaffner, S. F., et al. (2013). Genetic surveillance detects both clonal and epidemic transmission of malaria following enhanced intervention in Senegal. PLoS One 8, e60780. doi:10.1371/journal.pone.0060780

PubMed Abstract | CrossRef Full Text | Google Scholar

Do, C., Waples, R. S., Peel, D., Macbeth, G. M., Tillett, B. J., and Ovenden, J. R. (2014). NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol. Ecol. Resour. 14, 209–214. doi:10.1111/1755-0998.12157

PubMed Abstract | CrossRef Full Text | Google Scholar

Earl, D. A., and vonHoldt, B. M. (2012). Structure Harvester: a website and program for visualizing Structure output and implementing the evanno method. Conserv. Genet. Resour. 4, 359–361. doi:10.1007/s12686-011-9548-7

CrossRef Full Text | Google Scholar

Ewels, P., Magnusson, M., Lundin, S., and Käller, M. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047–3048. doi:10.1093/bioinformatics/btw354

PubMed Abstract | CrossRef Full Text | Google Scholar

Fadilah, I., Djaafara, B. A., Lestari, K. D., Fajariyani, S. B., Sunandar, E., Makamur, B. G., et al. (2022). Quantifying spatial heterogeneity of malaria in the endemic Papua region of Indonesia: analysis of epidemiological surveillance data. Lancet Regional Health - Southeast Asia 5, 100051. doi:10.1016/j.lansea.2022.100051

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, M. U., and de Oliveira, T. C. (2015). Challenges for Plasmodium vivax malaria elimination in the genomics era. Pathog. Glob. Health 109, 89–90. doi:10.1179/2047772415Z.000000000263

PubMed Abstract | CrossRef Full Text | Google Scholar

Fola, A. A., Harrison, G. L. A., Hazairin, M. H., Barnadas, C., Hetzel, M. W., Iga, J., et al. (2017). Higher complexity of infection and genetic diversity of Plasmodium vivax than Plasmodium falciparum across all malaria transmission zones of Papua New Guinea. Am. J. Trop. Med. Hyg. 96, 630–641. doi:10.4269/ajtmh.16-0716

PubMed Abstract | CrossRef Full Text | Google Scholar

Fola, A. A., Nate, E., Abby Harrison, G. L., Barnadas, C., Hetzel, M. W., Iga, J., et al. (2018). Nationwide genetic surveillance of Plasmodium vivax in Papua New Guinea reveals heterogeneous transmission dynamics and routes of migration amongst subdivided populations. Infect. Genet. Evol. 58, 83–95. doi:10.1016/j.meegid.2017.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Fola, A. A., Kattenberg, E., Razook, Z., Lautu-Gumal, D., Lee, S., Mehra, S., et al. (2020). SNP barcodes provide higher resolution than microsatellite markers to measure Plasmodium vivax population genetics. Malar. J. 19, 375. doi:10.1186/s12936-020-03440-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fola, A. A., Feleke, S. M., Mohammed, H., Brhane, B. G., Hennelly, C. M., Assefa, A., et al. (2023). Clonal spread ofPlasmodium falciparumcandidate artemisinin partial resistanceKelch13622I mutation and co-occurrence withpfhrp2/3deletions in Ethiopia. medRxiv Pre print. doi:10.1101/2023.03.02.23286711

CrossRef Full Text | Google Scholar

Francisco, A. P., Vaz, C., Monteiro, P. T., Melo-Cristino, J., Ramirez, M., and Carriço, J. A. (2012). PHYLOViZ: phylogenetic inference and data visualization for sequence based typing methods. BMC Bioinforma. 13, 87. doi:10.1186/1471-2105-13-87

PubMed Abstract | CrossRef Full Text | Google Scholar

García-Alcalde, F., Okonechnikov, K., Carbonell, J., Cruz, L. M., Götz, S., Tarazona, S., et al. (2012). Qualimap: evaluating next-generation sequencing alignment data. Bioinformatics 28, 2678–2679. doi:10.1093/bioinformatics/bts503

PubMed Abstract | CrossRef Full Text | Google Scholar

Gul, D., Rodríguez-Rodríguez, D., Nate, E., Auwan, A., Salib, M., Lorry, L., et al. (2021). Investigating differences in village-level heterogeneity of malaria infection and household risk factors in Papua New Guinea. Sci. Rep. 11, 16540. doi:10.1038/s41598-021-95959-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Gunawardena, S., Ferreira, M. U., Kapilananda, G. M. G., Wirth, D. F., and Karunaweera, N. D. (2014). The Sri Lankan paradox: high genetic diversity in Plasmodium vivax populations despite decreasing levels of malaria transmission. Parasitology 141, 880–890. doi:10.1017/S0031182013002278

PubMed Abstract | CrossRef Full Text | Google Scholar

Haubold, B., and Hudson, R. R. (2000). LIAN 3.0: detecting linkage disequilibrium in multilocus data. Bioinformatics 16, 847–849. doi:10.1093/bioinformatics/16.9.847

PubMed Abstract | CrossRef Full Text | Google Scholar

Henden, L., Lee, S., Mueller, I., Barry, A., and Bahlo, M. (2018). Identity-by-descent analyses for measuring population dynamics and selection in recombining pathogens. PLoS Genet. 14, e1007279. doi:10.1371/journal.pgen.1007279

PubMed Abstract | CrossRef Full Text | Google Scholar

Hetzel, M. W., Morris, H., Tarongka, N., Barnadas, C., Pulford, J., Makita, L., et al. (2015). Prevalence of malaria across Papua New Guinea after initial roll-out of insecticide-treated mosquito nets. Trop. Med. Int. Health 20, 1745–1755. doi:10.1111/tmi.12616

PubMed Abstract | CrossRef Full Text | Google Scholar

Hetzel, M. W., Pulford, J., Ura, Y., Jamea-Maiasa, S., Tandrapah, A., Tarongka, N., et al. (2017). Insecticide-treated nets and malaria prevalence, Papua New Guinea, 2008-2014. Bull. World Health Organ. 95, 695–705B. doi:10.2471/BLT.16.189902

PubMed Abstract | CrossRef Full Text | Google Scholar

Ippolito, M. M., Moser, K. A., Kabuya, J.-B. B., Cunningham, C., and Juliano, J. J. (2021). Antimalarial drug resistance and implications for the WHO global technical strategy. Curr. Epidemiol. Rep. 8, 46–62. doi:10.1007/s40471-021-00266-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Jennison, C., Arnott, A., Tessier, N., Tavul, L., Koepfli, C., Felger, I., et al. (2015). Plasmodium vivax populations are more genetically diverse and less structured than sympatric Plasmodium falciparum populations. PLoS Negl. Trop. Dis. 9, e0003634. doi:10.1371/journal.pntd.0003634

PubMed Abstract | CrossRef Full Text | Google Scholar

Juliano, J. J., Porter, K., Mwapasa, V., Sem, R., Rogers, W. O., Ariey, F., et al. (2010). Exposing malaria in-host diversity and estimating population diversity by capture-recapture using massively parallel pyrosequencing. Proc. Natl. Acad. Sci. U. S. A. 107, 20138–20143. doi:10.1073/pnas.1007068107

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaeuffer, R., Réale, D., Coltman, D. W., and Pontier, D. (2007). Detecting population structure using STRUCTURE software: effect of background linkage disequilibrium. Hered. (Edinb.) 99, 374–380. doi:10.1038/sj.hdy.6801010

PubMed Abstract | CrossRef Full Text | Google Scholar

Karl, S., White, M. T., Milne, G. J., Gurarie, D., Hay, S. I., Barry, A. E., et al. (2016). Spatial effects on the multiplicity of Plasmodium falciparum infections. PLoS One 11, e0164054. doi:10.1371/journal.pone.0164054

PubMed Abstract | CrossRef Full Text | Google Scholar

Kattenberg, J. H., Gumal, D. L., Ome-Kaius, M., Kiniboro, B., Philip, M., Jally, S., et al. (2020a). The epidemiology of Plasmodium falciparum and Plasmodium vivax in east sepik province, Papua New Guinea, pre- and post-implementation of national malaria control efforts. Malar. J. 19, 198. doi:10.1186/s12936-020-03265-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kattenberg, J. H., Razook, Z., Keo, R., Koepfli, C., Jennison, C., Lautu-Gumal, D., et al. (2020b). Monitoring Plasmodium falciparum and Plasmodium vivax using microsatellite markers indicates limited changes in population structure after substantial transmission decline in Papua New Guinea. Mol. Ecol. 29, 4525–4541. doi:10.1111/mec.15654

PubMed Abstract | CrossRef Full Text | Google Scholar

Kittichai, V., Koepfli, C., Nguitragool, W., Sattabongkot, J., and Cui, L. (2017). Substantial population structure of Plasmodium vivax in Thailand facilitates identification of the sources of residual transmission. PLoS Negl. Trop. Dis. 11, e0005930. doi:10.1371/journal.pntd.0005930

PubMed Abstract | CrossRef Full Text | Google Scholar

Koepfli, C., Ross, A., Kiniboro, B., Smith, T. A., Zimmerman, P. A., Siba, P., et al. (2011). Multiplicity and diversity of Plasmodium vivax infections in a highly endemic region in Papua New Guinea. PLoS Negl. Trop. Dis. 5, e1424. doi:10.1371/journal.pntd.0001424

PubMed Abstract | CrossRef Full Text | Google Scholar

Koepfli, C., Timinao, L., Antao, T., Barry, A. E., Siba, P., Mueller, I., et al. (2013). A large Plasmodium vivax reservoir and little population structure in the south Pacific. PLoS One 8, e66041. doi:10.1371/journal.pone.0066041

PubMed Abstract | CrossRef Full Text | Google Scholar

Koepfli, C., Ome-Kaius, M., Jally, S., Malau, E., Maripal, S., Ginny, J., et al. (2017). Sustained malaria control over an 8-Year period in Papua New Guinea: the challenge of low-density asymptomatic plasmodium infections. J. Infect. Dis. 216, 1434–1443. doi:10.1093/infdis/jix507

PubMed Abstract | CrossRef Full Text | Google Scholar

Koepfli, C., Waltmann, A., Ome-Kaius, M., Robinson, L. J., and Mueller, I. (2018). Multiplicity of infection is a poor predictor of village-level Plasmodium vivax and P. falciparum population prevalence in the southwest Pacific. Open Forum Infect. Dis. 5, ofy240. doi:10.1093/ofid/ofy240

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., and Mayrose, I. (2015). Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15, 1179–1191. doi:10.1111/1755-0998.12387

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y.-C., Korol, A. B., Fahima, T., Beiles, A., and Nevo, E. (2002). Microsatellites: genomic distribution, putative functions and mutational mechanisms: a review. Mol. Ecol. 11, 2453–2465. doi:10.1046/j.1365-294x.2002.01643.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Librado, P., and Rozas, J. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. doi:10.1093/bioinformatics/btp187

PubMed Abstract | CrossRef Full Text | Google Scholar

Lischer, H. E. L., and Excoffier, L. (2012). PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics 28, 298–299. doi:10.1093/bioinformatics/btr642

PubMed Abstract | CrossRef Full Text | Google Scholar

Miotto, O., Sekihara, M., Tachibana, S.-I., Yamauchi, M., Pearson, R. D., Amato, R., et al. (2020). Emergence of artemisinin-resistant Plasmodium falciparum with kelch13 C580Y mutations on the island of new Guinea. PLoS Pathog. 16, e1009133. doi:10.1371/journal.ppat.1009133

PubMed Abstract | CrossRef Full Text | Google Scholar

Mueller, I., Widmer, S., Michel, D., Maraga, S., McNamara, D., Kiniboro, B., et al. (2009a). High sensitivity detection of plasmodium species reveals positive correlations between infections of different species, shifts in age distribution and reduced local variation in Papua New Guinea. Malar. J. 8, 41. doi:10.1186/1475-2875-8-41

PubMed Abstract | CrossRef Full Text | Google Scholar

Mueller, I., Galinski, M. R., Baird, J. K., Carlton, J. M., Kochar, D. K., Alonso, P. L., et al. (2009b). Key gaps in the knowledge of Plasmodium vivax, a neglected human malaria parasite. Lancet Infect. Dis. 9, 555–566. doi:10.1016/S1473-3099(09)70177-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Mueller, I., Vantaux, A., Karl, S., Laman, M., Witkowski, B., Pepey, A., et al. (2022). Asia-Pacific ICEMR: understanding malaria transmission to accelerate malaria elimination in the Asia Pacific region. Am. J. Trop. Med. Hyg. 107, 131–137. doi:10.4269/ajtmh.21-1336

PubMed Abstract | CrossRef Full Text | Google Scholar

Neafsey, D. E., Galinsky, K., Jiang, R. H. Y., Young, L., Sykes, S. M., Saif, S., et al. (2012). The malaria parasite Plasmodium vivax exhibits greater genetic diversity than Plasmodium falciparum. Nat. Genet. 44, 1046–1050. doi:10.1038/ng.2373

PubMed Abstract | CrossRef Full Text | Google Scholar

Neafsey, D. E., Taylor, A. R., and MacInnis, B. L. (2021). Advances and opportunities in malaria population genomics. Nat. Rev. Genet. 22, 502–517. doi:10.1038/s41576-021-00349-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Nkhoma, S. C., Nair, S., Al-Saai, S., Ashley, E., McGready, R., Phyo, A. P., et al. (2013). Population genetic correlates of declining transmission in a human pathogen. Mol. Ecol. 22, 273–285. doi:10.1111/mec.12099

PubMed Abstract | CrossRef Full Text | Google Scholar

Noviyanti, R., Coutrier, F., Utami, R. A. S., Trimarsanto, H., Tirta, Y. K., Trianty, L., et al. (2015). Contrasting transmission dynamics of Co-endemic Plasmodium vivax and P. falciparum: implications for malaria control and elimination. PLoS Negl. Trop. Dis. 9, e0003739. doi:10.1371/journal.pntd.0003739

PubMed Abstract | CrossRef Full Text | Google Scholar

Omedo, I., Mogeni, P., Bousema, T., Rockett, K., Amambua-Ngwa, A., Oyier, I., et al. (2017). Micro-epidemiological structuring of Plasmodium falciparum parasite populations in regions with varying transmission intensities in Africa. Wellcome Open Res. 2, 10. doi:10.12688/wellcomeopenres.10784.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. doi:10.1093/genetics/155.2.945

PubMed Abstract | CrossRef Full Text | Google Scholar

Rabinovich, R. N., Drakeley, C., Djimde, A. A., Hall, B. F., Hay, S. I., Hemingway, J., et al. (2017). malERA: an updated research agenda for malaria elimination and eradication. PLoS Med. 14, e1002456. doi:10.1371/journal.pmed.1002456

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, L. J., Wampfler, R., Betuela, I., Karl, S., White, M. T., Li Wai Suen, C. S. N., et al. (2015). Strategies for understanding and reducing the Plasmodium vivax and Plasmodium ovale hypnozoite reservoir in Papua New Guinean children: a randomised placebo-controlled trial and mathematical model. PLoS Med. 12, e1001891. doi:10.1371/journal.pmed.1001891

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruybal-Pesántez, S., McCann, K., Vibin, J., Siegel, S., Auburn, S., and Barry, A. E. (2024). Molecular markers for malaria genetic epidemiology: progress and pitfalls. Trends Parasitol. 40, 147–163. doi:10.1016/j.pt.2023.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Senn, N., Rarau, P., Stanisic, D. I., Robinson, L., Barnadas, C., Manong, D., et al. (2012). Intermittent preventive treatment for malaria in Papua New Guinean infants exposed to Plasmodium falciparum and P. vivax: a randomized controlled trial. PLoS Med. 9 (3), e1001195. doi:10.1371/journal.pmed.1001195

PubMed Abstract | CrossRef Full Text | Google Scholar

Sutton, P. L. (2013). A call to arms: on refining Plasmodium vivax microsatellite marker panels for comparing global diversity. Malar. J. 12, 447. doi:10.1186/1475-2875-12-447

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamura, K., Stecher, G., and Kumar, S. (2021). MEGA11: molecular evolutionary genetics analysis version 11. Mol. Biol. Evol. 38, 3022–3027.

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, A. R., Schaffner, S. F., Cerqueira, G. C., Nkhoma, S. C., Anderson, T. J. C., Sriprawat, K., et al. (2017). Quantifying connectivity between local Plasmodium falciparum malaria parasite populations using identity by descent. PLoS Genet. 13, e1007065. doi:10.1371/journal.pgen.1007065

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, A. R., Jacob, P. E., Neafsey, D. E., and Buckee, C. O. (2019). Estimating relatedness between malaria parasites. Genetics 212, 1337–1351. doi:10.1534/genetics.119.302120

PubMed Abstract | CrossRef Full Text | Google Scholar

Tessema, S., Wesolowski, A., Chen, A., Murphy, M., Wilheim, J., Mupiri, A.-R., et al. (2019). Using parasite genetic and human mobility data to infer local and cross-border malaria connectivity in Southern Africa. Elife 8, e43510. doi:10.7554/eLife.43510

PubMed Abstract | CrossRef Full Text | Google Scholar

Tusting, L. S., Bousema, T., Smith, D. L., and Drakeley, C. (2014). Measuring changes in Plasmodium falciparum transmission: precision, accuracy and costs of metrics. Adv. Parasitol. 84, 151–208. doi:10.1016/B978-0-12-800099-1.00003-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Vasimuddin, M., Misra, S., Li, H., and Aluru, S. (2019). Efficient architecture-aware acceleration of BWA-MEM for multicore systems. In: 2019 IEEE international parallel and distributed processing symposium (IPDPS); 2019 May 20-24; Rio de Janeiro, Brazil: IEEE. p. 314–324.

CrossRef Full Text | Google Scholar

Volkman, S. K., Ndiaye, D., Diakite, M., Koita, O. A., Nwakanma, D., Daniels, R. F., et al. (2012). Application of genomics to field investigations of malaria by the international centers of excellence for malaria research. Acta Trop. 121, 324–332. doi:10.1016/j.actatropica.2011.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Waltmann, A., Koepfli, C., Tessier, N., Karl, S., Fola, A., Darcy, A. W., et al. (2018). Increasingly inbred and fragmented populations of Plasmodium vivax associated with the eastward decline in malaria transmission across the southwest Pacific. PLoS Negl. Trop. Dis. 12, e0006146. doi:10.1371/journal.pntd.0006146

PubMed Abstract | CrossRef Full Text | Google Scholar

WHO (2024). World malaria report. Geneva, Switzerland: WHO. Available online at: https://www.who.int/teams/global-malaria-programme/reports/world-malaria-report-2024 (Accessed April 25, 2025).

Google Scholar

Zhong, D., Koepfli, C., Cui, L., and Yan, G. (2018). Molecular approaches to determine the multiplicity of plasmodium infections. Malar. J. 17, 172. doi:10.1186/s12936-018-2322-5

PubMed Abstract | CrossRef Full Text | Google Scholar