All management and experimental procedures were approved by the Animal Care Committee, Institute of Subtropical Agriculture, the Chinese Academy of Sciences, Changsha, China.

Thirty-two Liuyang black goat kids were separated from their dams after birth, and maintained within individual pens for the duration of the study. From birth to day 20, kids consumed goat milk and 4 kids were slaughtered at day 0 and 7 (non-rumination phase), respectively. The remaining 24 kids were randomly allocated to two groups based on different feeding types, supplemental feeding (S) and grazing (G). Between day 20 and 40, kids in S group consumed milk, forage plus concentrate, while kids in G group consumed milk and forage. The animals were weaned at day 40. From day 40 to 70, S kids were provided with concentrate and forage, while G kids grazed 8 h on pasture daily, and returned to pens overnight without receiving any extra feed. In both groups, four kids were slaughtered at each of the following ages: day 28, 42 (transition phase) and day 70 (rumination phase). All kids had free access to water. Detailed feeding management, and ingredients of concentrate and forage have been described in our previous parallel study (Jiao et al., 2015b).

After the kids were slaughtered, ruminal digesta samples were collected from five locations (anterior dorsal, anterior ventral, medium ventral, posterior dorsal, and posterior ventral locations) within the rumen (∼10 g each). The samples were composited, and aliquots were stored at -80°C for subsequent DNA extraction and sequencing analysis. Ruminal mucosa samples were collected from the same sites where digesta samples were collected, and rinsed three times with sterile phosphate buffered saline (PBS, pH = 7.4). Afterward, the mucosa samples were separated from the underlying muscular layer, cut into 3 mm × 3 mm fragments, composited, frozen in liquid nitrogen and stored at -80°C for later molecular analysis.

Genomic DNA was extracted using the QIAamp DNA Stool Mini kit (Qiagen GmbH, Hilden, Germany) and quantified using NanoDrop ND1000 (NanoDrop Technologies, Wilmington, DE, United States). Conventional PCR was conducted to amplify bacterial V2–V3 and archaeal V6–V8 regions of the 16S rRNA genes, as well as fungal V3–V5 and protozoal V4–V5 regions of the 18S rRNA genes using primers detailed in Supplementary Table S1. The PCR conditions were as follows: pre-denaturation at 94°C for 3 min, 35 cycles of denaturation at 94°C for 30 s, annealing at 58°C for bacteria (62°C for archaea, 58°C for fungi and 54°C for protozoa) for 30 s, and elongation at 72°C for 45 s, and following with post-elongation at 72°C for 10 min. The PCR amplification products were purified using Wizard^® SV Gel and PCR Clean-Up System (Promega, Madison, WI, United States). Barcoded amplicons were mixed at equal molar rations, and submitted for 300 bp pair-end sequencing on an Illumina MiSeq PE300 instrument (Illumina, San Diego, CA, United States).

Raw data were filtered through Quantitative Insight into Microbial Ecology (QIIME) (Caporaso et al., 2010b) and Mothur (Schloss et al., 2009). Reads were deconvoluted into individual samples based on their barcodes. The pair-end reads were assembled based on their overlapped sequences using COPE (Connecting Overlapped Pair-End) (Liu et al., 2012) and trimmed of primers.

For bacteria and archaea, the assembled sequences were assigned to operational taxonomic units (OTUs) at a 97% identity threshold using UPARSE (Edgar, 2013). Chimeric sequences were removed using UCHIME. Taxonomy classifications were assigned against the latest Greengenes database (May 2013 release) for bacteria (DeSantis et al., 2006) and RIM-DB version13_11_13 for archaea (Seedorf et al., 2014), respectively. The taxonomic assignment was performed using Mothur based implementation of the RDP Bayesian classifier with a 0.80 confidence threshold. For fungi and protozoa, sequences were assigned into OTUs using UCLUST at a 97% identity threshold. The OTUs were classified using BLAST against SILVA database for fungi and the reference database constructed by Kittelmann and Janssen (2011) for protozoa. Sequences with identity scores, between 95 and 97% were resolved at the genus level, 90 and 95% at the family level, 85 and 90% at the order level, 80 and 85% at the class level, and 77 and 80% at the phyla level (Hristov et al., 2012). Sequences were aligned by PyNAST (Caporaso et al., 2010a), and phylogenetic trees were constructed using FastTree.

