Introduction

Plasmodium falciparum malaria remains a major global health burden, causing over half a million deaths annually, primarily in sub-Saharan Africa (WHO, 2023). The clinical spectrum of malaria ranges from non-severe to severe disease, with severe cases characterized by life-threatening complications such as impaired consciousness (IC), severe malarial anaemia (SMA), and respiratory distress (RD) (Marsh et al., 1995). A deeper understanding of host-parasite interactions at the molecular level that drive distinct clinical malaria syndromes is essential for the development of targeted therapies. A key unresolved question in malaria research is whether parasites causing severe disease biologically differ from those responsible for non-severe malaria, and whether specific severe manifestations can be attributed to particular parasite features.

Central to the pathogenesis of severe malaria is the adherence of infected red blood cells (iRBCs) to endothelial cells in the microvasculature, a process known as cytoadhesion. Cytoadhesion facilitates immune evasion by preventing spleen-mediated clearance of iRBCs (Rowe et al., 2009). The primary parasite ligands mediating this interaction are the Plasmodium falciparum erythrocyte membrane protein 1 (PfEMP1) family. PfEMP1 proteins mediate cytoadhesion by binding to endothelial receptors such as endothelial protein C receptor (EPCR), intercellular adhesion molecule-1 (ICAM-1), cluster of differentiation 36 (CD36), complement receptor 1 (CR1), and chondroitin sulfate A (CSA) (Bernabeu & Smith, 2016; Smith et al., 2013; Tuikue Ndam et al., 2017), with different binding specificities linked to distinct severe malaria phenotypes. For instance, the binding of iRBCs expressing VAR2CSA to CSA on placental syncytiotrophoblasts contributes significantly to placental malaria (Clausen et al., 2012). Additionally, rosetting, in which iRBCs bind uninfected RBCs, plays a role in severe malaria pathogenesis (Cockburn et al., 2003; Lee et al., 2022).

PfEMP1 is a large 200 to 400-kDa clonally variant antigen encoded by the highly polymorphic var gene family (Gardner et al., 2002; Smith et al., 1995). Each P. falciparum genome carries approximately 60 var genes that are transcribed in a mutually exclusive fashion. This allows the parasite to alter antigenic profiles and evade host immunity (Deitsch & Dzikowski, 2017). PfEMP1 typically includes an N-terminal “head structure” with an N-terminal sequence (NTS) and a Duffy Binding-Like (DBL) domain, followed by a Cysteine-rich InterDomain Region (CIDR) (Smith et al., 2000). Var genes are classified into three major groups (A, B, and C) based on their chromosomal location and promoter sequence (Smith, 2014). This classification is clinically relevant, as group A var genes are typically associated with severe malaria, whereas group B and C are more frequently expressed during mild or asymptomatic infections (Kyriacou et al., 2006). Additionally, var genes are grouped according to conserved domain arrangements known as domain cassettes (DC) (Rask et al., 2010), with parasites expressing DC8 and DC13 strongly correlated with severe disease presentations (Lavstsen et al., 2012). Further categorization is based on cysteine residue counts and conserved amino acid motifs (“MFK” and “REY”) in positions of limited sequence variation within DBL alpha domains referred to as cp grouping, with six distinct subgroups (cp1-cp6) (Table 1) (Bull et al., 2005, 2007). Of these, cp1 has specifically been linked to IC, whereas cp2 is predominantly associated with SMA (Warimwe et al., 2009).

Table 1:

These findings suggest that parasites expressing certain PfEMP1 variants are more virulent. Despite these advances, previous studies have been limited by small sample sizes and uneven representation of clinical syndromes, which constrain statistical power and generalizability. Kilifi, a coastal region in Kenya has experienced major shifts in transmission and in the clinical spectrum of malaria in the past two decades (Mogeni et al., 2016; Njuguna et al., 2019). This setting provides an opportunity to revisit these questions across well-represented syndromes using samples collected over 2 decades.

Furthermore, previous studies have predominantly relied on either classification of PfEMP1 based on segments derived through quantitative real-time polymerase chain reaction (RT-qPCR) or sequencing of expressed DBLα-tags in parasite isolates (Andisi & Abdi, 2022; Bull et al., 2005; Lavstsen et al., 2012; Mkumbaye et al., 2017). Each approach has strengths and limitations. RT-qPCR offers quantitative resolution for predefined targets but cannot capture the full var repertoire. DBLα-tag sequencing on the other hand, captures var repertoire breadth and diversity but lacks precise quantitation of individual transcripts. Reliance on a single method has likely contributed to inconsistent associations with clinical phenotypes reported in previous studies.

To overcome these limitations, this study used clinical data, and parasite isolates from 712 children from Kilifi recruited between 1994-2012 with various clinical malaria presentations. The var gene expression profiles of the parasite isolates was determined using both DBLα-tag sequencing and RT-qPCR. By using these complementary datasets, our aim was to refine and robustly test associations between PfEMP1 variants and severe malaria syndromes.

