Introduction

Malaria in high-transmission endemic areas of sub-Saharan Africa (SSA) is characterised by vast diversity of the Plasmodium falciparum parasites from the perspective of antigenic variation (Chen et al., 2011; Day et al., 2017; Otto et al., 2019; Ruybal-Pesántez et al., 2022, 2017). As with other hosts of hyper-variable pathogens (Futse et al., 2008), children experiencing clinical episodes of malaria, eventually become immune to disease but not to infection. This results in a large reservoir of chronic asymptomatic infections, in hosts of all ages, sustaining transmission to mosquitos. Given the goal of malaria eradication by 2050, it is therefore of interest to examine how the parasite population changes following perturbation by major intervention efforts, both in terms of its size and underlying population genetics.

So, what do we mean by the parasite population size in the case of P. falciparum and how do we measure it? Parasite prevalence, detected by microscopy or more sensitive molecular diagnostics (e.g., PCR), describes the proportion of infected human hosts. Studies of P. falciparum genetic diversity have shown that the majority of people in high-transmission endemic areas harbour diverse multiclonal infections measured as the complexity or multiplicity of infection (MOI) (e.g., (Anderson et al., 2000; Paul et al., 1995; Smith et al., 1999; Sumner et al., 2021)) with complex population dynamics (Bruce et al., 2000; Farnert et al., 1997). These genetic data indicate much larger parasite population sizes than observed by prevalence of infection alone. Thus, from an ecological perspective, we can consider a human host as a patch carrying a number of “antigenically-distinct infections” of P. falciparum. The sum of these antigenically-distinct infections over all hosts provides us with a census of the parasite count of relevance to monitoring and evaluating malaria interventions. We refer to this census population size hereafter simply as population size but make clear that this measure is distinct from effective population size (N_e) as measured by neutral variation.

Diversity of P. falciparum single copy surface antigen genes such as circumsporozoite protein (csp), merozoite surface protein 1 (msp1) or 2 (msp2), and apical membrane antigen 1 (ama1) have each been widely used to measure MOI (e.g., (Falk et al., 2006; Lerch et al., 2017; Nelson et al., 2019)). They have become part of newer genetic panels (e.g., Paragon v1 (Tessema et al., 2020) and AMPLseq v1 (LaVerriere et al., 2022)) specifically for MOI determination. Typically, MOI is reported as the maximum number of alleles or single locus haplotypes present at the most diverse of these antigen-encoding loci rather than the number of unique multilocus haplotypes of these genes combined, as it is challenging to accurately reconstruct or phase these haplotypes in hosts with an MOI > 3 (Lerch et al., 2019). Each of these genes is under balancing selection with a few geographically-common haplotypes and many very rare haplotypes in moderate-to high-transmission settings (Markwalter et al., 2022; Sumner et al., 2021). Where there is a high probability of co-occurrence of two or more common single locus haplotypes in a host, genotyping each of these single copy antigen genes alone will underestimate MOI. Single nucleotide polymorphism (SNP) panels have been used to define the presence of multiclonal infections with limited reliability to estimate MOI for highly complex infections, typical in high transmission, even with the use of computational methods (Labbé et al., 2023).

As an alternative to genotyping single copy antigen genes and biallelic SNP panels to estimate MOI, we have proposed the use of a fingerprinting methodology known as varcoding to genotype the hyper-diverse var multigene family (∼50-60 var genes per haploid genome). This method employs a ∼450bp region of a var gene, known as a DBLα tag encoding the immunogenic Duffy-binding-like alpha (DBLα) domain of PfEMP1 (Plasmodium falciparum erythrocyte membrane protein 1), the major surface antigen of the blood stages (Zhang and Deitsch, 2022). Bioinformatic analyses of a large database of exon 1 sequences of var genes showed a predominantly 1-to-1 DBLα-var relationship, such that each DBLα tag typically represents a unique var gene, especially in high transmission (Tan et al., 2023). The extensive diversity of DBLα tags together with the very low percentage of var genes shared between parasites (Chen et al., 2011; Day et al., 2017; Ruybal-Pesántez et al., 2022, 2017) facilitate measuring MOI by amplifying, pooling, sequencing, and counting the number of unique DBLα tags (or DBLα types) in a host (Ruybal-Pesántez et al., 2022; Tiedje et al., 2022). From a single PCR with degenerate primers and amplicon sequencing, the method specifically counts the most diverse DBLα types, designated non-upsA, per infection to arrive at a metric we call MOI_var. It is not based on assigning haplotypes but exploits the fact that var repertoires are non-overlapping, especially in high transmission. Instead, it assumes a set number of non-upsA types per genome based on repeated sampling of 3D7 control isolates accounting for PCR sampling errors to calculate MOI_var (Ghansah et al., 2023; Ruybal-Pesántez et al., 2022; Tiedje et al., 2022). Consequently, rather than looking at the diversity of a single copy antigen-encoding gene like csp, msp2, or ama1 to calculate MOI, by varcoding we are looking at sets of 45 non-upsA DBLα types per genome. Prior work has shown that varcoding is more sensitive to measure MOI in high transmission where there is an extremely high prevalence of multiclonal infections that cannot be accurately phased with either biallelic SNP panels (Ghansah et al., 2023; Labbé et al., 2023; Tessema et al., 2020; Watson et al., 2021) or combinations of single copy antigen genes (Sumner et al., 2021).

