# Identification of Differentially Expressed Genes in Different Types of Broiler Skeletal Muscle Fibers Using the RNA-seq Technique

Authors: Han Wang; Zhonghao Shen; Xiaolong Zhou; Songbai Yang; Feifei Yan; Ke He; Ayong Zhao
Journal: BioMed Research International
 (2020)

---

## Abstract

The difference in muscle fiber types is very important to the muscle development and meat quality of broilers. At present, the molecular regulation mechanisms of skeletal muscle fiber-type transformation in broilers are still unclear. In this study, differentially expressed genes between breast and leg muscles in broilers were analyzed using RNA-seq. A total of 767 DEGs were identified. Compared with leg muscle, there were 429 upregulated genes and 338 downregulated genes in breast muscle. Gene Ontology (GO) enrichment indicated that these DEGs were mainly involved in cellular processes, single organism processes, cells, and cellular components, as well as binding and catalytic activity. KEGG analysis shows that a total of 230 DEGs were mapped to 126 KEGG pathways and significantly enriched in the four pathways of glycolysis/gluconeogenesis, starch and sucrose metabolism, insulin signalling pathways, and the biosynthesis of amino acids. Quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) was used to verify the differential expression of 7 selected DEGs, and the results were consistent with RNA-seq data. In addition, the expression profile ofMyHC isoforms in chicken skeletal muscle cells showed that with the extension of differentiation time, the expression of fast fiber subunits (types IIA and IIB) gradually increased, while slow muscle fiber subunits (type I) showed a downward trend after 4 days of differentiation. The differential genes screened in this study will provide some new ideas for further understanding the molecular mechanism of skeletal muscle fiber transformation in broilers.


---

## Body

## 1. Introduction

Chicken is widely welcomed by consumers because it is low in fat, is low in cholesterol, and has no religious restrictions. Therefore, chicken has become the second-largest consumed meat product after pork [1, 2]. In recent years, with the continuous improvement in people’s living standards, their requirements for meat flavor have also become stricter. Age, sex, heredity, environment, nutrition, and intramuscular fat content (IMF) are all important factors affecting chicken quality [3–6]. In addition, studies have shown that muscle fiber properties were an essential factor affecting the meat quality of chicken [7, 8]. Skeletal muscle is one of the most important components of meat-producing animals, accounting for approximately 45%~60% of the whole animal body [9]. According to the characteristics of contraction and metabolism, skeletal muscle fibers can be roughly divided into two types: slow-twitch (type I) and fast-twitch (type II) muscles. Type II muscle fibers can be categorized into a fast-twitch oxidation type (IIA), a fast-twitch glycolytic type (IIB), and a super-fast-twitch type (IIX) [10, 11]. Previous studies showed that muscle with a high proportion of fast-twitch muscle fibers appeared paler, while the muscle with more slow-twitch muscle fibers tended to be redder [8]. Therefore, skeletal muscle was usually divided into two kinds of muscle, red muscle and white muscle [12]. Moreover, it was found that muscles with a higher proportion of oxidized muscle fibers had better meat quality [13].Currently, it has become a research hotspot to study the molecular mechanism of muscle fiber-type transformation in order to improve meat quality. It has been shown thatPGC-1 can promote the transformation of muscle into type I fibers in the skeletal muscle [14]. Cofilin2b (CFL2b) can affect the transformation pattern of myosin heavy chain (MyHC) subunits in the muscle of piglets [15]. Furthermore, noncoding RNA was also proven to play a role in slow-twitch muscle fiber formation and skeletal myogenesis [16]. At present, preliminary progress has been made on the molecular mechanism of chicken muscle fiber-type transformation. One study showed that the expression of PGC-1α might be closely related to the slow muscle fiber content in chicken skeletal muscle [17]. In addition, miR-1611 was discovered to be highly expressed in slow-twitch muscle fibers and could drive the transformation from fast-twitch to slow-twitch muscle fibers in chicken [18]. However, the regulatory network of muscle development is complex, and studies on the formation and transformation of chicken skeletal muscle fiber types are still incomplete. Therefore, more relevant, important genes need to be further explored.At present, transcriptome sequencing technology has been widely used in biological research, medical diagnostics, and therapeutic studies [19]. Meanwhile, this technique also plays an important role in farm animal muscle development. Based on the transcriptome analysis of biceps between Small-tail Han sheep and Duper sheep, a total of 1300 differentially expressed genes were identified, which could help further elucidate the mechanisms of muscle development as they are affected by breed [20]. Forty-nine differentially expressed genes were screened from the longissimus dorsi of Nellore cattle among different grades of marbling by RNA-seq, in an effort to find strategies to select animals with greater marbling [21]. Moreover, another study used RNA-seq to identify genes related to pig muscle development and meat quality [22]. RNA-seq technology was also used in broiler chicken groups to uncover 68 genes that might participate in the tenderization process [23]. Thus, it can be seen that whole-genome sequencing technology exerts an immense influence on the discovery of new functional genes.To better explore the molecular mechanism of skeletal muscle fiber transformation in broilers, this study used RNA-seq technology to screen differentially expressed genes between breast and leg muscles of Ross 308 broiler chickens. Subsequently, through GO and KEGG pathway analysis, genes and pathways related to skeletal muscle fiber-type conversion were discovered. In addition, to further explore the growth characteristics of broiler muscle fibers, a skeletal muscle cell differentiation model of broilers was constructed in vitro. Subsequently, the expression ofMyHC subunits during the differentiation of chicken skeletal muscle satellite cells was studied. These results will provide a theoretical basis for the genetic improvement in muscle quality in broilers in the future.


## 2. Materials and Methods

### 2.1. Experimental Animals and Tissues

