Animal experiments were carried out in compliance with animal care protocols and all efforts were made to minimize suffering. The protocol was approved by the Animal Care Committee of South China Agricultural University (Guangzhou, China) with approval number SCAU#0014.

Xinghua (XH) chickens at E10 were obtained from the Chicken Breeding Farm of South China Agricultural University (Guangzhou, China) and incubated in Automatic Incubator (Oscilla, Shandong, China) at 37.8°C, 50–70% relative humidity. Chicken embryo sexes were identified by PCR amplification of the CHD1 (chromodomain helicase DNA binding protein 1) gene (Fridolfsson and Ellegren, 1999). Leg muscles of female Xinghua chickens in three different development stages (E11, E16, and D1) were used for iTRAQ analysis.

Two female chickens of each stage E11, E16, and D1 were used for iTRAQ assays (Applied Biosystems, Foster city, CA, USA). Total proteins were extracted by using a urea lysis buffer (7 M urea, 2 M thiourea, and 1% SDS) containing 1 mM PMSF. Proteins were quantified using a BCA Assay Kit (Pierce, Thermo, USA) and detected by SDS-electrophoresis. Total proteins were treated follow reduction, cysteine alkylation, and trypsin digestion to obtain peptides, and then equal amounts of peptides from six samples were labeled individually with different iTRAQ reagents (E11, labeled by 113 and 114; E16, labeled by 115 and 116; and D1, labeled by 119 and 121) using instructions provided by the manufacturer. After iTRAQ-labeling, peptides were desalted using a C18 solid-phase extraction and submitted for Nano Liquid Chromatography–Mass Spectrometry/Mass Spectrometry (LC-MS/MS) analysis.

The experiments were performed on a Nano Aquity UPLC system (Waters Corporation, Milford, MA) connected to a quadrupole-Orbitrap mass Spectrometer (Q-Exactive) (Thermo Fisher Scientific, Bremen, Germany) equipped with an online nano-electrospray ion source. The Q-Exactive mass spectrometer was operated in the data-dependent mode to switch automatically between MS and MS/MS acquisition. Survey full-scan MS spectra (m/z 350–1,200) were acquired with a mass resolution of 70K, followed by fifteen sequential high energy collisional dissociation MS/MS scans with a resolution of 17.5K. In all cases, one microscan was recorded using dynamic exclusion of 60 s. MS/MS fixed first mass was set at 100.

The MS/MS data were analyzed with Proteome Discoverer software v1.4 (Thermo Scientific), and search in the Uniprot database (Gallus gallus). The target-decoy based strategy was applied to control peptide level false discovery rates (FDR) lower than 1%. Only unique peptides were used for protein quantification and normalization on protein medians was used to correct experimental bias, the minimum number of proteins that must be observed to allow was set to 200. For protein quantitation, the ratios of each sample were weighted and normalized by comparing the control group (sample tagged as 113) as the denominator. For quantitative changes, we set a >1.5 or <0.66-fold change cutoff and p-value (t-test) <0.05 for differentially expressed proteins.

