Measuring changes in Plasmodium falciparum census population size in response to sequential malaria control interventions

Between 2013 and 2015, three rounds of IRS with non-pyrethroid insecticides were implemented across all of Bongo District (Figure 1A). 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–2014) and the third (2012–2015) round of IRS, respectively (Figure 1B, Appendix 1—table 1). These declines in parasite prevalence were observed across all ages, with the greatest impacts being observed on the younger children (1–5 years) who were ~3 times less likely to have an infection in 2015 (post-IRS) compared to 2012 (pre-IRS) (Figure 1C, Appendix 1—table 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, Appendix 1—table 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, Appendix 1—table 1). This difference by age group in 2017 can be attributed to SMC, which only targets children between 3 and 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 (pre-IRS) (Figure 1C, Appendix 1—table 1).

Study design and changes in the prevalence of microscopic P. falciparum infection following the indoor residual spraying (IRS) and seasonal malaria chemoprevention (SMC) interventions in Bongo, Ghana.

(A) Four age-stratified cross-sectional surveys of ~2000 participants per survey were conducted in Bongo, Ghana, at the end of the wet seasons in October 2012 (Survey 1, baseline pre-IRS, red), October 2014 (Survey 2, during IRS, orange), October 2015 (Survey 3, post-IRS, green), and October 2017 (Survey 4, SMC, purple) (see Materials and methods, Appendix 1—table 1). The three rounds of IRS (grey areas) were implemented between 2013 and 2015 (Tiedje et al., 2022). SMC was distributed to all children <5 years of age during the wet seasons in 2016 (two rounds between August and September 2016) and 2017 (four rounds between September and December 2017) (Gogue et al., 2020). Both IRS and SMC were implemented against a background of widespread long-lasting insecticidal net (LLIN) usage (Tiedje et al., 2022). This figure was adapted from Tiedje et al., 2022, Figure 1 (CC BY 4.0 licence). The copyright holder has granted permission to publish under a CC BY 4.0 licence. Prevalence of microscopic P. falciparum infections (%) in the (B) study population and (C) for all age groups (years) in each survey (Appendix 1—table 1). Error bars represent the upper and lower limits of the 95% confidence interval (CI) calculated using the Wald interval.