Alpha diversity was examined using the alpha rarefaction pipeline of QIIME. Principle coordinate analysis (PCoA) was performed using unweighted UniFrac distance (jackknifed beta diversity from resampling 100 times at a depth of 75% of least number of sequences). Transformation-based principal component analysis (tb-PCA) was used to ordinate the samples based on OTU abundance data using the Hellinger distance (Legendre and Gallagher, 2001), and only the top 10 OTUs were presented.

Total RNA was extracted from mucosa samples (200 mg) using TRIzol reagent (Invitrogen, Carlsbad, CA, United States). The RNA integrity was evaluated using 1% agarose gel electrophoresis, and its quantity was measured using ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States). The RNA was reversely transcribed using the First Strand cDNA Synthesis Kit (Toyobo, Osaka, Japan), and the cDNA samples were stored at -20°C until use. The RT-PCR analysis was used to evaluate expression of genes encoding TLRs, inflammatory cytokines, toll-IL-1 receptor (TIR)-domain containing adaptors, and tight junction proteins (TJs), with GADPH and β-actin as the housekeeping genes. The primers are listed in Supplementary Table S2, and procedures are detailed in Supplementary Text 1.

The effect of feeding type was examined from day 28 to 70. Data were analyzed using the MIXED procedures of SAS 8.2 (SAS Institute Inc., Cary, North Carolina, United States). Feeding type, age, and feeding type × age interaction were regarded as fixed effects, animal nested within feeding type × age was regarded as the random effect, and each individual animal was the experimental unit. The following model was used:

where Y_ijk = dependent variable, μ = overall mean, A_i = fixed effect of age, F_j = fixed effect of feeding type, Ani_k = random effect of animal nested within feeding type × age, and e_ijk = the random residual error.

The effect of age on ruminal microbial diversity and expression of genes involved in host mucosal innate immune function was evaluated from day 0 to 70. The MIXED procedure of SAS and Tukey’s test to compare least squares means were used, with age as the fixed effect, animal within age as random effect, and individual animal as the experimental unit. The following model was used:

where Y_ij = dependent variable, μ = overall mean, A_i = fixed effect of age, Ani_j = random effect of animal within age, and e_ij = the random residual error. Orthogonal contrasts were performed to test linear and quadratic effects of age. If feeding type × age interaction was not significant within day 28 to 70, effect of age was assessed using combined data of both groups; otherwise, they were presented separately. All data were presented as least squares means, and significance was declared with P < 0.05.

The PROC CORR procedure of SAS was used to examine the Spearman’s rank correlations between ruminal microbial abundance and expression of mucosal innate immune function related genes. Only those correlation coefficients above 0.50, and P < 0.05 were defined as correlated with each other and plotted using corrplot package in the R studio.

Sequencing data in this study were deposited in the NCBI Sequence Read Archive (SRA) under accession numbers PRJNA362795.

As reflected by rarefaction curve, sequencing depth was adequate to cover the community of each microbial domain, with Good’s coverage ranged from 0.980 to 0.999 (Supplementary Figure S1 and Table S3). Alpha diversity of ruminal bacteria, archaea and fungi was more abundant (P < 0.05) in G kids compared to that in S kids. Alpha diversity indices (OTU number, Chao, Ace and Shannon) of ruminal bacteria and archaea in both groups increased with age (P < 0.01). For G kids, alpha diversity indices of ruminal fungi increased (P < 0.01), while those indices of ruminal protozoa decreased quadratically (P < 0.05) with age. PCoA analysis revealed that microbial communities could be clearly discriminated by both age and feeding type (Figure 1). The tb-PCA results revealed the top 10 OTUs within each microbial domain contributing to separation at different ages and for feeding types (Figure 2). Details of these are presented in Supplementary Text 2.

