In present research, 26 Mongolian subjects (22–35 years old) started the three-week experiments on day 1 upon arrival in China. All subjects had no gastrointestinal disorder and did not take any antibiotics for three months before the experiment started and until the end of the experiment. During the experiment, dietary staple food was replaced each week (the first week: wheat, the second week: rice and third week: oat), and the amount and species of food intake for each individual was the same every day. The information was shown in the Supplementary Table S5. And the diet was prepared by a restaurant according to the designed recipes. All volunteers were instructed to eat only the provided foods or allowable beverages (water or unsweetened tea). They also confirmed that no remarkable changes occurred in their diet intake and no medication was taken during the experiment. The fecal samples of each individual in 4 time points (from week 0 to week 3) were collected for further sequencing. The study protocol was approved by the Ethical Committee of the Inner Mongolia Agriculture University (Hohhot, China). Sampling and all subsequent steps described in the Materials and Methods have been conducted in accordance with the approved guidelines. After obtaining the written informed consent, we collected habitual long-term dietary information from all participants using a food frequency questionnaire (Supplementary Table S6). All subjects gave written informed consent in accordance with the Declaration of Helsinki.

The QIAamp^® DNA Stool Mini Kit (Qiagen, Hilden, Germany) was used for DNA extraction from the fecal samples. The quality of the metagenomic DNA was assessed by 0.8% agarose gel electrophoresis. All of the DNA samples were stored at –80°C until further processing.

The Illumina HiSeq2500 platform was used for shotgun metagenomic sequencing. Paired-end reads were generated with 100 bp in the forward and reverse directions. The length of each read was trimmed with Sickle. This set of high-quality reads was then used for a further analysis. An average of 10.34 gigabases (Gb) of paired-end reads were obtained for each sample, totalling more than 1 Tb of high-quality data free of host genomic and adaptor contaminants.

The Illumina reads were assembled into contigs using IDBA-UD (Peng et al., 2012) with default parameters. Genes were predicted on the contigs with MetaGeneMark (Zhu et al., 2010). A non-redundant gene catalog was constructed with CD-HIT (Li and Godzik, 2006) using a sequence identity cut-off of 0.95, with a minimum coverage cut-off of 0.9 for the shorter sequences. This catalog contained 1,617,412 microbial genes.

To assess the abundance of genes profile, reads were aligned to the gene catalog with Bowtie2 (Langmead et al., 2009) using the following parameters: -p 12 -x nt -1 R1.fastq -2 R2.fastq -S R.sam. Then, for any sample N, we calculated the abundance as follows:

Step 1: Calculation of the copy number of each gene:

Step 2: Calculation of the relative abundance of gene i:

ai: the relative abundance of gene i

bi: the copy number of gene i from sample N

Li: the length of gene i

xi: the number of mapped reads

We use MetaPhlAn2 (Segata et al., 2012) to produce organism abundance profiling with default parameters, which relied on about 1 million unique clade-specific marker genes identified from about 17,000 reference genomes. Meanwhile, the metagenomic reads were annotated by the NCBI virus database to build the virome profile.

We aligned the amino acid sequences that were translated from the gene catalog against the proteins/domains to the Kyoto Encyclopedia of Genes and Genomes (KEGG) and antibiotic resistance genes (ARGs) databases using BLASTP (e-value ≤ 1e-5 with a bit-score higher than 60). Each protein was assigned to the KEGG ortholog group (KO) or ARG number by the highest scoring annotated hit. CAZymes (Carbohydrate-Active Enzymes) were predicted from amino acid sequences by against to family-specific HMM of CAZymes in dbCAN database (Yin et al., 2012) using Hmmscan program in HMMER 3.0 package (Eddy, 2009).

Differentially enriched KO modules and pathways were identified according to the reporterscores (Patil and Nielsen, 2005) from the Z-scores of individual KOs. Accordingly, the Z adjusted pathway of each KEGG module and pathway was calculated as previously described (Patil and Nielsen, 2005). Then, the Z adjusted pathway was used as the final reporter score for evaluating the enrichment of specific pathways or modules. A reporter score of >2.3 (90% confidence according to the normal distribution) could be used as a detection threshold for significantly differentiating pathways.