Next, we wanted to explore changes in population size measured by MOIvar. 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 pairwise type sharing (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). Somewhat unexpectedly, the change was in the direction of reduced similarity (i.e. less overlap) 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 the measurement of MOIvar, the median PTS scores for non-upsA DBLα types were lower following the IRS intervention (PTSnon-upsA: 2014=0.013 and 2015=0.013 vs. PTSnon-upsA: 2012=0.020). In 2017, the non-upsA PTS distributions shifted back toward higher median PTS scores (PTSnon-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 (MOIvar>1), we also examined isolates with monoclonal infections (MOIvar=1) and found that this non-overlapping structure persisted regardless of infection complexity, particularly for the non-upsA DBLα types (Figure 2—figure supplement 1). 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 PTSnon-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 PTSnon-upsA≥0.5 (i.e. siblings or recent recombinants) (Figure 2).

Sharing of upsA and non-upsA DBLα types among the DBLα isolate repertoires in 2012 (pre-indoor residual spraying [IRS], red), 2014 (during IRS, orange), 2015 (post-IRS, green), and 2017 (seasonal malaria chemoprevention [SMC], purple).

The overlapping density and violin plots (upper right-hand corners) show the distribution of pairwise type sharing (PTS) scores (i.e. DBLα isolate repertoire similarity) between the (A) upsA and (B) non-upsA DBLα isolate repertoires for those isolates with DBLα sequencing data (Appendix 1—tables 2 and 3) in each survey. The PTS scales in the density plots have been zoomed in to provide better visualisation of the upsA and non-upsA DBLα type PTS distributions. The coloured dashed lines in the density plots indicate the median PTS scores in each survey for the upsA (2012 [red]=0.078, 2014 [orange]=0.063, 2015 [green]=0.054, and 2017 [purple]=0.064) and non-upsA (2012 [red]=0.020, 2014 [orange]=0.013, 2015 [green]=0.013, and 2017 [purple]=0.016) DBLα types. Note: The non-upsA median PTS values in 2014 (orange) and 2015 (green) were both 0.013 and overlap in the figure. In the PTS violin plots, the central box plots indicate the medians (centre line), interquartile range (IQR, upper and lower quartiles), whiskers (1.5x IQR), and outliers (points).

The raw data of non-upsA DBLα isolate repertoire sizes were used to estimate MOIvar as adjusted using the Bayesian approach based on pooling the maximum a posteriori MOI estimates (Figure 3, Figure 3—figure supplement 1) and the mixture distribution (Figure 3—figure supplement 2). We observed that at baseline in 2012, the majority (89.2%) of the population across all ages carried multiclonal infections (median MOIvar=4 [interquartile range (IQR): 2–6]) (Figure 3A). Following the IRS intervention, the estimated MOIvar distributions became more positively skewed, indicating that a lower proportion of participants harboured multiclonal infections with a lower median MOIvar in 2014 (64.5%; median MOIvar=2 [IQR: 1–3]) and 2015 (71.4%; median MOIvar=2 [IQR: 1–3]) compared to 2012 (Figure 3A). These reductions in median MOIvar and 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 the median MOIvar (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) predominantly targeted by SMC (Figure 3). While the prevalence of infection in 2017 remained low for the younger children (1–5 years), those infected still carried multiclonal infections (84.1% of those infected) (Figure 3B). Although the MOIvar 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 MOIvar values equal to one or two, while a smaller proportion had MOIvar values ≥5 (Figure 3B).

MOIvar distributions in 2012 (pre-indoor residual spraying [IRS], red), 2014 (during IRS, orange), 2015 (post-IRS, green), and 2017 (seasonal malaria chemoprevention [SMC], purple) based on pooling the maximum a posteriori multiplicity of infection (MOI) estimates.

Estimated MOIvar distributions for the (A) study population and (B) for all age groups (years) in each survey for those isolates with DBLα sequencing data (Appendix 1—tables 2 and 3). The median MOIvar values are indicated with the black dashed lines and have been provided in the top right corner (median MOIvar value [interquartile range (IQR), upper and lower quartiles]) along with the percentage of P. falciparum infections that were multiclonal (MOIvar>1) in each survey and age group (years).

Census population size, measured as the number of P. falciparum var repertoires circulating in the population during each survey, was estimated by summation of isolate MOIvar (see Materials and methods; Figure 4, Appendix 1—table 2). In 2014 during IRS, this number decreased by 71.4% relative to the 2012 baseline survey (pre-IRS) (Figure 4C and E), whereas prevalence decreased by 54.5% (Figure 4C and E). Although census population size increased slightly in 2015 relative to 2014 (Figure 4A), there were still 64.4% fewer var repertoires in the population compared to 2012 (Figure 4C and E) in comparison to a 42.6% decrease in prevalence (Figure 4C and E). Importantly, this loss of var repertoires in 2014 and 2015 following the IRS intervention was seen for all age groups (Figure 4B), with the greatest overall reductions (≥83.8%) being observed for the younger children (1–5 years) (Figure 4D and F). However, in 2017, the number of diverse var repertoires in the population rebounded, more than doubling between 2015 and 2017 (Figure 4C and E). This increase in the number of var repertoires was seen for all age groups in 2017, except for the younger children (1–5 years) where those up to 59 months were targeted by SMC (Figure 4B and D). In fact, the greatest overall increase was observed for the adolescents (11–20 years) and adults (>20 years), where the number of var repertoires in 2017 was ~1.2 times higher compared to 2012 (Figure 4F). Similar trends in the number of var repertoires were also observed for the older children (6–10 years) in 2017, although the rebound was not as striking as that detected for the adolescents and adults.

Estimated number and relative change in the number of P. falciparum var repertoires in 2012 (pre-indoor residual spraying [IRS], red), 2014 (during IRS, orange), 2015 (post-IRS, green), and 2017 (seasonal malaria chemoprevention [SMC], purple).

The estimated number of var repertoires (i.e. census population size) for those isolates with DBLα sequencing data (Appendix 1—tables 2 and 3) in the (A) study population and (B) for all age groups (years). The estimated number of var repertoires vs. P. falciparum prevalence for (C) study population and (D) for all age groups (years) (Appendix 1—table 2). The percentage change in P. falciparum prevalence (black dotted line) and the estimated number of var repertoires (black solid line) in 2014, 2015, and 2017 compared to the 2012 baseline survey (red dashed horizontal line at 0% change) for the (E) study population and (F) for all age groups (years). Error bars in (A–D) represent the upper and lower limits of the 95% confidence intervals (95% CIs). To account for differences in sampling depth across age groups and surveys, we performed subsampling with replacement by selecting the minimum number of individuals in each age group across all surveys. We then calculated the total number of var repertoires from these subsampled individuals within each age group in each survey. This approach ensures consistent sample sizes within each age group across all surveys. Finally, we summed the var repertoires across age groups to obtain the total var repertoire count for each survey. The mean (coloured solid points) and 95% CIs for the number of var repertoires were estimated by repeating the subsampling procedure 10,000 times. The CIs were then derived from the distribution of these repeated subsampling replicates. The 95% CIs for P. falciparum prevalence (%) were calculated using the Wald interval.

As census population size changed considerably during the sequential IRS and SMC interventions, we investigated how the removal or loss of P. falciparum var repertoires and subsequent rebound in 2017 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 = 2218; non-upsA=33,159) (Figure 5, Appendix 1—table 3) and limited overlap of var repertoires (i.e. median PTSnon-upsA≤0.020) seen in a relatively small study population of 685 microscopically positive individuals (Figure 2). In 2014, as P. falciparum prevalence and census population size declined (Figure 4), so too did the number of DBLα types, resulting in a 32.2% and 55.3% reduction in richness for the upsA and non-upsA DBLα types, respectively, compared to 2012 (Figure 5, Appendix 1—table 3). 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 46.0% reduction for upsA and non-upsA DBLα types, respectively) (Figure 5, Appendix 1—table 3). 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 census population size (Figures 4 and 5).

