Rectal swabs from six rhesus macaques (n=6) were analyzed at baseline and longitudinally while infected with P. cynomolgi. The monkeys assigned to the study were malaria-naive (animal codes RAd14, RBg14, RIb13, RJn13, ROc14, and ROh14), of Indian origin, male, 7-13 kg and 5-6 years of age. The overall experimental design, P. cynomolgi M/B strain sporozoite inoculations, clinical monitoring, collection of samples, and treatment regimens have been described previously (Joyner et al. (2019); DeBarry et al. (2022)) and the scheme is depicted here in Figure 1A . Rectal swabs were also obtained once from four malaria-naive control rhesus macaques to perform whole metagenomics shotgun sequencing.

Rectal swabs/feces were collected at predefined timepoints (TPs), chosen based on the infection dynamic (baseline, patent parasitemia, peak infection, post peak, relapse, and sub-patent phases) as shown in Figure 1A (Joyner et al. (2019)). Genomic DNA was extracted from the samples using the Power Soil DNA isolation kit (MoBio Laboratories Inc., Carlsbad, CA) according to the manufacturer’s protocols, and stored at -80°C for further downstream processing.

Bacterial 16S rRNA genes were amplified using bacteria-specific PCR primers for the hypervariable V3-V4 region (Klindworth et al. (2013)). The amplifications, barcoding, and sequencing experiments were performed following established protocols (Caporaso et al. (2012); Kozich et al. (2013)). Sequencing was performed on the Illumina MiSeq platform at the Yerkes Genomics Core at Emory University, resulting in 2 x 275 bp paired-end reads.

Quality filtering of reads and overall sequence analysis from 42 samples (six subjects at eight TPs, Table 1 ) were performed using the Quantitative Insights into Microbial Ecology, QIIME 2 software (version 2019.4). Fastq files were merged, de-multiplexed and the sequence depth for each sample was quantified (Bolyen et al. (2019)). Defined Operational Taxonomic units (OTU’s) were picked for taxonomic identity assignment using the Greengenes reference database (version 13.8) at 97% similarity. The samples were then subjected to QIIME’s core diversity analysis at 20,387 rarefaction depth for beta diversity analysis and TP comparisons. Alpha diversity indexes and their comparisons were computed in R, using the Phyloseq package (McMurdie and Holmes (2013)). The abundance of 16S reads was first normalized to proportional abundances and then the average proportional abundance for each phyla was calculated from each of the TPs ( Table S1 ). The relative percentage abundances of taxa at the family level for each of the subjects, at various TPs during the P. cynomolgi infections, were plotted in R ( Figure 1B ). Alpha diversity was calculated within qiime2 using estimates standard in the field: Chao-1 Richness, Shannon, and Simpson indices (Finotello et al. (2018)). One-way repeated measures analysis of variance (ANOVA) was used to determine whether there were any statistically significant differences between the mean relative abundances of bacterial families across the rhesus macaque subjects at various time points. In order to identify which taxa significantly increased or decreased in abundance between timepoints, a linear discriminant analysis effect size (LEfSe) analyses was utilized with default criteria (p <0.05 by Kruskal-Wallis test; linear discriminant analysis (LDA) score >2) and subject pairing and plotted in cladograms based on phylogenetic relationship (Segata et al. (2011)). Each TP was compared to its preceding TP for LEfSe analysis. For comparisons with previously published studies, raw sequencing data were obtained from publicly available sources. Experimental data and reference sequences from Yasuda and colleagues (Yasuda et al. (2016)) were merged into a single sequence file, run through the SEPP Fragment Insertion QIIME2 pipeline, and inserted into the Green Genes 13.8 tree. Unifrac distances based on the tree generated were used to run a Principal Coordinate Analysis (PCoA). PICRUSt2 was used to estimate gene functionality in the OTUs. The QIIME2-2021.2 version of the QIIME2 PICRUSt2 plug-in was utilized with EPA-NG sequence placement and “mp” hidden-state prediction. The data were analyzed using ALDEx2, which creates a centered log-ratio data distribution. T-tests with Benjamini-Hochberg correction were used for paired TP analyses.