Here we report an investigation of changes in parasite population size and structure through two sequential malaria control interventions between 2012 and 2017 in Bongo District located in the Upper East Region of Ghana, one of the 12 highest burden countries in Africa (World Health Organization, 2022). We present a novel Bayesian modification to the published varcoding approach (Ruybal-Pesántez et al., 2022; Tiedje et al., 2022) that takes into account under-sampling of non-upsA DBLα types in an isolate to estimate MOI_var (Ruybal-Pesántez et al., 2022; Tiedje et al., 2022) and therefore population size. We document P. falciparum prevalence as well as var diversity and population structure from baseline in 2012 through a major perturbation by a short-term indoor residual spraying (IRS) campaign managed under operational conditions, which reduced transmission intensity by > 90% as measured by the entomological inoculation rate (EIR) and decreased parasite prevalence by ∼40-50% (Tiedje et al., 2022). Next, we followed what happened to parasite population size more than two years after the IRS intervention was discontinued and seasonal malaria chemoprevention (SMC) was introduced for children between the ages of 3-59 months (i.e., < 5 years) (Wagman et al., 2018). Detectable changes in parasite population size were seen as a consequence of the IRS intervention but this quantity rapidly rebounded 32-months after the intervention ceased. Overall, throughout the IRS, SMC, and subsequent rebound, the parasite population in humans remained large in size and retained the var population genetic characteristics of high transmission (i.e., high var diversity, low var repertoire overlap), demonstrating the overall resilience of the species to survive significant short-term perturbations.

Materials and methods

Human subject ethical approval

The study was reviewed/approved by the ethics committees at the Navrongo Health Research Centre (Ghana), Noguchi Memorial Institute for Medical Research (Ghana), The University of Melbourne (Australia), and The University of Chicago (United States). Details on the study area, study population, inclusion/exclusion criteria, and data collection procedures have been previously described (Tiedje et al., 2022, 2017).

Study design and sample collection

Using an interrupted time-series study design, four age-stratified cross-sectional surveys of ∼2,000 participants per survey were undertaken to investigate the impacts of IRS and SMC in combination with long-lasting insecticide-treated nets (LLINs) impregnated with pyrethroids under operational conditions on the asymptomatic P. falciparum reservoir from two proximal catchment areas (i.e., Vea/Gowrie and Soe, with a sampling area of ∼60 km²) in Bongo District, Ghana (Tiedje et al., 2022, 2017) (Table supplement 1). Bongo District, located in the Upper East Region, is categorized as high-seasonal malaria transmission based on the World Health Organization’s (WHO) “A Framework for Malaria Elimination” (WHO/GMP, 2017) where P. falciparum prevalence was ≥ 35% at baseline in 2012 (Tiedje et al., 2022, 2017). These ∼2,000 participants of all ages (1-97 years) represent ∼15% of the total population that reside in these two catchment areas in Bongo (Tiedje et al., 2017). The four cross-sectional surveys were completed at the end of the wet season (i.e., high-transmission season) and the study can be separated into four distinct study time points: (1) October 2012 (Survey 1) prior to the IRS and SMC (i.e., baseline), (2) October 2014 (Survey 2) two months after the second round of IRS (Actellic 50EC), (3) October 2015 (Survey 3) seven months after the third round of IRS using a long-acting non-pyrethroid insecticide (Actellic 300CS) (Gogue et al., 2020; US President’s Malaria Initiative Africa IRS (AIRS) Project, 2016), and finally (4) October 2017 (Survey 4) 32-months after the discontinuation of IRS, but during the deployment of SMC to all children 3-59 months (i.e., < 5 years) (Figure 1). Details on the study area, study population, and data collection procedures have been previously described (Tiedje et al., 2022, 2017).

Figure 1.

Details of the IRS and SMC interventions

Starting in 2013, the AngloGold Ashanti Malaria Control Programme (AGAMal) in a public-private partnership with the Global Fund, scaled up IRS across all of the Upper East Region of northern Ghana (Gogue et al., 2020). As part of this initiative, three rounds of IRS with organophosphate formulations (i.e., non-pyrethroid) were rolled out prior to the start of the wet season between 2013 and 2015 (Figure 1A) in Bongo District (Tiedje et al., 2022). Based on AGAMal’s operational reports, IRS coverage in Bongo District was 91.8% in Round 1, 95.6% in Round 2, and finally 96.6% in Round 3 (AGAMal, personal communication). To monitor the impact of the IRS on the local vector population, monthly entomological collections were undertaken between February 2013 and September 2015 (Tiedje et al., 2022). Using these surveys we observed that the monthly entomological inoculation rate (EIR) (infective bites/person/month (ib/p/m)), a measure of local transmission intensity, declined by > 90% at the peak of the wet season between August 2013 (pre-IRS) (EIR = 5.3 ib/p/m) and August 2015 (post-IRS) (EIR = 0.4 ib/p/m) (Tiedje et al., 2022). Following the IRS, SMC was rolled out in the Upper East Region by the National Malaria Elimination Programme (NMEP)/ Ghana Health Service (GHS) starting in 2016 (Gogue et al., 2020) (Figure 1A). SMC is the intermittent administration of full treatment courses of an antimalarial to children between the ages of 3-59 months (i.e., < 5 years) (WHO, 2012). Like other countries, the SMC drug combination of choice in Ghana is amodiaquine plus sulfadoxine-pyrimethamine, which is administered at monthly intervals during the high-transmission season (i.e., wet season) (WHO, 2012). The goal of this age-targeted intervention is to both clear current malaria infections and prevent malarial illness by maintaining a therapeutic concentration of an antimalarial in blood over the period of greatest risk (i.e., high-transmission season). Reported SMC coverage in Bongo District was 92.6% in 2016 (two rounds between August – September 2016) and 94.6% in 2017 (four rounds between September – December 2017) (Gogue et al., 2020).

