All animal management and research procedures were approved by the Animal Welfare and Ethical Committee of China Agricultural University (Permit No. DK1318). Experiments were performed in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals (The State Science and Technology Commission of P. R. China, 1988). Fifteen yak bulls were used to analyze the composition of the ruminal bacterial and archaeal communities. The healthy yak bulls had a body weight of 118 ± 17 kg and were 45 ± 5 months old. They were randomly allocated to three feeding regimes: grazing (graze, n = 6), grazing and supplementary feeding (GSF, n = 5), or indoor feeding (feed, n = 4). The animals allocated to graze regime group grazed in their natural habitat (Mozhugongka county, Lhasa District, Tibet (29°50′13.86″ N, 91°43′46.81 E), about 4,000 m altitude) during summertime. Yaks grazed mixed pastures dominated by Kobresia pygmaea and Stipa purpurea for majority of the time. Grass samples were collected once every month for determining nutritive quality (from June to October, 2013). Samples were clipped from five 0.25-m² quadrats across the pasture to mimic forage selected by grazing yaks. Then samples were pooled, dried in a forced-air oven at 60°C for 48 h, grounded to pass through a 1-mm sieve, and subsequently analyzed for macronutrient composition. The macronutrient composition of the summer-season pasture was (DM basis): 19.5% crude protein (CP), 3.7% ether extract (EE), 58.6% neutral detergent fiber (NDF), 32.5% acid detergent fiber (ADF), and 8.2% ash. The animals allocated to the feed regime group were individually housed in tie stalls and fed ad libitum a mixed forage containing (DM basis): 8.6% CP, 2.9% EE, 68.9% NDF, 35.6% ADF, and 9.45% ash. In addition, the animals were fed 2 kg/day supplemental concentrate with the macronutrient composition (DM basis): 11.1% CP, 5.1% EE, 25.6% NDF, 9.4% ADF, and 11.62% ash. For GSF regime group, yaks were also group-fed daily to provide 2 kg of supplemental concentration when back from pastures at ~1830 h. All yaks had continuous access to high-quality water and trace mineral salt. The whole experiment lasted for 135 days, from June to October, 2013.

All yaks grazing on summer pasture and feed indoor were transported to a commercial slaughter house at the end of the experiment. The animals were stunned and samples (mixture of liquids and solids) were taken from the dorsal, central, and ventral regions of the rumen. Rumen samples were pooled and strained through four layers of cheesecloth to obtain rumen liquids. Aliquots of the pooled rumen liquid samples (~50 mL) were placed on ice, immediately transferred to the laboratory, and stored at −80°C until analysis.

Metagenomic DNA was extracted from 0.5 g of each homogenized ruminal semi-fluid sample using the repeated bead beating plus column purification (RBB+C) method (Yu and Morrison, 2004) and an oscillator (Precellys 24, Bertin Technology, France). The rotating speed of the oscillator was 5,500 rpm with two circulations (30 s per circulation). The DNA quality was assessed via agarose gel (1%) electrophoresis; metagenomic DNA concentrations were determined with a NanoDrop 2000 Spectrophotometer (Thermo Scientific, USA).

The PCR amplification of the bacterial 16S rRNA V3 region was performed with the primer set 341F/518R (forward primer 341F: 5′-CCT ACG GGA GGC AGC AG-3′; reverse primer 518R: 5′-ATT ACC GCG G CT GCT GG-3′; Zhang et al., 2010; Whiteley et al., 2012; Zhao et al., 2015; Paz et al., 2016). The 5′-end of the reverse primers were fused to an Ion A adaptor plus key sequence and a sample barcode sequence; the forward primers were fused to a truncated Ion P1 adapter sequence. PCR conditions and amplification reaction system have been previously described (Zhao et al., 2015). Amplicons were examined on 2% E-Gel Size SelectTM Agarose Gels and purified with Agencourt AMPure XP Reagent. Library sizes and molar concentrations were determined with a Agilent 2100 Bioanalyzer™ with the Agilent High Sensitivity DNA Kit (Agilent Technologies, Inc., Santa Clara, CA). Emulsion PCR was performed using the Ion OneTouch™ 200 Template Kit v2 DL (Life Technologies, Inc.) according to the manufacturer's instructions. Sequencing of the amplicon libraries was performed on a 318 chip with the Ion Torrent Personal Genome Machine (PGM) system using the Ion PGM™ Sequencing 300 kit (Life Technologies, Inc.).

