Transcriptome Analysis of Differentially Expressed mRNA Related to Pigeon Muscle Development

Simple Summary The growth and development of skeletal muscle determine the meat production performance of pigeons and are regulated by complex gene networks. To explore the genes involved in regulating the growth and development of pigeon skeletal muscle, RNA sequencing (RNA−seq) was used to characterise gene expression profiles during the development and growth of pigeon breast muscle and identify differentially expressed genes (DEGs) among different stages. This study expanded the diversity of pigeon mRNA, and it was helpful to understand the role of mRNA in pigeon muscle development and growth. Abstract The mechanisms behind the gene expression and regulation that modulate the development and growth of pigeon skeletal muscle remain largely unknown. In this study, we performed gene expression analysis on skeletal muscle samples at different developmental and growth stages using RNA sequencing (RNA−Seq). The differentially expressed genes (DEGs) were identified using edgeR software. Weighted gene co−expression network analysis (WGCNA) was used to identify the gene modules related to the growth and development of pigeon skeletal muscle based on DEGs. A total of 11,311 DEGs were identified. WGCNA aggregated 11,311 DEGs into 12 modules. Black and brown modules were significantly correlated with the 1st and 10th day of skeletal muscle growth, while turquoise and cyan modules were significantly correlated with the 8th and 13th days of skeletal muscle embryonic development. Four mRNA−mRNA regulatory networks corresponding to the four significant modules were constructed and visualised using Cytoscape software. Twenty candidate mRNAs were identified based on their connectivity degrees in the networks, including Abca8b, TCONS−00004461, VWF, OGDH, TGIF1, DKK3, Gfpt1 and RFC5, etc. A KEGG pathway enrichment analysis showed that many pathways were related to the growth and development of pigeon skeletal muscle, including PI3K/AKT/mTOR, AMPK, FAK, and thyroid hormone pathways. Five differentially expressed genes (LAST2, MYPN, DKK3, B4GALT6 and OGDH) in the network were selected, and their expression patterns were quantified by qRT−PCR. The results were consistent with our sequencing results. These findings could enhance our understanding of the gene expression and regulation in the development and growth of pigeon muscle.


Introduction
Pigeon meat is rich in nutrition, high in protein, low in fat, and high in medicinal value. In China, pigeon meat is called "animal ginseng" and is considered an advanced nourishment product that is increasingly favoured by consumers [1]. Meat production performance is an important index to measure the economic value of pigeons. However, the genetic improvement of the meat production performance of pigeons is relatively lagging in comparison with other poultry. The growth and development of skeletal muscle determine the meat production performance of pigeons [2]. Therefore, understanding the molecular regulation mechanism of pigeon skeletal muscle growth and development is an important prerequisite for improving meat production performance by molecular breeding technology [3].
The development of skeletal muscle is closely related to skeletal muscle cell differentiation. Myogenic regulatory factor family (MRF) members play essential roles in the process of skeletal muscle cell differentiation, including MyoD (myogenic determining factor), MyoG (myogenin), MRF4 (myogenic regulatory factor 4), and Myf5 (myogenic factor 5) [4]. However, the process of skeletal muscle growth and development involves multi−gene expression, signal transduction, and network regulation, and there are still a large number of regulatory factors to be identified.
RNA sequencing (RNA−seq) can explore the differences of gene types and expression levels at the overall level and directly link the changes of gene expression levels with phenotypic changes. In recent years, to explore the molecular mechanism of critical economic traits, considerable transcriptome studies on livestock and poultry were carried out using RNA−seq technology [5]. Xing et al. revealed essential genes related to muscle fat and abdominal fat deposition in chickens during the development process and identified 21 key genes in total, using RNA−Seq analysis [6]. Hu et al. identified several genes and pathways that may regulate skeletal muscle growth in the black Muscovy duck using RNA−seq [7]. RNA−seq technology has also been applied to pigeon transcriptome analysis. Wang et al. performed a comprehensive investigation into miRNA transcriptomes in livers across three pigeon developmental stages using RNA−seq and identified several vital target genes (e.g., TNRC6B, FRS2, PTCH1, etc.) of DE miRNA, which is closely linked to liver development [8]. Ye et al. compared the transcriptomes of muscle and liver tissues between squabs of two breeds to identify candidate genes associated with the differences in the capacity of fat deposition. A total of 27 genes were identified as a basis for further investigations to screen markers closely associated with intramuscular fat content and fatty acid composition in squabs [9].
However, the application of RNA−Seq technology in mining genes related to pigeon skeletal muscle growth and development has not been reported yet. The regulatory mechanisms of pigeon skeletal muscle growth and development remain poorly understood. Therefore, the present study aims to characterise gene expression profiles during pigeon muscle development and growth and identify key genes involved in this biological process. Our findings will contribute to a better understanding of the mechanisms by which genes regulate pigeon muscle development and growth.

