This study was implemented in line with the recommendations for animal care and ethics in China. All animal experiments were allowed by the Animal Ethics and Welfare Committee of Baotou Teacher’s College (Approval Number: AEWC-BTTC2023001).

The MSTN gene edited Luxi Yellow cattle used in this study were generated by CRISPR/Cas9 gene editing tools and somatic cell nuclear transfer as reported previously (8). Three 24-month-old wild-type female cattle and three female MSTN mutant cattle were chosen for sampling. The cattle were fed as common beef cattle in a local livestock farm.

The cattle were fasted for 24 h before slaughter. The slaughter process complied with the national standard operating procedures (GB/T 19477-2018, Cattle Slaughtering, China). During the slaughter process, jejunal tissues and their contents were collected from both MSTN gene-edited and wild-type cattle, and samples were quickly placed in pre-cooled tubes, then immediately stored in liquid nitrogen and continuously preserved at −80°C for successive analysis.

Samples for metabolite determination were extracted from the jejunal contents of three wild-type and three MSTN gene mutant cattle. A high-resolution mass spectrometer (Waters Xevo G2-XS QTOF) controlled by MassLynx V4.2 was used to collect primary and secondary mass spectral data in MSe mode (Waters). A series of processing operations including peak extraction, peak alignment were carried out on raw data collected by MassLynx V4.2 using Progenesis QI software. Metabolites (theoretical fragment identification deviation and mass deviation ≥100 ppm) were charactered according to the information provide by METLIN online database and a self-constructed Biomark database in Progenesis QI software.

After normalizing the raw peak area information to total peak area, further analysis was carried out as follows. PCA and Spearman correlation analysis were conducted to assess the reproducibility of intra-group samples and quality control samples. Classification and pathway information for the metabolites were analyzed by KEGG database. Fold change of each group was calculated and then compared based on the grouping information, the significance of the content differences of each metabolite was obtained by t-test calculations. R package ropls was used to create the OPLS-DA model, and to verify its reliability, 200 times of permutation tests were conduct. The VIP (variable importance in projection) values were calculated by multiple cross-validation. All of the above parameters including fold change, p-value and VIP obtained from OPLS-DA model, together determine the metabolites that differ between groups. The selection criteria were set as: fold change (FC) >1.5, p < 0.05, and VIP >1. Hypergeometric distribution test was employed to calculate the significance of KEGG pathway enrichment for differential metabolites.

Total RNA was extracted from the jejunal tissues of three wild-type and three MSTN gene-mutated cattle using the RNAiso Plus kit (9108; Takara), with purity and concentration assessed using a NanoDrop 2000 spectrophotometer. After confirming sample quality, eukaryotic mRNA was enriched by magnetic beads modified with Oligo(dT) to construct libraries; the first and second cDNA strands were synthesized and then the cDNA was purified; the purified double-stranded cDNA underwent end repair, A-tailing, and ligation of sequencing adapters, and then performed size selection using AMPure XP beads; finally, conducted PCR enrichment to obtain the cDNA library. After quality control, Illumina NovaSeq6000 was used for the library sequencing (in PE150 mode). The clean data obtained after filtering is subsequently aligned with the specified genome sequence to obtain mapped data. FPKM (fragments per kilobase of transcript per million fragments mapped) value was used as an indicator to normalize the transcript or gene expression levels. Principal component analysis (PCA) was performed to assess the sample dispersion. The DESeq2 was used as a tool to determine differentially expressed genes (DEGs) according to the number of genes in each sample, with a cutoff of fold change ≥1.5 and false discovery rate (FDR) value <0.05 as screening criteria. All screened DEGs were subsequently mapped to the KEGG database.

RNA from the tissue was isolated using an RNAiso Plus kit (9108; Takara), the RNA samples was subsequently reversed to cDNA by reverse transcription using a kit (R333-01; Vazyme). These cDNAs were used as templates for the next step, preparing a reaction system. In a 20 μL reaction system, the following components were mixed: cDNA, 8 μL; TB buffer, 10 μL; forward/reverse primer, 0.8/0.8 μL; sterile water, 0.4 μL. The reaction conditions were set as follows: initial denaturation (95°C for 30 s), 40 cycles of amplification (95°C for 10 s, 60°C for 30 s and 72°C for 30 s), with data collected during the elongation phase at 72°C. GADPH was served as the internal reference gene. qRT-PCR experiments were conducted, 2 − ΔΔCT method was used to calculate the differences in gene expression, with amplification conducted using a real-time PCR instrument. Primers used in this study were listed in Table 1.