The archaeal community composition was assessed by sequencing the hypervariable V6~V8 regions of the 16S rRNA gene. DNA was amplified using the universal archaeal primer set Ar915aF/Ar1386R (forward primer Ar915aF: 5′-AGG AAT TGG CGG GGG AGC AC-3′; reverse primer Ar1386R: 5′-GCG GTG TGT GCA AGG AGC-3′; Kittelmann et al., 2013). A unique 10-base error-correcting barcode added to the end of the reverse primer Ar1386R was used for each sample to allow sample multiplexing. The PCR reaction was performed using phusion high-fidelity PCR Mastermix [New England Biolabs (Beijing) LTD., China]. Cycling conditions were: 94°C for 2 min, then 30 cycles of: 94°C for 10 s, 68°C for 20 s, 72°C for 1 min. Amplicons were purified with a Qiaquick PCR purification kit (Qiagen, GmbH) according to manufacturer's instructions. The amplicons from each reaction mixture were pooled at equimolar ratios to generate amplicon libraries. Sequencing was conducted on an Illumina MiSeq 2 × 300 platform according to the protocols described as previous report (Caporaso et al., 2011).

The sequencing data were analyzed using the QIIME pipeline. Sequences were processed for quality control with Fast QC software, and only sequences without ambiguous characters were included in further analysis. FLASH 1.2.7v was used to merge paired-end reads from raw sequencing data (Magoc and Salzberg, 2011). Chimeric sequences were removed using USEARCH sbased on the UCHIME algorithm (Edgar et al., 2011). Samples were excluded from analyses if an insufficient number of sequencing reads were obtained. For bacteria sequencing analysis, one sample was removed from analysis since it had under 14,000 sequences. Thus, 14 samples were used for bacterial community analysis, including four samples from the feed group, six from the graze group, and four from the GSF group. For archaea sequencing analysis, three samples with <8,500 sequencing reads of archaeal 16S rRNA gene were excluded. Thus, 12 rumen liquid samples were used to assess the archaeal community, including four samples from the feed group, three from the graze group, and five from the GSF group. Prior to the calculation of downstream diversity characteristics (i.e., alpha and beta diversity), all samples were subsampled to equal size. The microbial diversity was analyzed using QIIME 1.7.0v (Caporaso et al., 2010) with Python scripts. The sequences were clustered into Operational Taxonomic Unit (OTUs) using the de novo OTU picking protocol with a 97% similarity threshold. Representative sequences of OTUs were aligned to the Greengenes database (version 13_8) for bacterial and archaeal 16S rRNA genes. Alpha diversity analysis (i.e., observed species, Chao1, Shannon-Weiner, and Simpson's indices) were generated and jackknifed beta diversity, including those based on both unweighted and weighted Unifrac distances, were visualized using principal coordinate analysis (PCoA; Lozupone and Knight, 2005). Sequence number, sample coverage, unique OTUs, sample richness, sample diversity, phylum relative abundance, and genus were evaluated using the generalized linear model (GLM) procedure in SAS (version 9.1.3; SAS Institute Inc., Cary, NC, USA). Means were separated by using the Student-Newman-Keuls test (SNK). After the multiple testing analyses, false discovery rate (FDR) adjusted p-values, also called q-values, were computed by using the QVALUE software (version 2.6.0) in R package (version 3.1.0). The q-value of FDR <0.05 represents significant difference in microbe relative abundances between the three feeding regimes.

We generated a total of 468,400 sequences from 14 samples that passed all filtering metrics; the average length across all samples was about 190 bp. The number of sequences per animal ranged from 14,082 (GSF group) to 46,450 (graze group), with mean of 33,457 sequences. A total of 43,525 OTUs were detected from 14 samples at 3% dissimilarity. Each sample had an average of 9,473 OTUs. Richness estimates and diversity indices were developed after normalization to 14,000 sequences (Table S1).

We found 19 distinct phyla across feeding regimes; the most abundant were Bacteroidetes (51.06%) and Firmicutes (32.73%). The top eight phyla, which had high relative abundances and were prevalent in all samples (Figure 1A), accounted for nearly 90% of all sequences. The less abundant phyla included Verrucomicrobia (1.73%), Tenericutes (1.43%), Actinobacteria (0.94%), Proteobacteria (0.76%), TM7 (0.67%), and SR1 (0.47%). The other known phyla with relative abundances <0.4% accounted for 0.66% of all sequence data, while the remaining sequences that were unclassified at the phylum level accounted for 9.54%. We detected 120 distinct genera across all feeding regimes, and 25 genera were present in all samples (Figure 1B); this demonstrates the variability in abundance across samples. The top nine genera had relative abundances above 0.04%. The predominant genera across all samples were Prevotella (15.81%), YRC22 (2.08%), Succiniclasticum (1.87%), Butyrivibrio (1.40%), and Ruminococcus (1.25%). Minor genera, such as CF231, BF311, Clostridium, and Bifidobacterium accounted for 0.81, 0.53, 0.43, and 0.25%, of relative abundance, respectively. The other known genera accounted for 2.74% of sequences, while sequences that were unclassified at the genus level accounted for 72.85% of sequences.

We generated 226,711 archaeal sequences from 12 samples with an average length of 470 bp (Table S2). This was an average of 18,892 sequences, or 851 OTUs per animal. Richness estimates and diversity indices were developed after normalization to 8,500 sequences. Four distinct phyla were found across feeding regimes. Euryarchaeota was the most abundant phylum in all samples, accounting for 99.67% of the all archaeal sequences (Figure 2A). Crenarchaeota accounted for 0.11% of archaeal sequences, while the remaining sequences (0.12% of sequences) could not be classified at the phylum level. We detected 10 distinct genera across feeding regimes, and there were seven genera that were shared by all samples (Figure 2B). The most abundant genus was Methanobrevibacter (91.60%). vadinCA11, Methanoplanus, Methanosphaera, Methanobacterium, Methanimicrococcus, and Methanosarcina accounted for 2.07, 1.65, 1.24, 0.19, 0.04, and 0.02% of archaeal sequences, respectively. The remaining 3.17% of sequences were unclassified at the genus level.

The Shannon-Weiner and observed species indices were similar (q > 0.05) between the three feeding regimens (Figures 3A,B). The similarity in archaeal OTUs between feeding regimes was visualized in a Venn diagram (Figure 3C). A thermal double dendrogram of the 80 most abundant bacterial OTUs illustrated that GSF regime samples could not grouped from other two treatments (Figure S1); while indoor feed regime samples could clearly separate from that of pasture graze regime. The PCoA analysis showed the relationships between the microbiome community structures (Figure 3D). There was no obvious dividing line and the two principal components covered 21.72% of the variation; furthermore, samples from the graze regime group could not be distinguished from samples from the GSF regime group. Interestingly, the microbiome of the indoor feed group was distinct from that of the graze group.

The relative abundance of bacterial taxa were used to describe the impact of feeding regimes on bacterial populations. At phylum level (Table 1), the ruminal microbiome from the graze group had significantly higher relative abundance of Bacteroidetes than the samples from the feed and GSF regime groups (q = 0.02). While there is a lower relative abundance in Firmicutes in graze regime group (q = 0.02). The graze regime group had higher relative abundance of Verrucomicrobia compared with the feed and GSF regime groups (q = 0.04). The feed regime group had higher relative abundances of Tenericutes and TM7 than the graze regime group (q < 0.01 and q = 0.04, respectively). We found no differences in the rest of the phyla (q > 0.05). At genus level (Table 2), the relative abundances of Prevotella in the graze and GSF regime groups were significantly higher than in the feed regime group (q = 0.03). The relative abundance of Ruminococcus in the graze regime group was lower than in the feed and GSF regime group (q = 0.02). No significant differences in community structure were observed in the rest of the genera (q > 0.05). These data suggest that the bacterial community structure differs between the three feeding regimes.

Analyses of archaeal 16S rRNA gene fragments revealed no differences in diversity between the feeding regimes (Figures 4A,B) by the Shannon-Weiner index or the observed species analyses. The similarity in archaeal OTUs between feeding regimes was visualized in a Venn diagram (Figure 4C). The number of archaeal OTUs in the feed regime group was similar to that of the graze and GSF regime groups (Table S2). A thermal double dendrogram of the 80 most abundant archaeal OTUs illustrated that samples in the same regime were not clearly grouped (Figure S2); this implies that the archaeal communities are similar between feeding regimes. Furthermore, the overall archaeal community structure analyzed by PCA analysis based on unweighted UniFrac method revealed no obvious differences between regime groups (Figure 4D). The relative abundance of the most detected genera were similar across feeding regime groups (q > 0.05; Table 3).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.