All statistical analyses were undertaken using the R software. PCoA and a kernel density distribution analysis were performed using the ade4 (Zapala and Schork, 2006) package of R software. Differential abundance of species, genes and KOs were tested by the Wilcoxon rank sum test, and the significantly different (p < 0.01). The heatmap was built in R using the “pheatmap” package. The correlation between the dietary components, metabolic pathways and relative genera were calculated by Spearman’s rank correlation coefficient and visualized by network in Cytoscape (Version 3.2.1).

The sequence data reported in this paper have been deposited in the NCBI database (Metagenomic data: SRP080787).

The fecal samples for all subjects were collected weekly on 4 time points, at week 0, 1, 2, and 3. Subsequently, microbial DNAs from a total of 104 fecal samples were extracted for metagenomic analysis. By employing the shotgun metagenomic sequencing, we obtained more than 1Tiga base (Tb) of pair-end reads (averagely 61,225,132 high-quality reads for each microbiota). The high-quality reads were assembled into contigs, and the genes were predicted onthe contigs. Then, a non-redundant gene catalog including 1,617,412 genes was constructed. Annotated by HMP database, thus taxonomy profile of each sample was obtained. Meanwhile, the reads were mapped to the collective gene catalog to reconstruct sample-specific gene profile, as well as the profiles for the associated KEGG database Orthologs (KO; Supplementary Table S1).

At the taxonomy level, Principal Coordinates Analysis (PCoA) was performed based on the unweighted UniFrac (Figure 1), weighted UniFrac (Supplementary Figure S1A) and Bray–Curtis distances (Supplementary Figure S1B) at the species level. No significant apparent clustering pattern was identified among samples in different time points, and the fluctuation in gut microbiota was limited by individuals. It is indicated the staple food-induced change was smaller than the inter-individual difference. This conclusion held true when the data was analyzed at microbial functional feature level (Supplementary Figure S2). By extracting and comparing the first principal component of Bray–Curtis distance matrices representing taxonomic and functional features in gut community structure among subjects at each time point, we found that changed staple food rapidly altered microbial community structure and metabolic pathway (Figure 2) for each individual. Then we concluded that the intestinal microbiota could be rapidly altered by staple food, however, the changes was limited at individual’s level. Additionally, in comparing the influence of the three-staple food on gut microbiome, it can be observed from the Weighted UniFrac distance between groups at different time points that the impact of wheat on intestinal microbiota was the largest, followed by rice and oat (Figure 3).

We identified differences in specific species at taxonomy level and metabolic pathway at functional level after the consumption of the respective staple food. As it is shown in Table 1, the species Lactobacillus delbrueckii, Ruminococcus gnavus, Bacteroides vulgatus, and Bacteroides massiliensis decreased significantly but the Bifidobacterium catenulatum, Bifidobacterium bifidum, and Alistipes indistinctus increased significantly when subjects consumed the wheat as staple food for one week. When rice and oat was assigned as the staple food in the second and third week, respectively, Bifidobacterium adolescentis, Bifidobacterium longum, Weissella cibaria, and Rothia mucilaginosa declined sharply in the second week whereas the Bifidobacterium adolescentis risen sharply in the third week. At the functional level, we calculated the reporter Z scores of each pathway or module by using the reporter feature algorithm (Figure 4). A reporter score = 2.3 (90% confidence according to the normal distribution) was used as the detection threshold for significantly differentiating modules or pathways. In particular, we found a significant decline in the microbial biosynthesis of amino acids (including valine, leucine, and isoleucine) in the first week. In contrast, the microbial metabolism of carbohydrate kept growing during experiment. Interestingly, it exhibited a high concordance of changes in microbial capacity for fructose metabolism and glycolysis and the changes in relative content of fructose and glucose in wheat, rice and oat. And these results could be confirmed by changes in metabolic models among different staple food intake group.

Using Wilcoxon sum t-test, we identified species which showed significant difference in each staple food group. We focused on the correlation between functionality and changes in the staple food. The species which showed correlations included Bacteroidesovatus, Faecalibacteriumprausnitzii, Lactobacillus ruminis, Bacteroidesstercoris, and Phascolarctobacterium succinatutens. We determined the correlation among the species and their typical metabolic pathway/model according to Spearman’s rank correlation coefficient (R > 0.4, Figure 5). Interestingly, a general positive correlation was observed between the species above and the phosphotransferase system (PTS) and the Aminoacyl-tRNA biosynthesis, and a negative correlation with the FoxO signaling pathway and Cyanoamino acid metabolism.

