12 Indian Rhesus Macaques (RM), consisting of 8 geriatric animals (age 19-24 years) and 4 young adult animals (age 3-4 years) were included in this study. Animals were housed and blood/tissue/stool samples were collected at the New Iberia Research Center (NIRC) at the University of Louisiana at Lafayette. Detailed characteristics of animals used including birth date, gender, relatedness and housing information are presented in the Supplementary Table 1 . All animals used in this study are specific pathogen free (free of STLV, SIV, SRV, Herpes B, and Foamy virus) except for RRV and RhCMV.

All animal experiments were conducted following guidelines established by the Animal Welfare Act and the NIH for housing and care of laboratory animals and performed in accordance with institutional regulations after review and approval by the Institutional Animal Care and Usage Committees (IACUC, # 2016-8798-055) of the University of Louisiana at Lafayette. All efforts were made to minimize suffering and stress. Appropriate procedures were performed to ensure that potential distress, pain, discomfort, and/or injury were limited to that unavoidable in the conduct of the research plan. Ketamine (10 mg/kg) and/or telazol (4 mg/kg) were used for collection of blood and stool samples, and analgesics were used when determined appropriate by veterinary medical staff. This study was carried out in compliance with the ARRIVE guidelines (https://arriveguidelines.org/resources/author-checklists).

Plasma biomarkers of immune activation and inflammation were determined using either an ultra-sensitive, multiplexed Milliplex magnetic bead-based assay panel (EMD Millipore) acquired on a MAGPIX instrument (Luminex Corporation) or commercially available. Multiplex assay was used in the measurement of plasma levels of IL-2, IL-6, IL-8, IL-10, IFN-γ, and TNF-α. Commercial ELISA kits were used in the measurement of plasma levels of, hsCRP (Life Diagnostics), and Neopterin (IBL international).

PBMCs were thawed and rested overnight and stained for markers specific for the identification of CD4 and CD8 T cell subsets using Rhesus specific Abs against CD3, CD4, CD8. CD28, CD95 and CCR7 along with live dead aqua dye and acquired on BD Fortessa flow cytometer and analyzed using FlowJo software.

Fresh fecal samples were collected from anesthetized animals by insertion of a fecal loop into the rectum and swiping several times and stored at -80°C. The stool samples from both the study groups were collected at the same time in July 2017 and microbiome sequencing was performed for all the samples together. Microbial DNA was isolated from fecal samples using DNeasy Powersoil Kit (Qiagen) and quantified using Qubit DNA analyzer (Invitrogen/Life technologies). DNA samples were sent to the University of Minnesota Genomic Center for whole genome sequencing of the microbial profile. MetAMOS or PartekFlow softwares were used for microbial profiling.

Shotgun metagenomics library was constructed with the Nextera DNA sample preparation kit (Illumina, San Diego, CA), as per manufacturer’s specification. A 94 dual-indexed Nextera XT libraries were prepared for this study. All libraries combined into a single pool and sequenced in a single lane of a HiSeq 2500 HO, 125-bp PE sequencing run using v4 chemistry.Generated ≥ 220 M reads for the lane (an average of 3.5-4 million reads per sample). Mean quality score for all high yield libraries is ≥Q30. Barcoding indices were inserted using Nextera indexing kit (Illumina). Products were purified using Agencourt AMpure XP kit (Beckman Coulter, Brea, CA) and pooled for sequencing. Samples were sequenced using MiSeq reagent kit V2 (Illumina).

Raw sequences were sorted using assigned barcodes and cleaned up before analysis. For assembly and annotation of sequences, MetAMOS pipeline or Partek Flow software (Partek Flow Partek Inc., St. Louis, MO) were used. Alpha and Beta diversity calculations were done using embedded programs within the metagenomics pipeline, or using Stata15 (StataCorp LLC, College Station, TX) or EXPLICET software. Functional profiling was performed using HUMAnN2-0.11.1 (29) with Uniref50 database to implement KEGG orthologies. Sequence files for the samples in this manuscript have been submitted to ArrayExpress with the accession number E-MTAB-10160.

We utilized bacterial-matrix based dimensionality reduction and clustering algorithms. Compared to Qiime-derived PCA, which is done sample-wise, these algorithms are identity agnostic and decipher qualitative association (and disassociations) between experimental groups. We used Graphia software (Edinburgh, UK) and its implementation of Markov Clustering algorithm (MCL) (30). Using our group-wise bacterial matrix, MCL looked for cluster structures using mathematical bootstrapping. While this method is agnostic to phylogenetic hierarchy in the matrix, we used the hierarchy as identifying markers to understand the bacterial clusters and changes in those clusters within the 2 study groups. MCL used the stochastic flow of the matrix to decipher the distances between the bacteria at equilibrium, thereby generating a cluster map by using correlation scores as distance. For generating the cluster, nodes scoring above a Pearson correlation value of 0.85 were used.

Metabolites from stool and plasma fractions were extracted using chloroform, methanol and water mixture to obtain phase separation. The extracted lipids corresponding to 100 μg proteins were used for IROA (described below) and untargeted liquid chromatography Q-Exactive Orbitrap tandem mass spectrometry (LCMS/MS) for profiling. All solvents were LC-MS/MS grade. Methanol (6 ml) and chloroform (3 ml) were added to each sample (corresponding to 100 μg of total protein from each fraction). After 2 min vigorous vortexing and 2 min sonication in an ultrasonic bath, the samples were incubated at 48°C overnight (in borosilicate glass vials, PTFE-lined caps). The following day, 3 ml of water (LC-MS grade) and 1.5 ml of chloroform were added, samples were vigorously vortexed for 2 min and centrifuged at 3000 RCF, 4°C for 15 min to obtain phase separation. Lower phases were collected. Remaining samples (upper and interphases) were reextracted by addition of 4.5 mL of chloroform and centrifugation at 3000 RCF, 4°C for 15 min to obtain phase separation. Lower phases from both extractions were combined and dried in a centrifugal vacuum concentrator. Samples were stored at -20°C and were reconstituted in 150μL of chloroform:methanol (1:1) before mass spectrometric analysis. To determine the impact of the microbiome on host metabolism in our RM model, we used mass-spectrometry based Isotope Ratio Outlier Analysis (IROA) from RM plasma and feces. This pipeline encompasses sample preparation, LC-MS based peak acquisition, proprietary software-based library creation, normalization and quantification of metabolites and we have standardized and validated IROA method in our laboratory (31). For LC-MS/MS, The Q Exactive (Thermo) mass spectrometer was operated under heated electrospray ionization (HESI) in positive and negative mode separately. The spray voltage was 4.4 kV, the heated capillary was held at 310°C (negative mode) or 350°C (positive mode) and heater at 275°C (positive mode). The S-lens radio frequency (RF) level was 70. The sheath gas flow rate was 30 (negative mode) or 45 units (positive mode), and auxiliary gas was 14 (negative mode) or 15 units (positive mode). Full scan used resolution 70,000 with automatic gain control (AGC) target of 1x106 ions and maximum ion injection time (IT) of 100 ms. MS/MS spectra were acquired using the following parameters: resolution 17,500; AGC 1x105; maximum IT 75 ms; 1.3 m/z isolation window; underfill ratio 0.1%; intensity threshold 1x103; dynamic exclusion time 3 s. Normalized collision energy (NCE) settings were: 15, 30, 45, 60, 75, 90 (in positive and negative mode separately; total 12 runs per sample). Spectral file analysis was done with MZmine2 with integrated Metacyc database. Pathway and disease mapping was performed using MetaboAnalyst 4.0 at https://www.metaboanalyst.ca/MetaboAnalyst/home.xhtml. Metacyc analysis classifies the involved pathways as biosynthesis and degradation pathways, which are essentially anabolic and catabolic pathways respectively.

A Mann-Whitney U test was conducted to determine whether there was a difference in the plasma cytokine and immune cell phenotypes of the young adult and geriatric groups. The Benjamini-Hochberg correction for multiple comparisons was implemented to calculate adjusted p-values. The alpha level of .05 was used to determine significance. Bray-Curtis dissimilarity matrix was statistically evaluated using PERMANOVA test with Monte-Carlo permutation (999 tests) and Bonferroni correction for FDR. For correlation analysis, microbial data was reordered according to difference in prevalence between young adult and geriatric animals of a given microbe in the microbiome. The 25 microbes most prevalent in geriatric or young adult animals were then correlated with plasma cytokines using prism software. The list of bacteria and cytokines were derived from a Benjamini-Hochberg FDR corrected list. In order to perform the selection, we lined up the list of taxa to look for significantly changing species in geriatric individuals, compared to young adults and vice-versa (with appropriate FDR correction), we picked taxa which had q<0.01 or below. For one comparison, we had 25 taxa that fit the criteria and we picked the top 25 from the other comparison as well. Correlation analyses were performed using Spearman rank order correlation using prism software. Test of significance for this analysis was a simple 2-tailed t-test for the correlation. Heat maps were generated to show the results and significant correlations. Results with a p value <0.05 were considered significant.

12 RM with 8 geriatric and 4 young adult animals were investigated to elucidate the effects of age on the immune system and microbiome. Analysis of plasma biomarkers depicted in Figure 1 showed significantly higher levels of acute inflammatory protein CRP ( Figure 1A ), neopterin, a systemic marker of immune activation ( Figure 1B ), and various inflammatory cytokines including TNF ( Figure 1C ), IL-10 ( Figure 1D ), IFN-γ ( Figure 1E ), IL-8 ( Figure 1F ), IL-6 ( Figure 1G ), and IL-2 ( Figure 1H ) in geriatric animals. Analysis of absolute numbers and frequencies of circulating immune cell subsets ( Supplementary Figure 1 ) did not show statistically significant differences between groups. However, geriatric animals exhibited a non-significant trend of lower mean absolute numbers of lymphocytes ( Supplementary Figure S1A ) and non-significant trend of higher mean absolute monocytes ( Supplementary Figure S1B ) compared to young adult animals. Within the lymphocytes, geriatric animals showed a non-significant trend of lower mean absolute number of CD4 T cells ( Supplementary Figure S1C ) while mean absolute number of CD8 T cells ( Supplementary Figure S1D ) were similar between young adult and geriatric animals.

In agreement with the trend of lower numbers of CD4 T cells, we also found a non-significant trend of lower mean frequencies of CD4 T cells ( Supplementary Figure S1E ) in geriatric animals. CD4 T cell maturation subsets showed a non-significant trend for lower mean frequencies of naïve (defined as CD4⁺CD28⁺CD95^-) ( Supplementary Figure S1F ) along with a non-significant trend of higher mean frequency of central memory (defined as CD4⁺CD28⁺CD95⁺) CD4 T cells ( Supplementary Figure S1G ) in geriatric animals. CD4 effector memory (defined as CD4⁺CD28^-CD95⁺) subset ( Supplementary Figure S1H ) did not differ between the study groups. Frequencies of CD8 T cells ( Supplementary Figure S1I ) did not differ between the groups with the geriatric animals exhibiting a non-significant trend of lower mean frequencies of naïve ( Supplementary Figure S1J ) and central memory ( Supplementary Figure S1K ) and higher mean frequencies of effector memory ( Supplementary Figure S1L ) cells. Taken together, our data support an age associated increase in inflammation and a trend for alterations in immune cell subsets in aging NHP that is similar to humans.

Stool samples from young adult and geriatric rhesus macaques were analyzed by shotgun metagenomics sequencing ( Figures 2 , 3 ). Group-wise compositional differences were determined by Principal coordinate analysis (PCA) using the Bray-Curtis dissimilarity matrix for β -diversity. PCA analysis showed a clear difference in microbiome composition between young adult and geriatric animals without any overlap as shown in the 2D rendering of the Bray-Curtis dissimilarity matrix (( Figure 2A ). These differences were further visible in the Euclidean and Jansen-Shannon PCA matrices ( Supplementary Figures 2A, B ). While there was a distinct clustering of the young macaque microbiota, the differences were not significant by the PERMANOVA analysis with FDR correction (‘q’ value of 0.57). Diversity changes between the two groups (α -diversity) using three different indices also showed a trend towards age associated decrease in the Shannon index ( Figure 2B ) and increases in Chao1 ( Figure 2C ) and Simpson ( Figure 2D ) indices. An unsupervised hierarchical clustering plot representing each animal is shown in Figure S2C . In the young group, 3 animals clustered well while in the old group, 7 animals clustered together.

Investigation of the bacterial phyla and classes showed that the gut microbiome in the animals was dominated by Proteobacteria and Firmicutes ( Figure 3 ), followed by phylum Bacteroidetes. Young adult animals showed a Firmicutes dominated microbiota while Proteobacteria dominated in geriatric animals together with relative expansion of Bacteroidetes, Tenericutes and Verrucomicrobia ( Figure 3A ). Despite a decrease with age, Firmicutes still remained the second largest phylum in the microbiome of geriatric animals. Classes within Proteobacteria and Firmicutes consisted of Bacilli (Firmicutes) dominating in young adult animals, followed closely by Gammaproteobacteria (Proteobacteria), Clostridia (Firmicutes), and Bacteroidia (Bacteroidetes). In contrast, in geriatric animals, the most dominant microbiome class was Gammaproteobacteria with greatly diminished representation from Bacilli which constituted the second most abundant class ( Figure 3B ). Clostridia and Bacteroidia occupied the third and fourth most dominant position within the geriatric group and their relative abundance was unchanged compared to young adult animals.

At the species level, 55 species were significantly different in the gut microbiome of the young adult and geriatric animals ( Figure 4 ). As expected, most of the species showing diminished relative abundance in geriatric animals belonged to phylum Firmicutes, while members of phylum Proteobacteria exhibited elevated relative abundance. It was clear from the statistical analysis that there were subtle compositional changes in the microbiome between young adult and geriatric macaques with some significant genus-level changes. While the microbiome changes were modest, we set out to understand the overall physiological role of these changes.

To better understand the association between the microbiome and inflammation, a correlation analysis between specific microbial species and plasma inflammatory cytokines was performed. The list of over 600 microbes identified was arranged according to difference in relative abundance between young adult and geriatric animals. Twenty-five microbes representing most abundant microbes in young adult or geriatric animals respectively, were correlated with plasma cytokines ( Figure 5 ). In geriatric animals almost all the microbes correlated positively with plasma inflammatory cytokines, with many microbes showing significant correlations with neopterin, CRP, TNF, IL-2, IL-6, IL-8 and IFN-γ. In young adult animals a trend for negative correlation of microbes was noted with the same cytokines, with neopterin showing a significant negative correlation.

Next, we performed correlation modeling of 100 significantly changing microbial species between the two groups and the cytokine levels in young adult and geriatric macaques. For this, we used a Spearman r correlation matrix (presuming a non-gaussian distribution) with 95% confidence interval. Since we had already selected correlates that were corrected for FDR, we performed a simple two-tailed Mann-Whitney U test for significance for the correlation matrix. As shown in Table 1 , many bacterial species were strongly correlated, positively (red) or negatively (blue) with more than one cytokine in geriatric animals. Interestingly, a majority of these bacterial species were from the phylum Firmicutes that are highly metabolically active and proinflammatory for the host. Supplementary Figure 3 shows representative matrices for CRP, IL6 and Neopterin and Table 1 shows the significant correlations with cytokines and bacterial species in the geriatric animals. The correlation matrix showed a recurring theme for cytokines CRP, IL6 and Neopterin, where increased levels correlated positively with various bacterial species in geriatric macaques (seen as warmer palette), whereas in young adults, less correlation was seen for the same cytokines ( Figure 5 ). Table 1 further shows the significant correlations with cytokines and bacterial species in old animals.

Bacterial associations were calculated using Markov Clustering Algorithm (MCL), which evaluates correlations of vector abundance of bacterial species within each group. The rendering for MCL only includes those with pearson’s r > 0.85 and unidentifiable feature-level sequences were not included for this analysis. This calculation results in clustering of species that are most likely to co-occur together in a network, using correlation values as the distance matrix. This approach also allows for each species to participate in multiple clusters, emulating real life participation of individual microbes in different networks and biofilms. As shown in Figure 6A , the microbiome of young adult animals exhibited a single large super-cluster at microbial homeostasis. However, in geriatric animals, major constraints have been introduced into the network structure with the emergence of satellite clusters and scattering of the super cluster ( Figure 6B ). In geriatric animals, various newly emerging breakout clusters could be seen with known bacterial composition (Representative clusters numbered, and composition shown in Supplementary Figure 4 ). The investigation of age-induced relative changes in the microbial composition to changes in bacterial associations are depicted in Supplement Figure 4 showing that bacteria from cluster 1 in Figure 6B belong to phyla Archaea, Thermotogae, Proteobacteria and Bacteroidetes. A common property of constituents of cluster 1 was that they were all extreme thermophiles. Similarly, in cluster 2, with the exception of Lactobacillus casei, all constituents were known pathobionts (i.e. potential mammalian pathogen). The putative separation/isolation of functionally distinct clusters of microbes to form distinct niches is associated with changes in metabolism and virulence. Hence, in aging, constraints are introduced into the homeostatic microbial networks and bacteria of known traits seemingly cluster together, forming their own niche; this could also result in change in their metabolic profiles.

One of the ways the microbiome can affect the host is through their own metabolism, producing specific metabolites that make their way into the host circulation. Metabolites in the serum of young adult and geriatric animals were analyzed by mass-spectrometry based Isotope Ratio Outlier Analysis (IROA) as previously described (31). Serum from young adult and geriatric animals showed a few metabolites that were significantly different between the two groups ( Figure 7A ). Geriatric animals showed higher levels of hypoxanthine, L-arginine, and an unidentified compound with a chemical composition of C7H11N6O4. Young adult animals showed higher levels of L-nicotinamine, fructoselysine, N-acetylneuraminate, leucyl phenylalanine, and 2’,3’-cyclic uridine monophosphate (UMP). Most of the biosynthetic pathways were found to be upregulated in the geriatric animals ( Figure 7B ). All degradation pathways except for carbohydrate and secondary metabolite catabolism, were higher in geriatric animals ( Figure 7C ).

To understand the metabolic landscape of the microbial compartment in the gut microbiome of the two age groups, we isolated metabolites from the microbial and microbe-free fecal matter. LC-MS/MS analysis of the pure microbial fraction and pathway enrichment analysis showed that glycerolipid, amino acid and secondary metabolite metabolism predominated in both age groups ( Figure 7D ) but with significantly different metabolites as shown in Figure 7E . In the microbe-free fecal compartment, secondary metabolite synthesis, nucleotide synthesis and amino acid degradation remained at synchrony with bacterial and host metabolic status ( Figure 8 ). Several metabolites were significantly different between geriatric and young adult RM ( Supplementary Table 2 ), which were mapped to various human pathways using KEGG orthologies. Macaque chow was used as a control for this experiment and all metabolites present in the chow were excluded from the differential analysis. Significantly different metabolites were mapped on human pathways using MetaboAnalyst 4.0, and significantly impacted pathways were identified for metabolites that were abolished or down-regulated in geriatric macaques ( Figures 8A, B , respectively). Next, we mapped the altered metabolites to Metabolite-disease interaction networks. Several human diseases were identified, centered around alterations in Phosphate ( Figure 8C ), L-Argenine/Oleic acid ( Figure 8D ) and L-Methionine/L-Phenylalanine ( Figure 8E ) metabolism. Many of these diseases are co-morbidities associated with aging.