Animal Ethics Statement
This study was performed following the Chinese guidelines for animal welfare, and the animal protocol was approved by the Animal Welfare Committee of Yangzhou University (permit number SYXK [Su] 2016-0020).

Experimental Animals
The Tarim pigeon is a local breed of semi−wild pigeons in Yeerqiang River Basin, which has strong disease resistance, high reproductive capacity and good stress resistance. It is mainly distributed in Hotan, Kashgar and Aksu. Tarim pigeons are a kind of meat and egg dual−purpose pigeon, which primarily feeds on various grains and can go out in groups.

Sample Preparation
Twelve pigeons were obtained from Wuxi Sanxiangan Agricultural Technology Development Co., Ltd. (Jiangsu, China). The fresh left breast muscle tissues were collected, immediately frozen in liquid nitrogen, and stored in liquid nitrogen until use. A total of 12 samples of four stages (3 replicates for each stage, 8 and 13 embryonic age, 1 and 10 day-olds) were collected.

Library Construction and RNA Sequencing
A total of twelve cDNA libraries were constructed with muscle tissues, and 3 µg total RNA per sample was used as the input material for a cDNA library. After total RNAs were extracted, rRNAs were removed, and then the enriched RNAs were fragmented into short fragments and reverse transcribed into cDNAs. Double−stranded cDNAs were synthesised by replacing dTTPs with dUTPs in the reaction buffer used in second−strand cDNA synthesis. The resulting double−stranded cDNAs were ligated to adaptors after being end−repaired and A−tailed. Then, uracil−N−glycosylase (UNG) was used to digest the second−strand cDNAs. The digested products were size selected by agarose gel electrophoresis, PCR amplified and sequenced by Gene Denovo Biotechnology Co. (Guangzhou, China) using Illumina HiSeq 4000 (150 bp paired−end).

mRNA Identification and Quantification
High−quality clean reads were obtained by removing reads containing adapters, reads consisting of all A bases, reads containing more than 10% of unknown nucleotides (N), and reads containing more than 50% of low quality (Q-value ≤ 20) bases. The high−quality clean reads were aligned against the ribosome RNA (rRNA) database using Bowtie2 [10] (version 2.2.8). The rRNA mapped reads then were removed. The remaining clean reads of each sample were then mapped to the Columba livia reference genome (https://www.ncbi. nlm.nih.gov/assembly/GCA_000337935.2/, 13 December 2017) by TopHat2 [11] (version 2.1.1) with default parameters, respectively. Fragments per kilobase of transcript per million mapped reads (FPKM) algorithm was used to quantify the expression levels of mRNAs.

Identification of Differentially Expressed Genes
The edgeR software was used to identify differentially expressed mRNAs. Genes with p-value < 0.05 and fold change ≥2 were considered significantly differentially expressed genes (DEGs).

Co−Expression Network Analysis and Visualisation
Co−expression network analysis was analysed using the weighted gene co−expression network analysis (WGCNA) R package [12] based on the DEGs. The co−expressed modules were detected using the automatic network construction function blockwiseModules with power = 5. The hierarchical clustering of identified modules was conducted by applying the Dynamic Tree Cut algorithm. Each of these modules was summarised by its first principal component referred to as its eigengene, providing a single value for a module's expression profile [13]. In order to identify modules associated with pigeon skeletal muscle development, eigengenes were correlated with the different stages of development and growth. Modules with an average eigengene−stage correlation of r > 0.4 were considered significantly associated with pigeon skeletal muscle development and growth. The top 1000 interactions of each significant module were exported and visualised using Gephi software (version 0.9.2) [14].

Functional Enrichment Analysis
To explore the functions of the significant modules, genes in each module were submitted to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The GO biological process terms and KEGG pathways with p value < 0.01 were considered significantly enriched GO terms and KEGG pathways.

qRT−PCR Confirmation of Differentially Expressed Genes
mRNA was reverse−transcribed respectively using HiScript Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme, Nanjing, China) following the manufacturer's recommendations. The mRNA forward primers were obtained commercially from Tsingke Biotechnology (Beijing, China) ( Table 1). The ChamQ SYBR qPCR Master Mix (Low ROX Premixed) (Vazyme, Beijing, China) was used to perform the qRT−PCR procedure following the manufacturer's recommendations. The reaction mixtures were added in a 96−well plate at 95 • C for 30 s, followed by 40 cycles of 95 • C for 10 s and 60 • C for 30 s, and then followed by a melting curve using an ABI 7500 Real−Time PCR System (Applied Biosystems, Foster City, CA, USA). All experiments were performed in triplicate for each biological repeat. Quantification of selected gene expression was performed using the comparative threshold cycle (2 −∆∆Ct ) method. The quantitative real−time PCR (RT−qPCR) results for all genes were statistically tested using the Student's t-test (version SPSS 20.0).

