Distinct gastrointestinal microbial signatures predict parasite levels in controlled Plasmodium infections in both rhesus macaques and humans

Abstract

Functions of the gastrointestinal (GI) microbiome include maintenance of immune homeostasis and protection against infectious disease. Current assessments of the role of the GI microbiome in Plasmodium infection have been primarily conducted using mouse models and observational human cohorts. Here, we experimentally assessed associations between pre-infection GI microbiome composition and acute Plasmodium parasitemia using 16S rRNA sequencing and samples from rhesus macaques (RMs) and adult humans enrolled in a previously conducted controlled human malaria infection (CHMI) trial (NCT04072302) originally designed to test the efficacy of KAF156, a novel imidazolopiperazine class of antimalarial drugs. We identified distinct pre-infection 16S microbial signatures that were associated with increased risk for above median parasitemia in RMs infected with P. fragile and CHMI participants infected with P. falciparum. Further, we identified a Bifidobacterium feature set that accurately stratified parasitemia risk and could therefore serve as a foundation for a potential biomarker panel to aid prevention efforts in malaria endemic regions. Together, our findings demonstrate that pre-infection GI microbiome composition is indicative of risk for Plasmodium parasitemia, and our observation that the pre-infection microbiome-P. fragile dynamic in RMs mirrors the pre-infection microbiome-P. falciparum interaction in CHMI participants supports the future use of this model in pre-clinical investigations of novel microbiome-targeting approaches to reduce malaria burden.

Fig. 1. Unsupervised hierarchical clustering of 16S rRNA sequences identifies association between GI microbiota composition, housing, diet, and risk of above median P. fragile parasitemia in RMs.

A Divisive hierarchical clustering to identify major clusters in baseline microbiota samples. RM microbiota profiles segregated according to housing location/diet using this unsupervised approach (Building A/high fiber=teal; Building B/high protein=purple; n = 8 RMs per group). B Percent infected RBC demonstrate significantly higher peak peripheral parasitemia in Building A/high fiber RMs as compared to Building B/high protein RMs (p = 0.00031; n = 8 RMs per group). Values along the left y axis indicate Z-scored %RBC, while right y axis values indicate %infected RBC. Individual data points represent Z-scored %RBC at peak parasitemia for each macaque; box limits represent the upper and lower quartiles, while box whiskers represent 1.5x the interquartile range. Horizontal lines within each box represent the median. Points overlaying the box and whisker bars represent data from individual RMs. Statistical significance between Building A/high fiber and Building B/high protein RM parasitemia was determined using a 2-sided Wilcoxon rank-sum test, with the p value shown above the horizontal line at the top of the plot. C %RBC data for each macaque across the entire study. Plots are ordered according to peak parasite load for each RM and colored according to cohort. W=Week.

Fig. 2. Aitchison distance demonstrates durable microbial distinctions between housing/diet cohorts.

A Principal Coordinate Analysis (PCoA) of Aitchison distance for animals. Colors indicate hierarchically determined microbiota clusters for each cohort (Building A/high fiber=teal; Building B/high protein=purple; n = 8 RMs per group). Shapes represent unique time-points, with large circles representing the earliest sample time point (Week [W]-2), small circles representing Week 0 (W0) through Week 4 (W4), and large squares representing the final observation (Week [W]6). Lines connect timepoints for each RM. B Total sample-to-sample distance across time for individual animals demonstrate higher microbiota variability over time in Building B/high protein cohort as compared to Building A/high fiber cohort (p = 0.028; n = 8 RMs per group). Box limits represent the upper and lower quartiles, while box whiskers represent 1.5x the interquartile range. Horizontal lines within each box represent the median. Points overlaying the box and whisker bars represent data from individual RMs. Statistical significance between Building A/High Fiber and Building B/High Protein RM total Aitchison distance traveled was determined using a 2-sided Wilcoxon rank-sum test, with the p value shown above the horizontal line at the top of the plot.

Fig. 3. Shannon diversity and relative abundance of core microbiota do not differ across housing/diet cohorts.