All identified proteins were annotated and classified by GO (Gene Ontology, http://www.geneontology.org/) and KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/) pathway. The differentially expressed proteins were further processed by DAVID 6.8 Functional Annotation Tool (http://david.abcc.ncifcrf.gov/) for term enrichment analysis (Huang et al., 2009a). The results were filtered based on a Fisher Exact statistic methodology as previously described (Huang et al., 2009b). The GO biological network was assessed using the ClueGO of Cytoscape software (http://www.cytoscape.org/). Protein-protein interaction analysis was performance by String v10.0 (http://www.string-db.org/), and a high coefficient value of 0.7 was used as a cutoff (Szklarczyk et al., 2015). Cluster analysis was performed to identify the expression patterns of differentially expressed proteins (fold change ≥ 1.2 or ≤ 0.8 and p < 0.05) using hcluster (https://pypi.python.org/pypi/hcluster/0.2.0).

Transcriptomes data (including mRNAs and lncRNAs) obtained from our previous study were using for integrative analysis with the proteome data (Li et al., 2017). Differentially expressed proteins were compared with the differentially expressed mRNAs, to identify proteins that were consistently expressed at the RNA and protein levels. The correlation between these consistently expressed proteins was analyzed by String v10.0 with a coefficient value of 0.7. The differentially expressed proteins were also compared with the predicted potential target genes of differentially expressed lncRNAs. The interaction network of the corresponding proteins and lncRNAs were constructed by String v10.0 and Cytoscape software.

Proteins were extracted from muscle tissues by using a urea lysis buffer (7 M urea, 2 M thiourea and 1% SDS) containing 1 mM PMSF. Total proteins (50 μg) were separated on a 12% SDS-PAGE, and transferred to a polyvinylidene difluoride membrane (Millipore, Bedford, MA). The membrane was blocked in 5% BSA blocking solution for 1 h at room temperature and incubated overnight at 4°C with primary antibodies (Abcam, Cambridge, UK) as follows: MYL1 (ab97427, diluted 1:500), MYL3 (ab137767, diluted 1:1,000), RPL4 (ab154907, diluted 1:1,000), STMN1(ab194670, diluted 1:1,000), RPS3A (ab101690, diluted 1:500), and TNNT3 (ab82784, diluted 1:500). After that, the membrane was washed with PBS-T and then developed with anti-mouse or rabbit horseradish peroxidase conjugated secondary antibodies (Sigma-Aldrich, diluted 1:5,000). Protein bands were visualized using enhanced chemiluminescence (ECL) system (GE Healthcare, USA) and quantified with an ImageQuant LAS4000 system (Fujifilm, Tokyo, Japan).

Embryonic growth and development of skeletal muscle has a significant effect on muscle mass in chickens (Halevy et al., 2004, 2006). To assess the protein expression profiles during embryonic skeletal muscle development, we used the same samples as our previous RNA sequencing (leg muscle tissues of female XH chicken at E11, E16, and D1; each stage two individuals) for iTRAQ-based proteomics analysis (Li et al., 2017).

The LC-MS/MS data analysis generated a total of 58,815 match spectra, 19,659 peptides, and 15,495 unique peptides. We could identify 3,240 proteins possessing at least one unique peptide with a confidence level above 95% (Table S1). All proteins were grouped according to biological process, cellular component, and molecular function by GO analysis (Table S2). The cellular process, metabolic process and biological regulation were mainly categories for these proteins in biological processes, and the binding and catalytic activity were the two most abundant categories in molecular function (Figure 1A). The KEGG analysis indicated that these proteins were primarily involved in pathway of signal transduction, translation, and transport and catabolism (Figure 1B).

To access proteins dynamic changes during embryonic development of skeletal muscle, we identified 491 differentially expressed proteins (fold change ≥ 1.5 or ≤ 0.666 and p < 0.05) from the three different developmental groups (Figure 2A). There were 19 up- and 32 down-regulated proteins in E11 vs. E16 (Table S3) group, 238 up- and 227 down-regulated proteins in E11 vs. D1 (Table S4) group, and 13 up- and 5 downregulated proteins in E16 vs. D1 (Table S5) group (Figures 2B,C). According to the p-values, the top 10 differentially expressed proteins of each comparison group were listed in Table 1. These differentially expressed proteins could play an important role on development of embryonic muscle, and served as candidate genes for further study.

In order to gain insights into the functions of these differentially expressed proteins, GO functional enrichment analysis was performed using DAVID tool. In total, 420 proteins were enriched into130 GO term. These were involved in epithelial cell differentiation, collagen fibril organization, muscle contraction, cytoplasmic translation and protein folding (Figure 2D). We also found significant enrichments for phagosome, ribosome, pyruvate metabolism, and Glycolysis/Gluconeogenesis by KEGG pathway analysis (Figure S1). We further analyzed 57 proteins taken from the top 10 biological processes and formed a functional interaction network. This analysis indicated that these significantly altered proteins were involved in structural constituent of ribosome as well as cardiac and striated muscle contraction (Figure 2E).

To characterize the expression patterns of differentially expressed proteins (fold change ≥ 1.2 or ≤ 0.8 and p < 0.05), cluster analysis was performed using hcluster. The heatmap of these differentially expressed proteins indicated that their expression patterns are similar in E16 and D1 (Figure S2). According to the expression pattern, these differentially expressed genes could be mainly grouped into four distinct clusters (Figure 3). The largest cluster (cluster 1, including 476 proteins) had a consistently overall pattern of down-regulation from E11 to D1 (Figure 3A). However, t cluster 2 (358 proteins) were up-regulated steadily from E11 to D1 and the 52 proteins of cluster 3 were up-regulated sharply from E11 to D1. The nine proteins of cluster 4 were up-regulated sharply from E11 to E16 but showed no significant change at D1 (Figure 3D).

To confirm the iTRAQ results, we performed western blotting analysis for several candidate proteins, including MYL1 (myosin light chain 1), MYL3 (myosin light chain 3), RPS3A (ribosomal protein S3A), STMN1 (stathmin 1) and TNNT3 (troponin T3) (Figure 4A). The no obvious changes protein of RPL4 (ribosomal protein L4) was used as reference gene. The expression levels of MYL1, MYL3, and TNNT3 all increased from E11 to D1, while the RPS3A and STMN1 decreased during this same period (Figure 4B). The expression patterns of these proteins were consistent with our iTRAQ data.

Proteins are usually not work alone, but rather interacted with each other to perform various functions. To explore the protein interaction networks altered in development of chicken embryonic skeletal muscle, differentially expressed proteins identified in this study were analyzed using STRING software. In E11 vs. D1 group, 396 of 465 differentially expressed proteins were detected in STRING, and 77 proteins were connected to nodes in the network (Figure 5A). Interestingly, we found interaction networks for ribosomal protein (RPS20, RPS25, RPLP1, and RPL9 etc.), muscle contraction (TNNC1, TNNC2, TNNI1, TNNI2, TPM2, TPM3, MYL1, MYL2, MYL3, and DMD), pyruvate metabolism (PDHA1, PDK3, and ACAC) and oxidative phosphorylation (NADH: ubiquinone oxidoreductase family, including NDUFA5, NDUFS6, and NDUFB9 etc.). In E11 vs. E16 group, 43 differentially expressed proteins were analyzed and 22 of them constituted an interaction network. This network was related to pathway of GTP and myosin binding and included ADSS, ARF4, PSMC1, ACTA1, and ATP2A2 (Figure 5B). In E16 vs. D1 group, 16 of 18 differentially expressed proteins were detected and eight of them constituted a network including NDUFA5, CYCS, PRDX3, SMYD1, and PCNA (Figure 5C).

Skeletal muscle development in the chicken is a complex physiological process, and gene expression is regulated at multiple levels. Thus, we performed a proteome and transcriptomes integrative analysis to give us clues as to the levels of regulation for the differentially expressed proteins. When we compared our differentially expressed proteins with the differentially expressed mRNAs, 189 differentially expressed proteins were regulated in the same manner as their mRNAs (Table S6). The interaction between these proteins was also involved in the pathway of muscle contraction and oxidative phosphorylation (Figure 6).

LncRNAs can regulate gene expression through cis-acting or trans-acting mechanisms (Bu et al., 2012; Han et al., 2012). We next integrated our lncRNA data with the proteomic data to construct an interaction network. We found that only a small fraction of differentially expressed proteins were from the target gene of corresponding differentially expressed lncRNAs. In E11 vs. E16 group and E16 vs. D1 group, only the MARCKS gene was identified as the target gene of lnc00057929, and KLHL15 gene was found as target gene of lnc00037615 and lnc00066198. For the E11 vs. D1 group, there were 21 differentially expressed proteins we found as target genes of 22 lncRNAs (Figure 7 and Table S7). Interestingly, four proteins (DMD, MYL3, TNNI2, and TNNT3) involved in muscle contraction were identified in this network, indicating that they may regulated by their corresponding lncRNAs (Figure 7).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.