The animal protocol in this study was approved by the Animal Care and Use Committee of Inner Mongolia University (Approval No. IMU-sheep-2020-041).

Healthy male East Friesian × Small-tail Han lambs were purchased from Inner Mongolia Lark Biotechnology Co., Ltd. Fourteen lambs with an average age of 60 days and with body weight of about 22.03 ± 1.08 kg were selected. The feeding trials were conducted at the Inner Mongolia University Meat Sheep Nutritional Base, where each lamb was housed individually in the same pens. Fourteen lambs were randomly allocated to two dietary treatment groups: one group was feed with concentrate and Leymus chinensis hay (Lc group); another was fed with concentrate and alfalfa (Medicago sativa) hay (Ms group) with free access to water. The Lc group consumed about 325 g of concentrate and 803 g of (Leymus chinensis hay each sheep per day, while the Ms group consumed about 325 g of concentrate and 925 g of alfalfa hay each sheep per day. Ingredient and chemical compositions of diets are shown in Supplementary Table S1. Before the experiment, lambs were acclimatized to the environment for about 1 week. The experiment lasted for 60 days. At the end of the experiment, lambs were selected for slaughter after fasting for about 12 h, and the LD and rumen contents samples were collected. Muscle samples were collected using a disposable scalpel at the same location in the LD of each lamb into 50 mL RNAase-free centrifuge tubes and stored at −80°C for metabolite analysis. After the rumen was opened, the rumen contents were filtered using 4 layers of sterile gauze, and the rumen contents were collected in a 50 mL RNAase-free centrifuge tube and stored at −80°C for 16S rDNA amplicon analysis.

Dry matter and crushed ash of Leymus chinensis hay and alfalfa hay samples were determined according to Ran et al. (23). Crude protein, acid detergent fiber, neutral detergent fiber, and crude fat in the Leymus chinensis hay and alfalfa hay samples were determined using an Automatic Kjeldahl Protein/Nitrogen Analyzer (K1160, Hanon Advanced Technology Group, China), Automatic Fiber Analyser (F2000, Hanon Advanced Technology Group, China) and Automatic Soxhlet Extractor (SOX606, Hanon Advanced Technology Group, China) respectively with reference to Shi et al. (24). The content of mineral elements was determined with using microwave (REVO, LabTech, China) and inductively coupled plasma-optical emission spectrometer (PQ 9000, analytikjena, German) (25).

Body weight and carcass weight were measured following overnight fasting and evisceration, respectively. Loin-eye area between the 12th and 13th ribs was traced on sulphate papers and calculated (26). Dressing percentage = carcass weight/body weight. Organ index = organ weight/body weight (27). Total weight = final body weight (FBW) − initial body weight (IBW). Average daily gain (ADG) = total weight gain/total days. Feed to gain ratio (F/G) = average daily feed intake (ADFI)/ADG.

The rumen contents were collected and stored at −80°C. Microbial DNA was extracted using a TIANGEN kit. The V3–V4 region of the bacterial 16S rDNA gene was amplified by PCR (95°C for 3 min, followed by 30 cycles at 95°C for 30 s, 50°C for 30 s, 72°C for 45 s, and an extension at 72°C for 10 min) using the F3 (ACTCCTACGGGAGGCAGCAG) and R4 (GGACTACHVGGGTWTCTAAT) primer pair (28). DNA Library Prep Kit for Illumina following manufacturer’s recommendations and index codes were added. High throughput sequencing was performed utilizing the Illumina MiSeq PE300 platform to detect the 16S rDNA amplicons according to standard protocols.

After sequencing, paired-end reads from sequencing were merged by FLASH, and low-quality reads were filtered by Trimmomatic. UPARSE was used to align operational taxonomic units (OTUs) at 97% identity. Taxonomy was assigned to OTUs by searching against the Greengenes database version 13.8. α- and β- diversities were calculated using QIIME2. The biomarkers with statistical differences were screened by Linear discriminant Effect Size (LEfSe) analysis, and the screening criteria were LDA threshold ≥4.0 and p-value <0.05. The redundancy analysis (RDA) and Spearman’s rank correlation analysis were performed using the R packages: vegan and pheatmap, respectively.

Untargeted metabolomics analysis was referenced to the previous method (29, 30). Briefly, 200 mg LD samples were mixed with 2-chlorophenylalanine (4 ppm) methanol (−20°C) and ground by a high-throughput tissue grinder for 90 s at 60 Hz. The samples were centrifuged at 4°C for 10 min at 12,000 rpm, and the supernatant was filtered through 0.22 μm membrane to obtain the prepared samples for LC-MS. Twenty microlitre from each sample were taken to the quality control samples, and the rest of the samples were used for LC-MS detection.

LC-MS detection was referenced against the previous method (31). Chromatographic separation was accomplished in a Thermo Ultimate 3000 system equipped with an ACQUITY UPLC^® HSS T3 column maintained at 40°C. The temperature of the autosampler was 8°C. Gradient elution of analytes was carried out with 0.1% formic acid in water (C) and 0.1% formic acid in acetonitrile (D) or 5 mM ammonium formate in water (A) and acetonitrile (B) at a flow rate of 0.25 mL/min. Injection of 2 μL of each sample was performed after equilibration. An increasing linear gradient of solvent B (v/v) was used as follows: 0–1 min, 2% B/D; 1–9 min, 2%–50% B/D; 9–12 min, 50%–98% B/D; 12–13.5 min, 98% B/D; 13.5–14 min, 98% B/D; 14–20 min, 2% D-positive model (14–17 min, 2% B-negative model) (31).

The ESI-MSn experiments were executed on the Thermo Q Exactive Focus mass spectrometer with the spray voltage of 3.8 kV and −2.5 kV in positive and negative modes, respectively. Sheath gas and auxiliary gas were set at 30 and 10 arbitrary units, respectively. The capillary column temperature was 325°C. The analyzer scanned over a mass range of m/z 81–1,000 for a full scan at a mass resolution of 70,000. Data-dependent acquisition (DDA) MS/MS experiments were performed with an HCD scan. The normalized collision energy was 30 eV. Raw data were converted to the mzXML format by Proteowizard (v3.0.8789), and then processed by the MetaboAnalystR package (32).

GraphPad Prism was used for statistical analysis. The results of growth performance, carcass performance, and organ index are presented as means ± standard error (SE). All data were considered statistically significant at p < 0.05. Statistical significance of organ index, growth and carcass performance, and metabolites were determined by t-test. The statistical significance of α-diversity was performed with Wilcoxon tests.

To investigate the effects of Leymus chinensis hay and alfalfa hay on growth and carcass performance, we first compared IBW, FBW, ADG, carcass weight, and loin-eye area between the two groups. IBW was not significantly distinct, while total weight gain, FBW, and ADG were significantly greater in the Ms. group. Carcass weight, dressing percentage, and loin-eye area were also obviously increased in the Ms. group. However, lambs in the Ms. group had higher ADFI. The F/G was similar between the two groups (Table 1). In addition, we analyzed the organ index in lambs. Kidney index was increased in the Ms. group, while others including spleen, lung, liver, and heart index were not significantly affected (Table 2). These data indicated that compared to Leymus chinensis hay, alfalfa hay improved growth and carcass traits in lambs.

To investigate the effects of Lemus chinensis hay and alfalfa hay on rumen microbiota, we performed and analyzed bacterial 16S rDNA sequencing of rumen contents between two groups. The observed OTU and Shannon curves reached the saturation phase, suggesting the sufficient and reliable sequence depth captured in this study (Figure 1A). A total of 2,140 OTUs were identified in both groups, among which 3,491 and 3,861 specific OTUs were observed in the Lc and Ms. groups, respectively (Figure 1B). Shannon and Simpson indexes reflected the bacterial community richness and diversity. No significant differences were observed in the Shannon and Simpson index between the two groups (Figure 1C). The β-diversity presenting as Principal Co-ordinates Analysis (PCoA), Non-metric multidimensional scaling (NMDS), Principal Components Analysis (PCA), and Partial Least Squares Discrimination Analysis (PLS-DA) was further applied to analyze the variation of bacterial structure. The NMDS, PCA, and PLS-DA are in Supplementary Figure S3. The Lc group rumen contents (LcRc) formed a different rumen bacterial community clustered separately from the Ms group rumen contents (MsRc) (Figure 1D).

To further define the composition of rumen microbiota in the LcRc and MsRc, the percentages of bacterial community abundance were analyzed at phylum and genus levels. At the phylum level, Bacteroidetes, Firmicutes, Spirochaetes, and Fibrobacteres were the dominant phyla (Figure 2A). The relative abundance of Fibrobacteres, Bacteroidetes, and Spirochaetes were obviously greater in the Lc group (Supplementary Table S2); the relative abundance of Firmicutes, Actinobacteria, Fusobacteria, Verrucomicrobia, and Proteobacteria were obviously greater in the MsRc (Supplementary Table S2). At the genus level, the top 20 genera were present with relatively high abundance (Figure 2B). The 17 genera significantly altered between two groups, among which only 5 genera, i.e., Clostridiu, Fibrobacter, Selenomonas, BF311, and Treponema, were greater in the LcRc. Most of the altered bacteria in the MsRc, including Moryella, Blautia, p_75_a5, Faecalibacterium, Mogibacteriaceae, Anaerostipes, Bifidobacterium, Eggerthella, and Collinsella belong to Firmicutes and Actinobacteria (Supplementary Table S3).

The rumen bacterial composition was further compared through LEfSe analysis with 4.0 as the threshold on the LDA score to identify specific species in each group. The constitution of rumen bacteria changed between LcRc and MsRc. A total of 17 taxa sequences were enriched in the LcRc and mainly belonged to Fibrobacteres, Bacteroidetes, and Spirochaetes (Figure 3). A total of 7 rumen bacteria mainly belonging to Firmicute were determined as enriched in the MsRc (Figure 3). Taken together, lambs in Lc and Ms. groups had distinct compositions of rumen bacterial community, which may lead to different growth performance, carcass performance, and metabolic profiles.

To further investigate the potential relationship between rumen microbiota variation and ADG, carcass weight, and body weight, we performed the RDA and Spearman’s correlation between differential rumen bacteria and host parameters. According to RDA, the Lemus chinensis hay and alfalfa hay treatments were distinct and two different clusters were observed. ADG, carcass weight, and body weight showed a significant correlation with bacterial community (Figures 4A,B). The Spearman’s correlation data showed that the relative abundance of BF311 and Pseudobutyrivibrio were negatively related to ADG, carcass weight, and body weight (Figure 4C). Additionally, most genus-level bacteria such as Enterococcus, Lactobacillus, Acidaminococcus, and Moryella belonging to Firmicutes showed a strongly positive correlation with ADG, carcass weight, and body weight (Figures 4C). Alistipes belonging to Bacteroidetes, Eggerthella and Bifidobacterium belonging to Bacteroidetes, and Anaeroplasma belonging to Tenericutes, were also positively associated with ADG, carcass weight, and body weight (Figure 4C).

To define similarities and differences between lamb meat fed with Leymus chinensis hay and alfalfa hay, we performed and analyzed the untargeted metabolome profiles of LD samples. PCA score plots with positive and negative ionization modes revealed no significant clustering of the two treatment groups (Figure 5A), but further analysis revealed significant differences in metabolite composition. These metabolites mainly consisted of peptides, lipids, organic acids, steroids, carbohydrates, nucleic acids, vitamins and cofactors, hormones, and transmitters (Figure 5B). These two comparisons identified 96 significant differential metabolites (67 in the positive mode and 29 in the negative mode), among which 52 and 44 differential metabolites were elevated in the Lc and Ms. groups, respectively. In detail, organic acid and derivatives including N-Acetyl-L-leucine, Pantetheine, Cysteinyl glycine, Alanyl tyrosine, and Ala-Ile were greater in the Ms group longissimus dorsi (MsLD), while Cinnamoyl glycine, L-arginosuccinate, and Glu-Thr were greater in the Lc group longissimus dorsi (LcLD). Lipids and lipid-like molecules such as 2-methylbutyroylcarnitine, Glycoursodeoxycholic acid, and 2-Hydroxymyristic acid were higher in the LcLD, while Hexadecanediouc acid, prostaglandin F2alpha, and tetradecanedioic acid were higher in the MsLD. Interestingly, LcLD contained higher organoheterocyclic compounds including Momocrotaline, (+)-alpha-lipoic acid, Methyl indole-3-acetate, benzenoids (i.e., hippuric acid), phenylpropanoids, and polyketides (i.e., 2-phenylpropionic acid). MsLD had a higher relative abundance of alkaloids and derivates including ecgonine methyl ester (Figures 5C,D). Between Leymus chinensis hay and alfalfa hay treatment, peptides including Ala-Ile, Glu-Thr, N-Acetyl-Asp-Glu, gamma-glutamyl-cysteine, and gamma-glutamyl-L-leucine, significantly changed (Figures 5C,D). Furthermore, based on the pathway analysis, we found that these differential metabolites were mainly enriched in alanine, aspartate and glutamate metabolism, D-glutamine and D-glutamate metabolism, phenylalanine metabolism, nitrogen metabolism, and tyrosine metabolism (Figures 6A,B). Changed metabolites in these pathways were shown in Table 3, indicating that different forages affected amino acid composition and contents, and its metabolism, leading to various meat nutritional qualities.

To facilitate the identification of potential biomarkers between the two groups, PLS-DA and OPLS-DA were applied to identify the major differences in metabolites. PLS-DA score plots showed the best separation of LcLD and MsLD samples (Supplementary Figure S2A). Several potential biomarkers with VIP >2 including Monocrotaline, Hippuric acid, 6-Hydroxytestosterone, N-Acetyl-L-leucine, LysoPE 18:0, Ala-Ile, Methadone-d9, Cysteinylglycine, Cinnamoylglycine, Pilocarpine, and 2-Phenylpropionic acid and were highlighted in the volcano plots (Supplementary Figure S2B). To minimize overfitting in the PLS-DA model, we further explored potential markers using OPLS-DA between groups. The analysis presented a cross-validated score plot (R2X = 0.311, R2Y = 0.993, Q2 = 0.737, p-value = 0.01) of the discriminating model between LcLD and MsLD in positive ionization mode. The cross-validated score plot (R2X = 0.454, R2Y = 0.995, Q2 = 0.526, p-value = 0.05) was shown in negative ionization mode. Identified biomarkers were the same as those identified by PLS-DA as having a trend toward being distinct between the two groups (Supplementary Figure S3). Additionally, random forest showed top differential metabolites including N-Acetyl-L-leucine, Hippuric acid, (+)-alpha-lipoic acid, and 2-Phenylpropionic acid (Supplementary Figure S4).

Differential metabolites and potential biomarkers have been clarified as mentioned above. The correlation of potential biomarkers with growth performance was further analyzed. A total of 30 differential metabolites were significantly associated with ADG, carcass weight, and body weight. Among them, 14 metabolites including N-Acetyl-L-leucine, Methadone-d9, Hexadecanedioic acid, and alanyltyrosine were positively related to ADG, carcass weight, and body weight. Sixteen metabolites such as LysoPE 18:0, Cinnamoylglycine, and 2-Methylbutyroylcarnitine exhibited a negative relationship with growth performance (Figure 7).

The 16S rDNA sequencing revealed that lambs with the consumption of dietary Leymus chinensis hay and alfalfa hay exhibited distinct rumen microbiota communities. The metabolomics analysis indicated that exposure to Leymus chinensis hay and alfalfa hay resulted in different metabolic profiles of LD. Therefore, we explored the correlations between the rumen microbiota and the altered LD metabolites through Spearman’s correlation analysis, revealing high correlations between key rumen microbiota with differential metabolites (Figure 8). In detail, Fibrobacter, Selenomonas, BF311, Treponema, and Clostridium exhibited a strong positive relationship with organic acid and derivatives such as Glu-Thr, benzenoids including hippuric acid, phenylpropanoids and polyketides including 2-phenylpropionic acid, and organoheterocyclic compounds including Methyl indole-3-acetate, Momocrotaline, and (+)-alpha-lipoic acid. And they were negatively associated with other differential metabolites including Ala-Ile, N-Acetyl-L-leucine, Tetradecanedioic acid, and Methadone-d9. Other genera including Bifidobacterium, Faecalibacterium, Blautia, Roseburia, Moryella, Streptococcus, [Ruminococcus], and Odoribacter, showed an obviously opposite correlation trend with the five genera mentioned above (Figure 8). These results indicated that the rumen bacterial compositions obviously and specifically impacted host LD metabolites.