A Shannon diversity assessment of aggregated baseline samples (Weeks -2 and 0) shows a trend toward increased baseline diversity in the Building B/high protein cohort (purple) as compared to Building A/high fiber cohort (teal; p = 0.083; n = 8 RMs per group). B The relative abundance of the core microbiota does not significantly differ between Building A/high fiber and Building B/high protein RMs (p = 0.150; n = 8 RMs per group). In both A and B, box limits represent the upper and lower quartiles, while box whiskers represent 1.5x the interquartile range. Horizontal lines within each box represent the median. Points overlaying the box and whisker bars represent data from individual RMs. Statistical significance between Building A/high fiber and Building B/high protein RM Shannon diversity and core relative abundance was determined using a 2-sided Wilcoxon rank-sum test, with the p value shown above the horizontal lines at the top of each plot.

Fig. 4. Univariate and multivariate assessments identify ASVs and higher taxonomic levels that differ and distinguish microbiota cohorts.

A Differential abundance of amplicon sequence variants (ASVs) across cohorts (Building A/high fiber=teal; Building B/high protein=purple; n = 8 RMs per group). Each lolliplot represents the log2fc for ASVs or higher taxonomic levels that differ significantly, as determined by DESeq2 analysis, with points to the left being more prevalent in Building A/high fiber and those to the right more prevalent in Building B/high protein. Taxonomic levels include the class, order, family, and genus levels. The color intensity indicates the degree of log2fc. B Bar graph depicting the prevalence levels of differentially abundant ASVs in Building A/high fiber and Building B/high protein RMs (n = 8 RMs per group). The height of each bar indicates the proportion of samples for each cohort that were positive for detection of a particular ASV. Significant differences in prevalence (p < 0.05) as determined using a 2-sided Wilcoxon rank-sum test with Benjamini-Hochberg correction are indicated by orange stars (Prevotella_copri_18, p = 0.318; Lactobacillus_salivarius_33, p = 0.0008; unresolved_Ruminococcaceae_77, p = 0193; unresolved_Lactobacillus_23,p = 1.85e-05; unresolved_[Mogibacteriaceae]_124, p = 0.0015; unresolved_Lachnospiraceae_11, p = 0.0374; unresolved_Ruminococcaceae_89, p = 0.0168; Lactobacillus_mucosae_63, p = 0.0004; unresolved_Christensenellaceae_185, p = 0.0077; unresolved_Mogibacterium_235, p = 0.0008; unresolved_Prevotella_105, p = 0.0008; unresolved_Treponema_147, p = 6.16e-07; unresolved_Moryella_61, p = 0.0286; unresolved_Treponema_138, p = 0.0124; unresolved_Lachnospiraceae_83, p = 0.0124; unresolved_Porphyromonadaceae_136, p = 2.49e-06; unresolved_Treponema_312, p = 0.0043; unresolved_Lachnospiraceae_71, p = 9.83e-10; unresolved_Parabacteroides_200, p = 0.0315; Lactobacillus_coleohominis_152, p = 2.49e-09; unresolved_Succinivibrio_144, p = 0.0012; unresolved_Lachnospiraceae_246, p = 1.60e-09; unresolved_Desulfovibrio_322, p = 3.12e-06; unresolved_RF16_353, p = 0.0124; unresolved_Ruminococcus_146, p = 7.73e-07; unresolved_Fibrobacter_231, p = 1.63e-05; Prevotella_copri_244, p = 0.0023; Dorea_formicigenerans_330, p = 0.0077; unresolved_Parabacteroides_427, p = 0.0124; unresolved_Lachnospiraceae_166,p = 8.49e-05; unresolved_Pasteurellaceae_385, p = 0.0124; unresolved_Prevotella_306,p = 0.0079; unresolved_Coriobacteriaceae_314, p = 0.0329; unresolved_Rosburia_262, p = 0.0124; unresolved_Lachnospiraceae_193, p = 0.0215; unresolved_S24-7_229, p = 0.0181; Prevotella_stercorea_412, p = 0.0318).

Fig. 5. Enrichment in Lactobacillus and lower levels of Prevotella observed in RMs with elevated P. fragile parasitemia.