Var genotyping and sequence analysis

Genomic DNA was extracted from the dried blood spots for all participants with a confirmed microscopic asymptomatic P. falciparum infection (i.e., isolate) (2,572 isolates, Table supplement 2) using the QIAamp DNA mini kit (QIAGEN, USA) as previously described (Tiedje et al., 2017). For var genotyping or varcoding, the sequences encoding the DBLα domains of the antigen-encoding P. falciparum var genes were amplified by PCR, pooled, and sequenced on a MiSeq 2×300 paired-end cycle protocol as previously published (New York University Genome Technology Center, New York, NY, USA; Australian Genome Research Facility, Melbourne, Australia) (Figure supplement 1) (Day et al., 2017; He et al., 2018; Ruybal-Pesántez et al., 2017). The raw sequence data was then cleaned and processed using our previously published customized bioinformatic pipeline (https://github.com/UniMelb-Day-Lab/DBLaCleaner) (He et al., 2018). This pipeline was used to de-multiplex and merge the paired-end reads as well as remove low-quality sequences and chimeras using several filtering parameters (see Figure supplement 2 for additional details). These steps resulted in a total of 291,498 cleaned DBLα sequences for all surveys (Figure supplement 2). To identify the unique DBLα types, we then clustered these cleaned DBLα sequences with 237,910 DBLα sequences previously published from Bongo (He et al., 2018; Pilosof et al., 2019; Rorick et al., 2018) at the standard 96% sequence identity (https://github.com/UniMelb-Day-Lab/clusterDBLalpha) (Barry et al., 2007; Day et al., 2017; Ruybal-Pesántez et al., 2017). Our dataset was then further curated by translating the DBLα types (N = 68,507) into amino acid sequences and removing any DBLα types that could not be translated (i.e., contained a stop codon) (N = 137; 0.2%). The remaining DBLα types were then assigned to their most likely DBLα domain class (i.e., DBLα0, DBLα1, or DBLα2) using a hidden Markov model, and further classified based on the association of specific domain classes with semi-conserved upstream promoter sequences (ups) as either upsA or non-upsA DBLα types (https://github.com/UniMelb-Day-Lab/classifyDBLalpha) (Ruybal-Pesántez et al., 2017). For additional information on the use of these bioinformatic pipelines a detailed tutorial is available (https://github.com/UniMelb-Day-Lab/tutorialDBLalpha). Using a cut off of ≥ 20 DBLα types, DBLα sequencing data was obtained from 2,099 P. falciparum isolates (81.6%) with a total of 289,064 DBLα sequences and 53,241 unique DBLα types being identified in the study population (Table supplement 2). This genotyping success was acceptable given that we were working with low-density asymptomatic infections (Table supplement 1). The median P. falciparum density was ∼4 times higher for isolates with ≥ 20 DBLα types compared to those that gave no or < 20 DBLα types (520 [IQR: 200-1880] parasites/μL vs. 120 [IQR: 40-200] parasites/μL, respectively).

DBLα type diversity

We monitored the impacts of the sequential interventions (i.e., IRS and SMC) on diversity by measuring changes in the population genetics of DBLα types at the population level (i.e., P. falciparum reservoir). Diversity was monitored using two measures, DBLα type richness and DBLα type frequency. Richness was defined as the number of unique DBLα types observed (i.e., DBLα type pool size) in each survey or study time point (i.e., 2012, 2014, 2015, and 2017). DBLα type richness, however, does not provide any information about the relative frequencies of the different DBLα types in the population, as they are all weighted equally whether they are observed once or more frequently (e.g., observed in > 20 isolates per survey). To further examine the impacts of the interventions on DBLα type diversity we also assessed the frequency of each unique DBLα type in 2012, 2014, 2015, and 2017. Here we defined DBLα type frequency as the number of times (i.e., number of isolates) a DBLα type was observed in each survey. Both upsA and non-upsA DBLα type diversities were measured due to their different biological features, chromosomal positions (i.e., subtelomeric regions vs. internal or central regions), as well as population genetics (Falk et al., 2009; Gardner et al., 2002; Jensen et al., 2004; Kaestli et al., 2006; Kalmbach et al., 2010; Kraemer et al., 2007; Kraemer and Smith, 2006; Kyriacou et al., 2006; Lavstsen et al., 2003; Normark et al., 2007; Rottmann et al., 2006; Warimwe et al., 2012, 2009; Zhang and Deitsch, 2022). The proportions of upsA and non-upsA var genes in a repertoire or single genome have been defined as ∼15-20% and ∼80-85%, respectively, based on whole genome sequencing (Rask et al., 2010). The upsA and non-upsA DBLα type proportions were partitioned as expected in our analyses, with the median proportions at the repertoire level being comparable in 2012 (19% upsA and 81% non-upsA), 2014 (22% upsA and 78% non-upsA), 2015 (21% upsA and 79% non-upsA), and 2017 (20% upsA and 80% non-upsA).

Repertoire similarity as defined by pairwise type sharing

To estimate genetic similarity between the DBLα repertoires (i.e., unique DBLα types identified in each isolate) identified from two isolates, pairwise type sharing (PTS) was calculated between all pairs of isolates in each survey as previously described (Barry et al., 2007). PTS, analogous to the Sørensen Index, is a similarity statistic to evaluate the proportion of DBLα types shared between two isolate repertoires (i.e., DBLα repertoire overlap) and ranges from 0 (i.e., no DBLα repertoire overlap) to 1 (i.e., identical DBLα isolate repertoires), where < 0.50 = unrelated, 0.5 = recent recombinants/siblings, > 0.5 = related, and 1 = clones. PTS is a measure of identity-by-state (IBS) used to assess repertoire similarity (or relatedness) between isolates and is not used to infer inheritance from a recent common ancestor (i.e., identity-by-decent (IBD)) (Speed and Balding, 2015).

DBLα isolate repertoire size and number of P. falciparum var repertories

For this study we have exploited the unique population structure of non-overlapping DBLα isolate repertoires to estimate isolate MOI_var. To calculate MOI_var, the non-upsA DBLα types were chosen since not only are they more diverse and less conserved between isolate repertoires (i.e., low median PTS_{non- upsA} ≤ 0.020) compared to the upsA DBLα types, but they have also been shown to have a more specific 1-to-1 relationship with a single var gene than upsA (Tan et al., 2023). The low to non-existent overlap of repertoires enables an estimation of MOI that relies on the number of non-upsA DBLα types sequenced from an individual’s isolate (Ruybal-Pesántez et al., 2022; Tiedje et al., 2022). A constant repertoire size or number of DBLα types in a parasite genome can be used to convert the number of types sequenced in an isolate to its estimated MOI (Ruybal-Pesántez et al., 2022; Tiedje et al., 2022). As such, the approach neglects the measurement error in this size introduced by targeted PCR and amplicon sequencing of var genes in an infection.

Bayesian estimation of MOI_var and associated population size

Here, we extend the method to a Bayesian formulation and estimate the posterior distribution for each sampled individual for the probability of different MOI values. From individual posterior distributions, we can then obtain the estimated MOI frequency distribution for the population as a whole. The two pieces of information required for our approach are the measurement error and the prior distribution of MOI. The measurement error is simply the repertoire size distribution, that is, the distribution of the number of non-upsA DBLα types sequenced given MOI = 1, which is empirically available (Figure supplement 3). We refer to it as P(s | MOI=1) where s denotes repertoire size. In the following, s denotes the number of non-upsA DBLα types sequenced when MOI > 1. The prior distribution of MOI refers to the belief we have for what the actual MOI distribution might look like at the population level before empirical evidence is taken into consideration. For example, the prior distribution of MOI is likely to centre around a higher value in high-transmission endemic areas than in low-transmission ones.

We can obtain P(s | MOI = m) from the serial convolution of the repertoire size distribution P(s | MOI = 1) and P(s | MOI = m - 1). Starting with the repertoire size distribution given a single infection, we can derive P(s | MOI = m) for m equal to 2,3,…, up until a maximum value of 20 (empirically determined), as follows:

where L and U are the lower and upper limit for the repertoire size, 10 and 45 respectively from the empirical repertoire size distribution.

For simplicity, we begin with a uniform prior. We use Bayes’ rule to derive a posterior distribution of MOI given a certain number of non-upsA DBLα types sequenced from an individual:

where k is the maximum value of MOI, here 20, as empirically determined.

To obtain the MOI distribution at the population level, we could either simply pool the maximum a posteriori MOI estimate for each sampled individual, or use a technique called mixture distribution. For the latter, we weighed each posterior MOI distribution for each sampled individual equally and sum over all posterior distributions at the individual level to derive the MOI distribution at the population level:

where n is the number of sampled individuals. These two approaches gave similar results for our empirical survey data as determined by the Kolmogorov-Smirnov Test. The obtained distance statistic is close to 0 and the corresponding p-value is non-significant across all surveys, indicating that the two estimates were drawn from the same distribution (Table supplement 3). Given this similarity, we present hereafter only the result of pooling the maximum a posteriori MOI estimates. Note that we focused on individuals who had a confirmed microscopic asymptomatic P. falciparum infection for our MOI estimation.

To examine alternative priors, we considered empirical MOI distributions described in the literature including the Poisson, hyper-Poisson, and negative binomial distributions (Dietz, 1988; Henry, 2020). The hyper-Poisson and negative binomial distributions can capture the overdispersion seen in the empirical distribution of MOI for certain areas and caused by factors such as heterogeneous biting. We investigated changing the parameters of a negative binomial distribution to generate priors with different means spanning a wide range of MOI values (mean MOI within [∼3, ∼9]), including those seen in high-transmission endemic areas. We found that the final MOI distribution at the population level is fairly robust against the specification of those parameters. In other words, the estimated MOI distributions at the population level for different parameters do not significantly diverge from each other. The pairwise distance between these distributions measured by KL-divergence centres around 0, and the difference between their mean MOI values across different priors is close to 0 (Figure supplement 3 and 4). Additionally, we performed Pearson correlation tests. We found that the correlation coefficients between the estimated MOI distributions at the population level across different parameters were close to 1 with a highly significant p-value (Figure supplement 5 and 6). Thus, we provide in our analyses the estimated population MOI distribution using a uniform prior.

Statistical analysis

We used the R v4.0.2 for all data analyses with the collection of R packages in tidyverse being used for data curation along with base, stats, gtsummary, and epiR for the statistical analyses (R Core Team, 2018; Sjoberg, D et al., 2021; Stevenson, 2020; Wickham et al., 2019). Continuous variables are presented as medians with interquartile ranges (IQRs) and discrete variables are presented using the observed proportions or calculated values with 95% confidence intervals (CIs). Kaplan-Meier survival curves were generated for the time (i.e., number of surveys) to first event (i.e., when the DBLα type was no longer observed/detected) comparing the upsA and non-upsA DBLα types; p-values were determined using the log-rank test using the R packages survival and survminer (Kassambara et al., 2021; Therneau, 2023). The time interval to first event considered for all survival curves was the number of surveys or year (i.e., 2012, 2014, 2015, and 2017) that each DBLα type was observed and only includes those upsA (N = 2,215) and non-upsA (N = 33,151) DBLα types that were seen at baseline in 2012 (i.e., those DBLα types observed prior to the IRS intervention).

Data availability

The sequences utilized in this study are publicly available in GenBank under BioProject Number: PRJNA 396962. All data associated with this study are available in the main text or the supplemental information.

Results

Between 2013 and 2015, three-rounds of IRS with non-pyrethroid insecticides were implemented across all of Bongo District. Coincident with the > 90% decrease in transmission following IRS (Tiedje et al., 2022), the prevalence of microscopic P. falciparum infections compared to the 2012 baseline survey (pre-IRS) declined by 45.2% and 35.7% following the second (2012 to 2014) and the third (2012 to 2015) round of IRS (Figure 1B, Table supplement 1). These declines in parasite prevalence were observed across all ages, with the greatest impacts being observed on the younger children (i.e., 1-5 years) who were ∼3 times less likely to have an infection in 2015 (post-IRS) compared to 2012 (pre-IRS) (Figure 1C, Table supplement 1). These reductions were however short-lived and in 2017, 32-months after the discontinuation of IRS, but during SMC, overall P. falciparum prevalence rebounded to 41.2%, or near pre-IRS levels (Figure 1B, Table supplement 1). Importantly, this increase in the prevalence of infection in 2017 was only observed among the older age groups (i.e., ≥ 6 years) (Figure 1C, Table supplement 1). This difference by age group in 2017 can be attributed to SMC, which only targets children between 3-59 months (i.e., < 5 years). A notable increase in parasite prevalence for adolescents (11-20 years) and adults (>20 years) was found in 2017 relative to 2012 (Figure 1C, Table supplement 1).

Next, we wanted to explore changes in population size measured by MOI_var. As this metric is based on non-overlap of var repertoire diversity of individual isolates, specifically non-upsA DBLα types, we investigated whether DBLα isolate repertoire similarity (or overlap), as measured by PTS, increased following the sequential interventions (i.e., IRS and SMC). Figure 2 shows that median PTS values for both upsA and non-upsA DBLα types remained low in all surveys although the PTS distributions for both groups changed significantly at each of the study time points relative to the 2012 baseline survey (pre-IRS) (p-values < 0.001, Kruskal-Wallis test) (Figure 2). However, the change was in the direction of reduced similarity with lower median PTS scores and a larger number of isolates sharing no DBLα types (i.e., PTS = 0) in 2014, 2015, and 2017 compared to 2012. Relevant to measurement of MOI_var, the median PTS scores for non-upsA DBLα types were lower following the IRS intervention (PTS_non-upsA: 2014 = 0.013 and 2015 = 0.013 vs. PTS_non-upsA: 2012 = 0.020). In 2017, the non-upsA PTS distributions shifted back towards higher median PTS scores (PTS_non-upsA = 0.016) and fewer isolates shared no DBLα types relative to 2014 and 2015 (Figure 2). To verify this pattern was not influenced by multiclonal infections (MOI_var > 1) we also examined isolates with monoclonal infections (MOI_var = 1) and found that this non-overlapping structure persisted regardless of infection complexity, particularly for the non-upsA DBLα types (Figure supplement 7). These PTS data make clear that we were dealing with a large, highly diverse parasite population where 99.9% of the isolate comparisons in all surveys had PTS_non-upsA scores ≤ 0.1 (i.e., shared ≤ 10% of their non-upsA DBLα types), indicating that DBLα isolate repertoires were highly unrelated (Figure 2). In fact, throughout the IRS, SMC, and subsequent rebound, very few DBLα isolate repertoires were observed to be related, with < 0.003% isolate comparisons in each survey having a PTS_non-upsA ≥ 0.5 (i.e., siblings or recent recombinants) (Figure 2)

Figure 2.

The raw data of non-upsA DBLα isolate repertoire sizes (Figure supplement 8) were used to estimate MOI_var as adjusted using the Bayesian approach (Figure 3). We observed that at baseline in 2012, the majority (89.2%) of the population across all ages carried multiclonal infections (median MOI_var= 4 [IQR: 2 – 6]) (Figure 3A). Following the IRS intervention, the estimated MOI_vardistributions became more positively skewed, indicating that a lower proportion of participants harboured multiclonal infections with a lower median MOI_var in 2014 (64.5%; median MOI_var = 2 [IQR: 1 – 3]) and 2015 (71.4%; median MOI_var = 2 [IQR: 1 – 3]) compared to 2012 (Figure 3A). These reductions in median MOI_varand the proportion of multiclonal infections, which were observed across all age groups (Figure 3B), are consistent with the > 90% decrease in transmission intensity following the IRS in turn reducing exposure to new parasite genomes. However, in 2017, both median MOI_var (3 [IQR: 2 – 4]) and the proportion of multiclonal infections (78.9%) rebounded in all age groups, even among the younger children (1-5 years) targeted by SMC (Figure 3). While the prevalence of infection in 2017 remained low for the younger children (1-5 years) targeted by SMC, those infected still carried multiclonal infections (84.1% of those infected) (Figure 3B). Although the MOI_var distributions across all age groups started to rebound in 2017 (i.e., less positively skewed compared to 2014 and 2015) they had not fully recovered to the 2012 baseline patterns (Figure 3). This was most apparent among the younger children (1-5 years), as a larger proportion of isolates in 2017, compared to 2012, had MOI_var values equal to one or two, while a smaller proportion had MOI_var values ≥ 5 (Figure 3B).

Figure 3.

Population size measured as the number of P. falciparum var repertoires circulating in the population during each survey was estimated by summation of isolate MOI_var. In 2014 during IRS, this number decreased by 72.4% relative to the 2012 baseline survey (pre-IRS) (Figure 4A), whereas prevalence decreased by 54.5% (Figure 4B). Although population size increased slightly in 2015 relative to 2014, there were still 64.5% fewer var repertoires in the population compared to 2012 (Figure 4A) in comparison to a 42.6% decrease in prevalence (Figure 4B). Importantly this loss of var repertoires in 2014 and 2015 following the IRS intervention was seen for all age groups, with the greatest overall reductions (> 85%) being observed for the younger children (1-5 years) (Figure 4CD). However, in 2017, the number of diverse var repertoires in the population rebounded, more than doubling between 2015 and 2017 (Figure 4AB). This increase in the number of var repertoires was seen for all age groups in 2017, except for the younger children (1-5 years) who were directly targeted by SMC (Figure 4CD). In fact, the greatest overall increase was observed for the adolescents (11-20 years), where the number of var repertoires in 2017 was ∼1.4 times higher compared to 2012 (Figure 4CD). Similar trends in the number var repertoires were also observed for the older children (6-10 years) and adults (>20 years) in 2017, although the rebound was not as striking as that detected for the adolescents.

Figure 4.

As parasite population size changed considerably during the sequential interventions, we investigated how the removal or loss of P. falciparum var repertoires and subsequent rebound altered DBLα type richness, measured as the number of unique upsA and non-upsA DBLα types in the parasite population in each survey. Richness at baseline in 2012 (pre-IRS) was high with a large number of unique DBLα types (upsA = 2,215; non-upsA = 33,151) (Figure 5, Table supplement 2) and limited overlap of var repertoires (i.e., median PTS_non-upsA ≤ 0.020) seen in a relatively small study population of 685 microscopically positive individuals (Figure 2). In 2014, as P. falciparum prevalence and population size declined (Figure 4) so too did the number of DBLα types, resulting in a 32.1% and 55.3% reduction in richness for the upsA and non-upsA DBLα types, respectively compared to 2012 (Figure 5, Table supplement 2). Again in 2015, as P. falciparum prevalence and population size remained low (Figure 4) DBLα type richness was still less than that observed in 2012 (24.6% and 45.9% reduction for upsA and non-upsA DBLα types, respectively) (Figure 5, Table supplement 2). Finally, in 2017, we found that upsA and non-upsA DBLα type richness rebounded relative to 2014 and 2015, coincident with the increase in P. falciparum prevalence and parasite population size (Figure 4, Figure 5).

Figure 5.

Given this reduction in DBLα type richness following the IRS intervention and subsequent rebound in 2017, we wanted to explore whether the loss of richness was influenced by the frequency of individual DBLα types in the parasite population within and among surveys. To answer this, we defined the relative frequency of individual DBLα types in all isolates in each survey (Figure 6, Figure supplement 9). We discovered that individual upsA and non-upsA DBLα types were not all at equal frequencies within a survey and among surveys. They could be classified as frequent (i.e., observed in 11-20 or > 20 isolates), less frequent (i.e., observed in 2-10 isolates) or only seen once, at baseline in 2012 (Figure 6AB). In 2014 and 2015, following the IRS, there was a significant increase in the proportion of upsA and non-upsA DBLα types in the lower frequency categories (p-value < 0.001, Mann-Whitney U test), with all DBLα types becoming rarer in the population (Figure 6CD). This change can be attributed to the removal of P. falciparum var repertoires (Figure 4) with associated loss of upsA and non-upsA DBLα type richness (Figure 5), which disproportionally affected those DBLα types seen once. This shift to all DBLα types becoming rarer following IRS changed in 2017, where the proportion of DBLα types in the more frequent categories (i.e., 2-10, 11-20, or > 20 isolates) significantly increased while the proportion seen once decreased (p-values < 0.001, Mann-Whitney U tests) (Figure 6CD). Data in Figure 6 (A-D) pointed to a differential effect of the IRS intervention and subsequent rebound on the less frequent upsA and non-upsA DBLα types vs. those that were classified as frequent, where those DBLα types that were most frequent persisted longitudinally.

Figure 6.

To explore this observation further we restricted the longitudinal analysis to those DBLα types from the baseline survey in 2012 (pre-IRS). We compared the probability of survival for the DBLα types identified at baseline and found that the upsA DBLα types persisted significantly longer in the population relative to the non-upsA DBLα types (p < 0.001, log-rank test), despite the IRS intervention. The simple explanation being that although the upsA DBLα types had lower richness (Figure 5), a larger proportion was classified as frequent indicating that multiple copies existed in the population compared to the non-upsA DBLα types (Figure 6CD). Furthermore, when we examined survival using the frequency categories, the upsA and non-upsA DBLα types that persisted and were observed at multiple study time points (i.e., 2012, 2014, 2015, and 2017) were those that were most frequent (i.e., observed in 11-20 and >20 isolates) in the population at baseline in 2012 (Figure 6EF). As expected, the DBLα types that were only observed once in 2012, were significantly less likely to be seen longitudinally (p-value < 0.001, log-rank test) (Figure 6EF). These differential changes in DBLα type richness with respect to rare vs frequent DBLα types are a consequence of changes in census population size with interventions where each isolate repertoire is composed of many rare DBLα types as defined by PTS (Figure 2).

Discussion

P. falciparum populations in high-transmission endemic areas in SSA, are characterised by extensive diversity, high rates of recombination, as well as frequent multiclonal infections. Here we defined census population size of P. falciparum to understand total parasite diversity in a human population and explore the age-specific efficacy of malaria interventions to reduce this metric in such areas as typified by Bongo, Ghana. Census population size proved more informative than parasite prevalence alone because it captures “within” host parasite population size as MOI, rather than using the infected host per se as a unit of population size. Whilst the concept of census population size is agnostic of how you measure MOI, the extensive DBLα isolate repertoire diversity presented makes a strong case for fingerprinting parasite isolates by varcoding in high transmission. This is opposed to looking at allelic diversity of a single copy antigen gene such as csp that is under balancing selection.

By varcoding we identified a very large census parasite population size in a relatively small human population of ∼2,000 individuals at baseline and captured age-specific changes in this metric in response to sequential malaria control interventions. IRS reduced the MOI_var parasite population size substantially with the greatest reductions (85%) seen in the younger children (i.e., 1-5 years). More than two years after cessation of IRS, rebound in 2017 was rapid in all age groups with the exception of younger children targeted by SMC. Population sizes in adolescents and adults showed they carried more infections in 2017 than at baseline in 2012. This is indicative of a loss of immunity during IRS which may relate to the observed loss of var richness, especially the many rare types. This warrants further investigation of changes in variant-specific immunity. During and following the IRS and SMC interventions, var diversity remained high and var repertoire overlap remained low, reflecting characteristic properties of high transmission and demonstrating the overall resilience of the species to survive significant short-term perturbations. Combining interventions and targeting older age groups or the whole community with chemoprevention would no doubt have a much greater impact on reducing the diversity of the reservoir of infection.

What was striking about the Bongo study was the speed with which rebound in MOI_var per person and census population size occurred, once the short-term IRS was discontinued. We looked for a potential explanation in our genetic data. PTS and population frequency data showed that many of the DBLα types occurred in multiple repertoires or genomes. This enabled the survival of these more frequent DBLα types through the interventions, facilitating rebound by maintenance of this diversity. The other notable population genetic result of our study was the failure to increase similarity (or relatedness by state) of var repertoires by reducing transmission by > 90% via IRS. From a baseline of a very large population size with very low overlap in repertoires you need outcrossing to create relatedness. However, this was less likely to happen due to reduced transmission as a result of IRS. Rebound, with associated increases in transmission, led to a small increase in var repertoire similarity. This is the opposite to what has been observed in areas of lower transmission under intense malaria control where the intensity of interventions led to increased genome similarity, as assessed by IBD, from a starting point of much lower genome diversity and greater relatedness (Daniels et al., 2015).

Our molecular approach to measure population size has been to sum MOI_var in individual hosts with microscopically-detectable infections. Like any diagnostic method there are limits to sensitivity and specificity, which can be more or less tolerated dependent upon the purpose of the study. Here we have looked at relative changes in population size with sequential interventions using an interrupted time series study design and observed changes by measuring MOI_var. We have accounted for missing DBLα type data where complete var repertoires may not have been sequenced using a Bayesian method based on empirical knowledge of the measurement error. This approach has conceptual relation to the Bayesian approach by Johnson and Larremore (Johnson and Larremore, 2022) to estimate complete repertoire size of, and overlap between, monoclonal infections from incomplete sampling of DBLα types. Our measurement of population size based on MOI_var will be subject to other sampling errors which may in the end be more significant (discussed in detail in Labbé et al. (Labbé et al., 2023)). For example, low parasitemia typical of asymptomatic infections, small blood volumes, clinical status, and/or within host dynamics, including synchronicity, will all create sampling problems but these are common to all measures of MOI.

Previously, we have drawn attention to underestimation of the number of DBLα types of related parasites from a cross by var genotyping (Labbé et al., 2023). Such related parasites must be created frequently in high transmission due to extensive outcrossing (Babiker et al., 1994; Paul et al., 1995). Single clone genomics experiments using biallelic SNPs from whole genome sequencing data have also detected related parasites using IBD in clinical infections from humans in a high-transmission area of Malawi (Nkhoma et al., 2020). Here we have analysed only low density, chronic asymptomatic infections under strong immune selection in semi-immune hosts whom we have shown select against related parasites (He et al., 2018; Ruybal-Pesántez et al., 2022). When considering the importance of possible exclusion of parasites related by descent, the only sure way to detect such parasites in high transmission is by single cell genomics, a methodology of limited application to malaria surveillance due to practicality and cost of scale up. Again, the error from failure to sample related infections must be weighed up against the issues of under-sampling as described above.

In summary, our findings provide parasite population insights into why rebound is the inevitable consequence of such short-term IRS interventions unless you simultaneously target the highly diverse, long-lived parasite population in humans, not just children < 5 years by SMC. Of potential translational significance for malaria molecular surveillance, we identify new metrics, especially MOI_var and census population size as well as var frequency category, informative to monitor and evaluate interventions in high-transmission areas with multiclonal infections and high rates of outcrossing. Such metrics could be used longitudinally to detect incremental gains of transmission-reducing interventions, including IRS, LLINs, and vaccines to perturb the high-transmission characteristics of the parasite population in humans in high-burden countries in SSA.

Supporting information

Supplemental Information

Acknowledgements

This research was supported by Fogarty International Center at the National Institutes of Health through the joint NIH-NSF-NIFA Ecology and Evolution of Infectious Diseases award R01-TW009670 to K.A.K, M.P, and K.P.D; and the National Institute of Allergy and Infectious Diseases, National Institutes of Health through the joint NIH-NSF-NIFA Ecology and Evolution of Infectious Diseases award R01-AI149779 to A.R.O, K.A.K, M.P, and K.P.D. We wish to thank the participants, communities, and the Ghana Health Service in Bongo District, Ghana for their willingness to participate in this study. We would like to thank the field teams in Bongo for their technical assistance in the field, as well as the laboratory personnel at the Navrongo Health Research Centre for their expertise and for undertaking the sample collections and parasitological assessments.

Data availability statement

The sequences utilized in this study are publicly available in GenBank under BioProject Number: PRJNA 396962. All data associated with this study are available in the main text or in the Supporting Information. All custom code is available in an open source repository: (1) DBLαCleaner pipeline is available athttps://github.com/UniMelb-Day-Lab/DBLaCleaner, (2) clusterDBLalpha pipeline is available at https://github.com/Unimelb-Day-Lab/clusterDBLalpha, and the (3) classifyDBLalpha pipeline is available at https://github.com/Unimelb-Day-Lab/classifyDBLalpha. A dataset to demo this custom code is available at https://github.com/UniMelb-Day-Lab/tutorialDBLalpha. For additional information on the use of the Bayesian approach to estimate MOI_var please see https://github.com/qzhan321/Bayesian-formulation-varcoding-MOI-estimation.

Benefit-sharing statement

A research collaboration was developed with scientists from Ghana based at the Navrongo Health Research Centre and the Noguchi Memorial Institute for Medical Research. All collaborators are included as co-authors and the relevant results from the research has been shared with key stakeholders and the local community (i.e., Paramount Chief of Bongo, divisional Chiefs, Queen Mothers, and community members) as well as the Bongo District Health Directorate and the Regional Health Directorate (Upper East Region). Before this research was undertaken, informed consent was sought and obtained from the key stakeholders and the local community in Bongo District. In addition, members of the local community were trained as field workers and were directly involved in liaising with the local community and in the collection of the study data. The contribution of these individuals to this research is described in the Acknowledgements. This research addresses a priority concern regarding malaria control and the impact of interventions. These concerns are relevant for both the local community as well as for the National Malaria Elimination Programme in Ghana.

Author contributions

K.P.D and M.P conceptualized the study of population size using MOI_var. K.A.K and K.P.D conceived and designed the field study with input from A.R.O and K.E.T, who together with O.B coordinated the field studies. K.E.T, S.R-P, and S.L.D processed the samples and performed the genotyping experiments. G.T-H developed and validated the bioinformatic pipelines. K.E.T, S.R-P, and D.C.A processed, cleaned, and curated the datasets for analysis. Q.Z developed and implemented the Bayesian estimations. K.E.T, Q.Z, M.P, and K.P.D analysed the data with contributions from S.R-P, Q.H, and M-H.T. K.E.T and Q.Z visualized the data. K.P.D, M.P, and K.E.T wrote the original draft of the manuscript. Q.Z., S.R-P, G.T-H, Q.H, M-H.T, D.C.A, S.L.D, O. B, A.R.O, and K.AK reviewed and edited manuscript. M.P and K.P.D supervised the research. A.R.O, K.A.K, M.P, and K.P.D acquired the funding.

Competing interests

The authors declare no competing interests.