A total of 20 phyla and 48 genera were identified in all samples, with 23 bacterial genera each accounting for >0.5% of all sequences. A great proportion of sequences (31.2–69.2%) remained unclassified at the genus level (Figure 3A). Relative abundances of Prevotella, Asteroleplasma, Coprococcus were greater (P < 0.05), while relative abundances of Butyrivibrio and Pseudobutyrivibrio were lower (P = 0.011) in the rumen of G kids compared to those in S kids. In both groups, Arcobacter (30.1%), Escherichia (13.3%), Porphyromonas (5.7%) and Fusobacterium (2.1%) were predominant ruminal genera at day 0, whist their relative abundances declined to minimal values (<0.2%) afterward. Conversely, Prevotella relative abundance increased drastically (P < 0.01) with age, being the most dominant genus after day 28 (>20%), while Bacteroides relative abundance was high (12.1%) at only day 7. Moreover, relative abundances of ruminal BF311, CF231, Butyrivibrio, Pseudobutyrivibrio, and Fibrobacter in both groups increased quadratically (P < 0.05) with age. After day 28, relative abundances of ruminal Asteroleplasma and Anaeroplasma ranged from 0.3 to 1.2% in G kids, while their ruminal relative abundances were low (<0.2%) in S kids. Relative abundances of ruminal Akkermansia, Corprococcus, Ruminococcus, Selenomonas, Succiniclasticum and Treponema ranged from 0.1 to 3.6% in both groups. And relative abundances of Enterococcus and Lactobacillus were low (<1.0%) in the rumen of both groups.

Archaea were not detectable at day 0, and began to be detected at day 7. The relative abundance of Mbb. gottschalkii within the Methanobrevibacter genus was greater (P = 0.022), while relative abundances of genus Methanosphaera, as well as Group 11 and 12 within the Methanomassilliicoccales order were lower (P < 0.05) in the rumen of S kids compared to those in G kids (Figure 3B). Within the Methanobrevibacter genus, relative abundances of Mbb. gottschalkii (16.1–66.7%) and Mbb. ruminantium (2.9–9.7%) clades increased drastically from day 28, constituting the top two abundant clades. The highest relative abundance of Mbb.acididurans in the rumen of S and G kids was 15.9% (day 70) and 3.7% (day 70), respectively. Relative abundances of Methanimicrococcus and Methanosphaera genera were low (<2.0%) in the rumen of both kids, while Methanomicrobium genus was only detected in the rumen of S kids at day 28 (6.4%). Furthermore, a significant proportion of Methanomassilliicoccales (from 12.1 to 34.8%) was found in the rumen of both groups.

Two fungal phyla, Ascomycota and Neocallimastigomycota, and 103 fungal genera were identified, with 16 fungal genera each accounting for >0.5% of all sequences (Figure 3C). The relative abundance of Neocallimastix was greater (P < 0.001), while relative abundances of Microcyclosporella, Knufia, Parastagonospora and Uwebraunia were lower (P < 0.05) in the rumen of S kids compared to those in G kids. Ruminal Boeremia (25.6%), Neurospora (11.3%), Arthrinium (7.6%), Aspergillus (4.5%), Purpureocillium (4.3%) were the top five abundant genera at day 0, while their relative abundances declined significantly (P < 0.01) afterward. The highest relative abundance of Candida in the rumen was 56.4% (day 7). The relative abundances of Knufia, Microcyclospora, Microcyclosporella, Parastagonospora, Penidiella, Phoma and Uwebaunia in the rumen of S kids were low (<2%) at all ages, while their relative abundances in the rumen of G kids increased with age (P < 0.05). Relative abundances of ruminal Neocallimastix and Piromyces in S kids increased quadratically with age (P < 0.05), peaking at 94.2% (day 28) and 61.3% (day 42), respectively, while their ruminal relative abundances in G kids increased linearly with age (P < 0.05), peaking at 7.1% (day 70) and 14.3% (day 70), respectively.

Protozoa began to colonize after day 28, with >99% sequences classified into six genera (Figure 3D). The relative abundances of Anoplodinium–Diplodinium and Polyplastron were greater (P < 0.05), while Polyplastron relative abundance was lower (P < 0.001) in the rumen of S kids compared to those in G kids. Entodinium was the most abundant genera (39.9–96.8%) in the rumen, with its relative abundance declined linearly (P = 0.001) in S kids, while increased linearly (P = 0.002) in G kids from day 28 to 70. Relative abundances of Anoplodinium–Diplodinium and Polyplastron were high (8.2–87.4%) in the rumen of S kids, but low (<3.0%) in the rumen of G kids. In contrast, the relative abundances of ruminal Epidinium and Eremoplstron-Diploplastron were low (<1.0%) in S kids, but high in G kids (7.3–32.4%).