A Line graph with line height along the y axis illustrating the CLR-transformed levels of Lactobacillus (circles) and Prevotella (squares) and time represented linearly on the x axis for RM cohorts (Building A/high fiber=teal; Building B/high protein=purple; n = 8 RMs per group). Data points are faceted by cohort. Week=W. B Partial Least Squares Discriminant Analysis (PLSDA) showing the separation of microbial communities between the two cohorts based on mean aggregated CLR-transformed counts of ASVs. The x axis represents the first discriminant component (comp1) which captures the maximum variance between the two cohorts. Points represent individual ASVs, with the color indicating the family and the size reflecting the relative prevalence of each ASV. The direction of the axis indicates increased abundance in Building A/high fiber to the right and Building B/high protein to the left.

Fig. 6. Unsupervised hierarchical clustering of 16S rRNA sequences identifies association between GI microbiota composition and risk of elevated P. falciparum parasitemia.

A Divisive hierarchical clustering identifies 3 major clusters in baseline microbiota samples. Major clusters were defined as human cluster 1 (HC1; gold; n = 10 participants), human cluster 2 (HC2; gray; n = 13 participants) and human cluster 3 (HC3; dark cyan; n = 9 participants). Branch tips indicate the parasitemia classification of each participant, which was based on detection of peak parasite levels (above median parasitemia=red or below median parasitemia=blue) via qRT-PCR. B Log2-transformed parasitemia values (qRT-PCR) for participants in each microbiota cluster show significantly higher peak peripheral parasitemia in HC1 participants (HC1 [n = 10] vs HC3 [n = 9], p = 0.043; HC1 vs HC2 [n = 13], p = 0.17; HC2 vs HC3, p = 0.23). Box limits represent the upper and lower quartiles, while box whiskers represent 1.5x the interquartile range. Horizontal lines within each box represent the median. Points overlaying the box and whisker bars represent data from individual CHMI participants. The left side of the figure magnifies the plot on right, which was compressed by outliers. Significant differences in parasitemia between the three groups were assessed using a 2-sided Wilcoxon rank-sum test. C Odds ratio (OR) assessment to determine likelihood of participants within each microbiota cluster progressing to above median parasitemia (HC1 [n = 10] vs HC3 [n = 9]: OR = 8.17, 95% confidence interval [CI] = 1.17-83.38, p = 0.047; HC2 [n = 13] vs HC3: OR = 4.08, 95% CI = 0.67–34.22, p = 0.149).

Fig. 7. Higher GI microbiota diversity and lower relative abundance of core microbes is linked to likelihood of developing elevated peripheral parasitemia.

A Principal Coordinate Analysis (PCoA) of Aitchison distance for CHMI participants. Colors indicate hierarchically determined microbiota clusters (HC1=gold, n = 10 participants; HC=gray, n = 13 participants; HC3=dark cyan, n = 9 participants). Shapes represent unique time-points, with large circles representing baseline samples, small circles representing pre-treatment and post-treatment timepoints, and large squares representing follow-up; lines connect timepoints for each participant. B Boxplots showing the total (Aitchison) distance travelled for participants within each microbiota cluster (HC1 [n = 10] vs HC3 [n = 9], p = 0.095; HC1 vs HC2 [n = 13], p = 0.088; HC2 vs HC3, p = 0.600). C Boxplots showing that the distance between baseline and follow-up samples within each microbiota cluster (HC1 [n = 10] vs HC3 [n = 9], p = 0.0076; HC1 vs HC2 [n = 13], p = 0.0015; HC2 vs HC3, p = 0.36). D Boxplots showing Shannon diversity levels at baseline and at aggregated points post-challenge (baseline: HC1 [n = 10] vs HC3 [n = 9], p = 0.0021; HC1 vs HC2 [n = 13], p = 0.0011; HC2 vs HC3, p = 0.32; post-challenge: HC1 [n = 10] vs HC3 [n = 9], p = 0.022; HC1 vs HC2 [n = 13], p = 0.03; HC2 vs HC3, p = 0.43). E Boxplots showing the relative abundance of the core microbiota at baseline and at aggregated points post-challenge (baseline: HC1 [n = 10] vs HC3 [n = 9], p = 4.3e-05; HC1 vs HC2 [n = 13], p = 1.7e-06; HC2 vs HC3, p = 0.39); post challenge: HC1 [n = 10] vs HC3 [n = 9], p = 8.7e-05; HC1 vs HC2 [n = 13], p = 2.1e-05; HC2 vs HC3, p = 0.69). For B–E box limits represent the upper and lower quartiles, while box whiskers represent 1.5x the interquartile range. Horizontal lines within each box represent the median. Points overlaying the box and whisker bars represent data from individual CHMI participants. Statistical significance between CHMI participant groups was determined using a 2-sided Wilcoxon rank-sum test, with the p value shown above the horizontal line at the top of each plot. F Scatterplot visualizing the relationship between the relative abundance of core microbes and the log2-transformed peak peripheral parasitemia. Each point represents data from individual CHMI participants (HC1, n = 10 participants; HC2, n = 13 participants; HC3, n = 9 participants). The blue line represents the linear model fit (least-squares regression line) and the shaded error bands represent the 95% confidence interval of the estimated mean response. The line indicates a negative association between peripheral parasitemia and core abundance (r² = 0.01541, p = 0.026).