The joint analysis of the metabolome and transcriptome was conducted using two-way orthogonal PLS (O2PLS) and Pearson correlation to integrate and analyze the associations between the transcriptome and metabolome (platform BMKCloud¹). A Venn diagram was applied to display the KEGG signaling pathways jointly enriched by the transcriptome and metabolome, while a bubble chart showed the significantly enriched signaling pathways (p ≤ 0.05).² A Sankey diagram illustrated the associations between the signaling pathways, genes, and metabolites (see text footnote 2).

LC-MS/MS was used to analysis the metabolites in the jejunal contents of WT and MT cattle detection platform. PCA was conducted to analyze the reproducibility between samples, and the results were displayed in Figure 1A. The MSTN gene-edited group samples clustered together, while the wild-type group (WT) samples formed another cluster, indicating good reproducibility within the groups and significant differences between them. PC1 and PC2 represented 35.6 and 23.2% of the total variations (Figure 1A). An orthogonal partial least squares discriminant analysis (OPLS-DA) model was constructed to analyze the differential metabolites between the MT and WT groups, with permutation tests conducted on the OPLS-DA model, as displayed in Figures 1B,C. A significant trend of separation between the two groups. The OPLS-DA model parameters R²Y and Q²Y were 1 and 0.8, respectively, proving that the OPLS-DA model was effective. These results indicate that there are significantly different metabolites in the jejunal contents of MT and WT cattle.

Metabolites with FC ≥1 were selected, and those with VIP ≥1 were identified as differential metabolites. Totally, 2,432 metabolites were detected, among them, the abundance of 304 metabolites showed significant differences between MT and WT samples (Figures 1D,E), of which 142 were upregulated and 162 were downregulated. A hierarchical clustering heat map of the selected differential metabolites suggested significant differences between the two groups (Figure 1F).

Differential metabolites were screened based on VIP values, and the top 10 differential metabolites in the MT and WT groups were listed in Supplementary Table S1, including acetylphosphate, 21,22-diprenylpaxilline, avermectin A1b aglycone, pyochelin, cholic acid, 11-dehydro-thromboxane B2, benzyl gentiobioside, N8,N′8-citryl-bis(spermidine), aucubin, and cobalt-precorrin 7. These metabolites are primarily involved in signaling pathways such as biosynthesis of metabolites (ko01110), biosynthesis of 12-, 14- and 16-membered macrolides (ko00522), biosynthesis of siderophore group nonribosomal peptides (ko01053), secondary bile acid biosynthesis (ko00121), arachidonic acid metabolism (ko00590), porphyrin metabolism (ko00860), and biosynthesis of cofactors (ko01240).

The metabolites with abundance differences in the top 20 in the MT and WT groups mainly include, carboxylic acids and derivatives, organooxygen compounds, steroids and steroid derivatives, prenol lipids, benzene and substituted derivatives, glycerolipids, glycerophospholipids, phenols, pteridines and derivatives, purine nucleosides, pyridines and derivatives, pyrimidine nucleosides, flavonoids, indoles and derivatives, naphthalenes, organonitrogen compounds, 2-arylbenzofuran flavonoids, azoles, and azolines (Figure 2A). Based on the differential metabolite scores, KEGG analysis was then performed. As was shown in Figure 2B. The KEGG-enriched metabolic pathways mainly included fatty acid biosynthesis, biosynthesis of type II polyketide products, fatty acid elongation, biosynthesis of various alkaloids, secondary bile acid biosynthesis, fatty acid degradation, nicotinate and nicotinamide metabolism, limonene and pinene degradation, atrazine degradation, biosynthesis of siderophore group nonribosomal peptides, dioxin degradation, nitrotoluene degradation, glycine, serine and threonine metabolism, phenylalanine, tyrosine and tryptophan biosynthesis, biosynthesis of unsaturated fatty acids, biosynthesis of various antibiotics, benzoate degradation, peptidoglycan biosynthesis, ether lipid metabolism, and pyruvate metabolism (Figure 2B).

