In this study, we have used a total of 18 cattle belonging to Simmental (S), Simmental x Mongolian (SM) and Simmental x Holstein (SH), taking 6 animals from each group. All healthy and of similar age (480.33 ± 43.27 days) and weight (573.00 ± 30.60 kg). All groups were fattened at the Tianmuzhen Farm in Xing’an League, Inner Mongolia. The experimental cattle were housed in mixed pens and fed a total mixed ration twice daily at 08:00 and 18:00 h, with free access to feed and water during the 180 d fattening period. They were fed using the same feeding strategy (see Table 1) during fattening, the three head in each group closest to the average weight were selected for slaughter and the cattle were slaughtered at the slaughterhouse of Inner Mongolia Tianmuzhen Meat Industry Co.

Before slaughter, blood was collected from each cattle after fasting in non-anticoagulated and anticoagulated blood collection tubes from nine cattle and allowed to stand for 30 min, Centrifugation in a centrifuge (Centrifuge 80-2, Jiangsu Kangjie Medical Devices Co., Ltd. China) at 1076.7 rcf for 10 min to separate the serum and plasma, which were then divided into centrifuge tubes and stored at −80°C for the determination of serum biochemical indices (three repetitions per cattle). The levels of glucose (GLU), triglycerides (TG), total cholesterol (TC), high-density lipoproteins (HDL), very low-density lipoproteins (VLDL) and beta-hydroxybutyric acid (BHBA) were determined using Bovine ELISA kits (HDL:CK-E150097B;BHBA:CK-E92679B;VLDL:CK-E150082-48;GLU:CK-E95613B-48;TG:CK-E90126B;TC:CK-E93389B) purchased from Quanzhou Ruixin Biotechnology Co., Ltd. (China) according to the manufacturer’s instructions.

Nine 50 mg samples from each group were placed in a 2 mL grinding tube with steel beads and 1 mL of dichloromethane: methanol (1:1),and the sample was then ground in a freezer grinder. The sample was then sonicated for 15 min at 4°C, followed by 15 min of standing at −20°C. The supernatant was removed by centrifugation (Eppendorf, Centrifuge 5,424 R,MJPRO-YQ-015) at 13,000 rcf for 10 min and dried under nitrogen atmosphere. Then, 0.5 mL of sodium hydroxide in methanol (0.5 mol/L) was added, vortexed for 30 s, and sonicated at 60°C for 30 min in a water bath. After cooling, 0.5 mL of hexane was added, vortexed for 30 s, centrifuged at 4°C for 10 min, and 100 μL of the upper layer (hexane layer) was transferred to the injection vial and detected by GC–MS (23) using an 8890B-5977B (Agilent Technologies Inc. CA, United States) gas chromatograph. The default parameters of the Masshunter quantification software (Agilent, United States, v10.0.707.0) were used for the automatic identification and integration of each ionic fragment of the target fatty acids with the aid of manual inspection. A linear regression standard curve was plotted, with the analyte concentration as the horizontal coordinate and the mass spectral peak area as the vertical coordinate. The concentration of each sample was calculated by substituting the value into a linear equation.