Quality Control
A total of 1,003,374,810 clean reads were obtained. The GC content of the clean reads across the 12 samples varied from 43.49% to 45.84%. The percentage of Q20 bases and Q30 bases for the clean reads were above 97.85% and 93.30%, respectively. After quality control, the percentage of high−quality clean reads of the 12 samples were above 98.89%, indicating that the sequencing data can be used for subsequent data analysis (Table 2).

Identification of Differentially Expressed Genes
A total of 30,900 transcripts (including 16,522 known transcripts and 14,378 new transcripts) were identified (Table S1). Figure 1B showed the correlation between the samples based on gene expression measured using Pearson's correlation coefficient. The results showed that the genes in pigeon skeletal muscle have developmental stage− and time−specific patterns. Samples of the embryonic stage (E8 and E13) are tightly clustered into the same subgroup, while the 1-day-old and 10-day-old samples were clustered into distinct subgroups, indicating varying gene expression of pigeon skeletal muscle between the embryonic and the growth stage.

Identification of Differentially Expressed Genes
A total of 30,900 transcripts (including 16,522 known transcripts and 14,378 new transcripts) were identified (Table S1). Figure 1B showed the correlation between the samples based on gene expression measured using Pearson's correlation coefficient. The results showed that the genes in pigeon skeletal muscle have developmental stage− and time−specific patterns. Samples of the embryonic stage (E8 and E13) are tightly clustered into the same subgroup, while the 1−day−old and 10−day−old samples were clustered into distinct subgroups, indicating varying gene expression of pigeon skeletal muscle between the embryonic and the growth stage.
A total of 11,311 genes were identified as DEGs among the four stages (p < 0.05, FC > 2) ( Table S2)   A total of 11,311 genes were identified as DEGs among the four stages (p < 0.05, FC > 2) (Table S2). The E8−D10 comparison group had the largest DEGs of 8242, including 5874 downregulated DEGs and 2368 upregulated DEGs. In contrast, the E8−E13 comparison group had the least number of DEGs of 1959, of which 1057 were downregulated, and 902 were upregulated ( Figure 1A,C).

Weighted Gene Co−Expression Network Analysis
To investigate the relationship between DEGs and muscle development and growth, we performed a WGCNA analysis based on the expression of 11,311 DEGs. WGCNA categorised 11,311 DEGs into 12 modules with a soft power of 5 (Figure 2A,B). Among the 12 modules, black and brown were significantly correlated with day 1 and 10 stages, while turquoise and cyan were significantly correlated with embryonic development stages of days 8 and 13 (r > 0.4) ( Figure 2C). blue to red represent low to high correlation. (C) Volcano plot of differential gene expression analysis. Upregulated genes are marked red, downregulated genes in yellow, and nonsignificant genes in blue, respectively.