Fig. 8. Univariate and multivariate assessments identify ASVs and higher taxonomic levels that differ and distinguish microbiota cohorts.

A Differential abundance of amplicon sequence variants (ASVs) across microbiota clusters. Baseline samples were compared between HC1 (n = 10 participants) and the combined participants from HC2 and HC3 (n = 22 participants). Each lolliplot represents the log2fc for ASVs or higher taxonomic levels that differ significantly, as determined by DESeq2 analysis, with points to the left being more prevalent in HC1 and those to the right more prevalent in HC2 and HC3. Taxonomic levels include the phylum, class, order, family, and genus levels. The color intensity indicates the degree of log2fc. B Bar graph depicting the prevalence levels of select differentially abundant genera in HC1 (n = 10 participants) and HC3 (n = 9 participants). The height of each bar indicates the proportion of samples for each cohort that were positive for detection of a particular ASV. Significant differences in prevalence as determined using a 2-sided Wilcoxon rank-sum test with Benjamini-Hochberg correction were assessed only for differentially abundant ASVs and are indicated by pink stars (unresolved_Aldercreutzia_109, p = 0.0057; unresolved_Aldercreutzia_182, p = 0.0304; Bifidobacterium_adolescentis_4, p = 0.0099; unresolved_Finegoldia 47, p = 0.0006; unresolved_Lactococcus_163, p = 0.0464; unresolved_Peptoniphilus_122, p = 0.0043).

Fig. 9. Bifidobacteriaceae feature set could serve as biomarker panel to predict risk of elevated parasite burden.

A Select results from Partial Least Squares Discriminant Analysis (PLS-DA) showing the separation of microbial communities between HC1 (n = 10 participants) and HC2/HC3 (n = 22 participants), based on mean aggregated CLR-transformed counts of ASVs. The x axis represents the first discriminant component (comp1) which captures the maximum variance between the microbial clusters. Points represent individual ASVs, with the color indicating the family and the size reflecting the relative prevalence of each ASV. The direction of the axis indicates increased abundance, with taxa associating with lower peripheral parasitemia (HC2/HC3; n = 22 participants) to the left, and taxa associating with higher peripheral parasitemia (HC1, n = 10 participants) to the right. B Receiver Operating Characteristic (ROC) curve for Partial Least Squares-Discriminant Analysis (PLS-DA). The ROC curve plots the true positive rate (Sensitivity) against the false positive rate (1-Specificity) to evaluate the diagnostic ability of the PLS-DA model at various threshold settings. The solid blue line represents the performance of the PLS-DA model, while the dashed red line indicates the line of no-discrimination, which represents a model with no diagnostic ability (random guessing). Bifidobacteriaceae levels achieved an AUROC of 0.801, suggesting a utility in clinical settings.