Age-related changes in microbial composition and function in cynomolgus macaques

General characteristics of the recruited monkeys

Sixteen healthy monkeys were included in this study under the same feeding and living conditions. According to age, these monkeys were divided into three groups: young, adult and old. The detailed characteristics of the recruited monkeys are shown in Supplementary Table 1.

Similar within-sample microbial diversity among the three groups

To characterize the gut composition of monkeys among the three age groups, 16S rRNA gene sequencing was initially used. In total, we identified 905,921 high-quality reads, ranging from 47,691 to 66,725 per sample, with an average read length of 439.06 bp (439.06 ± 2.88 bp). A Venn diagram showed that three age groups shared 677 of the 969 operational taxonomic units (OTUs, defined based on 97% sequence similarity), whereas 33, 42 and 35 OTUs were unique to young, adult and old monkeys, respectively (Figure 1A). At the phylum and family levels, the relative abundance of microbes in the three groups is presented in Figure 1B, 1C. Our findings revealed that the microbial composition of cynomolgus monkeys was dynamically changed with age. The phyla Firmicutes, Bacteroidetes and Proteobacteria dominated in the gut microbiota of monkeys. Prevotellaceae and Ruminococcaceae represented the most abundant families among the three groups, but showed a downward trend with age.

The α-diversity, including microbial community richness (Chao, Ace) and diversity (Shannon, Invsimpson), was compared among the young, adult and old groups. A downward trend was observed in these indexes with increased age, although no statistical difference was detected (Wilcoxon rank-sum test, all p values>0.05) (Supplementary Figure 1A–1D).

Significant differences in the microbial composition among the three groups

To further explore whether the microbial composition of the young, adult and old groups differed significantly, β-diversity analysis was carried out. Firstly, at the OTU level, a 3-D principal component analysis (PCA) plot displayed that there was a discriminative trend among the three groups, but no statistical difference was detected (PERMANOVA, all p values>0.05) (Figure 2A). Interestingly, we found that samples from the adult group were distributed in the central region, whereas samples from the young and old groups were distributed on both sides. This finding suggested that gut microbial composition may change dynamically with age. Moreover, we found that the gut microbiota of the adult group showed more variation than the young group, which was consistent with previous reports in humans [12]. To further display the differences of the bacterial communities among the three groups, the partial least squares discriminant analysis (PLS-DA) was performed, the bacterial communities of the three groups clustered separately (Figure 2B), indicating distinctive fecal microbial communities among three groups.

To identify the profiles of gut microbiota at different age phases, LEfSe analysis was performed. This analysis identified 148 different OTUs responsible for this discrimination (Figure 3A, Supplementary Figure 2A). With increased age, the relative abundances of Veillonellaceae (p=0.014), Coriobacteriaceae (p=0.022) and Succinivibrionaceae (p=0.009) were significantly increased, whereas Ruminococcaceae (p=0.003) and Rikenellaceae (p=0.009) were significantly decreased with age at the family level (Figure 3C, Supplementary Figure 3). Moreover, among these differentially represented microbes, we found that Ruminiclostridium 9 was the only OTU that dynamically changed in all three groups (Figure 3B).

Alternations of arginine and purine metabolism as well as changes in microbial peptidoglycan with age

To characterize the functions encoded by the gut microbial DNA, we performed whole-genome shotgun sequencing of stool samples obtained from the three groups. A total of 1,602,244,870 filtered reads and 45,088,293 open reading frames (ORFs) were used for functional annotation in the KEGG and CAZy databases. Using linear discriminant analysis (LDA) effect size (LEfSe) analysis, we identified 44 KEGG enzymes responsible for discrimination among the three groups (LDA>2.0). Pathway-enrichment analysis revealed that these differential enzymes were mainly involved in amino acid metabolism (especially arginine metabolism) and nucleotide metabolism (purine metabolism), as well as microbial polysaccharides including lipopolysaccharide (LPS), glycosaminoglycan (GAG) and peptidoglycan (PGN) (Figure 4A, Supplementary Figure 2B). The marked enzymes involved in arginine and purine metabolism, as well as microbial polysaccharide are presented in Figure 4B. Our data showed that arginine biosynthesis and purine metabolism were enriched in the adult and old groups, respectively. For polysaccharide metabolism, GAG and PGN were significantly enriched in the young group, whereas LPS was enhanced in the old group.