The chickens used in this study were Ross 308 broiler chickens from a farm in Fuyang District, Hangzhou, Zhejiang Province, China. Five 42-day-old Ross 308 broilers with similar body weights were selected for slaughter. The leg and breast muscles were collected and frozen at -80°C for RNA extraction.At the same time, we bought 15 fertilized eggs of the same breed and incubated them for 10 days to collect chicken skeletal muscle satellite cells. All animal procedures used in this study were approved by the Ethics Committee for Animal Experiments of Zhejiang A&F University and were performed in accordance with the Guidelines for Animal Experimentation of University (Hangzhou, China).


### 2.2. Total RNA Extraction, Library Construction, and Sequencing

Total RNA was extracted from the gastrocnemius of legs and the pectoralis major of breasts using the TRIzol method [24]. The quality inspection was carried out through the Sample Testing Center of Biomarker Technologies Corporation (Beijing, China) to ensure that the samples were up to the standard, and the test quality report of the RNA samples is presented in Table S1. All RNA samples with an OD 260/280 1.8-2.0 qualified for further analysis. Sequencing libraries were generated using a NEBNext Ultra™ RNA Library Prep Kit (NEB, USA). Briefly, mRNA was purified from 1 μg of total RNA per sample using poly-T oligoattached magnetic beads. Then, first-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The cDNA fragments of 240 bp were preferentially selected with an AMPure XP system (Beckman Coulter, USA) for cDNA adaptor ligating. At last, polymerase chain reaction (PCR) amplification was carried out and the products were purified (AMPure XP system). The library quality was assessed on the Agilent Bioanalyzer 2100 system. The libraries were sequenced by paired-end (PE 150 bp) sequencing on an Illumina HiSeq™ platform (Illumina, USA), and the raw reads were obtained. The leg muscle and breast muscle tissues of broilers were numbered PT1-PT5 and PX1-PX5, respectively, for subsequent analysis.


### 2.3. Sequencing Quality Assessment and DEG Screening

Clean data (clean reads) was obtained by removing reads containing adapters or poly-N; low-quality reads were also removed from the raw data. At the same time, Q20, Q30, GC-content, and sequence duplication levels of the clean data were calculated. Ultimately, all downstream analyses were based on high-quality clean data. The chicken genome sequence (GRCg6a) was downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/genome/?term=chicken). Then, Hisat2 software was used to map these clean reads to the chicken reference genome [25]. Transcript assembly was performed using StringTie software [26]. StringTie is an algorithm based on optimization theory, which uses comparison information to construct a multivariable clipping graph and traffic network to assemble and evaluate the expression of reads according to the maximum flow algorithm. Compared with other software such as Cufflinks, it can build a more complete transcript and better evaluate expression [26]. Quantification of gene expression levels was estimated by the fragments per kilobase of transcript per million fragments (FPKM) method [27]. Based on raw counts, the reads mapped to the reference genome were quantified by StringTie software to obtain FPKM [28]. Then, DESeq R package (1.39.0), a statistical program for determining differential expression in digital gene expression data based on a model of negative binomial distribution, was used to analyze the differential gene expression between samples with biological repetition [29]. Differentially expressed genes (DEGs) were identified using the following filter criteria: foldchange≥2 and FDR<0.05. Fold change indicates the ratio of expression quantity between two samples (groups). The false discovery rate (FDR) was obtained by correcting the significant difference P value.


### 2.4. Gene Ontology (GO) and KEGG Pathway Enrichment Analysis

