Samples and epidemiological data from two cohorts were used, being the Kilifi malaria longitudinal cohort [KMLC study (Muthui et al., 2019b)] and the assessment of the infectious reservoir of malaria [AFIRM study (Gonçalves et al., 2017)]. The KMLC cohort comprised three sub-cohorts of children followed up longitudinally and sampled at cross-sectional surveys to assess asymptomatic P. falciparum infections. The AFIRM cohort was a cross-section sampling carried out in the wet and dry seasons and comprised children and adults. A breakdown of the cohorts is provided in Table 1 .

Candidate antigens for the study were identified from a published dataset of the gametocyte proteome (Lasonder et al., 2016). From an initial list of 2,241 proteins, we shortlisted 24 proteins for further analysis. These proteins were shortlisted based on features suggestive of surface localization (signal peptides, transmembrane domains and glycosylphosphatidylinositol anchors). An additional antigen with potential association with naturally acquired transmission reducing immunity was identified from a conference abstract (Stone et al., 2015) to give a total of 25 proteins ( Supplementary Table 1 ). At the time of the search, the candidate antigens were predominantly uncharacterized as targets of naturally acquired immunity to mature gametocyte antigens. Sequence variation between the reference lab strain P. falciparum 3D7 (genome version 3.0) and a fully sequenced Kilifi clinical isolate – PfKE04 – was also analyzed by pairwise alignment using Geneious bioinformatics software (version 11.1.2). Where variation was identified [either an insertion/deletion or non-synonymous single nucleotide polymorphism (SNP)], both variants were prioritized for construct design.

Protein production was attempted in the wheat germ cell-free expression system (WGCFS, WEPRO7240H, Cell-Free Sciences) and, where unsuccessful in the WGCFS, in the mammalian HEK293E system. Of the 25 antigens targeted for recombinant protein production, seven antigens (including one variant antigen) were successfully expressed in sufficient quantities for this analysis. These were Pfs230-C (PF3D7_0209000), PF3D7_0208800, MDV1 (PF3D7_1216500), G377 (PF3D7_1250100, 3D7 and PfKE04 variants expressed), PF3D7_1314500, PF3D7_0303900 and PSOP1 (PF3D7_0721700)) giving a total of 7 candidate antigens (8 proteins) for immunoprofiling analysis. Notably, owing to the size of G377 protein that would be a challenge for protein production, a domain spanning amino acids 666-1146 termed ‘domain B’ was produced for analysis. This domain was identified from work by Alano et al. (1995), who were able to produce protein from this domain successfully (Alano et al., 1995). An extensive description of the protein identification and the recombinant protein production process is provided in the Supplementary Methods , and results from the protein production are provided in Supplementary Figures 2–5 . Additionally, we evaluated antibody responses to apical membrane antigen-1 protein (AMA1, PF3D7_1133400), a highly immunogenic asexual stage antigen (Drakeley et al., 2005; Jones et al., 2015; Idris et al., 2017) widely studied in the context of seroprevalence to place our results in context, and to a crude stage V gametocyte extract (GE).

The measurement of antibody responses to the candidate antigens, AMA1 and GE, was carried out using a previously described three-day ELISA protocol (Murungi et al., 2019; Omondi et al., 2021) with a few modifications. Additionally, to determine the optimal antigen coating concentration and serum dilution to use for the assay, checkerboard ELISAs were first carried out ( Supplementary Figures 6, 7 ). On day one of the assay, 100 µL of the purified recombinant protein was prepared in coating buffer (at a concentration of 1 μg/ml for the gametocyte antigens, 0.5 μg/ml for AMA1 and 1 in 250 dilution for GE) and coated onto a 96-well plate (Immulon 4 HBX, ThermoFisher Scientific). The plate was then incubated at 4°C overnight, and on the following day, the plate was washed four times with phosphate-buffered saline containing 0.05% tween-20 (Pierce™ 20X PBS Tween™ 20 Buffer, ThermoFisher Scientific). The plate was then blocked with 200 µL of blocking buffer for five hours at room temperature and subsequently washed three times with PBS/T. One hundred microlitres of test sera were then added at dilutions of 1 in 200 for the gametocyte antigens (or 1 in 100 for lowly reactive antigens - i.e., if average optical densities were below 0.15); 1 in 1000 (or 1 in 2000 for adult sera tested) for AMA1, and 1 in 500 for GE). A standardized ELISA format was used where defined concentrations of malaria immunoglobulin (MIG) were serially diluted to generate a standard curve against which antibody units were extrapolated.

