In the present study, all animal procedures were approved by the Institutional Animal Care and Use Committee of Northwest A&F University.

Field experiment were performed at a Saanen goat farm in Baoji, Shaanxi (34°41′N, 109°09′E). Thirty ewes delivering in 1 week (third birth), including 10 ewes producing single, 10 ewes producing twins, and 10 ewes producing triplets, with each litter including at least one female newborn, were selected. After born, one female kid was chosen from each litter, 30 female kids in total, were randomly selected from these 30 litters for sampling. Briefly, 10 female kids were randomly selected from single, twins, and triplets, respectively, and named as LZ1, LZ2, and LZ3. And the birth weight of those kids was weighed within 12 h after birth. And then, all kids were transferred to the kid barn with bottle-feeding of mixed milk collected from ewes. During pre-weaning phase (~90 days old), all goats were feed milk, alfalfa hay, and concentration mixture. After weaning, the goats were fed TMR with a 60:40 forage to concentration ratio. The animals were fed three times daily at 0730, 1300, and 1900 h. Water was available ad libitum. The detailed feeding programs and ingredient compositions are shown in Supplementary Table S1.

The body weight, wither height, body length, and heart girth of weaning and 6-month-old (188.9 ± 0.4 days old, mean ± SE) of those three group goats were measured. And the rumen fluid sample of those young goats were collected via esophageal tube before morning feeding. Briefly, the first ~50 ml of rumen fluid was discarded to avoid saliva contamination, and the text 30 ml rumen fluid strained through four layers of sterile cheesecloth under a constant flux of CO₂. Then the rumen fluid was aliquoted and stored at −80°C for further analysis.

The feed intake of all goats was measured 1 week before sampling (171.9 ± 0.4 days old). In brief, feed offered to and refused by each goat was recorded continuously for 7 days. The feed samples were dried at 65°C for 48 h to obtain dry matter content of ration. Daily dry matter intake (DMI) per goat was calculated by multiplying daily as-fed intake by dry matter content of ration.

Ruminal VFA (acetate, propionate, butyrate, valerate, isobutyrate, and isovalerate) analysis were performed through separation and quantification by a GC (Agilent 7820A, Santa Clara, CA, United States) with a capillary column (AE-FFAP of 30 m × 0.25 mm × 0.33 μm; ATECH Technologies Co., Lanzhou, China; Li et al., 2014). Briefly, thawed rumen fluid samples were centrifuged for 10 min at 13,500 rpm at 4°C. To remove the protein, a 1.5 ml supernatant was mixed with 300 μl 25%w/v metaphosphoric solution. After standing for 4 h at 4°C, the mixture was centrifuged for 10 min at 13,500 rpm at 4°C. Then 1 ml supernatant was moved into 200 μl 25% crotonic acid, and the mixture was collected into an EP tube passed through a 0.45-μm filter.

The detection of ruminal free amino acids abundance were conducted using LC–MS/MS (Exion LC AC, QTRAP 5500, AB SCIEX, Framingham, MA, United States) according to method described as previously described (Chen et al., 2022). In brief, in order to precipitate proteins rumen fluid samples were mixed through with sulfosalicylic acid (10%, wt/vol) and centrifuged at 3,000 × g for 10 min at 4°C. Then the supernatant fluid was collected and filtered through a 0.45-μm filter.

Ruminal fatty acids composition was analyzed using the method of Sun and Gibbs (2012). The freeze-dried sample (0.3–0.4 g) was methylated with 4 ml of 0.5 mol/L NaOH/methanol for 15 min at 50°C. Then the mixture was added 4 ml 2% HCl/methanol (v/v) and vortexed for 5 min, followed by a water bath (50°C) for 60 min. The extract was dissolved in 2 ml of heptane and then introduced to the GC (Agilent Technologies 7820A GC system, Santa Clara, CA, United States) equipped with a fused silica capillary column (SP-2560, 100 m × 0.25 mm × 0.2 μm; Supelco Inc., Bellefonte, PA, United States).

The total DNA was extracted using the QIAamp DNA Stool Mini kit (QIAGEN, Germany) according to the manufacturer’s protocol. The DNA concentration was determined by using a Nanodrop-2000 (Thermo Fisher Scientific, United States), and the purity was monitored on 1% agarose gel electrophoresis. The DNA was stored at −80°C until further processing. The amplicon library preparation was performed by polymerase chain reaction amplification of the V3-V4 region of the 16S rRNA gene using the primer pairs: 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). The PCR program was carried as previously described (Xue et al., 2018). Paired-end (2 × 300 bp) sequencing was performed on Illumina platform (Illumina, United States) according to the standard protocols (Caporaso et al., 2012).