Results

Participants and sample characteristics

Samples from 712 children (49% females) aged between 1-167 months who were slide positive for P. falciparum infection were included in this study (Figure 1A). 229 (32%), 252 (35%), and 231(32%) children were recruited between 1994-2002, 2003-2008 and 2009-2012 respectively corresponding to periods of high (pre-decline), rapidly declining (decline) and low malaria transmission (post-decline) respectively in Kilifi County, Kenya (Figure S1) (Mogeni et al., 2016). 340 of the participants had non-severe malaria and 372 had severe malaria that included 186 impaired consciousness (IC), 29 Respiratory distress (RD), 59 severe malarial anaemia (SMA), and 98 with overlapping severe malaria presentations (Figure 1B).

Figure 1:

As expected, the mean age of children included in this study at the time of admission increased as malaria transmission declined (Figure S1). Among these children, those admitted with any form of severe malaria syndrome were significantly younger (median 35 vs 44 months, p < 0.001) and mortality was significantly higher among severe cases (p < 0.001) (Table S2). By syndrome, IC and SMA were significantly associated with younger age compared to NS (IC p = 0.016; SMA p < 0.001), whereas RD wasn’t (Table S2). Peripheral blood parasitemia was only significantly elevated in RD/IC cases (p=0.034). Mortality was highest in IC (p < 0.001) with no significant difference in RD or SMA relative to NS (Table S2).

Impaired consciousness and severe malarial anaemia but not respiratory distress were associated with cys2 and group A-like var gene expression

We fitted five independent multiple linear regression models predicting the expression of cys2 var genes and its subtypes (bs1_cys2, cys2 MFK+REY-, cys2 MFK-REY+, cys2 MFK-REY-) hereafter referred to as group A-like, cp1 cp2 and cp3 respectively (Table 1, Methods) as dependent variables. The explanatory variables included the three primary clinical syndromes IC, RD, and SMA, along with age and peripheral parasite density.

Compared to the non-IC clinical group, IC was significantly and positively associated with cys2 [beta (β) (95%CI): 0.09(0.04, 0.14), P= 0.0002], group A-like [β(95%CI): 0.07(0.03, 0.12), P=0.001], and cp1 [β(95%CI): 0.09(0.05, 0.12), P= 00001] (Figure 2A, Table S3). RD did not show significant association with expression of any of the DBLα-tag-defined var subtypes (Figure 2B, Table S3). SMA was significantly associated with higher expression cys2 [β(95%CI): 0.07(0.01, 0.14), P=0.03] and group A-like [β(95%CI): 0.07(0.01, 0.14), P=0.02]. Although the association of SMA with cp2 did not reach statistical significance when all samples were considered (Figure 2C, Table S3), however, it became significant when children with overlapping syndromes were excluded from the analysis (Figure S2, Table S4).

Figure 2.

IC was differentially associated with higher DC8 transcript levels compared with SMA and RD

Next, we explored the differential association of transcript quantity obtained with primers targeting expression of general group A (gpA1, gpA2), DC8 (dc8-1 to dc8-4), group B (dc9, gpB) and group C (gpC) (Table S1) with the severe malaria syndromes. IC showed significant positive associations with the transcript levels detected by the gpA1 primer [β (95%CI): 0.36(0.09, 0.62), P=0.008]] (Figure 2A, Table S3). Additionally, IC was positively associated with the transcript levels obtained with the DC8 targeting primers; dc8-2 [β (95%CI): 0.41(0.08, 0.74), P=0.02], and dc8-3 [β (95%CI): 0.35(0.04, 0.66), P=0.03] as well as the overall DC8 transcript quantity, as captured by the median transcript quantity of the four DC8 primers, dc8median [β (95%CI): 0.33(0.04, 0.62), P=0.03] (Figure 2D and Table S3). When children with overlapping syndromes were excluded, IC remained significantly associated with transcript levels of gpA1 [β (95%CI) 0.44(0.14, 0.74), P=0.004], as well as the transcript levels detected with the DC8 targeting primers, dc8-2, dc8-3, and dc8-median (Figure S2D and Table S4). Notably, in this restricted analysis, IC showed significant positive associations with transcript levels measured with gpA2 [β (95%CI): 0.41(0.11, 0.71), P=0.01], and gpA-median [β (95%CI): 0.40(0.10, 0.70), P=0.01] (Figure S2D, Table S4).

RD showed no significant associations with any var transcript quantities (Figure 2E). However, when children presenting with overlapping syndromes were excluded, pure RD cases were significantly associated with higher transcript levels detected by the gpA2 [β (95%CI): 0.69(0.15, 1.24), P=0.01], gpA-median [β (95%CI): 0.59(0.05, 1.13), P=0.03], and dc8-2 [β (95%CI): 0.78(0.13, 1.44), P=0.02] (Figure S2E, Table S4).

