The experiment was executed in strict accordance with the guidelines for the care and use of experimental animals at Jilin University. All of the experiments were approved by the Institutional Animal Care and Use Committee of Jilin University (license number: 201809041).

Kazakh cattle and Xinjiang brown cattle were provided by the Xinjiang Yili Yixin Cattle and Sheep Breeding Cooperative. Three 30-month-old adult Xinjiang brown cattle and three 30-month-old adult Kazakh cattle were selected for slaughter. These adult steers were slaughtered in according to the procedure of the slaughterhouse, and then longissimus dorsi muscle was collected. All of the samples were immediately snap-frozen in liquid nitrogen and stored at −80°C until RNA extraction.

Total RNA was extracted from each group (the Kazakh cattle group and the Xinjiang brown cattle group) using TRIzol (Invitrogen, NY, USA). A NanoDrop 2000 spectrophotometer (Thermo, Waltham, MA, USA) was used to evaluate the concentrations and quality of the RNA, and agarose gel electrophoresis was used to evaluate the integrity of the RNA. Use epicentre Ribo-ZeroTM kit to remove sample rRNA, then add Fragmentation Buffer to randomly interrupt rRNA-depleted RNA. Use rRNA-depleted RNA as a template to synthesize the first cDNA strand with six-base random hexamers, and then add buffer, dATP, dUTP, dCTP, dGTP, RNase H, and DNA polymerase I to synthesize the second cDNA strand. Use AMPure XP beads to purify cDNA. The purified double-stranded cDNA is then repaired, added with A and connected to the sequencing adapter, and then used AMPure XP beads for fragment size selection. Finally, the U chain is degraded, and finally a cDNA library is obtained by polymerase chain reaction enrichment. The libraries were sequenced by BioMarker Technologies (Beijing, China) using an Illumina HiSeq X Ten platform.

StringTie uses FPKM (fragments per kilobase of transcript per million fragments mapped) as an indicator to measure the level of transcripts or gene expression. edgeR is suitable for differential expression analysis between samples (groups). edgeR is a Bioconductor software package that studies the differential expression of repeated count data. In the process of detecting differentially expressed genes, the fold change (FC) represents the ratio of the expression levels between two samples (groups), and the false discovery rate (FDR) is measured by the p-value obtained by calibration. We finally chose fold change ≥1.5 and p value <0.05 as the screening criteria [3,6].

Before the usage, removing each sample out of storage first, the applicable amount of tissue sample was weighed; then put each sample into a container that it was prepared with precooled liquid nitrogen and ground them into powder. Four volumes of lysis buffer were added into each sample. Then, a high-intensity ultrasonic processor (Scientz, Ningbo, China) was used to sonicate each sample on ice through three sound waves. To remove the cell fragments, the mixture was centrifuged at 12,000×g at 4°C for 25 min. The supernatant was then transferred to a new centrifuge tube and the protein concentration was detected using the BCA kit (Beyotime, Shanghai, China). The concentration of protein is 0.25 μg/μL. The peptides were dissolved in liquid phase A for liquid chromatograph (LC) and separated by an EASY-NLC 1000 ultra-high performance liquid phase system. Using tandem mass spectrometry (MS/MS) under a Q Exactive with an applied electrospray at voltage of 2.1 kV, the peptides were subjected to a nanospray ionization (NSI) source. First, each protein mixture was digested with 5 mM dithiothreitol for 30 min at 56°C and treated with 11 mM iodoacetamide at 25°C in the absence of light for 16 min. The concentrations of urea put into the last protein samples were less than 2 M. Then 1:50 trypsin:protein was used for the samples digested overnight at 37°C, and digested at a mass ratio of 1:100 for 4 hours. Then, the peptides were desalt with vacuum freeze-dried and using Strata X c18 column (Phenomenex, Los Angeles, CA, USA). Use TMT kit (Thermo, USA) to reconstitute the peptides in 0.5 MTEAB according to the manufacturer’s protocol. Thaw the labeling reagent and dissolve it in acetonitrile. Incubate for two hours at room temperature, then concentrate, desalinate and dry under vacuum.

The secondary mass spectrum data were collected using MaxQuant (v1.5.2.8). The database was Bos taurus_ensembl_1906 (37538 sequences). The anti-library was added to calculate the FDR caused by random matching, and a common pollution library was added to the database to eliminate contaminated proteins in the identification results. The decoy checkbox in Mascot was selected to search the decoy database automatically and generate a random sequence database. The random sequence database and the real database were tested, and the raw spectra were obtained. Peptides with a significance score of 20 in the 99% confidence interval and peptides greater than “identity” in Mascot probability analysis were regarded as identified, which reduces the probability of false peptides. Every reliable protein identification was associated with at least one unique peptide. Tandem mass tags (TMT)-labeled proteomics analysis was performed twice, and the fold change of differentially expressed proteins (DEPs) was >1.5 (p<0.05) or <0.67 (p<0.05).

First, the bovine protein was re-annotated with the aim of analyzing the potential function of the protein. In brief, bovine proteins were mapped to many public databases, covering the NCBI nonredundant, Swiss-Prot/UniProt, gene ontology (GO), and Kyoto encyclopedia of genes and genomes (KEGG) databases. With all proteins as the background, we used the numbers of DEPs to calculate the p-value and Q-value, respectively, representing the significance of the enriched GO terms/KEGG pathways and the FDR. The p-values were calculated by Fisher’s exact test, and the Q-values were calculated by an R package called “q-value”. The threshold of significance was defined as an FDR of 0.05.