All high-quality sequences were annotated by a virus protein database, and then we selected and constructed the phage or bacterial virus profile. A total of 144 species of phages were identified (Supplementary Table S2), and the quantity of phage of Enterobacteria, Bacteriophage, Shigella, Klebsiella, Escherichia, Streptococcus, Salmonella, Bacteroides, Lactococcus, Lactobacillus, Clostridium, Haemophilus, Staphylococcus, Brochothrix, and Bacillus accounted for more than 0.001% of the total microbiota (Supplementary Figure S3A). Meanwhile the correlations among the above phages were determined based on Spearman’s rank correlation (Supplementary Figure S3B), and a general positive correlation was found. Additionally, to test the degree of consistency between the observed clustering patterns for bacterial structure and phage structure, respectively, Procrustes analysis was performed based on the PCA matrix of bacterial genus-level organismal profile and that of genus level phage profile (Supplementary Figure S4). The results revealed a strong correspondence between the two profiles above (P < 0.001, using 10,000 Monte Carlo label permutations), which also indicated a potential synergistic interaction between the bacteria and bacterial virus in human gastrointestinal tract.

As the main component of the staple food in the present study was carbohydrate, we focused on the functional specificity of intestinal microbiota for carbohydrate utilization. The non-redundant gene catalog in our study was annotated by carbohydrate-active enzymes database (CAZy), and a total of 383 carbohydrate metabolism related enzymes were identified (Supplementary Table S3). By applying the Kruskal-Wallis test, we observed significant difference in enzymes among the different staple food groups, and these enzymes were mainly related to the metabolism of glucose (Supplementary Figure S5). Meanwhile, by clustering analysis, we exhibited the correlation profile between the microbial species and the carbohydrate-active enzymes (Supplementary Figure S6), and two apparent clustering patterns could be found in the profile which indicated two groups of intestinal microbiota with preference for specific carbohydrates utilization. Accordingly, we reconstruct the network between the species and their metabolic enzymes for each group (Figure 6). The first group including the species of Alistipe sputredinis, Bacteroides intestinalis, Bacteroides massiliensis, Bacteroides plebeius, Bacteroides salyersiae, Bacteroides vulgutas, Bacteroides uniformis, Butyrivibrio crossotus, Megaspha eraelsdenii, Parabacteroides merdae, and Paraprevotella clara, and they were closely related to the carbohydrate-active enzymes of arabinofuranosidase, pectin lyase, polygalacturonaseand xylanase. The second group including the species of Alistipesonderdonkii, Bacteroides fragilis, Bifidobacterium adolescentis, Bifidobacterium breve, Bifidobacterium pseudocatenulatum, Lactobacillus ruminis, Ruminococcus gnavus, and Streptococcus thermophilus, and they were closely related to the carbohydrate-active enzymes of beta-glucosidase, mannosidase, trehalose phosphorylase, fucosidase, and sialidase.

The results above based on the impacts of changes in staple food on intestinal microbiome offered us an excellent opportunity to reveal the linear relationship between intestinal microbial phylum and main nutrient elements, which was also the basic mechanism related to the interaction between the diet and gut microbiome. According to the fitted curves, we observed a declined in phylum Actinobacteria but an increase in phylum Firmicutes with the gradually increasing non-digestible carbohydrates (dietary fiber) content in diet (Figure 7A). At the same time, the phylum Bacteroidetes and Proteobacteria presented a positive correlation with the dietary fat and digestible carbohydrate, respectively (Figure 7A).

To further understand the dietary effects on human microbiome, we annotated our non-redundant gene catalog with the ARGs database, and constructed the ARGs profile of each individuals. By merging and classifying, 290 annotated ARGs were clustered into 53 antibiotic catalog (Supplementary Table S4). We found that the ARGs related to vancomycin were most abundant, followed by the ARGs related to macrolide, tetracycline, bacitracin, tigecycline, penicillin, chloramphenicol, lincosamide, and aminoglycoside, which gene abundance were more than 0.001%. To describe the relationship among the dietary nutrition, microbial phylum and ARGs above, a Circos on line diagram was constructed (Figure 7B).