The KEGG analysis of differential metabolites based on enrichment factors showed significant differences in the KEGG signaling pathways, including fatty acid biosynthesis, biosynthesis of type II polyketide products, fatty acid elongation, and biosynthesis of various alkaloids, as well as secondary bile acid biosynthesis (Figure 2C). The enrichment network diagram is shown in Figure 2D, where key differential metabolites in the fatty acid biosynthesis signaling pathway include lauric acid, decanoic acid, myristic acid, oleic acid, hexadecanoic acid, and palmitic acid. Key differential metabolites in the biosynthesis of type II polyketide products signaling pathway include urdamycin G, aquayamycin, landomycin D, C1′–C9-glycosylated UWM, and (S)-hemiketal. Key differential metabolites in the fatty acid elongation signaling pathway include palmitic acid and hexadecanoic acid. Key differential metabolites in the biosynthesis of various alkaloids signaling pathway include vanillylamine, L-tryptophan, benzoic acid, and shikimic acid. Key differential metabolites in the secondary bile acid biosynthesis signaling pathway include allocholic acid, alpha-phocaecholic acid, taurochenodeoxycholate, cholic acid, and bitocholate (Figure 2D).

RNA-seq was used to compare the gene expression differences in the jejunal smooth muscle tissues of MSTN gene-edited cattle (MT1, MT2, and MT3) and wild-type cattle (WT1, WT2, and WT3). The average reads of per individual was 42,005,317, with a range from 38,858,456 to 47,116,854. The GC content of the data was mostly around 50%, and the Q30 values were all above 93%, indicating a high quality of the sequencing data (Table 2).

PCA was used to analyze the reproducibility between samples, the results were showed in Figure 1A. The samples from the MSTN gene-edited group clustered together, while the wild-type group (WT) samples formed another cluster, indicates that samples within groups are well reproducible, while significant differences are manifested between groups. 32.31 and 28.94% of the total variations were attributed to PC1 and PC2, respectively (Figure 3A). Based on gene read count data, differential expression calculations were performed with significant criteria set as FDR <0.05, FC ≥1.50 or ≤0.67. The results showed that a total of 1,541 genes were differentially expressed in the MT group compared with the WT group, of which 536 genes were up-regulated and 1,005 genes were down-regulated (Figures 3B,C). Hierarchical clustering results indicated significant differences in the expression profiles of genes between the two groups (Figure 3D). The top 10 differentially expressed genes were FNDC1, NCOR2, LOC540321, LOC539009, LOC617219, LOC404103, PTI, NCKAP5L, ARHGEF33, and CIC. The signaling pathways they participated in primarily included TGF-beta signaling pathway (ko04350); notch signaling pathway (ko04330); Epstein–Barr virus infection (ko05169); autophagy-animal (ko04140); mTOR signaling pathway (ko04150); renin-angiotensin system (ko04614); salivary secretion (ko04970); spinocerebellar ataxia (ko05017) (Supplementary Table S2).

GO enrichment was performed to identify significantly enriched DEGs, revealing that they were significantly enriched in 44 GO terms (p < 0.05). The significantly enriched biological processes mainly related to cellular processes, biological regulation, response to stimulus, and metabolic processes. The significantly enriched cellular components were primarily related to cellular anatomical entities, intracellular components, and protein-containing complexes. The significantly enriched molecular functions were mainly associated with binding, catalytic activity, translation regulator activity, molecular function regulators, and transporter activity (Figure 3E). Compared to the wild-type jejunal tissue transcriptome data, the MSTN gene-edited cattle enriched 311 KEGG signaling pathways, including thermogenesis, ribosome, Huntington disease, thyroid hormone signaling pathway, lysine degradation, circadian entrainment, prion disease, focal adhesion, Parkinson disease, autophagy-animal, longevity regulating pathway-multiple species, Alzheimer disease, proteoglycans in cancer, aldosterone synthesis and secretion, protein digestion and absorption, ECM-receptor interaction, cardiac muscle contraction, adherens junction, circadian rhythm, and insulin secretion (Figure 3F).