Total RNA was extracted from the dorsal longissimus dorsi muscle tissues using TRIzol reagent (catalog number 15596026, Invitrogen, Carlsbad, CA, United States), according to the method described by Hao et al. (24). RNA quality was determined using a 5,300 Bioanalyzer (Agilent, Santa Clara, CA, USA) and quantified using an ND-2000 (NanoDrop Technologies). Only high-quality RNA sample (OD260/280 = 1.8 ~ 2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, >1 μg) was used to construct sequencing library. Using the eukaryotic A-tail structure, A-T base pairing was performed using oligomeric (dT) magnetic beads complementary to the A-tail (Illumina® Stranded mRNA Prep, Ligation, Illumina, United States). mRNA was isolated for transcriptome analysis, and DNA was digested with DNase I (TaKara) (Tiangen, Shanghai, China). Randomly interrupted mRNA, a small fragment of approximately 300 base pairs (bp) was isolated using magnetic bead screening (Illumina® Stranded mRNA Prep, Ligation, Illumina, United States). The mRNA was used as a template, and the cDNA strand was synthesized using a 6-base random primer. Double-stranded cDNA was purified using a double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina, United States); the purified double-stranded cDNA was repaired, attached to the a-tail and ligated to the sequencing connector. Library construction and RNA sequencing were performed by Shanghai Majorbio Bio-Pharm Biotechnology Co., Ltd. (Shanghai, China). cDNA target regions of 300 bp were selected from the libraries using 2% low-range ultra agarose, followed by PCR amplification using Phusion DNA polymerase (NEB) for 15 PCR cycles. The paired-end RNA-seq library was sequenced using a NovaSeq 6,000 sequencer (2 × 150 bp read length) after quantification using Qubit 4.0.

The raw sequencing data were subjected to quality control and the reads were mapped using fastp (25), The specific steps and sequences are as follows: (1) Remove the splice sequences from the reads, and remove the reads that do not have fragments inserted due to splice self-linkage and other reasons; (2) Trim off the low-quality (quality value less than 20: calculate the average quality value according to the sliding window size of 4) bases at the first end of the sequence (5″ end), and trim off the low-quality (quality value less than 3) bases at the end of the sequence (3″ end); (3) Removal of reads containing N (fuzzy bases): reads with more than 5 N are discarded; (4) Discard sequences less than 30 bp in length after de-adapter and quality trimming. After the completion of data quality control, the data after quality control were again subjected to statistics and quality assessment, which also included: base error rate distribution statistics and base content distribution statistics (Table 2).

As described by Man et al. (26). Subsequently, the mapped reads were individually aligned to the Bos taurus¹ genome using HISAT2 software (27). Finally, the reads were then assembled using StringTie using the method described in Pertea et al. (28) RSEM (29) software was used to quantify the gene and transcript levels. DESeq2 (30) was used to calculate differential expression with a threshold of (|log2FC| ≥ 1 and FDR < 0.05). GO and KEGG (31–33) were used for functional enrichment analyses to identify significantly enriched DEGs in metabolic pathways and GO terms.

Statistical analyses were performed using SPSS (34) software (SPSS v. 21, SPSS Inc.; Chicago, IL, United States). Significance of serum biochemical parameters and fatty acid contents. Statistical significance was determined by one-way analysis of variance (ANOVA), followed by the LSD multiple comparison test. The results in this study were presented as the means ± standard deviation (SD). p < 0.05 was considered statistically significant and indicated by “*,” whereas “**” indicated p < 0.01.

When comparing the blood serum indices between the different groups, as shown in Figure 1, the serum TG content in both the SM and S groups was significantly higher than that in the SH group (p < 0.05), the GLU content in the SM group was significantly higher than that in the S group (p < 0.05), and the differences in serum TC, BHBA, HDL, and VLDL contents between groups were not significant (p > 0.05).

The fatty acid contents of the dorsal longissimus dorsi muscles of the three groups of fattened cattle are shown in Table 3. There were 34 fatty acids detected in all three groups in total, which were classified by saturation. The contents of SFAs, unsaturated fatty acids (UFAs), monounsaturated fatty acids (MUFAs), and polyunsaturated fatty acids (PUFAs) were determined separately using GC–MS. The results showed that pentadecanoic acid and linoleic acid contents in the dorsal longissimus dorsi muscle of the SM group were found to be significantly higher than those of the S and SH groups (p < 0.05). The Docosanoic acid content in the S group was significantly higher than that in the two crossbreed groups (p < 0.01). The erucic acid contents in the dorsal longissimus dorsi muscles of the S and SM groups were found to be significantly higher than those in the SH group (p < 0.01). Eicosatrienoic acid levels in the S and SM groups were significantly higher than those in the SH group (p < 0.05). The UFA and MUFA contents in the SM group were significantly higher than those in the S and SH groups (p < 0.05).