Weighted Gene Co−Expression Network Analysis
To investigate the relationship between DEGs and muscle development and growth, we performed a WGCNA analysis based on the expression of 11,311 DEGs. WGCNA categorised 11,311 DEGs into 12 modules with a soft power of 5 (Figure 2A,B). Among the 12 modules, black and brown were significantly correlated with day 1 and 10 stages, while turquoise and cyan were significantly correlated with embryonic development stages of days 8 and 13 (r > 0.4) ( Figure 2C).

Visualisation of the Pigeon Muscle Development and Growth−Related Modules
To further explore candidate DEGs that regulate pigeon skeletal muscle growth and development, we visualised the four significant modules identified by WGCNA based on the top 1000 mRNA−mRNA interaction pairs in each module. Four mRNA−mRNA regulatory networks were constructed, corresponding to black, brown, cyan, and turquoise modules, respectively. The black, brown, cyan and turquoise modules contained 265, 308, 276, and 402 mRNAs, respectively (Figure 3).

Visualisation of the Pigeon Muscle Development and Growth−Related Modules
To further explore candidate DEGs that regulate pigeon skeletal muscle growth and development, we visualised the four significant modules identified by WGCNA based on the top 1000 mRNA−mRNA interaction pairs in each module. Four mRNA−mRNA regulatory networks were constructed, corresponding to black, brown, cyan, and turquoise modules, respectively. The black, brown, cyan and turquoise modules contained 265, 308, 276, and 402 mRNAs, respectively (Figure 3).

Identification of Candidate Genes Related to Pigeon Muscle Development and Growth
We used Cytoscape software to calculate the degree of connectivity of the genes in the mRNA−mRNA network, and DEGs with high connectivity were considered to play an essential role in the network. In this study, the top five DEGs ranked by the degree of connectivity in each network were considered as hub genes. A total of 20 hub genes were identified ( Table 3). The turquoise and cyan were development−specific modules, indicating that hub genes in turquoise and cyan modules might play potential roles in pigeon skeletal muscle development. The black and brown modules were significantly correlated to the 1−day−old and 10−day−old stages, suggesting that hub genes in the black and brown modules were involved in regulating the early growth stage of pigeon skeletal muscle.

Identification of Candidate Genes Related to Pigeon Muscle Development and Growth
We used Cytoscape software to calculate the degree of connectivity of the genes in the mRNA−mRNA network, and DEGs with high connectivity were considered to play an essential role in the network. In this study, the top five DEGs ranked by the degree of connectivity in each network were considered as hub genes. A total of 20 hub genes were identified ( Table 3). The turquoise and cyan were development−specific modules, indicating that hub genes in turquoise and cyan modules might play potential roles in pigeon skeletal muscle development. The black and brown modules were significantly correlated to the 1−day−old and 10−day−old stages, suggesting that hub genes in the black and brown modules were involved in regulating the early growth stage of pigeon skeletal muscle.

GO and KEGG Pathway Enrichment Analysis
To further explore the function of DEGs in the four significant modules and the possible molecular mechanism by which DEGs regulated pigeon skeletal development and growth, we performed GO and KEGG pathway enrichment analysis on the DEGs in the black, brown, cyan and turquoise network, respectively. GO enrichment analysis showed that DEGs in the black module were significantly enriched in 91 biological processes (p < 0.01), such as cellular protein localisation, cellular macromolecule localisation and protein import ( Figure 4A). DEGs in the brown module were significantly categorised into 47 biological processes (p < 0.01), including the generation of precursor metabolites and energy, the oxidation−reduction process and energy derivation by oxidation of organic compounds ( Figure 4B). DEGs in the cyan module were assigned to 14 significantly enriched biological processes, such as muscle organ development, response to osmotic stress and striated muscle tissue development ( Figure 4C). DEGs in the turquoise module were significantly categorised into 78 biological processes, including DNA−templated transcription, initiation and negative regulation of RNA metabolic process ( Figure 4D). The significantly enriched GO molecular function and cellular component terms are listed in Table S3. KEGG pathway enrichment showed that DEGs in the black module were enriched in signalling pathways such as PI3K−Akt signaling pathway and human papillomavirus infection ( Figure 5A). DEGs in the brown module were enriched in signalling pathways such as thermogenesis and metabolic pathways ( Figure 5B). DEGs in the cyan module were enriched in signalling pathways such as the rap1 signalling pathway and regulation of actin cytoskeleton ( Figure 5C). DEGs in the turquoise module were enriched in signalling pathways such as viral carcinogenesis and microRNAs in cancer and RNA transport ( Figure 5D).

