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.

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.

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).

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.

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).

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.

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).

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.

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.

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).

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 (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).

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).

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.