Conversely, SMA showed significant positive association only with the transcript levels measured by the gpA2 primer [β (95%CI): 0.39(0.02, 0.77), P=0.04] (Figure 2F, Table S3). Interestingly, SMA showed significant negative association with transcript quantity obtained with the primer targeting overall group B expression (gpB) and DC9 (dc9): gpB [β (95%CI): -0.49(-0.81, -0.16), P=0.003] and dc9 [β (95%CI): -0.44(-0.83, -0.04), P=0.03] (Figure 2F, Table S3). Upon exclusion of children with overlapping syndromes, SMA remained positively associated with transcript levels measured by the gpA2 primer [β (95%CI): 0.74(0.26, 1.22), P=0.002] but lost significance with gpB and dc9. Additional positive significant association emerged with transcripts detected by the primers gpA1 [β (95%CI): 0.51(0.03, 0.98), P=0.04], gpA-median [β (95%CI): 0.62(0.14, 1.09), P=0.01], dc8-1[β (95%CI): 0.79(0.18, 1.39), P=0.01],dc8-2 [β (95%CI): 0.59(0.02, 1.17), P=0.04], dc8-3[β (95%CI): 0.68(0.14, 1.22), P=0.01], and dc8-median [β (95%CI): 0.54(0.03, 1.04), P=0.04] (Figure S2F, Table S4).

SMA parasites have a higher homogenous var expression profile among admissions

To determine whether parasites associated with different severe malaria syndromes differ in the number of distinct var genes they express, we compared var expression homogeneity (VEH) across clinical malaria presentations. VEH was calculated from the expressed sequence tag data as previously described (Warimwe et al., 2013). Parasites from children with SMA exhibited the highest VEH, with significantly higher values than those observed in IC (p = 0.02), moderate (p = 0.0003), and mild (p = 0.03) malaria cases (Figure 3A). Consistent with this, children with SMA, either alone or in combination with other severe syndromes, showed higher values than non-SMA cases (p = 0.01) (Figure 3B).

Figure 3:

To explore factors underlying the elevated VEH observed in SMA, we examined associations with parasite and host variables. VEH was negatively associated with global var transcript abundance quantified by RT-qPCR. Notably, transcripts from group A var genes showed the weakest negative association, whereas group B and C transcripts showed the strongest negative associations (Figure 3C). This pattern suggests that increasingly homogeneous var expression profiles are more likely to be dominated by group A var genes, consistent with the negative association of group B and C expression with SMA and the positive association of VEH with this syndrome. In contrast, PfSir2a expression, HRP2 levels, and peripheral blood parasite density showed weak, non-significant negative correlations with VEH (Figure 3C), although PfSir2a expression was significantly lower in SMA compared with IC (p = 0.03) (Table S5).

Previous work has reported higher VEH in asymptomatic infections than in clinical malaria, a pattern attributed to antibody-mediated selection against infected erythrocytes (Warimwe et al., 2013). In line with this, antibody responses to infected erythrocytes were positively associated with VEH (ρ = 0.14, p = 0.00006) (Figure 3C). However, this mechanism is unlikely to account for the elevated VEH observed in SMA, as affected children were younger (Table S2) and had the lowest antibody levels (Figure 3D). Fever duration prior to admission was longest in SMA and shortest in IC (Figure 3E), consistent with SMA reflecting a more prolonged infection, although fever duration itself was not associated with VEH (ρ = 0.04, p = 0.45) (Figure 3C). In contrast, febrile temperature showed a significant negative association with VEH (ρ = −0.11, p = 0.01) (Figure 3C) and was lowest in SMA (median [IQR] 38°C [37.3-38.8]; Table S5), consistent with reduced PfSir2a expression in this group (Abdi et al., 2023).

Taken together, these findings suggest that SMA is associated with a more chronic infection state, in which sustained cytoadhesion-related selective pressures may favour parasites with a narrower, group A-dominated var expression profile compared with IC and respiratory distress.

Group A var genes and rosetting are differentially associated with severe malaria syndromes

Using 254 children with rosetting information, we fitted a binary logistic model to determine if rosetting is associated with severe malaria compared to non-severe. Rosetting frequency was significantly associated with severe malaria compared to non-severe malaria [OR: 3.16 (1.33, 7.95), p = 0.01] (figure 4A).

Figure 4:

To evaluate the relationship between rosetting and individual syndromes while controlling for overlapping presentations, we constructed a linear regression model with rosetting frequency as the outcome, and IC, RD, and SMA as predictors while adjusting for age and parasite density. Rosetting was significantly associated with RD [β(95%CI) 0.12(0.02,0.21), P=0.02] and SMA [β(95%CI) 0.15(0.02,0.23), P=0.01], but not with IC [β(95%CI) 0.05(-0.03,0.12), P=0.22] (Figure 4B). We further show that rosetting was positively correlated with the expression of group A and DC8 var genes, and negatively with expression of group B var genes (Figure 4C).