Gene Ontology (GO) enrichment analysis of the DEGs was implemented by GOseq R packages (http://www.geneontology.org/) based on Wallenius noncentral hypergeometric distribution [30], which could adjust for gene length bias in DEGs. Pathway significant enrichment analysis was run based on the hypergeometric test in the KEGG database (http://www.genome.jp/kegg/) to uncover the most significantly enriched pathways among the DEGs. KOBAS software [31] was used to test the statistical enrichment of DEGs in KEGG pathways.


### 2.5. Quantitative Real-Time RT-PCR Verification

Utilizing the NCBI gene bank, primers of 7 randomly selected DEGs, three different isoforms ofMyHC genes (MyHC I, IIA, and IIB) and the house-keeping gene GAPDH were designed by Primer Premier 5 software (Table S1). Primers were synthesized by Hangzhou Youkang Biotechnology Co., Ltd. (Zhejiang, China). Total RNA from cells or tissues was used to synthesize cDNA by a 5X All-In-One MasterMix Kit (abm, China). The cDNA was used as a template for quantitative PCR by using the EvaGreen 2X qPCR MasterMix Kit (abm, China). There were three replicates for each sample. The relative expression level of each related gene was calculated using the 2-△△CT statistical analysis method.


### 2.6. Isolation and Culture of Chicken Skeletal Muscle Satellite Cells In Vitro

The 10-day fertilized eggs were disinfected with ethanol. The eggshell on the upper part of the air chamber was removed using tweezers on the egg receptacle. A new pair of tweezers was used to clip out the chicken embryo. The leg muscle of the chicken embryo was isolated and placed in a new petri dish. Muscle samples were washed 3 times in phosphate buffer solution (PBS) (Logan, USA) containing penicillin-streptomycin liquid (Solarbio, China), and skin, blood vessels, adipose tissue, and connective tissue were removed. The sterile muscle tissue was cut into meat paste and digested with 0.25% trypsin-EDTA (Thermo Fisher, China) at 37°C for 20 min. DMEM/F12 (Logan, USA) containing 10% fetal bovine serum (FBS, Thermo Fisher, China) was added to terminate the digestion. The suspension was filtered through a 70μm mesh sieve and centrifuged at 1200 r/min for 8 min at room temperature. As the supernatant was discarded, cells were resuspended with DMEM/F12 containing 15% FBS and cultured in a 5% CO2 incubator at 37°C. One hour later, the cell suspension was transferred into a new cell petri dish; this was repeated twice to enrich for muscle satellite cells and eliminate fibroblasts. Cell differentiation was induced by replacing 15% FBS with 2% equine serum (Logan, USA). Once the isolated skeletal muscle satellite cells reached 70%-80% fusion, the growth medium was replaced with differentiation medium and the cells were cultured for 2, 4, and 6 days. Cells at 70%-80% fusion indicated day 0.


### 2.7. Immunofluorescence Staining

The isolated skeletal muscle satellite cells at 70%-80% fusion and the cells that were induced to differentiate after transfection were used for immunofluorescence detection following the method of Luo et al. [32]. Briefly, cells grown in 12-well plates were washed three times with precooled PBS and fixed with 4% paraformaldehyde for 15 min. Thereafter, the cells were permeabilized with 0.25% Triton X-100 per well for 10 min and blocked at 4°C overnight. Afterwards, cells were incubated with 1 : 100 diluted primary anti-Desmin (Abcam, China), anti-Pax7 (Abcam, China), or anti-MyHC (Santa, USA) antibodies for 1 hour at room temperature. Then, 1 : 100 diluted fluorescent secondary antibody (Thermo Fisher, China) was incubated with the cells for 1 hour at room temperature. After being washed by PBS three times, DAPI (Invitrogen, USA) was added to the cells and then incubated for 15 min at room temperature to stain the cell nuclei. Lastly, samples were captured using a fluorescence microscope (OLYMPUS, Japan).


## 2.1. Experimental Animals and Tissues

The chickens used in this study were Ross 308 broiler chickens from a farm in Fuyang District, Hangzhou, Zhejiang Province, China. Five 42-day-old Ross 308 broilers with similar body weights were selected for slaughter. The leg and breast muscles were collected and frozen at -80°C for RNA extraction.At the same time, we bought 15 fertilized eggs of the same breed and incubated them for 10 days to collect chicken skeletal muscle satellite cells. All animal procedures used in this study were approved by the Ethics Committee for Animal Experiments of Zhejiang A&F University and were performed in accordance with the Guidelines for Animal Experimentation of University (Hangzhou, China).


## 2.2. Total RNA Extraction, Library Construction, and Sequencing

Total RNA was extracted from the gastrocnemius of legs and the pectoralis major of breasts using the TRIzol method [24]. The quality inspection was carried out through the Sample Testing Center of Biomarker Technologies Corporation (Beijing, China) to ensure that the samples were up to the standard, and the test quality report of the RNA samples is presented in Table S1. All RNA samples with an OD 260/280 1.8-2.0 qualified for further analysis. Sequencing libraries were generated using a NEBNext Ultra™ RNA Library Prep Kit (NEB, USA). Briefly, mRNA was purified from 1 μg of total RNA per sample using poly-T oligoattached magnetic beads. Then, first-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase. Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The cDNA fragments of 240 bp were preferentially selected with an AMPure XP system (Beckman Coulter, USA) for cDNA adaptor ligating. At last, polymerase chain reaction (PCR) amplification was carried out and the products were purified (AMPure XP system). The library quality was assessed on the Agilent Bioanalyzer 2100 system. The libraries were sequenced by paired-end (PE 150 bp) sequencing on an Illumina HiSeq™ platform (Illumina, USA), and the raw reads were obtained. The leg muscle and breast muscle tissues of broilers were numbered PT1-PT5 and PX1-PX5, respectively, for subsequent analysis.


## 2.3. Sequencing Quality Assessment and DEG Screening

Clean data (clean reads) was obtained by removing reads containing adapters or poly-N; low-quality reads were also removed from the raw data. At the same time, Q20, Q30, GC-content, and sequence duplication levels of the clean data were calculated. Ultimately, all downstream analyses were based on high-quality clean data. The chicken genome sequence (GRCg6a) was downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/genome/?term=chicken). Then, Hisat2 software was used to map these clean reads to the chicken reference genome [25]. Transcript assembly was performed using StringTie software [26]. StringTie is an algorithm based on optimization theory, which uses comparison information to construct a multivariable clipping graph and traffic network to assemble and evaluate the expression of reads according to the maximum flow algorithm. Compared with other software such as Cufflinks, it can build a more complete transcript and better evaluate expression [26]. Quantification of gene expression levels was estimated by the fragments per kilobase of transcript per million fragments (FPKM) method [27]. Based on raw counts, the reads mapped to the reference genome were quantified by StringTie software to obtain FPKM [28]. Then, DESeq R package (1.39.0), a statistical program for determining differential expression in digital gene expression data based on a model of negative binomial distribution, was used to analyze the differential gene expression between samples with biological repetition [29]. Differentially expressed genes (DEGs) were identified using the following filter criteria: foldchange≥2 and FDR<0.05. Fold change indicates the ratio of expression quantity between two samples (groups). The false discovery rate (FDR) was obtained by correcting the significant difference P value.


## 2.4. Gene Ontology (GO) and KEGG Pathway Enrichment Analysis