The MIG was prepared from purified immunoglobulin G (98% IgG) obtained from a pool of 834 Malawian adult sera (Taylor et al., 1992). For positive controls, a characterized pool of hyper-immune sera (PHIS) was obtained from a random selection of adult residents of Junju location, Kilifi (Osier et al., 2014; Murungi et al., 2016), and a pool of sera from gametocyte positive individuals from the AFIRM cohort with high gametocytaemia were also included. Malaria naive sera from European adults (non-immune sera, NIS) were used as negative controls at the same dilution as the test serum for each respective antigen. The NIS were tested against the measles antigen to verify that they were indeed reactive ( Supplementary Figure 8 ). On day three, the plate was washed three times before 100 µL of secondary antibody (polyclonal rabbit anti-human IgG-HRP, Dako, UK) was added (1: 5,000 dilution) and the plate incubated at room temperature for three hours. The plate was washed four times before adding 100 µL of o-phenylenediamine dihydrochloride (OPD) substrate (Sigma, UK). Finally, the plate was incubated at room temperature for 15 minutes to allow colour development before the reaction was stopped using 25 µL of 2 M sulphuric acid (H₂SO₄). Each sample’s optical density (OD) was then determined by reading the absorbance at 492 nM. All samples were run in duplicate, with replicate samples run on a separate plate. As all samples were run in duplicate in the ELISA, the mean OD and co-efficient of variation were calculated before analysis. Samples with a coefficient of variation higher than 20% were repeated.

We summarised the demographic characteristics of participants for the two cohorts separately. We did this by sub-cohort (Ngerenya-early, Ngerenya-late and Junju) within the KMLC cohort and by season (wet and dry) within the AFIRM cohort. For continuous variables, we presented their means (standard deviation) or median (interquartile range, IQR) depending on the distribution. For categorical variables, we presented the percentages with comparisons carried out using the Chi-square test. We also summarised parasite density in the different cohorts by age group, sickle cell status, α–thalassaemia status, cohort and season. Comparisons between groups were performed using the Wilcoxon rank-sum test (where a categorical variable had only two groups) or Kruskal-Wallis test (for variables with more than two groups) with posthoc analysis done using the Dunn’s test with Bonferroni correction.

Next, we estimated the seroprevalence. We note that the analysis of seroprevalence was limited to the AFIRM cohort only, where sampling cut across all age groups. Antibody responses (IgG) were measured to the candidate antigens that were successfully produced as recombinant proteins. Responses were also measured to two markers of parasite exposure: Apical membrane antigen 1 (AMA1) (Drakeley et al., 2005; Jones et al., 2015; Idris et al., 2017) and gametocyte extract (GE) (Omondi et al., 2021). A four-parameter hyperbolic standard curve was generated from the titrated MIG, and the antibody concentration of each sample was then extrapolated. A cut-off was defined as the mean antibody concentration of the NIS plus three standard deviations (3SD) to determine seropositivity. Additionally, to compare the different methods of assigning seropositivity, we also calculated seropositivity using a mixed model with the distribution of OD values fitted as the sum of two Gaussian distributions (a narrow distribution of low responders and a broader distribution of high responders) (Bousema et al., 2010). We then estimated seropositivity from a cut-off defined as the mean OD of the low responders plus 3SD. These seroprevalences were compared across variant antigens and across sampling seasons by means of a Chi-square test. The Cochrane-Armitage test was used to test for a linear trend in seroprevalence when grouped by age. Additionally, we also evaluated the dependency of seropositivity to the different antigens using a correlation matrix, with correlation analyzed using Spearman’s rank correlation.

Furthermore, we examined how determinants of parasite exposure such as age, parasite prevalence, season (AFIRM cohort), transmission intensity (KMLC cohort), and haemoglobinopathies relate to the magnitude of antibody response to the candidate antigens. These covariates were chosen owing to their associations with gametocytaemia or anti-gametocyte antibody responses (Boudin et al., 1991; Lawaly et al., 2010; Grange et al., 2015; Skinner et al., 2015; Amoah et al., 2018). Factors associated with the magnitude and breadth of antibody response to the gametocyte antigens were identified using linear regression models, and we fitted both univariable and multivariable models. Cluster-robust standard errors were calculated with clustering specified to occur by participant ID, to account for instances of repeated sampling in the KMLC cohort. For the AFIRM cohort, parasitaemia and gametocytaemia were assessed by both microscopy and the highly sensitive nucleic acid sequence-based amplification; thus, it was also possible to compare the effect of patent versus sub-patent parasitaemia and gametocytaemia in a sub-analysis.