qRT−PCR
In order to confirm the differentially expressed genes among the modules obtained by RNA−Seq, five differentially expressed genes (LAST2, MYPN, DKK3, B4GALT6 and OGDH) in the network were selected, and their expression patterns were quantified by qRT−PCR. The results are consistent with our sequencing results, highlighting the reliability of our sequencing data ( Figure 6).

qRT−PCR
In order to confirm the differentially expressed genes among the modules obtained by RNA−Seq, five differentially expressed genes (LAST2, MYPN, DKK3, B4GALT6 and OGDH) in the network were selected, and their expression patterns were quantified by qRT−PCR. The results are consistent with our sequencing results, highlighting the reliability of our sequencing data ( Figure 6).

Discussion
Elucidating the molecular regulation mechanism of pigeon skeletal muscle growth and development is an essential prerequisite for using molecular breeding technology to improve pigeon meat production performance. However, no previous studies on identifying candidate genes and exploring mechanisms have been conducted in pigeon skeletal muscle development and growth to our best knowledge.
One of the aims of this study was to characterise gene expression patterns during pigeon skeletal muscle development and growth at the mRNA level using RNA−seq. We found a significant difference in gene expression patterns between the development stage and the growth stage of pigeon skeletal muscle. It suggests distinct mechanisms of genes to regulate the development and growth in pigeon skeletal muscle. In the study of Li et al., co−expression analysis was used to determine that four modules were related to the specific growth stage of chicken breast muscle development [15].
Weighted gene co−expression network analysis is a systems biology method for describing the correlation patterns among genes across microarray samples. WGCNA can be used for finding clusters (modules) of highly correlated genes, for summarising such clusters using the module eigengene or an intramodular hub gene, for relating modules to one another and to external sample traits (using eigengene network methodology), and for calculating module membership measures [16]. For example, Chen et al. used WGCNA to identify some key genes in abdominal aortic aneurysms [17]. Li et al. analysed