To verify the results of the TMT analysis, we used PRM to further analyze the results. The liquid gradient was set to increase from 6% to 25% in 40 minutes. After 40 minutes, the gradient increased from 25% to 35% over 12 minutes and from 35% to 80% over 4 minutes before holding at 80% for at least 4 minutes, with a constant flow of 500 nL/min. Using MS/MS at a 2.0 kV electrospray voltage, the peptides were placed in the NSI source. To detect the complete peptide, the full scan range of the primary MS in the Orbitrap was set to 350 to 1,000 m/z with a resolution of 70,000. The resolution of the secondary MS is 35,000. The peptide were then selected for MS/MS using NCE setting as 27. The AGC of the main MS was set to 3e6, and the maximum injection time was set to 50 ms. The AGC of the secondary MS was set to 1e5, and the maximum injection time was set to 200 MS. The isolation window was set to 1.6 m/z.

All transcription and proteomic samples were planned for two biological replications. The results were expressed as mean±standard error of the mean. The student’s t test was used to analyze the differences between the two groups. When p<0.05, difference was considered statistically significant.

To promote understanding of the entire experimental process, a flowchart was constructed. Three adult Xinjiang brown cattle and three adult Kazakh cattle were selected for slaughter. These cattle were slaughtered according to the procedure of the slaughterhouse, and then the longissimus dorsi muscle was collected. Later, we extracted the RNA from the sample, and the RNA was used to construct a cDNA library. The libraries were sequenced by BioMarker Technologies (Beijing, China) using an Illumina HiSeq X Ten platform. We performed differential expression analysis of RNA obtained from RNA-seq.Later, the protein was extracted from the sample, digested with trypsin and labeled with TMT. Then, high performance liquid chromatography classification, LC/MS, database searching and MS quality control analysis were performed, and differential expression analysis of the proteins obtained from protein sequencing was performed. Finally, we performed joint analysis of differentially expressed RNA and DEPs (Figure 1).

In this project, 22,677 transcripts and 1,874 proteins were identified through quantitative analysis of the transcriptome and proteome. By comparing the identified transcriptome and proteome, we found that 1,737 genes were identified at both the transcriptome and proteome levels (Figure 2A). The levels of transcriptome and proteome are significantly related (Figure 2B). When Log2 FC>0.585 and p value<0.05, there was a significantly upregulated transcript; when Log2 FC <−0.585 and p value<0.05, there was a significantly downregulated transcript. When the ratio >1.2, there was a significantly upregulated protein; when the ratio <1/1.2, there was a significantly downregulated protein. Through the above screening conditions, we obtained 929 transcripts with significant differential expression, of which 471 were upregulated and 458 were downregulated. There were 119 proteins that were significantly differentially expressed, of which 75 were upregulated and 44 were downregulated (Figure 2C).

We selected differentially expressed RNA and protein from 929 transcripts and 119 proteins (Supplementary Table S1 and S2). We compared the relationship between RNA and protein between Xinjiang brown cattle and Kazakh cattle. The expression of RNA and protein was mostly the same, but there were also differences. For example, the RNA and protein expression levels of troponin I1 (TNNI1), crystallin alpha B (CRYAB), cysteine and glycine rich protein 3 (CSRP3), nidogen 2 (NID2), dihydropyrimidinase like 2 (DPYSL2), and glutamate-oxaloacetic transaminase 1 (GOT1) were all upregulated, but those of phosphotriesterase-related (PTER), myosin-binding protein H (MYBPH), and glutathione s-transferase mu 3 (GSTM3) were all downregulated. The RNA expression of myosin light chain 3 (MYL3) was upregulated, but their protein expression was unchanged. The RNA expression of aspartoacylase (ASPA) and receptor accessory protein 5 (REEP5) were downregulated, but its protein expression were unchanged (Figure 3).

The GO annotations are divided into 3 categories: biological processes, cell composition, and molecular functions. Fisher’s exact test was used to test the DEPs on the background of the identified proteins, and a GO enrichment test p-value less than 0.05 was considered significant. The DEPs and transcript were enriched in the collagen trimer, transferase activity transferring pentosyl groups, haem oxygenase activity and flavin adenine dinucleotide binding GO terms (Figure 4).

The KEGG database was used for enrichment analysis of pathways. The Fischer exact two-sided test method was used to test whether the DEP was based on the identified protein as the background, and a p-value less than 0.05 in the pathway enrichment test was significant. Finally, these channels were classified according to the KEGG website channel-level classification method. The DEPs and transcripts were enriched in adrenergic signaling in cardiomyocytes, hypertrophic cardiomyopathy, dilated cardiomyopathy, cardiac muscle contraction, PI3K-Akt signaling pathway, AMPK signaling pathway, peroxisome proliferators-activated receptor signaling pathway, fatty acid degradation and histidine metabolism (Figure 5).

In this experiment, we performed PRM quantification of the 20 selected target proteins in 6 samples. Due to the characteristics of protein and the abundance of their expression, we ultimately quantified 7 of them. The protein expression of TNNI1, CRYAB, CSRP3, ASPA, NID2, DPYSL2, and GOT1 was verified by PRM. The protein expression of TNNI1, CRYAB, CSRP3, NID2, DPYSL2, and GOT1 was upregulated, and the protein expression of ASPA was unchanged (Figure 6A). The PRM results were consistent with the protein sequencing results (Figure 6B).