The raw sequences were merged with FLASH (v1.2.11; Magoč and Salzberg, 2011) and quality filtered with fastp (0.19.6; Chen et al., 2018). Sequences were demultiplexing using QIIME2 v2021.8, and the construction of an amplicon sequence variant (ASV) table using DADA2 (Callahan et al., 2016). Bacterial 16S ASVs were assigned a taxonomy using the SILVA database v138 as the reference.

The alpha diversity diversities of the bacterial communities were determined using various diversity indices (Sobs, ACE, Chao 1, Shannon, Simpson) and calculated using the procedures within QIIME 2. The NMDS (non-metric multidimensional scaling) graphs were performed based on Bray-Curtis distance and statistical significance was determined using analysis of similarities (ANOSIM) with 999 permutations at the ASV level.

Differential DMI, growth performance traits (body weight, ADG, wither height, body length, and heart girth), and ruminal metabolites (VFA, AA, and FA) were tested by one-way ANOVA among three groups, and Student’s t-test was used between two groups. Orthogonal polynomials were used to determine the linear effects of increasing the litter size.

The taxa with different relative abundances among those three groups were identified by a Kruskal–Wallis test, and Wilcoxon rank-sum test was used between two groups with a false discovery rate (FDR) value <0.05 to correct the p values. The differences were statistically significant at p < 0.05 or a tendency of difference at p < 0.10. Spearman’s rank correlation coefficients were used to examine the correlations between bacterial abundance and rumen metabolites concentrations.

In the present study, although there was no significant difference in birth weight between these three groups (p > 0.05). Numerically, the LZ3 group (2.94 ± 0.20 kg) with lower birth weight compared to LZ1 (3.21 ± 0.14 kg) and LZ2 (3.20 ± 0.12 kg) group (Supplementary Figure S1A).

The growth performance traits were shown in Figure 1. With the increased litter size, the weaning weight, ADG, and wither height of goat kids were linearly decreased (p < 0.05). Moreover, compared to the LZ1 group, the weaning weight, ADG, and wither height were significantly lower in LZ3 group (t-test, p < 0.05).

And, there was no significant difference in DMI between the three groups (p = 0.729, Supplementary Figure S1B) at the youth period. With the increased litter size, the youth body weight, and ADG of young goats were linearly decreased (p < 0.05), the body length and heart girth had a tendency to decrease (p < 0.1). Furthermore, compare to the LZ1 group, LZ3 group had a lower body weight, ADG, heart girth (t-test, p < 0.05), and body length (t-test, p = 0.093) at the youth period.

As shown in Table 1, with the increased litter size, the concentration of ruminal propionate, total VFA of young goats were linearly decreased (p < 0.05), the ruminal butyrate concertation had a tendency to decrease (p = 0.098). In addition, compare to the LZ1 group, LZ3 group had a lower propionate, total VFA (t-test, p < 0.05), and butyrate concentration (t-test, p = 0.093).

Then the difference in ruminal AA and FA were identified between the three groups. We found that the concentration of some essential amino acids (EAA), such as His, Trp, Lys, Phe, Val and total EAA concentration, were linearly decreased with the increased litter size (Table 2, p < 0.05). In addition, compared to the LZ1 group, the concentration of His, Cys, Val, Trp, Phe, Lys, and total EAA were significantly decreased in LZ3 group (t-student, p < 0.05), and the concentration of some AA, such as Leu, Pro, and total branched chain amino acids (BCAA) of LZ3 group, tended to be decreased than LZ1 group (t-student, p < 0.1).

As shown in Table 3, the proportion of iso C14:0, anteiso C15:0, C15:0, saturated fatty acids (SFA), odd and branched chain fatty acids (OBCFA) and the ratio of SFA to UFA were linearly increased with the increased litter size (p < 0.05). And the proportion of C16:0, C18:1, C18:1t11, C18:1c9, unsaturated fatty acids (UFA) and monounsaturated fatty acids (MUFA) were linearly decreased with the increased litter size (p < 0.05). In addition, we found that the proportion of iso C14:0, anteiso C15:0, C15:0, OBCFA were significantly increased and C16:0, C18:1t11, C18:1c9, C18:1, and MUFA were significantly decreased in the LZ3 group compared to the LZ1 group.

