Proteomic analyses of sheep (ovis aries) embryonic skeletal muscle

The growth and development of embryonic skeletal muscle plays a crucial role in sheep muscle mass. But proteomic analyses for embryonic skeletal development in sheep had been little involved in the past research. In this study, we explored differential abundance proteins during embryonic skeletal muscle development by the tandem mass tags (TMT) and performed a protein profile analyses in the longissimus dorsi of Chinese merino sheep at embryonic ages Day85 (D85N), Day105 (D105N) and Day135 (D135N). 5,520 proteins in sheep embryonic skeletal muscle were identified, and 1,316 of them were differential abundance (fold change ≥1.5 and p-value < 0.05). After the KEGG enrichment analyses, these differential abundance proteins were significant enriched in the protein binding, muscle contraction and energy metabolism pathways. After validation of the protein quantification with the parallel reaction monitoring (PRM), 41% (16/39) significant abundance proteins were validated, which was similar to the results of protein quantification with TMT. All results indicated that D85N to D105N was the stage of embryonic muscle fibers proliferation, while D105N to D135N was the stage of their hypertrophy. These findings provided a deeper understanding of the function and rules of proteins in different phases of sheep embryonic skeletal muscle growth and development.

). Results of comparable analyses showed that, in differential abundance proteins, 33 proteins were up-regulated and 21 were down-regulated in D105N vs. D85N, 272 proteins up-regulated and 220 down-regulated in D135N vs. D105N, 419 proteins up-regulated and 351 down-regulated in D135N vs. D85N. The subcellular localization results of differential abundance proteins were significant enriched in cytoplasm (Fig. 2). Function analyses of differential abundance proteins. After Go analyses, in three comparable groups, the differential abundance proteins were mainly involved in metabolic and cellular biological processes, binding and catalytic activation molecular functions in cells (Fig. 3). To deeply analyze the differential abundance proteins, the functional classification of clusters of orthologous group protein (COG) (Fig. 4) and KEGG pathways enrichment (Fig. 5) were used to regard the function and characteristic. Interestingly, the results of COG functional classification in differential abundance proteins showed that the cytoskeletal proteins were mainly expressed in D105N vs. D85N. And they played an important role in regulating the contraction and maintaining the morphology of skeletal muscle cells at the early development phase 17 . While, the energy production and conversion proteins were significant expressed in D135N vs. D105N, the signal transduction and the energy production and conversion were expressed in D135N vs. D85N. And all results of function analyses indicated that the rules of the differential abundance proteins were consistent with the growth and development of skeletal muscle. And all analysing results indicated that the number of embryonic skeletal muscle fibres were increased mainly in D85N to D105N. But in D105N to D135N, the embryonic skeletal muscle fibres were hypertrophic with a higher energy metabolism. While the regulatory proteins of muscle contraction were up-expressed in D135N vs. D85N. PRM validation of TMT data. The results of the validation protein fragment ion peptides indicated that 16 proteins were validated and 11 of them were quantified with two unique peptides (Table 1). After analyses, the results demonstrated that the quantitation ratios of TMT and those of PRM were consistent each other (Table 2). So the reliability of the TMT quantitation data and the advantages of the PRM technology were notable. (the detail validation of the PRM as Supplementary Table S3)

Discussion
The embryonic phase is important for skeletal muscle development and growth 18 . During this phase, the muscle fibres are formed and proliferated at early embryonic stage and hypertrophic at late stage 19,20 , which has great influence on production of meat mass 21 . In this paper, the protein profile of the embryonic skeletal muscle in sheep was established by using the TMT 22,23 and validated by PRM technology 24,25 . TMT is widely implemented in the higher throughput proteomics, and PRM can be used to quantify protein abundance without specific antibodies and detect specific predetermined with known fragmentation characteristics more accurately in complex backgrounds than the traditional immunoassays [26][27][28][29][30] . And the results of TMT quantitation were entirely validated by PRM (Table 1). Finally, a total number of 5,520 proteins were identified and 1,316 of them were differential abundance proteins, which will help us to unveil the mechanisms of embryonic skeletal muscle development and growth of sheep.
After deep analyses of these data, we found that the cytoskeletal proteins and the energy production and conversion proteins were mainly expressed in D105N vs. D85N and D135N vs. D105N respectively (Fig. 4). And the results of GO and KEGG analyses revealed that these proteins were enriched in the protein binding, muscle contraction and energy metabolism pathways which took part in the regulation in the process of embryonic skeletal muscle development and growth of sheep [31][32][33][34][35] . And the cytoskeletal proteins with higher abundance in stage of muscle fibers proliferation and the energy metabolism played a major role in the embryonic muscle hypertrophy 13,36 , this indicated that D85N to D105N was the stage of embryonic muscle fibers proliferation and D105N to D135N was hypertrophy.
The metabolic pathway was very significant in the process of embryonic skeletal muscle development. And the metabolic pathway decided the muscle fibers proliferation and hypertrophy and regulated the protein differentiation and cell cycle 37 . Meanwhile, a few of reports showed that energy metabolism significantly regulated the myoblast to be in hypertrophy, but not in the cycle of proliferation withdrew from the cycle of cell proliferation and turned into muscle hypertrophy by precisely regulating the protein activity of myoblast cycle. The myoblast cycle was regulated by changing the phosphorylation levels of their cycle regulation proteins under the different metabolic condition 38,39 .
In both of our study (Fig. 5) and other reports, all results showed that the metabolic and oxidative phosphorylation pathways were significantly associated with the development of muscle [40][41][42] . In a word, the metabolism provided not only energy in whole embryonic skeletal muscle development but also fine modulation of embryonic muscle development styles [43][44][45] . As to the protein profile constructed in this study, it will be certainly an important basis for deep elucidation of the mechanisms of sheep embryonic skeletal muscle development and growth.