UpsA and non-upsA DBLα type richness in 2012 (pre-indoor residual spraying [IRS], red), 2014 (during IRS, orange), 2015 (post-IRS, green), and 2017 (seasonal malaria chemoprevention [SMC], purple).

Number of unique (A) upsA and (B) non-upsA DBLα types (i.e. richness) observed in each survey vs. P. falciparum prevalence based on those isolates with DBLα sequencing data (Appendix 1—tables 2 and 3). Error bars represent the upper and lower limits of the 95% confidence intervals (95% CIs) for the P. falciparum prevalence (%; x-axis) and ±2 standard deviations (±2 SD) for the number of unique upsA and non-upsA DBLα types (y-axis). The 95% CIs for P. falciparum prevalence (%) were calculated using the Wald interval. The ±2 SD for the number of unique upsA and non-upsA DBLα types was calculated based on a bootstrap approach. We resampled 10,000 replicates from the original population-level distribution with replacement. Each resampled replicate has the same size as the original sample. We then derive the standard deviation (SD) based on the distribution of the resampled replicates.

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 6—figure supplement 1). 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 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 6C and D). Data in Figure 6A–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.

UpsA and non-upsA DBLα type frequencies and survival in 2012 (pre-indoor residual spraying [IRS], red), 2014 (during IRS, orange), 2015 (post-IRS, green), and 2017 (seasonal malaria chemoprevention [SMC], purple).

Heatmaps showing the patterns of diversity for the (A) upsA and (B) non-upsA DBLα types. The columns represent all the upsA and non-upsA DBLα types observed in the four surveys, and the rows represent each of the 2802 upsA DBLα types and the 50,436 non-upsA DBLα types (Appendix 1—table 3). White rows are used to denote the absence of a specific DBLα type, while the presence of a DBLα type is indicated in colour and further categorised (colour gradations) based on the frequency or the number of times (i.e. number of isolates) a DBLα type was observed in each survey (frequency categories: 1, 2–10, 11–20, >20 isolates). Note the frequency category cut-offs were chosen based on the frequency distributions in Figure 6—figure supplement 1. The proportions of (C) upsA and (D) non-upsA DBLα types in each survey based on the number of times (i.e. number of isolates) they were observed in each survey. Kaplan-Meier survival curves for the (E) upsA and (F) non-upsA DBLα types across time (2012–2017) categorised based on their frequency at baseline in 2012 (pre-IRS, red). The coloured shaded areas represent the upper and lower limits of the 95% confidence intervals (95% CIs), with the number (N) of upsA and non-upsA DBLα types in each frequency category provided in parenthesis. These survival curves include only those upsA (N=2218) and non-upsA (N=33,159) DBLα types that were seen at baseline in 2012 (pre-IRS) as indicated in red (Appendix 1—table 3). The x-axis indicates time where time ‘0’ denotes 2012 (pre-IRS), ‘1’ denotes 2014 (during IRS), ‘2’ denotes 2015 (post-IRS), and finally ‘3’ denotes 2017 (SMC). Note: In the survival curves, the 11–20 and >20 frequency categories for both the (E) upsA and (F) non-upsA DBLα types overlap in the figure.

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 in 2012 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 6C and D). Furthermore, when we examined survival using the frequency categories, the upsA and non-upsA DBLα types that were observed at multiple study time points (i.e. 2012, 2014, 2015, and 2017), albeit in different isolate repertoires, were those that were most frequent (i.e. observed in 11–20 and >20 isolates) in the population at baseline in 2012 (Figure 6E and F). 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 6E and F). 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 (Figure 4) where each isolate repertoire is composed of many rare DBLα types as defined by PTS (Figure 2).