Gene Ontology (GO) enrichment analysis of the DEGs was implemented by GOseq R packages (http://www.geneontology.org/) based on Wallenius noncentral hypergeometric distribution [30], which could adjust for gene length bias in DEGs. Pathway significant enrichment analysis was run based on the hypergeometric test in the KEGG database (http://www.genome.jp/kegg/) to uncover the most significantly enriched pathways among the DEGs. KOBAS software [31] was used to test the statistical enrichment of DEGs in KEGG pathways.


## 2.5. Quantitative Real-Time RT-PCR Verification

Utilizing the NCBI gene bank, primers of 7 randomly selected DEGs, three different isoforms ofMyHC genes (MyHC I, IIA, and IIB) and the house-keeping gene GAPDH were designed by Primer Premier 5 software (Table S1). Primers were synthesized by Hangzhou Youkang Biotechnology Co., Ltd. (Zhejiang, China). Total RNA from cells or tissues was used to synthesize cDNA by a 5X All-In-One MasterMix Kit (abm, China). The cDNA was used as a template for quantitative PCR by using the EvaGreen 2X qPCR MasterMix Kit (abm, China). There were three replicates for each sample. The relative expression level of each related gene was calculated using the 2-△△CT statistical analysis method.


## 2.6. Isolation and Culture of Chicken Skeletal Muscle Satellite Cells In Vitro

The 10-day fertilized eggs were disinfected with ethanol. The eggshell on the upper part of the air chamber was removed using tweezers on the egg receptacle. A new pair of tweezers was used to clip out the chicken embryo. The leg muscle of the chicken embryo was isolated and placed in a new petri dish. Muscle samples were washed 3 times in phosphate buffer solution (PBS) (Logan, USA) containing penicillin-streptomycin liquid (Solarbio, China), and skin, blood vessels, adipose tissue, and connective tissue were removed. The sterile muscle tissue was cut into meat paste and digested with 0.25% trypsin-EDTA (Thermo Fisher, China) at 37°C for 20 min. DMEM/F12 (Logan, USA) containing 10% fetal bovine serum (FBS, Thermo Fisher, China) was added to terminate the digestion. The suspension was filtered through a 70μm mesh sieve and centrifuged at 1200 r/min for 8 min at room temperature. As the supernatant was discarded, cells were resuspended with DMEM/F12 containing 15% FBS and cultured in a 5% CO2 incubator at 37°C. One hour later, the cell suspension was transferred into a new cell petri dish; this was repeated twice to enrich for muscle satellite cells and eliminate fibroblasts. Cell differentiation was induced by replacing 15% FBS with 2% equine serum (Logan, USA). Once the isolated skeletal muscle satellite cells reached 70%-80% fusion, the growth medium was replaced with differentiation medium and the cells were cultured for 2, 4, and 6 days. Cells at 70%-80% fusion indicated day 0.


## 2.7. Immunofluorescence Staining

The isolated skeletal muscle satellite cells at 70%-80% fusion and the cells that were induced to differentiate after transfection were used for immunofluorescence detection following the method of Luo et al. [32]. Briefly, cells grown in 12-well plates were washed three times with precooled PBS and fixed with 4% paraformaldehyde for 15 min. Thereafter, the cells were permeabilized with 0.25% Triton X-100 per well for 10 min and blocked at 4°C overnight. Afterwards, cells were incubated with 1 : 100 diluted primary anti-Desmin (Abcam, China), anti-Pax7 (Abcam, China), or anti-MyHC (Santa, USA) antibodies for 1 hour at room temperature. Then, 1 : 100 diluted fluorescent secondary antibody (Thermo Fisher, China) was incubated with the cells for 1 hour at room temperature. After being washed by PBS three times, DAPI (Invitrogen, USA) was added to the cells and then incubated for 15 min at room temperature to stain the cell nuclei. Lastly, samples were captured using a fluorescence microscope (OLYMPUS, Japan).


## 3. Results

### 3.1. Expression of MyHC Genes in Different Muscle Tissues of Broilers

In some previous studies, compared with histochemical methods, qRT-PCR was a fast and sensitive approach to distinguish different muscle types by identifying the expression ofMyHC subunits [33, 34]. This method was also applied to identify the muscle types of chickens [35, 36]. Therefore, in order to compare the fiber types of different muscles in broilers, the expression of MYH3, a fast muscle fiber marker gene [37], and MYH1C, a slow muscle fiber marker gene [38], were detected by qRT-PCR between breast and leg muscles of 10 randomly selected broilers. The results showed that the expression of MYH3 in the breast muscle was significantly higher than that in the leg muscle, while the expression of MYH1C showed the opposite result (Figure 1). The results indicated that the breast and leg muscles of broilers were two different types of skeletal muscle. The former was mainly composed of slow muscle fibers, while the latter was mainly composed of fast muscle fibers.Figure 1
Expression ofMyHC genes in different muscle tissues of broilers. Fast muscle marker gene MYH3 and slow muscle marker gene MYH1C. ∗∗An extreme significant difference between the two groups (P<0.01); ∗a significant difference between the two groups (P<0.05); PT represents leg muscle; PX represents breast muscle, n=5.


### 3.2. Sequencing Quality Evaluation

To verify the authenticity of the sequencing data, samples were tested for quality. In Table1, it could be found that the Q20 (a basequality>20 and errorrate<0.01) and Q30 (a basequality>30 and errorrate<0.001) levels were both above 90%. Meanwhile, the GC content was between 49% and 53%. The results showed that the sequencing data had a small probability of error.Table 1
Quality assessment of sequencing data.
Sample nameClean readsClean basesGC (%)Q20 (%)Q30 (%)PT129135911870571662251.0397.6593.63PT232518744969986324451.1997.7493.87PT332658211975214774451.2797.7393.85PT4377530321127887772052.4097.9294.24PT5341460171015990375652.0597.8694.13PX125454197759069852651.3697.1192.33PX2337092651002674241251.1997.6193.31PX329558343879437972251.2197.2992.71PX425035478747856905451.8997.6493.46PX529461523880713331851.8997.5793.37It was also found that more than 75% of reads aligned to the reference genome and that the number of reads aligned to a unique location of the reference genome was between 74% and 92%, which indicated that these samples were highly available (Table2). To summarize, the sequencing data were good and the subsequent analysis could be carried out.Table 2
Statistical results of the comparison with the reference genome.
Sample nameTotal readsMapped readsMapped ratio (%)Uniquely mapUniquely mapped ratio (%)Reads mapped to “+” strandReads mapped to “+” strand ratio (%)Reads mapped to “-” strandReads mapped to “-” strand ratio (%)PT1582718225129009388.02%4872015783.61%2547017943.71%2543745443.65%PT2650374885544594585.25%5273856481.09%2756363642.38%2747627943.65%PT3653164225695756987.2%5408938882.81%2829360943.32%2822888443.22%PT4755060646426915285.12%6102129280.82%3190414542.25%3188806442.23%PT5682920345346505078.29%5082631074.42%2655233342.25%2649356838.79%PX1509083944302472684.51%4084578580.23%2142425342.08%2131694641.87%PX2674185305898436387.49%5570266882.62%2932217243.49%2925832743.40%PX3591166865137431686.90%4859343882.20%2553013143.19%2547489843.09%PX4500709564296591485.81%4078688981.46%2135165342.64%2134206142.62%PX5589230465104640286.63%4848996382.29%2537516243.06%2531836342.97%


### 3.3. Screening of DEGs

To show the distribution of gene expression, a volcano plot of the differentially expressed genes was constructed, with -log10(FDR) as the ordinate and log2(FC) as the abscissa. Because the differential expression analysis of transcriptome sequencing was an independent statistical hypothesis test for a large number of gene expression values, there was a false-positive problem. Therefore, in the process of differential expression analysis, the recognized Benjamini-Hochberg correction method was used to correct the significantP values obtained by the original hypothesis test. Finally, FDR was used as the key index for differentially expressed gene screening. In our study, FDR<0.05 was used as a screening condition for DEGs. The results showed that there were a large number of differentially expressed genes. A total of 767 DEGs were identified. Compared with leg muscle, there were 429 upregulated genes and 338 downregulated genes in the breast muscle (Figure 2 and Table S2).Figure 2
Volcano plot for differential gene expression.Meanwhile, in order to ensure the accuracy of the analysis, the screened differentially expressed genes were analyzed by hierarchical cluster analysis where genes with the same or similar expression patterns clustered. The abscissa represented the sample clustering. Among them, one column represents a sample, and clustering is based on the similarity of gene expression between samples. The more similar the gene expression between samples, the closer the distance. This result suggested that these DEGs could significantly separate the samples into the leg muscle group and breast muscle group (Figure3), indicating the reliability of our analysis.Figure 3
The heat map of differential gene expression (PT represents broiler leg muscle; PX represents breast muscle).


### 3.4. Gene Ontology (GO) Annotation and Enrichment Analysis

To further reveal the molecular characterization of the DEGs, GO annotation was used to analyze the DEGs. The GO database contains three main branches: biological process, molecular function, and cellular component. The results showed that a total of 468 DEGs were annotated into 61 GO terms, including 24 biological processes, 19 cellular components, and 18 molecular function annotations. In biological processes, the significantly enriched GO terms mainly contained cellular process, single organism process, biological regulation, metabolic process, response to stimulus, multicellular organismal process, localization, developmental process, cellular component organization, biogenesis signalling, etc. Cells, cellular parts, and organelles were significantly enriched in cellular components. Simultaneously, the main functional terms of molecular function were binding and catalytic activity, which were enriched with 273 and 183 genes, respectively (Figure4).Figure 4
GO analysis of the DEGs.In addition, some important genes (FDR<0.01) that may be involved in the transformation of muscle fiber types are enriched in these GO terms, such as mitotic cell cycle process (FHL1, PDS5A) and calcium ion binding (NKD1, RYR3) (Table S3).


### 3.5. KEGG Pathway Enrichment Analysis

To know the signal pathways where the genes were distributed, a KEGG pathway analysis was performed. A total of 230 DEGs were mapped to 126 KEGG pathways, of which 20 pathways were significantly enriched (P≤0.05) (Table S4). There were several extremely significantly enriched pathways (Figure 5), including glycolysis/gluconeogenesis, starch and sucrose metabolism, insulin signalling, and biosynthesis of amino acids. Some important genes were found to be enriched in these pathways: PKM, LDHA, LDHA, PYGM, and HRAS (FDR<0.01) (Table S4).Figure 5
Top 20 significantly enriched KEGG pathways.


### 3.6. Validation of DEGs by qRT-PCR

The expression of 7 randomly selected DEGs, including 4 downregulated genes (FHL1, MYBPC1, ADPRHL, and TPM2) and 3 upregulated genes (HOMER3, MYH3, and BPGM), was verified through quantitative real-time PCR. The results illustrated that the expression of FHL1, MYBPC1, ADPRHL, and TPM2 in the breast muscle was significantly lower than that in the leg muscle, while the expression of HOMER3, MYH3, and BPGM was exactly the opposite (Figure 6(a)). As shown in Figure 6, the expression trend of these 7 genes was consistent with the results of the RNA-seq experiment, which proved that the RNA-seq results were reliable (Figure 6(b)).Figure 6
Validation of DEGs by qRT-PCR: (a) q-PCR results and (b) RNA-seq results (∗∗P<0.01, ∗P<0.05; PT represents broiler leg muscle; PX represents breast muscle, n=5).
(a)(b)


### 3.7. Identification of Chicken Skeletal Muscle Satellite Cells

To facilitate the study of muscle fiber growth characteristics in broilers, we isolated skeletal muscle satellite cells from the leg muscle of 15 chicken embryos which were incubated for 10 days in vitro. The cells in the process of proliferation and differentiation were identified by immunofluorescence with the skeletal muscle satellite cell-specific markersPax7, Desmin, and MyHC [39]. The results showed that Pax7 was positive in the nucleus (Figure 7(a)), and Desmin and MyHC were positive in the cytoplasm (Figures 7(d) and 7(g)). Therefore, the results showed that the isolated cells are indeed skeletal muscle satellite cells and could be used in the following research.Figure 7
Identification of skeletal muscle satellite cells: (a)Pax7 was detected with immunofluorescence; (b) DAPI stain of the nuclei of skeletal muscle satellite cells; (c) merged image of (a) and (b); (d) Desmin was detected with immunofluorescence; (e) DAPI stain of the nuclei of skeletal muscle satellite cells; (f) merged image of (d) and (e); (g) MyHC was detected with immunofluorescence; (h) DAPI stain of the nuclei of skeletal muscle satellite cells; (i) merged image of (g) and (h).


### 3.8. Expression Profile of MyHC Genes

In order to further explore the growth of muscle fibers in broilers, the expression patterns ofMyHC subunits during cell differentiation were analyzed by qRT-PCR. MyHC IIA is a genetic marker of fast oxidative muscle fibers [40]. MyHC IIB is a genetic marker of fast glycolytic muscle fibers [41]. MyHC I is a genetic marker of slow oxidative muscle fibers [42]. The results suggested that the expression of MyHC IIA was significantly upregulated during cell differentiation (Figure 8(a)). However, the expression of MyHC IIB was significantly downregulated on the second day of differentiation and then gradually upregulated from day 4 to day 6 (Figure 8(b)). Additionally, the expression of MyHC I was significantly upregulated at the beginning and decreased after the 4th day of differentiation (Figure 8(c)). The results indicated different development of muscle fiber types in broilers.Figure 8
Expression ofMyHC genes in the differentiation of chicken skeletal muscle satellite cells (∗∗P<0.01, ∗P<0.05; n=3).
(a)(b)(c)


## 3.1. Expression of MyHC Genes in Different Muscle Tissues of Broilers

In some previous studies, compared with histochemical methods, qRT-PCR was a fast and sensitive approach to distinguish different muscle types by identifying the expression ofMyHC subunits [33, 34]. This method was also applied to identify the muscle types of chickens [35, 36]. Therefore, in order to compare the fiber types of different muscles in broilers, the expression of MYH3, a fast muscle fiber marker gene [37], and MYH1C, a slow muscle fiber marker gene [38], were detected by qRT-PCR between breast and leg muscles of 10 randomly selected broilers. The results showed that the expression of MYH3 in the breast muscle was significantly higher than that in the leg muscle, while the expression of MYH1C showed the opposite result (Figure 1). The results indicated that the breast and leg muscles of broilers were two different types of skeletal muscle. The former was mainly composed of slow muscle fibers, while the latter was mainly composed of fast muscle fibers.Figure 1
Expression ofMyHC genes in different muscle tissues of broilers. Fast muscle marker gene MYH3 and slow muscle marker gene MYH1C. ∗∗An extreme significant difference between the two groups (P<0.01); ∗a significant difference between the two groups (P<0.05); PT represents leg muscle; PX represents breast muscle, n=5.


## 3.2. Sequencing Quality Evaluation

To verify the authenticity of the sequencing data, samples were tested for quality. In Table1, it could be found that the Q20 (a basequality>20 and errorrate<0.01) and Q30 (a basequality>30 and errorrate<0.001) levels were both above 90%. Meanwhile, the GC content was between 49% and 53%. The results showed that the sequencing data had a small probability of error.Table 1
Quality assessment of sequencing data.
Sample nameClean readsClean basesGC (%)Q20 (%)Q30 (%)PT129135911870571662251.0397.6593.63PT232518744969986324451.1997.7493.87PT332658211975214774451.2797.7393.85PT4377530321127887772052.4097.9294.24PT5341460171015990375652.0597.8694.13PX125454197759069852651.3697.1192.33PX2337092651002674241251.1997.6193.31PX329558343879437972251.2197.2992.71PX425035478747856905451.8997.6493.46PX529461523880713331851.8997.5793.37It was also found that more than 75% of reads aligned to the reference genome and that the number of reads aligned to a unique location of the reference genome was between 74% and 92%, which indicated that these samples were highly available (Table2). To summarize, the sequencing data were good and the subsequent analysis could be carried out.Table 2
Statistical results of the comparison with the reference genome.
Sample nameTotal readsMapped readsMapped ratio (%)Uniquely mapUniquely mapped ratio (%)Reads mapped to “+” strandReads mapped to “+” strand ratio (%)Reads mapped to “-” strandReads mapped to “-” strand ratio (%)PT1582718225129009388.02%4872015783.61%2547017943.71%2543745443.65%PT2650374885544594585.25%5273856481.09%2756363642.38%2747627943.65%PT3653164225695756987.2%5408938882.81%2829360943.32%2822888443.22%PT4755060646426915285.12%6102129280.82%3190414542.25%3188806442.23%PT5682920345346505078.29%5082631074.42%2655233342.25%2649356838.79%PX1509083944302472684.51%4084578580.23%2142425342.08%2131694641.87%PX2674185305898436387.49%5570266882.62%2932217243.49%2925832743.40%PX3591166865137431686.90%4859343882.20%2553013143.19%2547489843.09%PX4500709564296591485.81%4078688981.46%2135165342.64%2134206142.62%PX5589230465104640286.63%4848996382.29%2537516243.06%2531836342.97%


## 3.3. Screening of DEGs

To show the distribution of gene expression, a volcano plot of the differentially expressed genes was constructed, with -log10(FDR) as the ordinate and log2(FC) as the abscissa. Because the differential expression analysis of transcriptome sequencing was an independent statistical hypothesis test for a large number of gene expression values, there was a false-positive problem. Therefore, in the process of differential expression analysis, the recognized Benjamini-Hochberg correction method was used to correct the significantP values obtained by the original hypothesis test. Finally, FDR was used as the key index for differentially expressed gene screening. In our study, FDR<0.05 was used as a screening condition for DEGs. The results showed that there were a large number of differentially expressed genes. A total of 767 DEGs were identified. Compared with leg muscle, there were 429 upregulated genes and 338 downregulated genes in the breast muscle (Figure 2 and Table S2).Figure 2
Volcano plot for differential gene expression.Meanwhile, in order to ensure the accuracy of the analysis, the screened differentially expressed genes were analyzed by hierarchical cluster analysis where genes with the same or similar expression patterns clustered. The abscissa represented the sample clustering. Among them, one column represents a sample, and clustering is based on the similarity of gene expression between samples. The more similar the gene expression between samples, the closer the distance. This result suggested that these DEGs could significantly separate the samples into the leg muscle group and breast muscle group (Figure3), indicating the reliability of our analysis.Figure 3
The heat map of differential gene expression (PT represents broiler leg muscle; PX represents breast muscle).


## 3.4. Gene Ontology (GO) Annotation and Enrichment Analysis

To further reveal the molecular characterization of the DEGs, GO annotation was used to analyze the DEGs. The GO database contains three main branches: biological process, molecular function, and cellular component. The results showed that a total of 468 DEGs were annotated into 61 GO terms, including 24 biological processes, 19 cellular components, and 18 molecular function annotations. In biological processes, the significantly enriched GO terms mainly contained cellular process, single organism process, biological regulation, metabolic process, response to stimulus, multicellular organismal process, localization, developmental process, cellular component organization, biogenesis signalling, etc. Cells, cellular parts, and organelles were significantly enriched in cellular components. Simultaneously, the main functional terms of molecular function were binding and catalytic activity, which were enriched with 273 and 183 genes, respectively (Figure4).Figure 4
GO analysis of the DEGs.In addition, some important genes (FDR<0.01) that may be involved in the transformation of muscle fiber types are enriched in these GO terms, such as mitotic cell cycle process (FHL1, PDS5A) and calcium ion binding (NKD1, RYR3) (Table S3).


## 3.5. KEGG Pathway Enrichment Analysis

To know the signal pathways where the genes were distributed, a KEGG pathway analysis was performed. A total of 230 DEGs were mapped to 126 KEGG pathways, of which 20 pathways were significantly enriched (P≤0.05) (Table S4). There were several extremely significantly enriched pathways (Figure 5), including glycolysis/gluconeogenesis, starch and sucrose metabolism, insulin signalling, and biosynthesis of amino acids. Some important genes were found to be enriched in these pathways: PKM, LDHA, LDHA, PYGM, and HRAS (FDR<0.01) (Table S4).Figure 5
Top 20 significantly enriched KEGG pathways.


## 3.6. Validation of DEGs by qRT-PCR

The expression of 7 randomly selected DEGs, including 4 downregulated genes (FHL1, MYBPC1, ADPRHL, and TPM2) and 3 upregulated genes (HOMER3, MYH3, and BPGM), was verified through quantitative real-time PCR. The results illustrated that the expression of FHL1, MYBPC1, ADPRHL, and TPM2 in the breast muscle was significantly lower than that in the leg muscle, while the expression of HOMER3, MYH3, and BPGM was exactly the opposite (Figure 6(a)). As shown in Figure 6, the expression trend of these 7 genes was consistent with the results of the RNA-seq experiment, which proved that the RNA-seq results were reliable (Figure 6(b)).Figure 6
Validation of DEGs by qRT-PCR: (a) q-PCR results and (b) RNA-seq results (∗∗P<0.01, ∗P<0.05; PT represents broiler leg muscle; PX represents breast muscle, n=5).
(a)(b)


## 3.7. Identification of Chicken Skeletal Muscle Satellite Cells

To facilitate the study of muscle fiber growth characteristics in broilers, we isolated skeletal muscle satellite cells from the leg muscle of 15 chicken embryos which were incubated for 10 days in vitro. The cells in the process of proliferation and differentiation were identified by immunofluorescence with the skeletal muscle satellite cell-specific markersPax7, Desmin, and MyHC [39]. The results showed that Pax7 was positive in the nucleus (Figure 7(a)), and Desmin and MyHC were positive in the cytoplasm (Figures 7(d) and 7(g)). Therefore, the results showed that the isolated cells are indeed skeletal muscle satellite cells and could be used in the following research.Figure 7
Identification of skeletal muscle satellite cells: (a)Pax7 was detected with immunofluorescence; (b) DAPI stain of the nuclei of skeletal muscle satellite cells; (c) merged image of (a) and (b); (d) Desmin was detected with immunofluorescence; (e) DAPI stain of the nuclei of skeletal muscle satellite cells; (f) merged image of (d) and (e); (g) MyHC was detected with immunofluorescence; (h) DAPI stain of the nuclei of skeletal muscle satellite cells; (i) merged image of (g) and (h).


## 3.8. Expression Profile of MyHC Genes

In order to further explore the growth of muscle fibers in broilers, the expression patterns ofMyHC subunits during cell differentiation were analyzed by qRT-PCR. MyHC IIA is a genetic marker of fast oxidative muscle fibers [40]. MyHC IIB is a genetic marker of fast glycolytic muscle fibers [41]. MyHC I is a genetic marker of slow oxidative muscle fibers [42]. The results suggested that the expression of MyHC IIA was significantly upregulated during cell differentiation (Figure 8(a)). However, the expression of MyHC IIB was significantly downregulated on the second day of differentiation and then gradually upregulated from day 4 to day 6 (Figure 8(b)). Additionally, the expression of MyHC I was significantly upregulated at the beginning and decreased after the 4th day of differentiation (Figure 8(c)). The results indicated different development of muscle fiber types in broilers.Figure 8
Expression ofMyHC genes in the differentiation of chicken skeletal muscle satellite cells (∗∗P<0.01, ∗P<0.05; n=3).
(a)(b)(c)


## 4. Discussion

The skeletal muscle of vertebrates is composed of muscle fibers with different morphological, contractile, and metabolic characteristics. Meanwhile, the properties of muscle fibers have an important effect on the meat quality of livestock [8]. In our study, we found that slow/fast muscle fiber marker genes were differentially expressed between the breast and leg muscles of broilers. The results were consistent with several previous studies demonstrating that the breast muscle of the chicken was a fast muscle, while the leg muscle was a slow muscle [43]. Therefore, RNA-seq was performed on these two skeletal muscles to screen candidate genes related to muscle fiber-type conversion in our research. A total of 767 DEGs were discovered between both investigated groups, of which 7 DEGs were randomly selected and verified by qRT-PCR analysis. Although the expression fold changes of the selected genes between the two methods were not the same, the expression trend was highly consistent, proving that the results obtained by chip analysis were credible.The basal processes of muscle cell growth and muscle fiber formation in broilers include cellular processes, single biological processes, and biological regulation. GO functional analysis showed that 11 DEGs were enriched in GO term mitotic cell cycle processes (TableS3). Among these DEGs, the FHL1 gene was found to be highly abundant in oxidative muscle fibers [44]. Overexpression of FHL1 promoted the transformation of muscle fibers to the oxidative type, while knockdown of FHL1 could cause muscle fiber hypertrophy [45]. Slow muscle fibers are usually small in diameter, whereas the fast muscle fibers are usually large in diameter [46]. In this study, the PDS5A gene was also screened as an upregulated DEG in the breast muscle. PDS5A was found to be highly expressed in the skeletal muscle but low in other tissues and to have an effect on the proliferation of skeletal muscle cells [47]. These genes might play a key role in regulating muscle cell proliferation and muscle fiber transformation. In addition, intracellular calcium signalling is also an important factor that influences muscle fiber types [48]. Our study showed that 15 DEGs were enriched in the GO term calcium ion binding (Table S3). Among these genes, NKD1 was upregulated in leg muscles. Some studies have shown that Wnt/β-catenin signalling could transform slow muscle fibers into fast muscle fibers and that NKD1 has a strong antagonistic effect on Wnt/β-catenin signalling [49, 50]. Another gene, RYR3, controls calcium ion release channels in muscle cells [51]. RYR3 has a high calcium sensitivity and regulates rapid muscle contraction together with the RYR1 gene [52]. Some studies have shown that calcium sensitivity is related to the generation of fast myofibril subtypes, which was consistent with the results of our study [53]. KEGG pathway analysis also showed that the expression of calcium signalling pathway-related genes was altered significantly (Table S4). The function of these aforementioned DEGs in skeletal muscle fiber transformation needs further investigation.One of the methods to classify fast or slow muscle is by their energy metabolism [54]. GO functional annotation analysis showed that a total of 198 DEGs were enriched in metabolic process-related GO terms (Table S3). Previous studies showed that slow-twitch muscles have lower glycolytic enzyme activity, while fast-twitch muscles act mainly through the glycolytic pathway to meet the energy requirements for rapid muscle contraction [55, 56]. In our study, the KEGG pathway enrichment analysis results indicated that the four most significant enrichment pathways were related to the energy metabolism of muscle cells. Most of the DEGs involved in the glycolysis/gluconeogenesis pathway were upregulated in the breast muscle (15 of 17) (Table S4). Protein synthesis during skeletal muscle development is accompanied by energy metabolism, and the glycolysis/gluconeogenesis pathway plays a significant role [57]. PKM is a gene that catalyses the last step of glycolysis [58] and was upregulated in the breast muscle. The muscle type of pyruvate kinase (PKM) is one of the key mediators of the Warburg effect and tumor metabolism [59] and was closely related to the development of PSE meat [60]. Another gene, LDHA, the predominant LDH isoform in the skeletal muscle, is also an important enzyme involved in glycolysis [61]. The mutation at the C.423A>G site in exon 6 of porcine LDHA had a significant effect on the muscle fiber type II ratio [62]. In addition, 13 DEGs were enriched in the biosynthesis of amino acid pathway, of which 12 DEGs were upregulated in the breast muscle (Table S4). Another study found that genes related to glycolysis/gluconeogenesis and the biosynthesis of amino acid pathway were gradually upregulated during the growth and development of chicken breast muscle [63]. Therefore, we inferred that these DEGs were activated during fast-twitch muscle fiber formation.Our KEGG results further discovered that 22 DEGs were enriched in the insulin signalling pathway (TableS4). Insulin plays a key role in a series of biological reactions, such as glucose uptake, lipid and protein synthesis, cell growth, and glycogen synthesis [64]. In skeletal muscle cells, the insulin signalling pathway regulates glycogen synthesis and decreases the rate of glycogen breakdown [65]. It has been shown that the content of glycogen in type I skeletal muscle fibers is significantly lower than that in type II muscle fibers [66]. We found that 5 DEGs in the insulin signalling pathway were also enriched in the GO term glycogen metabolic process (Table S3), including GYS1 and PYGM; these were all upregulated in the breast muscle. The insulin signalling pathway stimulates the expression of GYS1 mRNA and promotes the synthesis of glycogen in the skeletal muscle [67]. Skeletal muscle fiber type IIB predominantly uses glycogen and glucose as fuel, and the upregulation of PYGM accelerates glycogen phosphorylation [68]. In addition, 3 DEGs in the insulin signalling pathway were upregulated in the leg muscle, including HRAS (Table S4). The expression of HRAS was negatively correlated with metabolic traits related to the insulin signalling pathway [69]. Thus, HRAS may be a candidate gene involved in slow muscle fiber formation.Skeletal muscle fibers are generally divided into four types according to the myosin heavy chain (MyHC) isoforms, namely, types I, IIA, IIB, and IIX [70]. However, there are very few type IIX muscle fibers in poultry skeletal muscle [71]. A previous study showed that Pax7, Desmin, and MyHC were specific genetic markers of skeletal muscle satellite cells [72]. To further explore the mechanism of skeletal muscle fiber transformation in broilers, an in vitro cell model was successfully constructed in our study. Then, the expression profiles of MyHC isoforms (type I, type IIA, and type IIB) during broiler skeletal muscle satellite cell differentiation were detected. MyHC IIA mRNA expression increased throughout the time course. Expression of the MyHC IIB isoform was low on day 2 but rose until day 6. MyHC I expression first grew dramatically from day 0 to day 4 and then declined until day 6. Interestingly, these results coincided with those in a C2C12 mouse myoblast cell differentiation model [73]. Therefore, future studies will use these approximate time points of skeletal muscle fiber transformation in vitro to explore the function of the key DEGs in broilers.


## 5. Conclusion

In the current research, RNA-seq was carried out to compare the leg and breast muscles of broilers. A total of 767 DEGs were identified. Compared with leg muscle, there were 429 upregulated genes and 338 downregulated genes in the breast muscle. The expression profile ofMyHC isoforms during broiler skeletal muscle cell differentiation was detected. Therefore, our study provides useful information for understanding the molecular mechanisms of muscle fiber-type transition in broilers. However, the transformation of skeletal muscle fiber in vivo is a dynamic process and these transcriptome data are preliminary, so the function of the DEGs still requires further investigation.



---

*Source: 9478949-2020-07-06.xml*