Methods
Experimental design and samples preparation. Embryonic skeletal muscle samples were taken from different development phases of foetal (the detail can be found in Supplementary Information S1). And on the basis of reports before 46 , the pregnant of female sheep were selected on D85N, D105N and D135N for caesarean section. The foetuses with similar embryonic ages were taken out and the longissimus dorsi was collected and there were three biological repeats per time phases (no mixing between per time phase samples). Furthermore, three comparable groups such as the D105N vs. D85N, D135N vs. D105N and D135N vs. D85N were set to study. Finally, these samples were immediately washed clear with PBS (pH7.2), then frozen in liquid nitrogen. All of samples were kept at −80 °C until further analyses.
Animal welfare disclaimer. The works herein described was conducted according to relevant animal care protocols, which were approved by the relevant guidelines and regulations of the Ministry of Agriculture of the People's Republic of China (Beijing, China). Furthermore, in this study, all use and handing of experimental animals were made to minimize suffering. Protein extraction and TMT labeling. Samples used for extraction of total proteins in each stage D85N, D105N and D135N were prepared according to the protocol as follows. The samples were taken out at −80 °C, and an appropriate amount of tissue sample was weighed into a liquid nitrogen pre-cooled mortar, and liquid nitrogen was sufficiently ground to a powder. Each group of samples were separately added to 4 volumes of lysis buffer (8 mol urea, 1% protease inhibitor and 2 mmol EDTA) and sonicated. After centrifugation at 12,000 g for 10 min at 4 °C, the cell debris were removed, the supernatant was transferred to a new centrifuge tube, and the www.nature.com/scientificreports www.nature.com/scientificreports/ protein concentration was determined using the BCA kit (Beyotime Biotechnology, China). Post addition of the dithiothreitol into the protein solution to a final concentration of 5 mmol, digestion was performed at 56 °C for 30 min. Then the iodoacetamide was added to a final concentration of 11 mmol and incubated for 15 min at room temperature in dark condition. Finally, the urea concentration in sample was diluted to less than 2 mol. Trypsin (1 μg/μl) was added at a mass ratio of 1:50 (pancreatin: protein sample, V/V) and digested overnight at 37 °C. Trypsin was added at a mass ratio of 1:100 (pancreatin: protein sample, V/V) and continued to digest for 4 h.
After trypsin digestion, the peptides were desalted by Strata X C18 SPE column (Phenomenex) and vacuum-dried, all stage groups were reconstituted in 0.5 mol TEAB and labelled with three labels according to the manufacturer's direction for TMT kit (Thermo, USA). Briefly, one unit of the TMT reagent was thawed and reconstituted in acetonitrile. Then, the peptides mixture with an equal molar ratio were incubated at room temperature for 2 h. These were combined and desalted. Finally, the TMT labelling peptide segments were completed by vacuum centrifugal drying.