As shown in Table 4, after matching the clean reads of each sample to the reference genome, 4.39–5.19% and 91.45–92.48% of the clean reads for each sample were repeatedly and uniquely mapped, respectively, with an average match rate of 96.68%.

Each sample was analyzed individually in order to identify DEGs. As per an empirical study by Du et al. (23), genes with p-values less than 0.01 and a fold change ≥2 were identified as DEGs. In the present study, 1,031 DEGs were identified between the SH and S groups. Of these, 251 genes were upregulated and 780 were downregulated, as shown in the volcano plot in Figure 2A(FC > 2, and FDR < 0.05), Specific FDR values are shown in Table 4. A total of 49 DEGs were identified between the SM and S groups, of which 18 were upregulated and 31 were downregulated, as shown in Figure 2B. Each point in the graph represents a specific gene, with the red points indicating significantly upregulated genes, blue points indicating significantly downregulated genes, and grey points indicating non-significantly different genes. After mapping all of the genes, it can be seen that the points on the left are genes that are differentially downregulated, and the points on the right are genes that are differentially upregulated; the closer to the two sides and the top the points are, the more significant the difference in expression. The genes upregulated in the SH group were 5-8-S rRNA and COX3, and the genes downregulated were CFDP2, PDK4, SLC9A2, SNED1, STRIP2 and KLF9. The genes upregulated in the SH group in the SM group were PLEKHH3, 5-8-S rRNA and COX3 and the genes downregulated were CFDP2, SNED1, PRSS2, VCPIP1,CFDP2.

The KEGG database was used to analyze gene functions, and link genomic information and functional information. Genes in the gene sets were classified according to the pathways in which they were involved or their functions. In Figure 3, the horizontal coordinate is the number of genes or transcripts in each pathway, and the vertical coordinate is the name of the pathway in the KEGG database. As shown in Figure 3, DEGs in the SH and S groups were enriched in six KEGG pathways: metabolism, genetic information processing, cellular processes, organismal systems, human diseases, and environmental information processing. DEGs in the SM and S groups were enriched in five major KEGG metabolic pathways: cellular processes, organismal systems, metabolism, human diseases, and environmental information processing.

KEGG pathway enrichment analysis was performed on the DEG gene sets and two sets of gene enrichment bubbles were obtained. The vertical coordinate is the name of the pathway and the horizontal coordinate is the ratio of the number of genes enriched in the pathway to the number of annotated genes. The size of each dot is proportional to the number of genes in the pathway, and the colors of the dots correspond to different P-adjusted regions. As shown in Figure 4A, the top 20 enriched pathways in the SH and S groups were colorectal and insulin signaling pathways. The pathways that were found to be more highly enriched in the SM and S groups were oxidative phosphorylation, thermogenesis, cardiac muscle contraction, and non-alcoholic fatty liver disease (Figure 4B).

We selected nine genes with significant differences between the SH, SM, and S groups and performed a relevance analysis of the blood biochemical parameters and fatty acid content of the longest dorsal muscle, as shown in Figure 5. In group SH, we found that HDL content was negatively correlated with the expression of 5-8-S rRNA, COX1 and COX3 and positively correlated with the expression of GHR and SLC9A2; UFA content was negatively correlated with the expression of 5-8-S rRNA and positively correlated with the expression of SLC9A2, GHR and SEND1; MUFA content was negatively correlated with the expression of 5-8-S rRNA and positively correlated with the expression of SLC9A2, GHR and SEND1.In group SM, HDL content in SM and S groups was negatively correlated with the expression of 5-8-S rRNA, COX1 and COX3; GLU levels were positively correlated with the expression of 5-8-S rRNA and COX3; UFA levels were positively correlated with the expression of 5-8-S rRNA and COX1; and MUFA levels were positively correlated with the expression of 5-8-S rRNA and COX1, and negatively correlated with the expression of SNED1.