To determine whether RD and SMA associations with rosetting were independent of var gene expression, we conducted logistic regression models for each syndrome with rosetting frequency as a predictor while adjusting for expression levels of var gene subtypes that showed significant association with rosetting and RD or SMA individually (Table 2). Rosetting remained significantly associated with RD even after adjusting for group A-like [OR (95%CI): 3.33 (1.10, 9.79), P=0.03], gpA1 [OR (95%CI): 3.29 (1.06, 10.02), P=0.04], and gpA2 [OR (95%CI): 3.2 (1.03, 9.67), P=0.04] (Figure 4D, Table 2(models 5B, 6B&7B)). Similarly, Rosetting remained significantly associated with SMA after adjusting for group A-like [OR (95%CI): 3.74(1.13, 12.13), P=0.04] while adjusting for gpA1 and gpA2 resulted in loss of statistical significance in the association between SMA and rosetting (Figure 4E, Table 2 (models 6C&7C)).

Table 2:

Group A-like var expression is independently associated with malaria mortality

Among the 712 patients included in this study, 40 died (Table S2). Among the 40 patients that died, 14 had overlapping syndromes, 5 had moderate malaria, 20 had IC, and 1 had SMA while none died from either mild or RD (Figure 5B). We performed a logistic regression model with both group A-like var expression levels and parasite density as explanatory variables with death as an outcome. Parasite density was not associated with child death [OR (95%CI): 1.00 (0.99, 1.00), p= 0.22]. However, death was predicted by group A-like var gene expression levels [OR (95%CI): 5.61 (1.80, 17.54), p=0.003], parasite density adjusted) (Figure 5A). The association was maintained even when the analysis was restricted to only children with impaired consciousness [OR (95%CI): 14.40 (3.140, 74.5), P = 0.0008, adjusted for parasite density). cys2 and cp2 were also associated with increased risk of mortality while gpB and gpC were associated with reduced risk of mortality, while the rest of var subtypes showed no significant association with mortality outcome (Table S6).

Figure 5:

Discussion

This study provides a comprehensive analysis of the relationship between var gene transcription profiles and the severe malaria syndromes IC, RD and SMA (Marsh et al., 1995). By integrating var expression data from both DBLα-tag sequencing and RT-qPCR analysis of parasite samples collected over 18-year period of changing transmission intensity in Kilifi, Kenya, we provide important insights into the relationship between expression of var gene subtypes and clinical presentation of severe malaria. Our findings reveal distinct associations between expression of var gene subgroups and severe malaria phenotypes. Notably, transcripts belonging to the cys2 and DC8 were associated with an increased risk of IC while RD showed no association with var expression. SMA was characterised by parasites exhibiting a more homogeneous var expression profile, while rosetting was associated with RD and SMA, but not with IC.

The pathogenesis of cerebral malaria, the more severe form of IC, has been closely link to sequestration of iRBC in the brain microvasculature (MacPherson et al., 1985; Taylor et al., 2004). This process is mediated by PfEMP1 variants capable of binding host endothelial receptors such as endothelial protein C receptor (EPCR), a characteristic feature of DC8-expressing parasites (Avril et al., 2013; Storm et al., 2019; Turner et al., 2013). Our data is consistent with previous findings that demonstrated higher expression of EPCR-binding PfEMP1 in parasites from cerebral malaria than SMA (Shabani et al., 2017), reinforcing the role of DC8 variants in the pathogenesis of cerebral malaria. The elevated expression of group A and DC8 var genes in IC cases further underscores their importance in disease pathogenesis and highlights these PfEMP1 variants as potential targets for intervention. In this context, recent a recent study has demonstrated that, despite their high sequence variability, EPCR-binding CIDRα1 domains maintain a conserved structural motif (Rajan Raghavan et al., 2023). This structural conservation may allow recognition by broadly neutralizing monoclonal antibodies, but this remains to be tested on diverse parasite isolates (Reyes et al., 2024). Similarly, Tessema et al., reported that antibodies against cp1 recombinant proteins were significantly associated with protection against severe malaria (Tessema et al., 2019). Our findings together with previous studies (Kyriacou et al., 2006; Shabani et al., 2017; Warimwe et al., 2009), illustrate the potential value of targeting PfEMP1 subtypes associated with cerebral malaria, which is the severe syndrome that carries the highest case fatality rate.

Our results suggest that PfEMP1 subtypes targeted by gpA2 and dc8-1 primers, as well as those belonging to cp2 may play a central role in the pathogenesis of SMA. We also observed that parasites from SMA cases displayed a more homogeneous var expression profile compared to the other severe syndromes. Previously, higher VEH has been observed in asymptomatic infection compared to clinical cases (Warimwe et al., 2013), with this pattern attributed to selection pressure exerted by antibodies against iRBCs. In SMA, however, antibody levels are unlikely to drive VEH, as affected children are predominantly young with low antibody levels. Instead, the chronicity of infection, reflected by the reported long fever duration, and the lower febrile temperatures observed in SMA cases in this study may represent a shared feature underlying the observed VEH patterns in both SMA and asymptomatic infections. Overall, our data indicate that SMA is associated with parasites exhibiting a more homogenous var expression pattern, dominated by group A subtypes, including those containing the cp2 motif. While PfEMP1-mediated cytoadhesion likely contributes to SMA pathogenesis, chronicity of the infection as previously suggested in a review by White, (2018) may also lead to loss of erythrocytes due to continuous destruction of red cells (Moss et al., 2026), and inflammation induced dyserythropoiesis (Knuttgen, 1987; Liu & Li, 2025; Lyke et al., 2004; Phillips et al., 1986).