LC-MS/MS analyses.
The peptides in each fractions were dissolved in liquid phase A phase and analysed by using the EASY-nLC 1200 ultra-high performance liquid system interfaced with a reversed-phase analytical column (15 cm length, 75 μmol ID). The liquid phase A is an aqueous solution containing 0.1% formic acid and 2% acetonitrile. While the liquid phase B is an aqueous solution containing 0.1% formic acid and 90% acetonitrile. The liquid phase gradient settings were in 8-16% B for about 30 min, 16-30% B for 15 min, 30-80% B for 2 min and 80% B for 3 min step by step. The flow rate maintained at 400 nl/min. Then, the peptides were separated by an ultra-high performance liquid phase system and they were injected into an NSI ion source for ionization and analysed by Orbitrap Lumos ™ mass spectrometry. The ion source voltage was set to 2.0 kV and the peptide precursor and its secondary fragments were detected and analysed using the high-resolution Orbitrap.
The primary mass spectrometer scan range was set to 350~1,550 m/z, and the scan resolution was set to 60,000. The secondary mass spectrometry scan range had a fixed starting point of 100 m/z, and the secondary scan resolution was set to 15,000. In the end, the data acquisition mode used a data-dependent scanning (DDA) program. After the first-stage scanning, the first 20 peptides with highest signal intensity were selected to enter the high-energy collision-induced dissociation (HCD). Subsequently, 32% of the energy was used for fragmentation and also performed in peptides. In order to improve the effective utilization of mass spectrometry and to avoid parent ions. In this analyses, the automatic gain control (AGC) was set to 50,000, the signal threshold 50,000 ions/s, the maximum injection time 70 ms and the dynamic exclusion time of tandem mass spectrometry 30 s. Data analyses. All proteins in 9 fractions were analyzed and compared with Maxquant engine v.1.5.2.8 (http://www.maxquant.org/) and tandem mass spectra were searched against the NCBI Ovis aries Oar_v4.0 (https://www.ncbi.nlm.nih.gov/genome/? term = Ovis + aries 43,036 sequences) database. Meanwhile, in order to eliminate the influence of contaminating proteins in the identification results, the reverse decoy database to calculate the false discovery rate (FDR) caused by random matching and common pollution database were concatenated. And the Trypsin/P was specified as cleavage enzyme allowing up to 2 missing cleavages. The mass tolerance for precursor ions was set as 20ppm in the first search and 5ppm in the main search. But the mass tolerance for fragment ions was set as 0.02 Da. And the quantitative method was set as TMT-6plex, the FDR for protein identification was set as 1%. The differential abundance proteins were calculated by the mean of the three repeats to another mean of the three repeats. Then T-tests between two-samples were used to compare differential www.nature.com/scientificreports www.nature.com/scientificreports/ abundance proteins. A significance level (P-value) of 0.05 was used for statistical analyses any time. Post the statistical analyses, the fold change 1.2, 1.3, 1.5 and 2 were set for the quantitation proteins. Based on the standard which the ratio of quantitative protein to differential abundance protein is less than 30% (Table 3), the fold change ≥ 1.5 and p-value < 0.05 were considered to be differential abundance proteins. For TMT quantification, the ratios of the TMT reporter ion intensities in MS/MS spectra (126-131 m/z) from raw data sets were used to Bioinformatics analyses. All identified proteins were annotated and classified according to the protein sequence-based algorithm software InterProScan v.5.14-53.0 (http://www.ebi.ac.uk/interpro/), which predicted the GO function of proteins and classifies them based on cellular composition, molecular function and biological process. Then protein pathways were clustered by using the KEGG pathway database to annotate. At first, the submitted proteins were annotated using the KEGG online service tool KAAS v.2.0 (http://www.genome.jp/kaas-bin/ kaas_main), then the annotated proteins were matched to the corresponding pathway in the database by the KEGG mapper V2.5 (http://www.kegg.jp/kegg/mapper.html). Furthermore, the Fischer's exact double-end test (Fisher's exact test) was used to probe differential abundance proteins as a background for enrichment analyses by the Perl module v.1.31 (https://metacpan.org/pod/Text::NSP::Measures::2D::Fisher). The clustering relationships of the differential abundance proteins were visualized by using the Heatmap which can be drawn by the R package of heatmap2 and Gplots v.2.0.3 (https://cran.r-project.org/web/packages/cluster/). And the Wolfpsort v.0.2 (http://www.genscript.com/psort/wolf_psort.html) software was used to perform subcellular structural localization prediction and classification statistics for differential abundance proteins.
Quantitative validation based on targeted proteomics. The extraction, digestion and mass spectrometry of validation proteins were same as before description and each sample of 1.5 µg hydrolysed peptides was validated. The results of MS data were processed using Skyline v. 3.6. Indexes for Peptide were set as: Trypsin digestion [KR/P], the max missed cleavage 2 and the peptide length 8-25 aa that the max variable modification of carbamidomethyl on Cys and oxidation on Met was set to 3. For transition sets, precursor charges were set as 2 or 3, ion charges as 1 or 2 and ion types as b, y and p. But the ions of product were set as from ion 3 to the last ion, and the value of ion matching tolerance was set as 0.02 Da. The all data of analyses were visualized by Skyline 47 .
In our study, to validate the TMT quantitation results by using the PRM. 40 differential abundance proteins were randomly selected from three comparable groups (the detail validation of the PRM as Supplementary  Table S3).
*(The current research is based on a small amount of animals and may be considered as preliminary). Table 3. Differential fold change abundance proteins summary. Based on the standard which the ratio of quantitative protein to differential abundance protein is less than 30%, the fold change ≥ 1.5 and p-value < 0.05 were considered to be differential abundance proteins.