Overview of transcriptome data and miRNA-seq data
After quality control and data filtering, clean RNA-seq paired-end reads with a length of 90 bp were obtained from longissimus dorsi muscles, and all Q20 percentages were more than 95%. More than 99% of clean reads were uniquely mapped to donkey, cow, and goat reference genome (EquDZ1.0, ARS-UCD1.2, and ARS1, respectively). Detailed information of other groups is shown in Additional file 1 .
Small RNA libraries were constructed and sequenced, which revealed the characteristics of miRNAs in skeletal muscle of donkey. After quality filtering and trimming of contaminant and adaptor sequences in donkey muscle group, a total of 17,807,876, 16,685,260, and 16,601,535 clean reads were obtained (approximately 95% of total reads). Clean reads were mapped to the donkey genome to remove ncRNA and repeat sequences. Then, mature miRNAs were identified by mapping the final reads to miRBase. Results showed that a total of 312 known miRNAs were expressed in donkey skeletal muscles. Meanwhile, 181 novel miRNAs were identified (Additional file 2) and blasted to the 40 predicted novel miRNAs [23]. The result showed 171 novel miRNAs in donkey muscle.
Differential gene expression and functional annotation clustering
GO enrichment analysis and KEGG pathway annotation were conducted on highly expressed genes in donkey meat, beef, and mutton. A total of 539 highly expressed genes (FPKM >100) were found in donkey muscle, and 157 GO terms were considerably enriched (Additional file 3). Energy metabolism, sarcomere organization, muscle contraction, and skeletal muscle fiber development were involved in biological process terms. In addition to making up the main organelle components, many genes participated in the basic blocks of myofibrils (I band, M band, A band). Genes involved in molecular functions such as troponin T, calmodulin, and actin binding play an important role in muscle contraction and relaxation. Results of KEGG signaling pathway annotation (Additional file 4) showed that ribosome gene pathway and oxidative phosphorylation signaling pathway were substantially enriched.
To make gene expression among the three species comparable, we first identified 5160 putative single-copy orthologous genes in cow, donkey, and goat. A total of 1338 and 1780 DEGs were identified and further visualized in a volcanic map in cow vs donkey and goat vs donkey (Fig. 1 and Additional file 5, 6).
To explore the functional implications of these DEGs, we performed gene set enrichment analysis to identify GO terms. In cow vs donkey, a total of 196 significant GO terms (p < 0.05) were annotated, including 101 biological processes, 36 cellular components, and 59 molecular functions (Additional file 7). The significant GO annotation terms between goat and donkey muscles are presented in Additional file 8. Most DEGs were involved in biological processes, including small GTPase-mediated signal transduction, protein ubiquitination, protein glycosylation, protein transport, protein folding, and MAP kinase tyrosine/serine/threonine phosphatase activity. In the goat and donkey muscle comparison group, five genes were annotated in skeletal muscle contraction (GO: 0003009), including TCAP (Telethonin), CCDC78 (Coiled-coil domain-containing protein 78), JSRP1 (Junctional sarcoplasmic reticulum protein), TNNI2 (Troponin I, fast skeletal muscle), and STAC3 (SH3 and cysteine-rich domain-containing protein 3). To shed further light on putative functions of DEGs, we mapped them to terms in the KEGG database. The predicted KEGG pathways are shown in Additional file 9 and file 10. Most DEGs were significantly assigned to oxidative phosphorylation, pyrimidine metabolism, and purine metabolism pathways in the two comparison groups. In addition, some DEGs were enriched in amino acid and fatty acid metabolism pathways. TGF-β, Wnt, and AMPK signaling pathways and regulation of actin cytoskeleton were also annotated.
TFs and transcription cofactors play a role in gene expression. All DEGs in two compared groups were aligned to animal TFs database (Animal TFDB 2.0). Among the DEGs, at least 464 genes were identified to encode for 39 TF families (Additional file 11). A total of 156 TFs were up-regulated, and 308 TFs were down-regulated. To gain insights into the functional interactions between DEGs, we created protein-protein interaction networks, and detailed distinct clusters are shown in Additional file 12, file 13 and file 14. We presented the top six and seven clusters with more enriched genes in the two compared groups. Results showed larger interaction networks in the goat vs donkey group, and fewer genes interacted in the cow vs donkey group. There were same clusters with node genes: PCNA (Proliferating cell nuclear antigen), MRPL2 (Mitochondrial ribosomal protein L2), and NADH dehydrogenase compounds. These cluster genes were mainly related to cell proliferation, ribosomal synthesis, and energy generation.
AS and SNP analysis in donkey muscle
At the transcriptome level, AS and SNPs can increase the abundance of transcripts and thus express more complex proteins. A total of 25,022 AS events were detected in donkey muscle under junction count only mode (Table 1). Among the five events, skipped exon (SE) and alternative 3ʹ splice site (A3SS) events predominated, accounting for 76.72% of alternative transcripts. Intron retention (RI) was the least frequent, accounting for 0.92%. We also analyzed statistically significant AS events. As shown in Fig. 2 and Additional file 15, 396 different events were determined. A total of 118 AS events occurred in at least two compared groups, and nine different AS events were detected in three comparison groups. Among them, NEB (nebulin) gene occurred in six different events including two MXE and four SE. Other AS events occurred on genes including actin, actinin, tropomyosin, integrin, myosin, titin, tropomodulin, and protein phosphatase. In comparing the different AS events of homologous genes among the three kinds of muscle samples, 24 AS events occurred in 16 genes in donkey muscle samples. Six events were different among the three genes, including TMOD1 (Tropomodulin-1), TPM2 (Tropomyosin beta chain), and TMOD4 (Tropomodulin-4). Protein interactions were observed among these genes, and they participate in tropomyosin binding (GO: 0005523), actin filament binding (GO: 0051015), and striated muscle thin filament (GO: 0005865).
Table 1. All AS events were identified in three donkey muscles.
Event type
|
SE
|
A3SS
|
A5SS
|
MXE
|
RI
|
Events Num
|
12,300
|
6896
|
4235
|
1338
|
230
|
Using sequence reads aligned to the donkey reference genome, we identified a total of 57,201 putative SNPs. SNPs were then classified to several categories on the basis of the following criteria: position in intergenic, up-, or downstream regions of a gene and their effect on gene coding (synonymous or nonsynonymous polymorphisms) according to the donkey genome annotation (Table 2). We identified transitions and transversions, but most of the nucleotide changes could be classified as transitions with a transition/transversion ratio of 2.62 (Additional file 16). The majority of SNP variations were located in known genes, even if there were 12,866 SNPs falling in the intergenic position. Results showed 2306 SNPs in the exon region. In these mutations, only one synonymous SNP was annotated, whereas three nonsynonymous SNPs were annotated. The function of other SNPs could not be determined. The reason for this result may be that the donkey samples are closely related to the reference genome (two of the samples came from the parents of the reference donkey).
Table 2. SNP distribution in accordance with gene locations in three donkey muscles.
Upstream
|
Splicing
|
Exonic
|
Intronic
|
Downstream
|
Upstream/downstream
|
Intergenic
|
ncRNA
|
Total
|
1614
|
3
|
2306
|
618
|
7638
|
116
|
12,866
|
32,040
|
57,201
|
Overview of proteomic analysis
A total of 388,206 spectra were detected from iTRAQ experiment using donkey, cattle, and goat muscle as materials. The data were analyzed using Mascot software version 2.3.02. In this study, Mascot identified 62,140 known spectra. After deleting low-scoring spectra, 52,329 unique spectra that were matched to 9119 peptides, 8615 unique peptides, and 1768 proteins were identified with FDR of less than 1%. The annotation information of all proteins is shown in Additional file 17.
Identification and functional annotations of DEPs
Among the 764 DEPs in the cattle vs donkey group, 420 proteins were up-regulated, and 344 proteins were down-regulated with 99% confidence and a 1.2-fold change cutoff. In goat vs donkey muscle, we identified 1206 DE proteins (839 up-regulated and 367 down-regulated). Cluster analysis of DEPs were shown in Fig. 3. All DE protein informations were presented in Additional file 18 and file 19.
Functional annotation of DEPs
Functional analysis of all DEPs was carried out using three criteria of protein functional annotation according to the GO database. The results of GO classifications (p < 0.05) in the cow vs donkey group are presented in Additional file 20. A total of 80 categories of biological process were annotated, in which most proteins were related to small molecule metabolic process. Moreover, seven DEPs participated in aspartate family amino acid metabolic process, including CTH, ASNS, AHCYL1, MTHFD1, GOT1, AHCY, and ASS1. Ten DEPs were annotated to regulation of fatty acid metabolic process and cholesterol biosynthetic process (PDK4, FABP3, EIF6, PLIN5, Adipoq, ACLY, APOA1, G6PD, CYB5R3, ACAA2). In addition, there were four proteins annotated in Acetyl-CoA metabolic/biosynthetic process (PDHA1, ACLY, ACSS1, MPC2). Only 19 terms corresponded to molecular function, and 16 terms corresponded to cellular component. In the goat vs donkey group, 107 significant GO terms were assigned to DE proteins, which were classified into independent categories: 57 terms corresponding to biological processes, 20 terms corresponding to molecular functions, and 30 terms corresponding to cellular components (Additional file 21). A total of 21 DEPs were involved in fatty acid catabolic process and fatty acid metabolic process; fatty acid beta-oxidation, fatty acid oxidation, and lipid oxidation.
KEGG pathway analysis of DEPs in the two comparison groups (Fig. 4 and Additional file 22 and 23) showed that the enriched pathways were also mainly associated with metabolism, including Oxidative phosphorylation, Fatty acid metabolism, Biosynthesis of amino acids, Citrate cycle (TCA cycle), and several types of amino acid metabolism (p < 0.05).
Correlation of proteomics and transcriptomics
The gene expression was spatially and temporally regulated at the transcriptional, post-transcriptional, and translational levels. Thus, combined transcriptomic and proteomic analysis was performed. Integrating RNA-seq and iTRAQ data, the overall expression of mRNAs and proteins is shown in Fig. 5. We found 59 and 129 genes with significant differences at the mRNA and protein levels between the two comparison groups (Additional file 24 and file 25).
We analyzed four types of mRNA and protein expression patterns, including transcripts and proteins in concordance (up-regulated DEGs corresponding to up-regulated DEPs, down-regulated DEGs corresponding to down-regulated DEPs), transcripts and proteins with opposite expression (up-regulated DEGs corresponding to down-regulated DEPs, down-regulated DEGs corresponding to up-regulated DEPs). In cow vs donkey, the expression of 24 genes was consistent in mRNA and protein levels, whereas that of 35 genes showed the opposite effects on mRNA and protein levels (Fig. 6A). In goat vs donkey, 129 genes were differentially expressed, of which 36 and 11 genes were consistently up- or down-regulated, respectively. As shown in Fig. 6B, 82 genes had inconsistent expression in mRNA and protein levels. We conducted regulatory network analysis using STRING to explore protein-protein interactions for analyzing the regulatory relationship between proteins and TFs. The interaction of DEGs in cow vs donkey was shown in Fig. 7. We annotated nine TFs, including CRKL, NME4, OLFML3, HSF4, Ppplr27, Psmc3ip, NEM2, Rab2b, and PURB. The regulated genes of only NME4 and HSF4 were found. The significant KEGG signaling pathways included oxidative phosphorylation, metabolic pathways, fatty acid degradation, arginine, proline, and cholesterol metabolism (Additional file 26). In goat vs donkey, 14 TFs were annotated, including CALR, UBE2M, CIRBP, YWHAH, YWHAZ, NME4, RAN, UBE2D3, RRAS, SERBP1, RHOC, ECHDC3, RAD21, and NPM1. The PPI network is shown in Additional file 27. Five genes (Tnni2, TMOD1, CFL2, MYOZ3, and MYL12B) were associated with muscle fiber synthesis and differed at transcriptional and proteome levels. CFL2 (Cofilin-2), an actin-binding protein, is the only specific isoform that is expressed in mature skeletal muscles [36]. This protein was localized near the M-bands in myofibrils. It can enhance actin filament conversion and regulate the length of actin filaments, which play an important role in the development of muscle morphology [37, 38]. Tnni2 (Troponin I2, fast skeletal muscle) is a subunit of troponin complex and plays a role in calcium regulation during muscle contraction and relaxation. The SNP mutation (g. 1167 C>T) of Tnni2 gene was significantly associated with pH, meat color value, and intramuscular fat content in pigs [39]. The other three genes play important roles in maintaining the stability of myosin and cellular integrity. Significant KEGG signaling pathways of 129 genes were associated with energy metabolism, including oxidative phosphorylation; metabolic pathways; glutathione metabolism; arginine, proline, and tryptophan metabolism; purine metabolism; and glycolysis/gluconeogenesis (Additional file 26). In the two compared groups, more subunits of NADH dehydrogenase and cytochrome genes, which play a role in muscle energy metabolism were found.
qRT-PCR confirmation of DEGs and MRM validation of DEPs
To validate the expression profiles obtained using RNA-seq, quantitative real-time PCR was performed on 15 genes, which was normalized by GAPDH using mRNA samples from donkey, cow, and goat muscles. The 15 DEGs included six up-regulated genes (DUSP1, CA3, LAPTM4B, CARNS1, STAC3, COX5B) and nine down-regulated genes (PACSIN3, FGL2, TNNI2, GSTM3, HSPA13, ATF4, SUMO2, DGUOK, RBP1). qRT-PCR and RNA-seq data were significantly correlated with a Pearson correlation coefficient r of 0.72 in the cow vs donkey group and 0.84 in the goat vs donkey group (Fig. 8). Therefore, qRT-PCR analysis confirmed that the selected DEGs were differentially expressed in three muscles, demonstrating that the high-throughput sequencing was accurate. The primers used for these genes are shown in Additional file 28.
We chose MRM methods to verify the reliability of protein iTRAQ results. A total of 40 DEPs were selected to establish an MRM method. Results showed that 19 proteins have MS/MS spectra and unique peptides. Therefore, the MRM detection and analysis were performed for these 19 DEPs (Additional file 29). The corresponding correlation coefficient r was 0.75 and 0.78 for protein expression using MRM and iTRAQ methods, respectively (Fig. 9 and Additional file 30). The protein expression levels in MRM were basically consistent with the iTRAQ expression pattern. Therefore, MRM analysis confirmed that the proteomic analysis was reliable.