To identify the species of bacteria from the significantly perturbed clades in the 16S analysis of our six-monkey infected cohort, we performed whole metagenomics shotgun sequencing on rectal swab DNA from four malaria-naive control rhesus. DNA from rectal swabs was isolated using the QIAGEN DNeasy PowerSoil Pro Kit, according to the manufacturer’s protocol. Extracted DNA samples were quantified using Qubit 4 fluorometer and Qubit™ dsDNA HS Assay Kit (Thermofisher Scientific). Between 10 and 20 million high quality reads were generated for each sample. Unassembled sequencing reads were directly analyzed by CosmosID-HUB Microbiome Platform (CosmosID Inc., Germantown, MD) described elsewhere (Hasan et al. (2014); Lax et al. (2014); Ottesen et al. (2016); Ponnusamy et al. (2016)) for multi-kingdom microbiome analysis and profiling of antibiotic resistance and virulence genes and quantification of microorganisms relative abundance. Briefly, the system utilizes curated genome databases and a high performance data-mining algorithm that rapidly disambiguates hundreds of millions of metagenomic sequence reads into the discrete microorganisms engendering the particular sequences. To identify the most abundant bacterial species in the gut microbiome of these control animals, we determined the relative abundance of all bacteria in the samples and calculated the relative abundance of each one in each sample. We then took the average relative abundance of each microbe, averaged across the four control animals, and ranked them in order of highest abundance. The top twenty most abundant bacterial species were determined to be the species most highly represented in the uninfected healthy rhesus macaque gut microbiome.

Plasma extraction and liquid chromatography mass spectrometry (LC-MS) analysis procedures were performed as described previously (Cordy et al. (2019)). Briefly, 50 μ l of plasma was spiked with 2.5 μ l of stable-isotope-labeled internal standards and 100 μ l of acetonitrile to precipitate protein. The clean extract was collected after centrifuging the plasma mixture at 14,000 g for 10 minutes at 4°C. 10 μ l was injected in triplicate into a Thermo Fisher Scientific Q Exactive HF high-field mass spectrometer using an autosampler maintained at 4°C. The sample order was randomized prior to LC-MS sample preparation, and then the samples were prepared for analysis in sequential batches. The metabolites were chromatographically separated using HILIC columns (Thermo Fisher Scientific Accucore 50 × 2.1 mm) with a 5-minute formic/acetonitrile gradient. Electrospray ionization was used in the positive-ion mode. Data quality was monitored by injecting 10 μ l of quality control samples (NIST SRM 1950 and internal standards) after every 20 samples. Raw data were preprocessed using apLCMS (Yu et al. (2009)) and xMSanalyzer (Uppal et al. (2013)) to extract retention time, m/z, and intensity information. Throughout the manuscript, an m/z feature refers to a unique combination of m/z and retention time. The preprocessed data were further treated to correct for batch effects using ComBat (Johnson et al. (2007)). The metabolic features were annotated and identified using xMSannotator (Uppal et al. (2013)). The untargeted LC-MS host plasma metabolomics data were analyzed using the R package xmsPANDA (Uppal (2021)). xmsPANDA output was then run through an annotation python script using MS/MS information from standardized metabolites that were analyzed on the same LC-MS machinery as the experimental samples. This annotation program had an error range of 30 seconds for retention time and 0.02 for mass to charge ratio. The median of three LC-MS replicates was used to quantify metabolite data. The LC-MS intensity data were used for plotting the values of kynurenine and tryptophan metabolites. In the longitudinal analysis shown in values below the limit of detection were represented as 0. For the ratio analysis of tryptophan and kynurenine, values under the limit of detection were analyzed as half of the minimum value seen at that TP. A paired student’s t-test was utilized for the ratio analysis.