In addition to assessing the factors associated with the magnitude of anti-gametocyte responses, we also assessed the parameters associated with the recognition of a greater number of gametocyte antigens (breadth of response). For this, we considered how the pre-specified covariates related to the breadth of antibody responses to the candidate antigens using linear regression models. All analyses were performed using R version 4.0.3 (R Core Team, 2020).

Prior to evaluating the relationships between the magnitude of antibody responses to the candidate antigens and the covariates age, transmission intensity (sub-cohorts), season and sickle cell and α-thalassaemia genotypes, we first analyzed how these covariates relate to parasite densities. This was done to evaluate whether the influence of these covariates on parasite densities relates to their impact on the magnitude of antibody responses to the candidate antigens. For the KMLC, microscopically determined asexual parasite densities were higher in children aged between 0 – 5 years of age compared to those aged 6 – 5 years (p = 0.019, Figure 1A ). In contrast, gametocyte densities were similar across age groups. There were no differences in both asexual and gametocyte densities among the different sickle ( Figure 1B ) and α-thalassaemia genotypes ( Figure 1C ) or across the subcohorts ( Figure 1D ). As few participants had patent gametocytaemia in the AFIRM cohort ( Table 3 ), we compared sub-patent parasite densities by the covariates mentioned above. In this case, sub-patent gametocyte densities did not vary by age ( Figure 2A ) or by sickle trait ( Figure 2B ) and α-thalassaemia genotype ( Figure 2C ). However, parasite densities were higher in the wet season compared to the dry season, and this was statistically significant for both all parasites (p = 0.015) and gametocytes (p = 0.013, Figure 2D ).

For seropositivity estimated from a population of malaria naïve individuals, there was a broad range in the seroprevalence estimates, ranging from 37.0% (95% CI: 30.6%, 43.9%) for PSOP1 to 77.8% (95% CI: 71.6%, 83.1%) for G377B 3D7 variant ( Table 4 ). Seroprevalence to G377B PfKE04 was similarly high (70.4%, 95% CI: 63. 8%, 76.4%), and there was a significant association between participants who responded to G377B 3D7 and those who responded to G377B PfKE04 (χ² statistic = 110, p <0.0001; Chi-square test of independence with Yates’ continuity correction, 2 x 2 tables provided in Supplementary Table 2 ). Seroprevalence to Pfs230-C was estimated as 64.2% (95% CI: 55.9%, 71.9%), with responses to both AMA1 and GE being the most prevalent at 87.9% (95% CI: 82.8%, 91.9%) and 88.3% (95% CI: 83.1%, 92.3%) respectively.

When seroprevalence was stratified by age, with age categorized as 0 – 5 years, 6 – 15 years and >15 years, there was a general trend towards increasing seroprevalence with age for all antigens being statistically significant for all except PF3D7_0208800 (p = 0.064, Figure 3 ). There was no evidence for a difference in seroprevalence between the dry and wet seasons for any of the antigens ( Figure 4 ). When seroprevalence was estimated using a subset of low responders as the seronegative population, the estimated seroprevalence was lower, particularly for the gametocyte antigens, likely reflecting the low ODs in the seropositive samples ( Supplementary Table 3 ). Despite this, the trend towards increasing seroprevalence with age remained significant for all antigens, except for PEB-P ( Supplementary Figure 9 ). Seroprevalence did not associate with season same as before ( Supplementary Figure 10 ).

We also examined the dependency between antibody responses to the candidate gametocyte antigens, AMA1 and GE using a correlation matrix ( Supplementary Figure 11 and Supplementary Table 4 ). The strongest correlation was observed between seropositivity to the two G377B variants (Spearman’s correlation ρ = 0.75, p <0.0001). Furthermore, there was a trend towards a stronger correlation between seropositivity to the different gametocyte antigens than between seropositivity to the gametocyte antigens and positivity to either AMA1 or GE. While seropositivity to AMA1 and GE did not strongly correlate with seropositivity to the gametocyte antigens, their strongest correlation was with each other (Spearman’s correlation ρ = 0.51, p <0.0001).