The expression of genes encoding TLRs (1–10) was down-regulated in the ruminal mucosa of S kids when compared to that in G kids. For S kids, mRNA abundance of TLRs1, 2, 4, 7, and 10 increased slightly (P < 0.05) from day 0 to 28, and then decreased afterward (Figure 4A). The mRNA abundance of other TLRs increased during the first week (P < 0.05), and declined afterward. Conversely, for G kids, expression of TLRs increased from day 0 to 28 (P < 0.01), remaining relatively stable afterward.

When genes involved in cytokine production were compared, the expression of IL1α and IL1β was up-regulated (P < 0.01), while expression of IL10 and TNFα was down-regulated (P < 0.01) in the ruminal mucosa of G kids when compared to those in S kids (Figure 4B). The mRNA abundance of IL1α, IL1β and IL6 in the ruminal mucosa of S kids increased from day 0 to 28, declined afterward, while their expression in the ruminal mucosa of G kids increased (P < 0.05) from day 0 to d 70. In both groups, mRNA abundance of IL10, IL18 and TNFα in the ruminal mucosa declined with age (P < 0.05).

For genes encoding TIR domain-containing adaptors, MYD88 expression was age- and feeding type-independent (Figure 4C). The expression of TICAM2 (P = 0.003) was up-regulated in the ruminal mucosa of G kids when compared to that in S kids. Age exerted a quadratic effect (P = 0.006) on TICAM1 mRNA abundance in the ruminal mucosa of S kids. TICAM2 mRNA abundance in the ruminal mucosa of S kids increased during the first month and declined afterward, whist its value in G kids increased (P = 0.043) from day 0 to 70.

For genes encoding tight junction proteins (TJs), the expression of Occludin, Claudin1 and Claudin4 was elevated (P < 0.05) in the ruminal mucosa of S kids when compared to those in G kids (Figure 4D). Furthermore, expression of TJs declined (P < 0.05) from day 0 to 70 in the ruminal mucosa of both groups.

As presented in Figure 5, TLR1 mRNA abundance was positively correlated with (R > 0.50) relative abundances of Anaeroplasma, Knufia, Piromyces, Uwebraunia, and Entodinium, while negatively correlated with (R < -0.50) relative abundances of Mbb. gottschalkii and Anoplodinium–Diplodinium. TLR2 mRNA abundance was negatively correlated with (R = -0.52) Anoplodinium–Diplodinium relative abundance, TLR4 mRNA abundance was negatively correlated with (R = -0.50) Mbb. gottschalkii relative abundance, while TLR3 mRNA abundance was positively correlated with (R > 0.50) relative abundances of Knufia and Uwebraunia. TLR6 mRNA abundance was positively correlated with (R > 0.50) relative abundances of Anaeroplasma, Kufia, Neocallimastix, and Uwebraunia). TLR8 mRNA abundance was positively correlated with (R > 0.50) relative abundances of Kufia, Neocallimastix, and Uwebraunia, while negatively correlated with (R = -0.61) Mbb. gottschalkii relative abundance. TLR9 mRNA abundance was positively correlated with (R > 0.50) relative abundances of Fibrobacter, Anaeroplasma, Kufia, Neocallimastix, Uwebraunia, and Epidinium. TLR10 mRNA abundance was positively correlated with (R = 0.51) Neocallimastix relative abundance, while negatively correlated with (R < -0.50) relative abundances of Anoplodinium–Diplodinium and Polyplastron. The mRNA abundance of IL1α and IL1β was negatively correlated with (R < -0.50) relative abundances of Mbb. gottschalkii, Anoplodinium–Diplodinium and Polyplastron. IL6 mRNA abundance was positively correlated with (R = 0.60) Mmc.Group10 relative abundance, while IL18 mRNA abundance was negatively correlated with (R = -0.59) Epidinium relative abundance. The mRNA abundance of IL10 and TNFα was positively correlated with (R > 0.50) Bacteroides relative abundance, but negatively correlated with (R < -0.50) relative abundances of Bacteroides, Prevotella, Coprococcus, Fibrobacter, Anaeroplama, Methanosphaera, Phoma, Piromyces, and Uwebraunia. The mRNA abundance of Occludin, Claudin1 and Claudin4 was negatively correlated with (R < -0.50) relative abundances of most microbes. TICAM1 mRNA abundance was positively correlated with (R = 0.63) Mmc.Group10 relative abundance, while TICAM2 mRNA abundance was negatively correlated with (R < -0.50) relative abundances of Mbb. gottschalkii and Anoplodinium–Diplodinium.