Age-related microbial functions in CAZymes

To further uncover the function of the gut microbiota, the differential CAZymes among the three groups were analyzed. In total, we identified 21 CAZymes responsible for the discrimination among the three groups (LDA>2.0) (Figure 5A, Supplementary Figure 2C). These CAZymes were mainly involved in the utilization of plant and animal carbohydrates, and microbial polysaccharide metabolism (Figure 5A, Supplementary Figure 4). The gut microbiota of the adult group had a greater utilization capacity for mucin rather than plants compared with the old group (Figure 5B). Interestingly, we found that the CAZymes related to PGN were also significantly changed (Figure 5C), which was consistent with the altered KEGG pathways. Moreover, analysis of CAZyme profiles showed that the microbiota of the adult group had greater functional capacity for the utilization of starch than the old group, and PGN utilization occurred predominantly in the old group.

Site and sample collection

Our study was conducted in Zhongke Experimental Animal Co., Ltd. (hereafter referred to as "Zhongke"), which is a breeding base for cynomolgus monkeys in Suzhou, China (E 31°07'03" to 31°07'06", N 120°19'08" to 120°19'15"). The housing conditions and animal care procedures were detailed in a previous report [29] and were in accordance with Chinese regulatory requirements and accredited by the Association for the Assessment and Accreditation of Laboratory Animal Care International (AAALAC). In brief, monkeys in Zhongke were provided standard sanitation, an adequate and regular diet, and a stable social structure. All procedures involving non-human primates were approved by the Animal Care and Use Committee of Chongqing Medical University and were in compliance with the Guide for the Care and Use of Laboratory Animals [29].

All 16 cynomolgus macaques, as listed in supplemental data, were randomly selected from different enclosures. In accordance with previous reports [30], the 16 macaques were assigned into three age groups: young (2–4 years), adult (5–15 years) and old (17–20 years). All macaques were housed in free enclosures without manual intervention measuring 8 × 3 × 3 m (L × W × H), were provided with water ad libitum, and were fed daily with fresh fruit, vegetables and compound high-nutrition monkey food [29]. Collection of fecal samples occurred immediately after the first defecation using a germ-free device, and samples were stored in liquid nitrogen before being transported on ice to maintain a temperature chain of < −80°C.

16S rRNA gene sequencing and data processing

DNA extraction and amplification of the 16S rRNA gene were performed as previously described [31]. Briefly, DNA was extracted from stool samples using the E.Z.N.A® DNA kit (Omega Bio-Tek, USA). The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified with primers 338F (5ʹ- ACTCCTACGGG AGGCAGCAG-3ʹ) and 806R (5ʹ-GGACTACHVGGG TWTCTAAT-3ʹ) by the thermocycler PCR system (GeneAmp 9700, ABI, USA). The PCR conditions were 3 min of denaturation at 95°C, followed by 27 cycles of 30 s at 95°C for denaturation, 30 s for annealing at 55°C, and 45 s for elongation at 72°C, and a final extension at 72°C for 10 min. The PCR was performed in triplicate in a 20 μL mixture containing 4 μL of 5 × FastPfu buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu polymerase and 10 ng of template DNA. The resulting PCR products were separated on 2% agarose gels, purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified using QuantiFluor™-ST (Promega, USA) according to the manufacturer’s protocol. Purified amplicons were pooled in equimolar concentrations and paired-end sequenced (2 × 300) on an Illumina MiSeq platform (Illumina, San Diego, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (Accession Number: SRP218284).

Raw fastq files were demultiplexed and quality-filtered using QIIME (version 1.9.1, http://qiime.org/). The 250 bp reads were truncated at any site of more than three sequential bases receiving an average quality score of <20. Low-quality reads (reads shorter than 50 bp or barcode/primer errors or reads with a quality value <20) were removed.

Operational taxonomic units (OTUs) were clustered with 97% similarity cutoffs using UPARSE (version 7.1, http://drive5.com/uparse/) with a novel ‘greedy’ algorithm that performs chimera filtering and OTU clustering simultaneously. The taxonomy of each 16S rRNA gene sequence was analyzed by the RDP classifier algorithm (http://rdp.cme.msu.edu/) against the Silva (SSU128) 16S rRNA database using a confidence threshold of 70%.

Shotgun metagenome sequencing and data processing

The fecal samples were further investigated by metagenomic sequencing. The concentration of the extracted DNA was determined by the TBS-380 method, its purity was measured with a NanoDrop 2000 spectrophotometer, and its quality was confirmed by electrophoresis on a 1% agarose gel. The DNA was then fragmented to an average size of about 300 bp using Covaris M220 (Gene Company Limited, China) for paired-end library construction. A paired-end library was constructed using the TruSeq™ DNA Sample Prep Kit (Illumina, San Diego, CA, USA). Adapters containing the full complement of sequencing primer hybridization sites were ligated to the blunt-end of fragments. Paired-end sequencing was performed on the Illumina HiSeq4000 platform (Illumina Inc., San Diego, CA, USA) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using the HiSeq 3000/4000 PE Cluster Kit and the HiSeq 3000/4000 SBS Kit according to the manufacturer’s instructions (www.illumina.com). Sequence data associated with this project have been deposited in the NCBI Short Read Archive database (Accession Number: SRP218382).

Adapter sequences were stripped from the 3' and 5' ends of the paired end Illumina reads using SeqPrep (https://github.com/jstjohn/SeqPrep). Sickle (https://github.com/najoshi/sickle) was used to remove the low-quality reads (length <50 bp or with a quality value <20 or having N bases).

The metagenomics data obtained were assembled using MEGAHIT [32] (https://github.com/voutcn/megahit). The final assembly containing contigs of 300 bp or more, was used for further gene prediction and annotation.

The prediction of ORFs from each assembled contig was performed using MetaGene software [33] (http://metagene.cb.k.u-tokyo.ac.jp/). The predicted ORFs with lengths of 100 bp or over were retrieved and translated into amino acid sequences.

All genes predicted to have 95% sequence identity (90% coverage) were clustered using CD-HIT software [34] (http://www.bioinformatics.org/cd-hit/), and representative sequences containing the longest sequences from each cluster were used to construct non-redundant gene catalogs. High quality reads were mapped to the representative sequences with 95% identity using SOAPaligner [35] (http://soap.genomics.org.cn/) to evaluate gene abundance in each sample.

KEGG annotation was conducted using BLASTP (Version 2.2.28+) against the Kyoto Encyclopedia of Genes and Genomes database [36] (http://www. genome.jp/keeg/) with an e-value cutoff of 1e-5. A total of 362,971,704 enzymes were annotated by the KEGG database. Carbohydrate-active enzyme annotation was conducted using hmmscan (http://hmmer.janelia.org/search/hmmscan) against the CAZy database version 5.0 (http://www.cazy.org/) with an e-value cutoff of 1e-5. A total of 60,024,802 enzymes were annotated by the CAZy database.

Statistical analysis

The α-diversity indexes were assessed according to species richness (Ace and Chao), species evenness (Shannon) and species diversity (Simpson). Beta diversity was assessed with package ‘vegan’ in R (version R-3.3.1), and generated on the basis of principal component analysis (PCA). Descriptive modelling and discriminative variable selection were evaluated by partial least squares-discriminant analysis (PLS-DA). PERMANOVA was performed to identify differences in β-diversity among three groups. The key bacterial taxa, and the CAZymes and KEGG categories responsible for discrimination among the groups, were identified using linear discriminant analysis effective size (LEfSe; http://huttenhower.sph.harvard.edu/galaxy/root?tool_id=lefse_upload) [37]. Only LDA values >2.0 at a P value <0.05 were considered significantly enriched.

Statistical analyses were conducted using the software SPSS, R package, and plots were generated from R and GraphPad Prism version 8.0. The strategy of multiple cooperation was one-against-all when performing LEfSe. One-way ANOVA or the Kruskal–Wallis H test were performed to assess the alpha diversity and abundance of some species among the three age groups, and the two-tailed Student’s t-test or Wilcoxon rank-sum test were performed to determine differences between each set of two groups. The correlation between microbial abundance and age was tested by Pearson correlation analyses and adjusted by partial correlation analysis to exclude the confounders for controlling the false discovery rate (FDR). Adjustments were performed with SPSS (version 22.0).