The relationship between the transcriptome and metabolome was integrated using two-way orthogonal PLS (O2PLS), the results were displayed in Figure 4A, indicating significant associations between genes and metabolites. Correlation analysis and clustering of genes and metabolites from different samples revealed a significant clustering trend between the WT and MT group samples (Figure 4B). Pearson correlation analysis was conducted between the top 10 differential genes and the top 10 differential metabolites, revealing that FNDC1, LOC539009, NCOR2, NCKAP5L, ARHGEF33, and CIC were significantly positively correlated with acetylphosphate, cobalt-precorrin 7, benzyl gentiobioside, and aucubin; while they were significantly negatively correlated with 21,22-diprenylpaxilline, avermectin A1b aglycone, pyochelin, N8,N′8-citryl-bis(spermidine), cholic acid, and 11-dehydro-thromboxane B2. Conversely, LOC540321, PTI, and LOC404103 were significantly negatively correlated with acetylphosphate, cobalt-precorrin 7, benzyl gentiobioside, and aucubin; while they were significantly positively correlated with 21,22-diprenylpaxilline, avermectin A1b aglycone, pyochelin, N8,N′8-citryl-bis(spermidine), cholic acid, and 11-dehydro-thromboxane B2 (Figure 4C). Venn diagram analysis showed that there were 34 KEGG signaling pathways jointly enriched by the transcriptome and metabolome, with 277 KEGG signaling pathways unique to the transcriptome and 48 unique to the metabolome (Figure 4D). Among them, three pathways showed dramatic differences (p < 0.05) (Figure 4E), particularly pathways related to fatty acid metabolism including its biosynthesis, elongation and degradation. The differential genes and metabolites involved in these three signaling pathways were analyzed. Genes including ACACA and FASN which expressed differently between sample groups are involved in the fatty acid biosynthesis signaling pathway, while the main differential metabolites were hexadecanoic acid, oleic acid, lauric acid, myristic acid, decanoic acid, and palmitic acid. The differential genes involved in the fatty acid elongation signaling pathway included PTPLAD2 and THEM4, with the primary differential metabolite being palmitic acid. The differential genes involved in the fatty acid degradation signaling pathway included PTPLAD2, THEM4, ACOX3, and C7H5orf63, with the main differential metabolites being hexadecanoic acid and palmitic acid (Figure 4F).

Analysis of the genes and metabolites involved in the three significantly different KEGG signaling pathways jointly enriched by the transcriptome and metabolome showed that ACOX3, ACACA, and FASN were significantly downregulated in the jejunal tissues of MSTN mutant cattle (Figure 5A), and this results were proved by real-time PCR (Figure 5B). PTPLAD2, C7H5orf63, and THEM4 were significantly upregulated in the jejunal tissues of MSTN mutant cattle (Figure 5A). Hexadecanoic acid, lauric acid, and decanoic acid were significantly downregulated in the jejunal tissues of MSTN gene mutant cattle, while oleic acid, palmitic acid, and myristic acid were significantly upregulated (Figure 5C). Correlation analysis of metabolites and genes involved in fatty acid metabolism showed that myristic acid was significantly negatively correlated with ACOX3, ACACA, and FASN, while it was significantly positively correlated with PTPLAD2 and C7H5orf63; oleic acid was significantly negatively correlated with FASN, while it was significantly positively correlated with C7H5orf63 and THEM4; palmitic acid was significantly negatively correlated with ACOX3, ACACA, and FASN, while it was significantly positively correlated with PTPLAD2, C7H5orf63, and THEM4; hexadecanoic acid was significantly positively correlated with FASN and ACOX3, while it was significantly negatively correlated with C7H5orf63; lauric acid was significantly positively correlated with ACACA, while it was significantly negatively correlated with PTPLAD2, C7H5orf63, and THEM4; decanoic acid was significantly positively correlated with ACACA and FASN, while it was significantly negatively correlated with PTPLAD2, C7H5orf63, and THEM4 (Figure 5D).