To identify and catalog gut microbiota during malaria infections, we sequenced and analyzed 1,862,024 Illumina 16S rRNA reads originating from rectal swabs/feces samples from six na¨ıve rhesus macaques acquired at baseline and predefined TPs during a longitudinal P. cynomolgi infection study (Joyner et al., 2019). Sequence analysis at the phyla level showed that Firmicutes, Bacteroidetes, Proteobacteria and Spirochaetes were the four major phyla in all subjects at any given TP, accounting for >90% of the total bacterial communities in all samples ( Table S1 ). Comparing all the TPs ( Figure 1A ), the difference between baseline (TP1A) and peak infection (TP2) was noteworthy. Proteobacteria increased dramatically during peak infection by TP2 while Firmicutes, Bacteroidetes, Spirochaetes amongst others decreased compared to baseline and other TPs as well ( Figure 1B , Table S1 ). Members of the family Helicobacteraceae, which lies within the Epsilon subdivision of the Proteobacteria, were significantly greater in numbers during peak infection compared to the baseline microbial community structure ( Figure 1B , depicted in blue). Members of Ruminococcaceae belonging to phylum Firmicutes, class Clostridia decreased considerably by TP2 ( Table S2 ). Another notable family, Prevotellaceae belonging to the phylum Bacteroidetes, declined considerably by TP2 ( Figure 1B , Table S2 ). ANOVA at the family level revealed Helicobacteraceae were significantly increased by TP2, as shown in Figure 5 . Conversely, selected families such as Streptococca ceae, Lactobacillaceae, Prevotellaceae, and Lachnospiraceae were statistically lower in numbers during peak infection ( Figure 5 ). Alpha diversity metrics support the above results, affirming lowered diversity and richness, as well as evenness during the peak phase of infection ( Table 1 ; Figure 2 ).

We obtained a cumulative 37,727 OTU’s in all 42 samples from our six-monkey cohort, after clustering at a 97% similarity level ( Table 1 ). A good sequencing depth for all samples was indicated by the rarefaction curve (subsampled at the lowest sample size) reaching asymptote as the number of sequences increased ( Figure S1 ). The alpha diversity indices were calculated to evaluate the overall microbiota richness and diversity across the seven experimental TPs. The averages of Chao 1 richness indices for all subjects at each TP were calculated to estimate species richness in terms of the number of bacterial species present. To assess the evenness and relative abundance of the gut bacterial communities, we calculated the Shannon and Simpson diversity indices (Kim et al. (2017)). Based on the alpha diversity metrics, the overall fecal microbiota richness and diversity decreased considerably during peak infection at TP2 compared to baseline and relapse stages of infection ( Table 1 , Figure 2 ). This decrease suggests there is a loss of microbial diversity at peak malaria leading to host gut dysbiosis. To evaluate the similarity of gut microbial communities between the two most distinct TPs, 1A and 2, we examined the beta diversity using a PCoA plot based on weighted Unifrac distances as shown in Figure S2 . The PCoA plot corroborated our previous results confirming that the gut microbiota at peak malarial infection was distinct from the baseline and relapse microbiota. The alteration in gut microbiome between the TPs in this study was analyzed using LEfSe (LDA score >2, p <0.05). The LEfSe results support the findings of increased Helicobacteraceae and decreased Firmicutes at TP2 ( Figure 3B ). This analysis also revealed increased Lactobacillaceae at TP3A ( Figure 3D ), the first relapse TP. Lactobacillaceae is often associated with gut health Turroni et al., 2014. Changes in taxa abundance were observed overwhelmingly in the primary infection compared to the relapse infection. No significant changes were detected between TPs 3A and 4 (Early relapse and peak relapse) nor TPs 4 and 5 (peak relapse and relapse resolution). Interestingly, of the 28 taxa that decreased between TP1B (pre-peak) and TP2 (peak) 19 of those increased between TP2 (peak) and TP3 (post-peak). These increases indicate a partial restoration of the gut microbiome post-peak ( Figure 3C ). The increases of Helicobacteraceae, an anaerobic, mucosal family, and the family Tissierellaceae (of which most members are anaerobic or aerotolerant), at peak infection. The increases of Helicobacteraceae, an anaerobic, mucosal family, and the family Tissierellaceae (of which most members are anaerobic or aerotolerant), at peak infection may indicate damage to the endothelium of the gut ( Figure 3B ). The results of this analysis indicated that the larger changes in mucosal versus luminal bacteria may be of interest.

