Changes between longissimus dorsi muscle from donkey, cow, and goat assessed by transcriptomic and proteomic analyses

Background: Domestic donkeys (Equus asinus) are farm animals; they are used mainly for carrying loads in the past centuries. Recently, interests on donkey milk and meat production have increased in many countries. Donkey meat is extremely popular in China due to its high nutritional value and unique avor. Except the origin and domestication of donkeys, few genetic studies have been conducted. Moreover, data from transcriptional proling and microRNA (miRNA) regulation of skeletal muscle tissues in donkey are scarce. Recent developments in high-throughput sequencing techniques can offer large-scale analysis of gene expression in different species. This study aimed to explore the differences of the molecular and regulation mechanisms among donkey meat, beef, and mutton using genome-wide transcriptomic analysis and proteomic methods to provide more effective genetic information. Methods: RNA sequencing and proteomic analysis (on the basis of isobaric tags for relative and absolute quantitation) on donkey, cow, and goat muscles, and miRNAs in donkey muscles were detected. We performed a comprehensive research on total RNA, including miRNAs from donkey muscles. The mRNA expression proles were characterized and differences in single-copy homologous genes among species were analyzed. Results: Most differentially expressed genes were associated with pathways related to protein and fat synthesis and metabolism, muscle formation, and development. We identied single nucleotide polymorphisms (SNPs) in donkey muscle and alternative splicing (AS) events. A total of 57,201 putative SNPs were found, and the main SNP variations were located in known genes. Several AS events occurred in genes related to muscle ber. Different AS events were also noted among species. Muscle proteomic data were obtained, including all expressed proteins and differentially expressed proteins. Combined transcriptomic and proteomic analysis effectively revealed pivotal mRNAs and proteins in muscles. We found ve genes associated with thin and thick laments, which indirectly explained the characteristics of donkey muscle bers. Conclusion: RNA-seq and iTRAQ analysis revealed altered expression of genes and proteins in three kinds of muscle. In the highly expressed genes of donkey muscle, we discovered 31 expressed genes were involved in muscle contraction and skeletal muscle ber development. Compared with mutton, ve genes related to muscle ber synthesis were found and showed differences 3 (cid:0) splice site, RI: Intron retention, A5SS: alternative 5 (cid:0) splice site, MXE: Mutally exclusive exons, DEP: Differentially expressed protein, PCNA: Proliferating cell nuclear antigen, MRPL2: Mitochondrial ribosomal protein L2, NADH: Mitochondrially encoded DADH dehydrogenase 4, NEB: Nebulin, TMOD1: Tropomodulin 1, TPM2: Tropomyosin 2, TMOD4 : Tropomodulin 4, CFL2: Colin 2.