Consistent with a previous study (Warimwe et al 2012), our data showed that rosetting was significantly associated with RD, but not with IC. Additionally, we showed that rosetting was associated with SMA. However, when we adjusted for var gene expression, the association between rosetting and SMA dropped and lost statistical significance while that with RD remained significant. This pattern suggests that rosetting may be a driving factor for RD pathogenesis but not for SMA. In contrast, adjusting for rosetting did not eliminate the associations between group A-like and gpA1 var transcripts for IC, and adjusting for gpA2 did not alter its association with SMA, supporting the hypothesis that these var gene subgroups play a primary role in the pathogenesis of respective syndromes.

A previous smaller study (Warimwe et al 2012) whose samples were also included in this study, showed no var gene group remained associated with RD after adjusting for rosetting, suggesting that rosetting may be a more dominant determinant of RD pathogenesis than var gene expression per se. Alternatively, the lack of association could reflect the absence of specific var subtypes relevant to RD pathogenesis in our dataset. Supporting this, McLean et al., (2025) identified domain cassettes DC11 and DC15 as mediators of rosetting in laboratory isolates, subtypes not directly captured by our primers. Future studies should therefore incorporate primers targeting these domain cassettes to elucidate the molecular basis of RD in field isolates. It is also possible that other variant surface antigens (VSAs), beyond the var gene family contribute to RD pathogenesis. Members of the surface-associated interspersed gene family (SURFIN), repetitive interspersed family (RIFIN), and subtelomeric variable open reading frame (STEVOR) proteins have been implicated in malaria pathogenesis (Kanoi et al., 2020). Future investigations of RD should therefore consider these additional VSA families to achieve a more complete understanding of RD pathogenesis.

Mortality was highest among children with IC. This finding corroborates previous work showing that mortality is highest in children with cerebral malaria, a severe form of IC (Idro et al., 2007). In our analysis, group A-like var gene expression was significantly associated with both IC and mortality, aligning with previous reports (Warimwe et al., 2009, 2012). These findings suggest that group A-like expressions may serve as a marker for malaria-related mortality and highlight its potential clinical utility in prognostic models.

A key strength of this study is the use of DBLα-tag sequencing and RT-qPCR approaches, which have provided complementary understanding of the var transcriptome. Whereas DBLα-tag sequencing captures the overall breadth of expressed var repertoires, RT-qPCR allows precise quantification of clinically important PfEMP1 subsets. Considering both datasets together revealed subtype-syndrome associations that would likely have been missed using either method alone. Despite this advantage, the cross-sectional design limits causal inference, and undetermined host genetic or immunological factors may also shape var expression profiles. Therefore, future studies should incorporate these variables to refine the observed relationships.

In conclusion, the use of DBLα-tag sequencing and RT-qPCR has allowed us to build a clearer and more reliable picture of subtype-specific var gene expression across the major severe malaria syndromes. Our findings consistently point EPCR-binding observed in impaired consciousness may be driven by DC8, cp1 and group A-like PfEMP1 variants, reinforcing the central role of sequestration in the pathogenesis of impaired consciousness. Taken together, these results highlight the importance of focusing on PfEMP1 subsets mediating sequestration and provide a strong rationale for their consideration in the development of future therapeutics and vaccines.

Materials and Methods

Study population

This study was conducted in Kilifi County, a malaria-endemic region located along the Kenyan coast. The patient samples analyzed in this study have been previously described (Abdi et al., 2023). Briefly, the cohort comprised 521 children diagnosed with malaria who were admitted to Kilifi County Hospital (KCH) and 191 treated for mild malaria at outpatient clinics in Kilifi between 1994 and 2012. This period represents three distinct phases of malaria transmission in the region: pre-decline, decline, and post-decline (Mogeni et al., 2016). These children were enrolled through the hospital’s admission surveillance system, as well as from children presenting with mild malaria at an outpatient clinic.

Definition of malaria syndromes

Severe malaria was defined as a hospital admission with any of the following; impaired consciousness (IC) defined using age-appropriate Blantyre coma score (BCS) thresholds (BCS of less than 5 for children over 8 months old or a BCS of less than 4 for children younger than 8 months) (Berkley et al., 2009), or respiratory distress (RD) (deep breathing) (English et al., 1996) or severe malarial anaemia (SMA) (hemoglobin < 5g/dl) (Marsh et al., 1995) or overlapping syndromes with more than one severe malaria syndrome presentation.