Detailed results from the univariable and multivariable analysis of predictors of antibody levels are provided in Supplementary Tables 5, 6 , with the results of the multivariable analyses discussed below.

Statistically significant associations between increased age and increased magnitude of antibody responses to AMA1, Pfs230-C, G377B 3D7 and G377B PfKE04 were observed. Furthermore, patent asexual parasitaemia was significantly associated with an increased magnitude of antibody responses to all antigens. Gametocytaemia, on the other hand, was only significantly associated with increased responses to AMA1, GE, PF3D7_1314500, PF3D7_0303900, and PF3D7_0208800 (p <0.05). No significant associations were observed between sickle cell trait and antibody levels for any of the antigens. Conversely, being heterozygous (estimate -0.24, 95% CI: -0.45, -0.03, p = 0.022) or homozygous (estimate -0.30, 95% CI: -0.57, -0.04, p = 0.025) for α-thalassaemia was associated with reduced responses to GE. Generally, transmission intensity did not appear to predict the magnitude of antibody response to the gametocyte antigens independent of parasite prevalence, except in the case of GE and G377B 3D7. Relative to the Junju sub-cohort (moderate transmission intensity), residing in the low transmission sub-cohort (Ngerenya-late) was associated with reduced immune responses to GE (estimate -0.30, 95% CI: -0.54, -0.06, p = 0.015) and G377B 3D7 (estimate -0.12, 95% CI: -0.23, -0.004, p = 0.043).

Both G377B 3D7 and G377B PfKE04 had similar associations between increasing age and increased magnitude of antibody responses, and increased antibody levels with concurrent parasitaemia. Additionally, R² values of the predictive models for both variants were similar (G377B 3D7 R² = 0.36, G377B PfKE04 R² = 0.36).

Detailed results of the univariable and multivariable analyses are provided in Supplementary Tables 7, 8 , with a summary of the results of the multivariable analyses provided below.

There was a statistically significant trend towards increasing magnitude of antibody responses with increasing age for Pfs230-C, AMA1 and GE, with only older age (>15 years of age) associated with an increased antibody response (p <0.05) to all other antigens except for PF3D7_0208800. Similar to the KMLC, parasitaemia was an important predictor of the magnitude of antibody response to the gametocyte antigens, with the exception of PSOP1. Both patent and sub-patent parasitaemia were associated with an increased magnitude of antibody responses to PF3D7_1314500, PF3D7_0303900, MDV1, G377B (both variants) and GE. To explore whether patent parasitaemia better predicted the magnitude of antibody responses to these antigens, a second analysis limited to positive parasite individuals was carried out ( Supplementary Tables 9, 10 ). There was no indication that patent parasitaemia was a better predictor of the magnitude of antibody responses relative to sub-patent parasitaemia in this sub-analysis. Moreover, in the adjusted analysis, gametocytaemia did not independently predict the magnitude of antibody responses to the gametocyte antigens except for Pfs230-C (estimate 0.23, 95% CI: 0.02 – 0.44, p = 0.030).

Unlike in the KMLC, where no associations with sickle cell were seen, sickle cell trait was associated with a reduced magnitude of antibody response to PF3D7_0303900 in both the univariable and multivariable analysis (multivariable analysis: estimate -0.24, 95% CI: -0.42 – 0.05, p = 0.013). No associations among the α-thalassemia genotypes and antibody levels were seen. Similarly, no association was observed between sampling season and antibody levels for any of the antigens.

Similar to the KMLC, both variants of G377B had associations between increasing age and increased magnitude of antibody responses and between and concurrent parasitaemia and increased level of antibody responses. Additionally, the R² values of models predicting the magnitude of antibody responses to either antigen varied marginally (G377B 3D7 R² = 0.34, G377B PfKE04 R² = 0.36).

The mean (standard deviation) breadth of response to the eight gametocyte proteins was 4.5 (2.5) in the KMLC cohort and 5 (2.6) in the AFIRM cohort ( Figure 5 ).

Adjusted linear regression analysis of factors associated with increased breadth indicated an association between increasing age and increased antibody breadth for both cohorts ( Tables 5 and 6 ). Additionally, concurrent parasitaemia, whether asexual parasites or gametocytes, was also associated with increased breadth. No statistically significant associations with sickle and α-thalassemia genotypes, transmission intensity or season were observed.