We compared our results with previously published data (Yasuda et al. (2016)) on the biogeography of the intestinal mucosal and lumenal microbiome in rhesus macaques using PCoA plots. As shown in Figure 6 , we found that the gut microbiota samples at TP2 clustered together with the reference mucosal microbiome samples from the previous study. However, the other TPs did not seem to have a notable correlation with either the mucosal or lumenal microbiome samples. These results suggest that, at peak malaria, there is an increase of mucosal derived taxa such as Helicobacter spp.

In addition, we analyzed the oxygen requirement profiles of the major bacterial clades at all TPs. The three predominant categories observed were the microaerophiles, obligate anaerobes, and facultative anaerobes ( Figure S3 ). Interestingly, the relative percentage abundance of microaerophiles increases at TP2 (Solnick and Vandamme (2001); Gueneau and Loiseauxdegoer (2002)), coinciding with the peak levels of P. cynomolgi, which is also a microaerophile.

To determine the identities of the most abundant bacterial species perturbed in the gut microbiome, we performed whole metagenome sequencing on rectal swab samples collected from four control rhesus macaques. The bacterial species found in highest abundance in the Helicobacteraceae family was Helicobacter macacae, a known commensal microbe of rhesus macaques, previously identified in Yasuda et al. and other reports, but typically shown to be adherent to the mucosa and not present in high numbers in the luminal contents (Yasuda et al. (2016)). Additional bacterial species identified included members of the Prevotellaceae family (Prevotella copri) and the Lactobacillaceae family (Lactobacillus johnsonii,Ligilactobacillus animalis, and Limosilactobacillus reuteri) ( Figure S4 ).

Untargeted metabolomics was performed on plasma collected from venous blood taken from the same animals at the same TPs in which rectal swabs were collected. Using LC-MS and a high resolution metabolomics workflow (Cordy et al. (2019)), dynamic shifts in multiple metabolites were shown to occur during the peak infection, at TP2. Between TP1A and TP2, 95 metabolites were found to be significantly changed in abundance ( Figure 4A ). Using previously published methodology involving the comparison of mass:charge (m/z) ratios and retention times for reference compounds run on the same machine, we were able to annotate seven of them ( Table 2 ). Among the seven, two members of the kynurenine pathway, acetyltryptophan and hydroxykynurenine, were identified. Acetyltryptophan peaked at TP1B and dropped significantly at TP2 while hydroxykynurenine peaked at TP2 and dropped significantly at TP3A ( Figures 4B–D ; Table 2 ). Other metabolic changes measured included a decline in hypoxanthine and an increase in multiple fatty acyl carnitines during peak parasitemia.

A software program for predicting the functions of microbes within a microbial community based on 16S data alone was used to assess whether microbes could have any involvement in the metabolic changes occurring in the host. Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) (Douglas et al. (2020)) was used to estimate active metabolic pathways of the gut microbiome at each TP. This revealed 124 significant MetaCyc pathways between TP1A and TP2 (t-test, p ≤ 0.05), and 592 significant individual predicted enzymes classified by Enzyme Commission number (EC) between TP1A and TP2. ECs connected to tryptophan metabolism as described previously (Kaur et al. (2019)) were specifically analyzed to investigate ties to the host metabolic shift seen at TP2, but all enzymes were either not present in our data or were insignificant. However, the pathway analysis did reveal a significant increase in the L-tryptophan biosynthesis pathways at TP2 (MetaCyc ID : TRPSYN-PWY) and its superpathway (MetaCyc ID: COMPLETE-ARO-PWY), suggesting the possibility of increased capacity for producing tryptophan by the gut microbiota of the macaques while experiencing peak parasitemia. A Spearman correlation was used to evaluate a possible association between the prevalence of the tryptophan biosynthesis pathway in the gut and host tryptophan and kynurenine levels during primary infection (TP1A-TP3). This analysis revealed host kynurenine and tryptophan levels to be highly correlated with gut tryptophan production, although not significantly correlated with each other ( Figure 4E ).