Sample inclusion and size

Participants included in this study were those in the KEMRI database with complete information on clinical syndrome, age, parasite density and var gene expression of the iRBC determined by both DBLα-tag sequencing and RT-qPCR quantification approaches. 712 participants fulfilled these criteria including 340 with non-severe malaria and 372 children with severe malaria including pure syndromes (59 SMA, 29 RD, 186 IC) and 98 with overlapping syndromes in the form of IC/RD, IC/SMA, RD/SMA or IC/RD/SMA (Figure 1B). Non-severe cases included 191 malaria positive children managed at outpatient clinics referred to as mild and 141 children admitted to hospital who did not present with any of the three severe syndromes described above, referred to as moderate.

Peripheral Blood Parasitemia

Peripheral parasitemia was quantified from Giemsa-stained thick and thin blood films following standard microscopy procedures. Parasite density was initially recorded as the number of infected erythrocytes per 500, 200, or 100 red blood cells (RBCs) examined. Parasitemia per microlitre of blood was then calculated using World Health Organization (WHO) guidelines (WHO, 2010). Briefly, parasite density (parasites/µl) was derived using either:

or, for counts performed against white blood cells (WBCs):

In instances where direct measurements of RBC or WBC concentrations were unavailable, standard assumptions of 5 × 10⁶ RBCs/µl and 8,000 WBCs/µl were applied, in accordance with WHO recommendations.

PfHRP2 Assay

Plasma concentrations of P. falciparum histidine-rich protein 2 (PfHRP2) were measured using a sandwich ELISA as previously described (Abdi et al., 2023). In brief, Nunc MaxiSorp 96-well plates (Thermo Fisher Scientific) were coated overnight at 4°C with mouse anti-PfHRP2 monoclonal antibody (MPFM-55A; MyBioscience) diluted in PBS to a final concentration of 0.9 µg/ml. Plates were washed three times with PBS containing 0.05% Tween-20 and subsequently blocked with PBS supplemented with 3% skimmed milk for 2 hr at room temperature (RT).

After an additional wash step, plasma samples and recombinant PfHRP2 standards (MBS232321; MyBioscience), diluted in PBS with 2% BSA, were added in duplicate and incubated for 2 hr at RT. Plates were washed and incubated for 1 hr with HRP-conjugated mouse anti-PfHRP2 detection antibody (MPFG-55P; MyBioscience) diluted in PBS/2% BSA to a final concentration of 0.2 µg/ml. Following three washes, colour development was achieved using o-phenylenediamine dihydrochloride (OPD) substrate for 15 min, and reactions were stopped with 2M H₂SO₄. Optical density (OD) was read at 490 nm using a BioTek Synergy4 plate reader. The ODs where log-transform for further analysis.

RNA and cDNA preparation

The method for RNA extraction and cDNA preparation from the samples used in this study was described previously (Abdi et al., 2016; Warimwe et al., 2012). Briefly, samples stored in Trizol (Invitrogen, Paisley, United Kingdom) at −80°C were thawed, and then total RNA was extracted by chloroform phase separation and GlycoBlue/isoamyl alcohol precipitation. Pellets were washed with 75% ethanol, briefly air-dried, and re-suspended in RNAsecure. Residual genomic DNA was removed by DNase I treatment (DNase-Free, Ambion) and inactivated with the Ambion DNase Inactivation Reagent according to the manufacturer’s instructions. cDNA was prepared by reverse transcribing RNA using reverse transcriptase (Invitrogen Superscript III). For each cDNA reaction, a negative control reaction was performed in the absence of reverse transcriptase to ensure that all contaminating DNA had been removed by DNAse pre-treatment. No sample showed any amplification in these controls.

Var sequencing and classification

The expressed sequence tag (EST) sequencing approach used to generate this dataset has been previously described (Abdi et al., 2016; Warimwe et al., 2009). In summary, degenerate primers (Bull et al., 2005), were used to amplify a region within the DBLα-tag domain that is highly conserved across most var genes. Two sequencing methods were employed: capillary sequencing and the 454-pyrosequencing platform.

For the capillary sequencing approach, amplified DBLα-tags were cloned into the TOPO®-TA cloning vector and transformed into E. coli. From each sample, up to 96 colonies were picked and sequenced individually. While for the 454 platform, the amplicons were barcoded, multiplexed, and sequenced using high-throughput 454 sequencing (Kivisi et al., 2019).