human milk for infant nutrition [3,4]. Donkey meat consumption has substantially increased in China due to favourable fatty acid pro le, characterized by signi cant proportions of unsaturated fatty acid, particularly polyunsaturated fatty acids (PUFAs) [5]. The skeletal muscle makes up approximately 40% of the total body weight of donkeys. The composition of donkey meat was measured in terms of major chemical components such as proteins, fats, glycogen, and minerals [6]. In a previous study investigating the physical (color) and chemical (fatty and amino acids content) characteristics of two donkey muscles, the results showed a high content of PUFAs and essential amino acids, especially glutamine and lysine [7]. You et al. [8,9] compared amino acid and fatty acid compositions between donkey meat and meat of other livestock including goat, cattle, and pig. The results indicated that donkey meat had the highest content of glutamine; essential amino acids, especially lysine; and PUFAs compared with mutton, beef, and pork. Given its low fat content, high digestibility, and good proportion of unsaturated/saturated fatty acids, donkey meat can be considered as an optional red meat for the production of salami and other fermented meat products.
Meat quality traits, including meat color, muscle ber diameter and density, tenderness, avor, and IMF content, have been a major focus of genetic improvement programs. These features are closely associated with meat composition. Development of muscles and adipose tissues are complex physiological processes. At present, the main constituents of donkey meat are known and are different from other meats. However, few molecular biology studies have analyzed donkey meat. Skeletal muscle development is a complex process that requires the concerted expression and interaction of multiple factors [10].
A number of studies have employed high-throughput analysis, proteomics, and metabolomics to explore the molecular mechanism of controlling meat quality and understanding gene regulatory networks of important biological processes [11,12,13]. RNA sequencing (RNA-seq) based on deep sequencing for transcriptome pro ling has considerable advantages in examining transcriptome data. Bongiorni et al. [14] sequenced mRNAs of two cattle skeletal muscles to identify genes affecting meat tenderness and correlated metabolic pathways.
Proteomics has been applied in meat science to achieve a better understanding of molecular mechanisms of meat quality [15]. Isobaric tags for relative and absolute quantitation (iTRAQ) labeling is another e cient quantitative proteomic technology that uses a family of isobaric isotope tags to label tryptic peptides, followed by liquid chromatography-tandem mass spectrometry (LC-MS/MS), which is an important reference for biomarker discovery [16,17]. Integrative analyses using high-throughput technologies, including cDNA microarray, 454-sequencing, and iTRAQ-based proteomics, have revealed the regulation of transcript and protein levels across heart and muscle tissues [18]. Zhang et al. [19] investigated the proteomic alteration of longissimus muscle in growing pigs after supplementation with 0.6% non-starch polysaccharide enzymes (NSPEs) in the diet using iTRAQ. A total of 51 proteins were differentially expressed between the control and NSPE groups, and functional analysis showed a change in protein abundance of different pathways in NSPE-fed pigs. Bjarnadottir et al. [20] researched protein changes between tender and tough meat from bovine Longissimus thoracis muscle using iTRAQ and 2-DE. They found three novel differential proteins related to tenderness.
Over the past decade, microRNAs (miRNAs) are the most abundant type of noncoding regulatory RNAs, which degrade or inhibit the translation of their target genes by binding to the 3′-untranslated region of coding mRNAs [21]. Previous studies showed that miRNAs could regulate myogenesis or synthesis of various nutrients. Skeletal muscle-speci c miR-1, miR-133, and miR-206 play roles in myoblast proliferation and differentiation [22].
A total of 808, 715, and 267 known miRNA sequences from Bos taurus, Equus caballus, and Ovis aries, respectively, are deposited in Sanger miRNA database release 21.0 (miRBase). Huang et al. [23] identi ed 118 known miRNA families and 40 novel miRNAs speci c to donkeys from eight types of donkey tissues, including heart, liver, spleen, lung, kidney, brain, spinal cord, and muscle. Xie et al. [24] assembled the transcriptome of donkey white blood cells de novo and identi ed three ear size-associated proteins, namely, HIC1, PRKRA, and KMT2A (the outer ears of donkey are notably larger than those of horses).
In the present study, we successfully assembled and annotated chromosome-level donkey genome from one male Dezhou donkey, which was used as a reference for our analysis. To investigate molecular differences and interaction networks among donkey meat, beef, and mutton, we investigated the mRNA, miRNA, and protein pro les among these three kinds of meat. We brie y analyzed the miRNA-mRNAprotein regulation network to enhance our understanding on the effect of key genes on muscle characteristics in donkey. This association analysis will help improve our understanding of the phenotypic differences between donkey meat and other meats.