Together, our results revealed that with the increased litter size, the growth performance and rumen fermentation efficiency were linearly decreased.

Based on the above results, the difference in rumen microbiota between LZ3 group and LZ1, LZ2 group were identified, to evaluate whether alterations of ruminal fermentative capacity were influence by rumen microbiota.

After size filtering, quality control, and chimera removal using the QIIME 2 pipeline, an average of 45,662 reads were generated for bacteria from a total 1,369,847 quality reads from 30 rumen fluid samples, with an average length of 416 bp. The calculated Good’s coverage exceeded 99.9% for all samples.

To compare bacterial community diversity across different groups, alpha-diversity and beta-diversity were evaluated. There were no significant differences in Sobs index, Ace index, Chao 1 index, Shannon index, and Simpson index among the three groups (Table 4, p > 0.05). And there were no significant difference in alpha-diversity indices between LZ1 and LZ3 group (p > 0.05). In a beta-diversity analysis (NMDS based on Bray-Curtis), the LZ1 and LZ2 group were clustered together and could not be distinguished (p = 0.477), whereas the LZ3 group was clearly distinguished from the LZ1 and LZ2 group (Figure 2, p < 0.05).

In total, 21 phyla, 32 classes, 71 orders, 120 families, and 237 genera were identified. More specifically, Firmicutes (50.2%), Bacteraioidota (43.1%), and Patescibacteria (3.4%) were the predominant phyla. The abundance of Firmicutes in LZ3 group was significantly lower than LZ2 group, while the abundance of Bacteraioidota in LZ3 group was significantly higher compared to LZ1 and LZ2 group (p < 0.05). Meanwhile, the Firmicutes: Bacteraioidota ratio was significantly decreased in the LZ3 group (Supplementary Figure S2, p < 0.05).

As shown in Figure 3, Prevotellaceae (24.8%) was the most abundant family in LZ3 group. While, F082 was the most abundant family in LZ1 (13.7%) and LZ2 (11.5%) group. Prevotellaceae abundance in LZ3 group was significantly higher than that in LZ1 and LZ2 group (p < 0.05). Conversely, the relative abundance of Rikenellaceae in LZ3 group was significantly lower than that in LZ1 group (p < 0.05). At the genus level, norank f F082 (11.6%), Prevotella (10.9%), Ruminococcus (9.6%), Rikenellaceae RC9 gut group (8.5%), Christensenellaceae R-7 group (6.8%), Oscillospiraceae NK4A214 group (6.6%), norank f Muribaculaceae (4.6%), Candidatus Saccharimonas (3.4%), unclassified c Clostridia (2.6%), and unclassified f Selenomonadaceae (2.5%) were the dominant (Figure 4). The abundance of Prevotella, unclassified f Prevotella and Prevotellaceae UCG-001 in LZ3 group were significantly higher than that in LZ1 and LZ2 group (p < 0.05). Succiiniclasticum abundance in the LZ3 group were significantly higher than that in LZ1 group (p < 0.05). Conversely, the abundance of Rikenellaceae RC9 group and unclassified c Clostridia were significantly lower than that in LZ1 or LZ2 group (p < 0.05).

As shown in Figure 5A, those LZ3 enriched bacteria were negatively correlated with ruminal VFA concentration and growth performance traits of young goats. For example, the relative abundance of Prevotellaceae UCG-001 was negatively correlated with ADG, wither height, and ruminal butyrate concentration (p < 0.05). The relative abundance of Prevotella, and Succinniclasticum were negatively correlated with ruminal propionate, butyrate, valerate, and total VFA concentration (p < 0.05). And Rikenellaceae RC9 gut group had a significantly positively correlated with ruminal isobutyrate and butyrate concentration (p < 0.05).

And then, we analyzed the correlation between ruminal AA, FA, and litter size related bacteria genera (Figures 5B,C). The results showed that those LZ3 enriched bacteria were significantly negatively correlated with ruminal AA, UFA, and MUFA, such as Prevotella and ruminal BCAA, total EAA, anteiso C15:0, and OBCFA were negatively correlated but positively correlated with C18:1t11, UFA, and MUFA (p < 0.05). unclassified f Prevotellaceae were significantly negatively associated with almost all ruminal AA, such as Leu, His and total BCAA, and C18:1t11 but positively associated with OBCFA (p < 0.05).