Sequences obtained from each sample were then classified and counted following methods described previously (Bull et al., 2005; Warimwe et al., 2009). EST proportional expression scores were calculated as the proportion of sequences in each sample that belonged to specific sequence groups. These groups were defined based on the number of cysteine residues, the presence or absence of conserved amino acid sequence motifs (MFK and REY) located in positions of limited variation (PoLV) and presence of block sharing block within the DBLα-tag. Overall, according to number of cysteine residues DBLα-tag sequences were classified into three groups namely, cys2 for those containing two cysteine residues, cys4 for those with four and cysX for those containing one, three, five, or six cysteine residues. Combining the number of cysteine residues and PoLV (cp grouping) resulted into six groups (cp1–cp6) as shown in Table 1. All known cp1 sequences fell in group A var genes; however, not all group A sequences can be grouped as CP1 (Bull et al., 2007). A network of recombining sequence that consisted of sequence blocks that tent to recombine with each other was used to determine block sharing groups (BS). The BS groups were made up of polymorphic blocks together with the number of cysteine residues in the DBLα-tag. Sequences that fall into block-sharing Group 1 and have two cysteine residues (BS1_cys2) belong to group A-like var genes. Only cp1, cp2, cp3 and group A-like are used in the analysis because of the previous association with severe malaria (Warimwe et al., 2009).

RT-PCR var transcription profiles

The qRT-PCR primers used to determine var transcript profile data in this study was published by Abdi et al, (Abdi et al., 2016) and was originally designed by (Lavstsen et al., 2012; Rask et al., 2010; Rottmann et al., 2006). Data was generated using primers that showed differential expression between severe and non-severe malaria in previous studies (Lavstsen et al., 2012). These included four primer sets targeting DC8 (named dc8-1, dc8-2, dc8-3, dc8-4), one primer set for each of DC13 (dc13), DC4 (dc4), and DC9 (dc9), two primer sets targeting the majority of group A var genes (gpA1 and gpA2). Group B and C var gene transcription was quantified using primer sets b1 and c2 respectively. The names and sequences of the original primers are shown in Table S1. In addition, we included a primer targeting Pfsir2a, which plays a role in epigenetic control (Abdi et al., 2016).

Levels of var gene transcription were normalized relative to the housekeeping gene seryl tRNA synthetase and fructose bisphosphate aldolase. The arbitrary transcript unit (Tus) was calculated using the formula (Tus = 2^(5−Δct)) where Δct is the averaged change in cycle threshold of var gene relative to housekeeping genes (Abdi et al., 2016; Lavstsen et al., 2012).

Variant Expression Homogeneity (VEH) Calculation

To quantify the degree of dominance among var gene transcripts expressed by each parasite isolate in a sample, we assessed Variant Expression Homogeneity (VEH). VEH reflects the extent to which a small number of var gene sequences dominate an isolate’s expression profile (Warimwe et al., 2013).

We employed Simpson’s Index of Diversity to capture the overall homogeneity of the var expression repertoire (Simpson, 1949; Warimwe et al., 2013). For each sample, we calculated Simpson’s Index (λ) as the sum of the squares of proportions of all distinct var gene expressed in each sample:

Where; π_i is the proportion of var gene i among N total var genes counted in a sample. The index ranges from (indicating maximal diversity and minimal homogeneity) to 1 (indicating complete dominance by a single sequence and thus maximal VEH).

Rosetting Assay

Cryopreserved parasite isolates were thawed following previously described methods (Kinyanjui et al., 2004). A suspension containing trophozoite-infected erythrocytes at 4% hematocrit was prepared in RPMI 1640 medium (Gibco) supplemented with 5μg/ml acridine orange (Sigma-Aldrich). For each isolate, 9.5μl of this suspension was mixed with 2.5μl of non-immune AB serum in a U-bottomed 96-well plate. The mixture was rotated vertically at room temperature for 30 minutes. Subsequently, the sample was transferred onto a microscope slide, covered with an 18mm x 18mm coverslip, and examined under a fluorescence microscope (Nikon Eclipse 80i, Japan) at 400x magnification. Rosetting frequency was determined as the percentage of 200 trophozoite-infected erythrocytes binding to two or more uninfected erythrocytes. Scoring was conducted without knowledge of the clinical category of the patient from whom each isolate originated.

Antibody Recognition of trophozoite-infected erythrocytes (α-IE Assay)

Total IgG binding to the surface of trophozoite-infected erythrocytes (IEs) was quantified by flow cytometry, as described previously (Kivisi et al., 2019). Briefly, plasma samples were incubated with IEs from four heterologous clinical isolates, and antibody binding was detected using a fluorescein isothiocyanate (FITC)-conjugated anti-Fcγ antibody (AF004, The Binding Site). Infected and uninfected erythrocytes were differentiated by ethidium bromide staining, enabling correction for nonspecific IgG binding by subtracting the MFI of uninfected erythrocytes from that of matched IEs for each sample. To account for isolate-specific background reactivity, the highest MFI observed among plasma from four malaria-unexposed European donors was subtracted from the corresponding values obtained with each child’s plasma. The resulting background-adjusted MFI values were used in all downstream analyses.

Statistical analysis

We selected samples with both RT-qPCR and DBLα-tag sequence transcript data for analysis. RT-qPCR-derived var transcript expression values were log-transformed, while proportions of DBLα-tag-defined var subtypes were arcsine square root transformed (Warimwe et al., 2009). To identify associations between var gene expression and severe malaria syndromes (IC, RD, and SMA), binary variables of the severe syndromes were created and non-severe cases were considered as controls. We fitted multiple linear regression models using stats R core package (R Core Team, 2021). The model followed the general structure below in equation 1:

Where:

y represents specific var gene defined by RT-qPCR primer sets or DBLα-tag subtype, b represents binary representation where all children with IC syndrome were assigned 1 and those not containing it were assigned 0, c represents binary representation where children with RD syndrome were assigned 1 and those not containing were assigned 0; d; represents binary representation where all children with SMA syndrome were assigned 1 and those not containing where assigned 0 e; represents age of the participant in months while f represent log transformed values of peripheral blood parasite density. This analysis was repeated where children with overlapping syndromes were excluded with the aim of identifying syndrome-specific var gene variants.

The output of the model are p-values and beta coefficient of estimate with 95% confidence interval (β (95%CI). Positive and negative beta (β) values with p-values less than 0.05 show positive and negative significant associations with specific malaria syndrome respectively in comparison to the control group. The results are presented as forest plots or tables.

We compared VEH scores across pure clinical presentations using two-sided Wilcoxon rank-sum tests to determine the syndrome with more homogenous var expression. To explore determinants of VEH, we assessed associations with host factors (α-IE, fever duration before admission and febrile temperature) and parasite factors (Pfsir2a, peripheral blood parasite density and RT-qPCR var gene transcripts expression levels) using Spearman’s rank correlation. We report correlation coefficients (rho) and p values, and visualize group differences and associations using ggplot2 in R.

To dissect the relationship between var expression, severe malaria syndromes and rosetting, the frequency of rosetting was arcsine square root transformed. Binary logistic regression model on 253 children with rosetting data including 133 severe and 120 non-severe cases was fitted to determine the association between rosetting and severe malaria using logistic regression. Correlation between rosetting frequency and var gene expression was assessed using Spearman’s rank correlation coefficient. To estimate 95% confidence intervals (95% CIs) for the correlation coefficients, we applied nonparametric bootstrapping with 5000 replicates using boot R package (Canty & Ripley, 2024). This was followed by evaluating the relationship between rosetting and individual severe malaria syndromes. We constructed a linear regression model as in equation 1 and replaced y with rosetting frequency. To test whether the relationship between rosetting and the severe syndromes is independent of var expression, we conducted logistic regression models considering each of syndromes as outcome with rosetting frequency as a predictor while adjusting for expression levels of var gene subtypes.

To identify parasite factors associated with mortality, we fitted binary logistic regression models with death (yes/no) as the outcome. We first ran univariable models for each var transcript measure and for parasite density separately and then fitted multivariable models including the var measure and parasite density together to estimate the effect of each on death while adjusting for age in all models. Odds ratios (ORs) with 95% CIs are reported.

All analyses were conducted in R (version 4.4.1) (R Core Team, 2021). Forest plots of effect sizes from linear and logistic regression models were generated using ggplot2 (Wickham, 2016). Data manipulation was performed using dplyr and tidyverse packages (Wickham et al., 2023), and model tidying was implemented with broom.

Supplementary Figures

Figure S1:

Figure S2:

Supplementary Tables

Table S1:

Table S2:

Table S3:

Table S4:

Table S5:

Table S6:

Table S7:

Data availability

The code for the analyses presented in this paper and data is deposited at https://doi.org/10.7910/DVN/W4FFT9

Acknowledgements

This work was supported by Wellcome Trust Awards: 209289/Z/17/Z (to Abdirahman I. Abdi). Henry Ndugwa is fully funded by the Science for Africa Foundation to the Developing Excellence in Leadership, Training and Science in Africa (DELTAS Africa) programme [DEL-22-012] with support from Wellcome Trust and the UK Foreign, Commonwealth & Development Office and is part of the EDCPT2 programme supported by the European Union.

Additional information

Ethical approval

This study obtained ethical approval from the Scientific Ethics Review Unit of the Kenya Medical Research Institute (KEMRI/SERU/3149) and written informed consent was obtained from parents or guardians of the study participants.

Author contribution

Henry Ndugwa Conceptualization, Investigation, methodology, Formal analysis, Visualization, Writing – original draft, review and editing,

Michelle Muthui methodology, Writing–review and editing,

J. Alexandra Rowe Writing –review and editing

Kinyanjui SM, Project administration, Supervision, Conceptualization, analysis, Writing-original draft, review and editing

Cheryl Andisi KivisiProject administration Supervision, Conceptualization. Methodology, analysis, Writing – original draft, review and editing

Abdirahman I. Abdi Conceptualization, Supervision, formal analysis, Writing – original draft, review and editing, All authors contributed to the interpretation of the analyses and revised the draft manuscript.

Funding

Wellcome Trust (WT)

https://doi.org/10.35802/209289

• Abdirahman I Abdi

Science for Africa Foundation (SFA Foundation) (DEL-22-012)

• Henry Ndugwa