Sample preparation
The study was performed using nine adult individuals, including three Dezhou donkeys, three Holstein cows, and three Boer goats. Donkeys were obtained from the Black Donkey Research Institute (Liaocheng, Shandong, China). Cows and goats were reared in the same area. All procedures involving animals were carried out in accordance with the guiding principles for the care and use of laboratory animals and were approved by the Institutional Animal Care and Use Committee at Shandong Academy of Agricultural Sciences. After slaughtering the animals, approximately 300 g of muscle sample was taken from the longissimus dorsi between the 12th and 13th rib from the left side of the carcasses. All samples were trimmed and cut into small pieces, cleaned with RNase-free water, and immediately stored in liquid nitrogen for further processing. The samples were stored in tubes and sent to Beijing Genome Institute (BGI-Shenzhen) for library construction and RNA-seq, proteomic, and multiple reaction monitoring (MRM) quantitative analyses.
RNA extraction, library preparation, and deep transcriptome sequencing Total RNA extraction, mRNA puri cation, and cDNA library construction were conducted at BGI-Shenzhen. Brie y, total RNA was isolated separately from nine samples using TRIzol reagent according to the manufacturer's protocol. Following puri cation, Dynal oligo (dT) beads (Invitrogen) were added, mRNA was cut into small fragments, and the rst cDNA strand was synthesized using random hexamer primers. Double-stranded cDNA fragments were processed and ligated with Illumina's adaptor. The products were puri ed and enriched by PCR ampli cation. The average insert size for paired-end libraries was ~200 bp. The nine ampli ed libraries were sequenced using Illumina HiSeq™ 2000 sequencing platform (Illumina, CA, USA).
Small RNA library construction and sequencing from three donkey muscle samples Small RNA library was constructed as follows: rst, total RNA fractions that range from 18 nt to 30 nt in length were isolated by denaturing polyacrylamide gel electrophoresis. Second, small RNA fractions were ligated to a pair of adaptors at the 5 and 3 ends according to the manufacturer's instructions (Illumina). Third, reverse transcription was used to create cDNA, which was ampli ed using PCR ampli cation.
Fourth, the PCR products were puri ed, and the library was validated. The puri ed cDNA was used for sequencing on Illumina HiSeq 2000 (Illumina) following the manufacturer's instruction.
Raw sequencing reads and small RNA categories were obtained. Sequences were blasted to the NCBI GenBank and RFam database to annotate all known rRNA, scRNA, snoRNA, snRNA, tRNA, other noncoding RNAs, repeat sequences, and mRNAs. The unannotated sequences were aligned to miRBase database (Release 21, http://www.mirbase.org/ftp.shtml) to identify the conserved miRNAs using BowTie software (v. 2.0.2). Only small RNAs whose mature and precursor sequences perfectly matched all the species in miRBase were considered to be known miRNAs in donkey muscle. To discover novel miRNA precursor sequences, the anking sequence of unique sequences was used for secondary structure analysis and evaluated by Mireap (http://source forge.net/projects/mireap/). For miRNAs, the nal target genes were predicted using miRanda, RNAhybrid, and TargetScan.

Identi cation of differentially expressed genes (DEGs)
To improve data quality, high-quality clean reads were generated from raw reads by removing adapter sequences, duplicated sequences, and sequences that showed more than 10% unidenti ed nucleotides using lter_fq software. Then, the clean reads of each group were mapped to their respective genome (Equus asinus, Bos taurus, Capra hircus) using BWA with default parameters. The gene expression intensity was calculated using the FPKM method [25]. Orthologous genes were used for normalization to perform comparison of gene expression in the three species.
For each pair of species, an all-vs-all sequence similarity search was conducted on the protein sequence using BLASTP with an e-value cutoff of 10 −5 and an identity threshold of 40%. The BLASTP results were imported into OrthoMCL software for orthologous group clustering using default settings [26]. Finally, we identi ed DEGs using gene expression after normalization with fold change ratio >1.5 and p value <0.05 [27] in each comparable group.

Functional annotation of DEGs
Functional analysis of DEGs was performed using DAVID v6.7 (http://david.abcc.ncifcrf.gov/) [28]. Gene Ontology (GO) and KEGG metabolic pathway enrichments were estimated using modi ed Fisher's exact test. Differential gene expression analysis was integrated with protein-protein interactions using Cytoscape version 3.6.1. The interactions were taken from STRING database version 11.0 (http://stringdb.org/). The criterion for cluster genes was log2FC >±1.5, and all gene interactions with a connection reliability >700 were considered.

Alternative splicing (AS) and single nucleotide polymorphism (SNP) analyses
The rMATS software was used to summarize the splicing forms in three donkey muscles [29]. We examined ve splicing events and differences in AS events using junction count only quantitative method. False discovery rate (FDR) <0.05 was used as screening criterion for different AS events. SNPs of donkey muscles were identi ed using GATK v3.6 program [30], and the standard used was the same as that in a previous study: the value of SNP is more than 2.0, and single-base mismatch occurs from 3 to 35 bp.
Protein extraction, trypsin digestion, and iTRAQ labeling Frozen samples were ground and extracted with Lysis Buffer 3 (8 M Urea, 40 mM Tris-HCl or TEAB, pH 8.5) containing 1 mM PMSF and 2 mM EDTA ( nal concentration). After placing on ice for 5 min, 10 mM DTT was added to the solution. Then, the solution was centrifuged at 4 °C and 25,000 g for 20 min. The supernatant was incubated at 56 °C for 1 h. Subsequently, after cooling to room temperature, the sample was incubated with 55 mM IAM for 45 min in the dark for alkylation. Protein samples were collected after centrifugation at 4 °C and 25,000 g for 20 min. Protein concentrations were measured using Bradford assay (Bio-Rad, Berkeley, CA). Approximately 30 µg of total protein from each sample was subjected to SDS-PAGE analysis. Then, 100 µg of protein from each sample was digested with Trypsin Gold (Promega, Madison, WI, USA) at 37 °C overnight with protein-to-trypsin ratio of 40:1. After trypsin digestion, the peptides were desalted and vacuum dried. They were dissolved in 30 µL of 0.5 M TEAB. iTRAQ labeling reagents were transferred and combined with corresponding samples. Peptide labeling was performed using the iTRAQ Reagent 8-plex Kit according to the manufacturer's protocol. Mixtures of peptides were dried by vacuum centrifugation and then fractionated using strong cationic exchange (SCX) chromatography.
Separation SCX chromatography and LC-MS/MS analysis The peptides were fractionated using a Shimadzu LC-20AB HPLC Pump System. Peptide mixtures were reconstituted with buffer A. The peptides were eluted at a ow rate of 1 mL/min with a gradient of 5% buffer B for 10 min, 5%-35% buffer B for 40 min, and 35%-95% for 1 min. The system was then maintained in 95% buffer B for 3 min before equilibrating with 5% for 10 min. Elution was monitored by measuring the absorbance at 214 nm, and peptides were pooled into 20 fractions, desalted by Strata X column, and vacuum dried. Each fraction was resuspended in buffer A and centrifuged. Then, 5 µL of supernatant was loaded onto a Shimadzu LC-20 AD nano HPLC coupled with a C18 trap column, and peptides were eluted onto an analytical C18 column that was packed in-house. Then, the samples were analyzed using 8%-35% buffer B, which was increased to 60% and 80%, and nally returned to 5% in 1 min and equilibrated for 10 min. The peptides were subjected to tandem mass spectrometry using Q Exactive mass spectrometer (Thermo Fisher Scienti c, San Jose, CA) by nano-electrospray ionization.
Parameters for MS analysis were as follows: electrospray voltage: 1.6 kV; precursor scan range: 350-1600 m/z at a resolution of 70,000 in Orbitrap; MS/MS fragment scan range: >100 m/z at a resolution of 17,500 in HCD mode; normalized collision energy setting: 27%; dynamic exclusion time: 15 s; automatic gain control for full MS target and MS2 target: 3e6 and 1e5, respectively; number of MS/MS scans following one MS scan: 20; and most abundant precursor ions above a threshold ion count of 20,000.

Protein identi cation and quanti cation
Raw MS/MS data were converted into MGF le format with Proteome Discoverer 2.2 (Thermo Fisher Scienti c, Waltham, MA, USA). Protein identi cation was performed using Mascot software (Matrix Science, London, U.K.; version 2.3.02) against a combined database containing E. asinus, Bos taurus, and Ovis aries biotypes. Relative quanti cation of identi ed proteins was performed using IQuant, which provided reliable signi cance measures [31]. All proteins with FDR less than 1% proceeded with downstream analysis.
De nition of differentially expressed proteins (DEPs) Differential expression ratios for proteins were determined using Student's t-test (p value £ 0.05), and proteins with a 1.2-fold change or greater were considered to be differentially expressed. Proteins usually interact with each other to play roles in certain biological functions. We performed pathway enrichment analysis of DEPs based on KEGG database. Proteomics data was deposited to the ProteomeXchange Consortium via PRIDE [32] partner repository.
qRT-PCR analysis and validation of target proteins Total RNA was extracted according to the manufacturer's specifications. RNA yield was determined using a NanoDrop 2000 spectrophotometer (Thermo Scienti c, USA), and the integrity was evaluated using agarose gel electrophoresis stained with ethidium bromide. Total RNAs were reverse transcribed for cDNA synthesis with random primers and RT-PCR kit according to the manufacturer's protocol (Promega, Madison, WI). Quantification was performed via a two-step reaction. GAPDH was selected as internal reference. Relative expression data were analyzed using the 2 −ΔΔct method [33]. MRM is a type of mass spectrometry technology that can detect the reliability of iTRAQ result using QTRAP 5500 mass spectrometer (SCIEX, Framingham, MA, USA). MRM method was successfully established only if the protein had at least one unique peptide and met the criteria in a previous study [34].

Results
Overview of transcriptome data and miRNA-seq data After quality control and data ltering, 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 le 1 .
Small RNA libraries were constructed and sequenced, which revealed the characteristics of miRNAs in skeletal muscle of donkey. After quality ltering 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 identi ed by mapping the nal reads to miRBase. Results showed that a total of 312 known miRNAs were expressed in donkey skeletal muscles.
Meanwhile, 181 novel miRNAs were identi ed (Additional le 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 le 3). Energy metabolism, sarcomere organization, muscle contraction, and skeletal muscle ber development were involved in biological process terms. In addition to making up the main organelle components, many genes participated in the basic blocks of myo brils (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 le 4) showed that ribosome gene pathway and oxidative phosphorylation signaling pathway were substantially enriched.
To make gene expression among the three species comparable, we rst identi ed 5160 putative singlecopy orthologous genes in cow, donkey, and goat. A total of 1338 and 1780 DEGs were identi ed and further visualized in a volcanic map in cow vs donkey and goat vs donkey ( Fig. 1 and Additional le 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 signi cant GO terms (p < 0.05) were annotated, including 101 biological processes, 36 cellular components, and 59 molecular functions (Additional le 7). The signi cant GO annotation terms between goat and donkey muscles are presented in Additional le 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, ve 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 le 9 and le 10. Most DEGs were signi cantly 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 identi ed to encode for 39 TF families (Additional le 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 proteinprotein interaction networks, and detailed distinct clusters are shown in Additional le 12, le 13 and le 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 ve 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 signi cant AS events. As shown in Fig. 2 and Additional le 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 lament binding (GO: 0051015), and striated muscle thin lament (GO: 0005865). Using sequence reads aligned to the donkey reference genome, we identi ed a total of 57,201 putative SNPs. SNPs were then classi ed 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 (  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 identi ed 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 identi ed with FDR of less than 1%. The annotation information of all proteins is shown in Additional le 17.

Identi cation 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% con dence and a 1.2-fold change cutoff. In goat vs donkey muscle, we identi ed 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 le 18 and le 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 classi cations (p < 0.05) in the cow vs donkey group are presented in Additional le 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 signi cant GO terms were assigned to DE proteins, which were classi ed into independent categories: 57 terms corresponding to biological processes, 20 terms corresponding to molecular functions, and 30 terms corresponding to cellular components (Additional le 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 le 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 signi cant differences at the mRNA and protein levels between the two comparison groups (Additional le 24 and le 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 downregulated, 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 signi cant KEGG signaling pathways included oxidative phosphorylation, metabolic pathways, fatty acid degradation, arginine, proline, and cholesterol metabolism (Additional le 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 le 27. Five genes (Tnni2, TMOD1, CFL2, MYOZ3, and MYL12B) were associated with muscle ber synthesis and differed at transcriptional and proteome levels. CFL2 (Co lin-2), an actin-binding protein, is the only speci c isoform that is expressed in mature skeletal muscles [36]. This protein was localized near the M-bands in myo brils. It can enhance actin lament conversion and regulate the length of actin laments, 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 signi cantly 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. Signi cant 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 le 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 con rmation of DEGs and MRM validation of DEPs
To validate the expression pro les 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 signi cantly correlated with a Pearson correlation coe cient 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 con rmed that the selected DEGs were differentially expressed in three muscles, demonstrating that the highthroughput sequencing was accurate. The primers used for these genes are shown in Additional le 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 le 29). The corresponding correlation coe cient r was 0.75 and 0.78 for protein expression using MRM and iTRAQ methods, respectively ( Fig. 9 and Additional le 30). The protein expression levels in MRM were basically consistent with the iTRAQ expression pattern. Therefore, MRM analysis con rmed that the proteomic analysis was reliable.

Discussion
Donkey meat consumption is marginal when compared with other conventional meat such as pork, beef, mutton or chicken. In China, the price of raw donkey meat is comparable to beef and mutton. If it's processed, the price will be twice as much as before, and the cooking percentage is about 55%-60% after testing. Donkey meat has recently been recognized as a nutritive food for human consumption with good quality proteins [35]. In the past, donkey meat was tough and not very acceptable because donkey was slaughtered at the end of working live. Nowadays, consumers prefer lean meat with less fat and balanced nutrition. So donkey meat is mainly produced from young asses to avoid lack of tenderness. In the study of meat quality of horse meat, factors such as Equine breed, sex, and slaughter age have been shown to affect carcass and meat quality. The nutritional composition of horse meat is characterized by low levels of fat and cholesterol, and a high protein content by comparison with pork, beef or poultry [35,36]. Ropka-Molik et al. [37] investigated the differently expressed genes' patterns and molecular genetic mechanisms of skeletal muscle in two pig breeds using transcriptome data. In this study, we rstly analyze the highly expressed genes in donkey muscle, functional annotation discovered 31 expressed genes were involved in muscle contraction and skeletal muscle ber development, including 5 Troponin genes (TNNC1,   TNNC2 (Tm) and actin laments [38].
More DE genes were detected in two comparison groups. Compared to goat meat, there were ve DEGs annotated in GO term of skeletal muscle contraction in the donkey muscle and three genes (JSRP1, TNNI2, STAC3) were up-regulated. In another compared group, STAC3 gene was also highly expressed in the donkey muscle. STAC3 mRNA was predominantly expressed in skeletal muscle, and regulated myogenic differentiation and myo brillar assembly [39]. STAC3 knockout mouse was embryonically lethal and the skeletal muscles showed signi cant abnormalities [40,41]. In the analysis of protein-protein interactions for DEGs, PCNA, MRPL2, and NADH dehydrogenase compounds were the key nodes. Among them, PCNA is a protein participated in DNA replication and expressed at a high level in proliferating cells [42]. MRPL2 and NADH were involved in protein synthesis and energy metabolism. In the alternative splicing analysis, NEB gene was detected in the differential AS events, this is a large gene, encompassing 183 exons, which coexists with thick and thin laments within the sarcomeres of skeletal muscles. On the other hand, Different AS events of TMOD1, TPM2, and TMOD4 were detected in two comparison groups. The sarcomeric tropomodulin (TMOD) isoforms TMOD1 and TMOD2 cap thin lament pointed ends and functionally interact with the leiomodin isoforms to control myo bril organization, thin lament lengths, and actomyosin crossbridge formation in skeletal muscle bers.
The content of polyunsaturated fatty acids (PUFAs) was found 25.16g/100g total fatty acids in the Longissimus [9]. In beef and mutton, the content of PUFAs was only 5-10%. In differential protein analysis, we found that FADS2 protein was highly expressed in donkey meat. FADS2 gene encoded Δ-6 fatty acid desaturases which are rate-limiting enzyme involved in metabolic processes of fatty acids [43]. Zhu et al. [44] reported that two SNPs in the transcriptional regulatory region of FADS2 were associated with fatty acid content. In another study, FADS2 and FADS1 genes showed a greater expression in slow growing chickens and consequently a higher PUFA content in the breast meat than other two kinds chickens [45]. Zhang et al. [46] conducted GWAS for fatty acid metabolic traits in pig populations and found that one signi cant SNP located in the fth intron of the FADS2 gene.
In donkey meat, the content of lysine and glutamate were signi cantly higher than other meats. In differential protein analysis of two compared group, we found the differential proteins ACO2 was upregulated expression in donkey meat in the lysine synthesis pathway. ACO2 gene encodes an essential metabolic enzyme mitochondrial aconitase. The functional polymorphisms of this gene might be associated with genetic diseases (Idiopathic Parkinson's Disease, Obesity). The function of this protein in muscle has been little studied. The differential proteins IDH1, GSR, and PGD were also high expressed in donkey meat and Lap3 protein was down-regulated in Glutamate metabolism pathway. The function of these differentially expressed proteins in muscle needs to be further validated.

Conclusions
Using a combination of transcriptomic and proteomic analyses, we analyzed the molecular differences of expressed genes and protein levels in donkey muscle, beef, and mutton. Meanwhile, AS events and SNP variations of the expressed genes in donkey muscle were examined. Regulatory factors including miRNAs and transcription factors were identi ed. The above information was integrated as a database to provide more valuable molecular information for further analysis of donkey muscle characteristics. In addition to key genes related to metabolism, signi cant differences were observed in genes associated with muscle ber composition, indicating the characteristics of donkey muscle ber and indirectly re ecting the distinct avor of donkey muscle.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.