Discussion
Elucidating the molecular regulation mechanism of pigeon skeletal muscle growth and development is an essential prerequisite for using molecular breeding technology to improve pigeon meat production performance. However, no previous studies on identifying candidate genes and exploring mechanisms have been conducted in pigeon skeletal muscle development and growth to our best knowledge.
One of the aims of this study was to characterise gene expression patterns during pigeon skeletal muscle development and growth at the mRNA level using RNA−seq. We found a significant difference in gene expression patterns between the development stage and the growth stage of pigeon skeletal muscle. It suggests distinct mechanisms of genes to regulate the development and growth in pigeon skeletal muscle. In the study of Li et al., co−expression analysis was used to determine that four modules were related to the specific growth stage of chicken breast muscle development [15].
Weighted gene co−expression network analysis is a systems biology method for describing the correlation patterns among genes across microarray samples. WGCNA can be used for finding clusters (modules) of highly correlated genes, for summarising such clusters using the module eigengene or an intramodular hub gene, for relating modules to one another and to external sample traits (using eigengene network methodology), and for calculating module membership measures [16]. For example, Chen et al. used WGCNA to identify some key genes in abdominal aortic aneurysms [17]. Li et al. analysed the gene co−expression network and functional modules of chicken breast muscle at different developmental stages [15]. Zhang et al. used WGCNA to analyse microbial and metabolic data sets [18]. WGCNA is widely used in biological analysis.
We constructed the gene regulatory network in each module to screen further the key mRNA regulating the growth and development of pigeon skeletal muscle. In the four modules, according to the degree, a total of 20 candidate key mRNAs were identified, including the tumoural expression of the ABC transporter A8b (ABCA8b), TCONS−00004461, von Willebrand factor (VWF), 2−oxoglutarate dehydrogenase (OGDH), TGF−beta induced factor homeobox 1 (TGIF1), Dickkopf homolog 3 (DKK3), glutamine−fructose−6−phosphate transaminase 1 (GFPT1) and replication factor C 5 (RFC5). Many of these mRNAs are involved in the regulation of muscle growth and development. For example, the kinesin family member 1C (KIF1C) is known to regulate podosomes, actin−rich adhesion structures that remodel the extracellular matrix during physiological processes [19], Gache et al. found that KIF1C affects the differentiation of muscle cells [20]. Myopalladin (MYPN) is a striated muscle−specific, immunoglobulin−containing protein located in the Z−line and I−band of the sarcomere as well as the nucleus. MYPN promotes skeletal muscle growth by activating the serum response factor (SRF) pathway in muscle [21]. DKK3 is a stress−induced, renal tubular epithelia−derived, secreted glycoprotein (molecular mass, 38 kDa) [22] that induces tubulointerstitial fibrosis through its action on the canonical Wnt/β−catenin signalling pathway [23]. DKK3 is a divergent member of the DKK protein family. Wilde searched for new genes related to specific muscle types and found that the expression of DKK3 in mouse quadriceps was higher than that in other tissues [24]. Yin et al. confirmed that Dkk3 is the key secretory factor of muscle production [25]. The remaining key genes, such as Abca8b, TCONS−00004461, VWF, OGDH, TGIF1, Gfpt1 and RFC5 have not been reported in muscle growth and development.
We analysed the mRNA of black, brown, cyan and turquoise modules by GO and KEGG enrichment to explore the possible mechanism of mRNA regulation of pigeon skeletal muscle growth and development. Go enrichment showed that genes in the cyan module were enriched in many biological processes related to the growth and development of skeletal muscle, including muscle organ development, striated muscle tissue development, muscle tissue morphogenesis, muscle tissue development, muscle organ morphogenesis, muscle structure development and embryo development. It suggests that the mRNA in cyan may play an important regulatory role in the differentiation and development of pigeon skeletal muscle cells. The mRNA in black, brown and turquoise is significantly enriched in biological processes, such as biological process regulation, cell process regulation and tissue morphogenesis, indicating that the above biological processes may be related to the growth of pigeon skeletal muscle.
Pigeon muscle growth is a complex process affected by multiple genes and regulated by multiple pathways. Many pathways of mRNA enrichment in the four modules are related to muscle growth and development. GRP94 promotes muscle differentiation by inhibiting PI3K/AKT/mTOR signalling pathway [26]. AMPK plays an important role in controlling skeletal muscle development and growth due to its effects on anabolism and catalytic cellular processes [27]. Quach [28] found that FAK regulates the expression of a set of muscle−specific genes specifically involved in myoblast fusion during early myogenic differentiation, including β1D integrins. Cui et al. found that ISLR can relieve skeletal muscle atrophy and prevent muscle cell apoptosis through the IGF/PI3K/AKT−FOXO signalling pathway [29]. Ge et al. [30] studied the important effect of the interaction between bioactive nanomaterials and muscle cells on enhancing skeletal muscle tissue regeneration and found that bioactive nanomaterials can activate the p38+MAPK signalling pathway and enhance myogenic differentiation. Jaafar's research found that the size of muscle cells is positively regulated by phospholipase D (PLD). Phospholipase D regulates the size of skeletal muscle cells by activating mTOR signals, thereby affecting muscle growth and development. Ambrosio [31] found that the proliferation and differentiation of skeletal muscle stem cells are completely dependent on the role of the thyroid hormone.
Ribosome biogenesis in eukaryotes has become an important regulator of skeletal muscle growth and maintenance [32].

Conclusions
Our research results identified 11,311 DEGs and four co−expressed gene modules associated with the development and growth of pigeon skeletal muscle. Twenty candidate genes were identified, including Abca8b, TCONS−00004461, VWF, OGDH, TGIF1, DKK3, Gfpt1 and RFC5. Our results profiled gene expression of pigeon skeletal muscle samples at different stages and identified DEGs, which will contribute to our understanding of the mechanisms underlying the development and growth of pigeon skeletal muscle.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ani11082311/s1, Table S1. Statistics of the number of transcripts detected of each sample. Table S2. Differentially expressed mRNAs among the four groups. Table S3. Significantly enriched GO molecular function and cellular component terms of DEGs in the four significant modules.