Sequence processing and analysis

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.