Comparative Transcriptome Profiling of Skeletal Muscle from Black Muscovy Duck at Different Growth Stages Using RNA-seq

In China, the production for duck meat is second only to that of chicken, and the demand for duck meat is also increasing. However, there is still unclear on the internal mechanism of regulating skeletal muscle growth and development in duck. This study aimed to identity candidate genes related to growth of duck skeletal muscle and explore the potential regulatory mechanism. RNA-seq technology was used to compare the transcriptome of skeletal muscles in black Muscovy ducks at different developmental stages (day 17, 21, 27, 31, and 34 of embryos and postnatal 6-month-olds). The SNPs and InDels of black Muscovy ducks at different growth stages were mainly in “INTRON”, “SYNONYMOUS_CODING”, “UTR_3_PRIME”, and “DOWNSTREAM”. The average number of AS in each sample was 37,267, mainly concentrated in TSS and TTS. Besides, a total of 19 to 5377 DEGs were detected in each pairwise comparison. Functional analysis showed that the DEGs were mainly involved in the processes of cell growth, muscle development, and cellular activities (junction, migration, assembly, differentiation, and proliferation). Many of DEGs were well known to be related to growth of skeletal muscle in black Muscovy duck, such as MyoG, FBXO1, MEF2A, and FoxN2. KEGG pathway analysis identified that the DEGs were significantly enriched in the pathways related to the focal adhesion, MAPK signaling pathway and regulation of the actin cytoskeleton. Some DEGs assigned to these pathways were potential candidate genes inducing the difference in muscle growth among the developmental stages, such as FAF1, RGS8, GRB10, SMYD3, and TNNI2. Our study identified several genes and pathways that may participate in the regulation of skeletal muscle growth in black Muscovy duck. These results should serve as an important resource revealing the molecular basis of muscle growth and development in duck.


Introduction
Skeletal muscle is the largest and most important tissue in animals, accounting for about 50% of body weight [1,2]. It participates in body movement, protection, and metabolic regulation [3]. Because meat yield directly determines the level of economic benefits, the study of skeletal muscle development is important in animal husbandry production. The growth and development of duck skeletal muscle is an essential economic trait, which is determined by genetic, and influenced by nutritional and environmental factors. The development of embryonic skeletal muscle is a tightly regulated process that is critically modulated by genes and related signaling pathways [4,5]. During prenatal and very early postnatal development, muscle growth in vertebrates depends on an increasing number of muscle fibers (hyperplasia) [6,7]. After postnatal growth, the increase in skeletal muscle is mainly due to muscle hypertrophy, accompanied by the proliferation of satellite cells and then the new myonuclei is incorporated into existing myofibers [8]. Black Muscovy duck, an excellent meat duck breed in China, has the advantages of fast growth and high lean meat rate. Most notably, few studies have investigated the mechanisms regulating growth patterns of skeletal muscle in black Muscovy duck. Since the growth of skeletal muscle is controlled by multiple genes, a more systematic understanding of the genes expressed at different growth stages in black Muscovy duck is needed.
Transcriptome sequencing (RNA-seq) is a technology that analyzing the transcripts by deep sequencing technology, which detects the whole transcriptome at the single nucleotide level [9]. It analyzes the structure and expression level of transcripts, which is an important method of gene expression and transcription analysis [10]. In recent years, RNA-seq has been widely used in study on livestock and poultry transcriptomes. Compared with other gene expression profiling methods, RNA-Seq has advantages in detecting mRNA expression in different tissues or at different developmental stages, which is helpful to reveal new genes and splicing variations [11], as well as the pathways [12]. Several studies at mRNA level have been performed in poultry birds, with the aims to investigate and analyze the genes and pathways that influence the growth and development of skeletal muscle. Xue et al. (2018) compared the mRNA expression of leg muscles of Jinghai Yellow Chickens at early growth stages by RNA-seq, a total of 978 differentially expressed genes (DEGs) were identified, and five significantly enriched pathways were found: EMC-receptor interaction, focal adhesion, tight junction, insulin signaling pathway, and regulation of the actin cytoskeleton [13]. Transcriptome analysis of breast muscle in Commercial Layers of Roman, White Broiler, and Daheng chickens showed that genes related to positive cell proliferation, growth, cell differentiation and developmental processes were more enriched, and the pathways, including ECM-receptor interaction, MAPK signaling pathway and focal adhesion, were the most enriched DEGs [14]. Xu et al. (2017) analyzed the gene expression profiles of Pekin ducks during incubation period, and the results showed that the DEGs related to cell division (M phase of mitotic cell cycle, cell division, mitosis, and mitotic prometaphase), and the pathways, including DNA replication, Cell cycle, Gap junction, were significantly enriched at hatched day 19 [15]. Zhu et al. (2017) analyzed potential candidate genes and signaling pathways related to the development of breast muscle during late-term embryonic to neonatal development, a total of 393 DEGs were identified, and the DEGs involved in different metabolism pathways (such as metabolic pathways, citrate cycle, glycine, serine, and threonine metabolism, sulfur metabolism, carbon metabolism, and pyruvate metabolism) [16].
The purpose of this study was to investigate the major DEGs and their expression pathways in skeletal muscle of black Muscovy duck by using RNA-Seq technology and bioinformatic tools. We found the genes and molecular mechanisms involved in this developmental process by carrying out a comprehensive analysis of genes with expression levels that reflected the growth pattern of breast and leg muscles in black Muscovy duck. Our findings are useful for understanding the mechanisms regulating the development of skeletal muscle and the pattern of duck growth.

Animal and Muscle Tissue Collection
One hundred eggs of black Muscovy duck were incubated in a standard incubator after disinfection. A total of 8 eggs were sampled from the day 17 (BE17), day 21 (BE21), day 27 (BE27), day 31 (BE31), and day 34 (BE34) of the incubation period, and the breast and leg muscle were separated for DNA and RNA extraction. DNA of muscle tissue was extracted according to the phenol chloroform protocol. According to the sequence of chromatin helix DNA binding protein 1 (CHD1) gene on sex chromosomes of duck, sex identification primers (gCHD) were used [17] and female embryos were selected as the research objects (Table 1), because the female Muscovy duck accounted for the majority of duck farms due to the reason of laying eggs, and the same gender can avoid the error of sequencing data. Besides, 6-month-old female ducks (BM6), raising under the same environmental conditions and free access to feed and water (Table 2), were slaughtered quickly to collect the breast muscle and leg muscle. Muscles were excised and immediately frozen in liquid nitrogen and stored at −80 • C until further use. Animal care, slaughter and experimental procedures were approved by Institutional Animal Care and Institutional Ethic Committee of Northwest A&F University (ethic code: #0326/2019).

RNA Extraction, Library Construction and Sequencing
The total RNA was extracted from the breast and leg muscle separately using QIAzol Lysis Reagent according to the manufacturer's instructions (QIAGEN, BerlinCity Germany). The RNA integrity was measured using a 2100 Bioanalyzer (Agilent Technologies, San JoseCity USA), the RNA purity and concentration were verified by agarose gel electrophoresis and Nanodrop 2000 (Thermo, Waltham, CA, USA).
Thirty-six sequencing libraries were construct by the TruSeq PE Cluster Kit v4-cBot-HS (Illumina HiSeq X Ten platform, San Diego, CA, USA) according to the manufacturer's instructions, and the indexed samples were clustered by the Illumina's cBot cluster generation system following the manufacturer's instructions. After cluster generation, the libraries were sequenced on the Illumina platform (Illumina HiSeq X Ten platform, San Diego, CA, USA), and a paired end read of 150 bp was generated.

Sequencing Analysis
After sequencing, raw reads were filtered: (1) Removing reads containing adapters or poly-N; (2) removing reads containing more than 10% of unknown nucleotides; and (3) removing low-quality reads containing more than 50% of low-quality (Q-value ≤ 10) bases. Besides, quality parameters for filtered data including Q30, GC content, and sequence-duplication level were used for data filtering. All downstream analyses were based on high-quality clean data. The filtered reads were mapped to the Anas platyrhynchos (Ap) genome sequence (https://www.ncbi.nlm.nih.gov/genome/?term=DUCK) and annotated transcripts (https://www.ncbi.nlm.nih.gov/assembly/GCF_003850225.1). Based on the Ap genome, further analyses and annotation have been carried out for a perfect matched reads or one mismatched reads. Then, the HISAT 2 tool software were used to map with the Ap genome.

Analysis of SNP/InDel
According to the results of the HISAT 2 comparison between the reads and the Ap genome sequence of each sample, the GATK software was used to find the single base mismatch between the sequenced sample and the Ap genome, and identify the potential single nucleotide polymorphism (SNP) site. The insertion-deletion (InDel) of the sample was also detected by the GATK software. According to the position of heterotopia in the Ap genome, the region where the mutation occurred (intergenic region, gene region or CDS region, etc.) and the effect of mutation (synonymous or non-synonymous mutation, etc.) were got by the SNPEFF software.

Prediction of Variable Splices
The results of HISAT 2 comparison were spliced by the StringTie software. The variable splicing types and corresponding expressions of each sample were obtained by the ASprofile software. The differentially expressed isoforms were estimated by the Cufflink.

Analysis of DEGs
To analyze DEGs, gene coverage and gene expression level were calculated. Gene coverage is the percentage of a gene covered by reads. The unigene expression was calculated and normalized to FPKM (fragments per kilobase per transcript per million mapped reads) based on FPKM (A) = cDNA Fragments/[Mapped Fragments (Millions) Transcript Length (kb)], where cDNA Fragments referred to the number of fragments compared to a transcript; Mapped Fragments (Millions) referred to the total number of fragments compared to a transcript, in 1 × 10 6 units; Script Length (kb) referred to the length of the transcript, in 1 × 10 3 bases units.
Differential expression was analyzed using the DESeq 2, and the false discovery rate (FDR) was calculated to calibrate the significant level and eliminate the influence of random fluctuations and errors. The unigenes with Fold change ≥ 2 and FDR < 0.01 were considered DEGs. The Benjamin-Hochberg correction method was used to correct the significant p-value of the original hypothesis test, and FDR was used as the key indicator for screening DEGs.

RNA-seq Validation by qPCR
In order to confirm the reliability of RNA-seq results, 26 representative responsive genes of breast and leg muscle in black Muscovy duck were selected for qPCR, respectively. Duck β-actin (accession number: NC_040060.1) was measured as the housekeeping gene. RNA with an A260/280 ratio between 1.9 and 2.1, and A260/230 ratio > 2.0 were used for cDNA synthesis. The list of primers were described in Table 1. Each reaction mixture contained 5 µL of TransStart Tip Green qPCR SuperMix (Transgen, Beijing, China), 0.8 µL of cDNA (400 ng/µL), 0.3 µL of each primer (10 µM) and 3.6 µL ddH2O in 10 µL final volume using EcoRT48 (OSA, London, UK). The optimal reaction procedure included 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s, 60 • C for 30 s, then 95 • C for 15 s, 55 • C for 15 s, 95 • C for 15 s. To derive the relative expression value, 2 − Ct method was adopted. The results were expressed as mean ± SD of at least three independent biological replicates. Statistical analyses of the data were conducted in SPSS 20.0. Significant differences in duck muscle gene expression in different growth periods were determined by one-way ANOVA followed by Tukey's test and Duncan's test.

Transcriptome Profiles
In this study, a total of 36 libraries were established from the breast and leg muscles of duck in BE17, BE21, BE27, BE31, BE34, and BM6. The gene expression of black Muscovy duck skeletal muscle transcriptome at different growth stages was systematically analyzed by high-throughput RNA sequencing. All samples had an RNA integrity number (RIN) > 7.5, and a RNA concentration ≥ 125.2 ng/µL (Table S1). After quality control, clean reads of samples ranging from 19,700,766 to 33,105,790 (24,559,860 on average), GC contents of the samples were between 49.19% and 56.80%, and ≥Q30 (%) were 91.36% to 93.80% (Table 3). In total reads of the samples, over 26,201,608 high-quality reads per sample were mapped to the Ap genome and used for gene expression analysis, where "Uniq Mapped Reads" ranging from 23,077,005 (51.13%) to 30,289,226 (69.98%) and "Multiple Map Reads" ranging from 1,297,423 (2.92%) to 11,929,431 (23.55%). Besides, the number of Reads aligned to the positive strand in the Ap genome were 9,228,720 to 20,801,712 and the number of reads aligned to the negative strand in the Ap genome were between 13,086,584 and 21,955,456 (Table S2).

Annotation and Classification of SNP/InDel
The number of SNP loci, the proportion of transition type and transversion type, as well as the ratio of heterozygous SNP loci from each sample were counted, the results were shown in Table 4. There were 23,381 to 634,028 SNPs in skeletal muscle of black Muscovy duck at different growth stages, where the total numbers of SNPs in the genic region were 20,572 to 574,484, the numbers of SNPs between genes were 2809 to 59,544. Besides, the percentage that the transition-type SNP accounts for all SNP locis were over 71.12%, the percentage that the transversion-type SNP loci accounts for all SNP sites were between 23.66% and 28.88%, and the percentage that the heterozygous SNPs account for all SNPs were 4.27% to 48.09%. The annotation results of SNP and InDel were shown in Figure 1. The annotations of SNP were mainly distributed in "INTRON", "SYNONYMOUS_CODING", and "UTR_3_PRIME". The annotations of InDel were mainly distributed in "INTRON", "UTR_3_PRIME", and "DOWNSTREAM".

Prediction of Alternative Splice (AS)
The chromosomal position of each transcript was obtained by aligning the sequence to the Ap genome. Twelve different splice patterns in the transcriptome data of duck skeletal muscle were detected:

Analysis of DEGs
The relative expression of gene was normalized as fragments per kilobase of exon model per million mapped reads (FPKM), which was proportional to the number of cDNA fragments originated by gene transcription. Statistical analysis of DEGs in breast and leg muscle of black Muscovy ducks at different stages were shown in Table 5. In breast muscle, differential gene expression analysis of BE17B_vs_BE21B showed that 410 genes were significantly differentially expressed (Fold change ≥ 2 and FDR < 0.01 at p < 0.05, the same below), including 218 up-regulated and 192 down-regulated genes. Moreover, there were 1958 significantly expressed genes from BE17B_vs_BE27B, among which 1162 were up-regulated genes and 796 were down-regulated genes. A number of 1517 DEGs were detected from BE17B_vs_BE31B, the number of up-regulated genes were 925 and the down-regulated genes was 592 respectively. There were 1460 DEGs from BE17B_vs_BE34B, including 852 up-regulated genes and 608 down-regulated genes. Besides, 5377 DEGs were found in BE17B_vs_BM6B, of which 2580 were up-regulated genes and 2797 were down regulated genes ( Figure S1).
In leg muscle, the results showed that there were 655 DEGs, including 371 up-regulated genes and 284 down-regulated genes from BE17L_vs_BE21L. Moreover, 2866 DEGs were discovered in BE17L_vs_BE27L, among which 1606 genes were up-regulated and 1260 genes were down-regulated. From BE17L_vs_BE31L, there were 4413 significantly expressed genes, among which 2440 were up-regulated genes and 1973 were down-regulated genes. 4326 DEGs were detected from BE17B_vs_BE34B, out of which 2374 up-regulated and 1952 down-regulated genes were identified as significantly differentially expressed. Besides, 4560 DEGs were discovered in BE17L_vs_BM6L, and 2303 DEGs were up-regulated genes and 2257 DEGs were down-regulated genes ( Figure S2).

Analysis of GO Annotation and KEGG Pathway
Functional annotations of DEGs in the database were performed, and the number of annotated genes were shown in Table S4. The total number of DEGs annotated were between 17 and 5253, where the COG were 3 to 1775, the GO were 14 to 4191 and the KEGG were 5 to 3542. Besides, the number of KOG ranged from 10 to 3788, the NR ranged from 17 to 5231, the Pfam ranged from 14 to 4816. Moreover, the Swiss-Prot were between 12 and 3788, the eggNOG were between 15 and 5118.
To gain valuable insight into the molecular functions of the genes potentially associated with muscle development, the identified DEGs were categorize into three functional groups depending on gene ontology: Biological Process, Molecular Function and Cellular Component. In breast muscle, GO analysis was performed on the common DEGs in BE17B_vs_BE21B indicated that biological processes, including regulation of calcium ion import, regulation of muscle filament sliding speed, dorsal root ganglion development, and negative regulation of fibroblast growth factor receptor signaling pathway were significantly enriched. The common DEGs in BE17B_vs_BE27B were primarily enriched in biological processes of actin binding, positive regulation of myoblast proliferation, cell division, positive regulation of cell proliferation and regulation of actin cytoskeleton organization. Besides, most genes that were involved in the biological processes of regulation of muscle filament sliding, skeletal muscle fiber development, positive regulation of myoblast differentiation and muscle cell cellular homeostasis were mainly enriched in BE31B compared to that of BE17B. The terms of regulation of cell cycle, positive regulation of fibroblast proliferation and positive regulation of substrate-dependent cell migration, cell attachment to substrate were enriched in BE34B compared to BE17B. Additionally, several fundamental biological processes were found to be notably enriched in BE17B_vs_BM6B, such as translation, immune response, regulation of cell size and regulation of cell growth (Figures 3-5).  In leg muscle, DEGs were mainly enriched in biological processes of cell proliferation, skeletal muscle fiber development, skeletal muscle tissue development, regulation of muscle filament sliding, positive regulation of cell division in BE17L_vs_BE21L, and regulation of transcription involved in cell fate commitment, calcium-mediated signaling, glucose transport in BE17L_vs_B27L, respectively. Compared to BE17L, most genes that were involved in the biological processes of cell maturation, embryonic limb morphogenesis, chordate embryonic development and immune response were mainly enriched in BE31L. DEGs in BE17L_vs_BE34L and BE17L_vs_BM6L were primarily enriched in biological processes of embryonic hindlimb morphogenesis, positive regulation of protein process, regulation of actin cytoskeleton organization, positive regulation of cell proliferation and Wnt receptor catabolic process, embryonic hindlimb morphogenesis, glucose transport, protein folding, and immune response, respectively (Figures 6-8). In the comparison of breast and leg muscle, there were some terms mainly reached in biological processes of myoblast migration involved in skeletal muscle regeneration, positive regulation of glucocorticoid receptor signaling pathways, regulation of multicellular organism growth in BE17B_vs_BE17L, and positive regulation of cellular process, metabolic process, fibroblast migration in BE21B_vs_BE21L, respectively. In BE27B_vs_BE27L and BE31B_vs_B31L, DEGs were primarily enriched in biological processes, including positive regulation of MHC class I biosynthetic process, immune response, metabotic process and regulation of cell shape, embryonic organ development, translation, and muscle structure morphogenesis, respectively. DEGs in BE34B_vs_BE34L were primarily enriched in biological processes of skeletal muscle cell differentiation, skeletal muscle fiber adaptation, myotube differentiation involved in skeletal muscle regeneration, positive regulation of skeletal muscle tissue regeneration. Additionally, several biological processes were found to be notably enriched in BM6B_vs_BM6L, such as muscle structure development, embryonic skeletal joint morphogenesis, myoblast fate commitment, and negative regulation of skeletal muscle tissue development (Figures 9-11, Table 6).   Besides, Cluster of Orthologous Groups (COG) database was constructed using phylogeny of bacteria, algae, and eukaryotic. The products of genes could be orthologously classified using the COG database ( Figures S7-S9). KEGG analysis of DEGs revealed that Focal adhesion, Regulation of actin cytoskeleton, MAPK signaling pathway, Wnt signaling pathway, ECM-receptor interaction, neuroactive ligand-receptor interaction, purine metabolism, calcium signaling pathway, endocytosis, ErbB signaling pathway, glucagon signaling pathway, RIG-I-like receptor signaling pathway, Insulin signaling pathway, cell cycle, apoptosis, oxidative phosphorylation phosphatidylinositol signaling system, cell adhesion molecules (CAMs), biosynthesis of amino acids and Ribosome were the most enriched for the DEGs at different developmental stages (Figures 12-20, Table 7).    adhesion, Regulation of actin cytoskeleton, MAPK signaling pathway, Wnt signaling pathway, ECM-receptor interaction, neuroactive ligand-receptor interaction, purine metabolism, calcium signaling pathway, endocytosis, ErbB signaling pathway, glucagon signaling pathway, RIG-I-like receptor signaling pathway, Insulin signaling pathway, cell cycle, apoptosis, oxidative phosphorylation phosphatidylinositol signaling system, cell adhesion molecules (CAMs), biosynthesis of amino acids and Ribosome were the most enriched for the DEGs at different developmental stages (Figures 12-20, Table 7).  The bigger the Rich factor is, the more significant the pathway is. The color of circle represented q value which is adjusted p value by multiple hypothesis testing. Thus, the smaller the q value is, the more significant the pathway is; the circle size represented the number of differentially expressed genes annotated with the pathway, the bigger circle size is, the higher number of genes is. The same below.
pathways, (b) was the ration of genes in the pathway with all genes in all pathways. The bigger the Rich factor is, the more significant the pathway is. The color of circle represented q value which is adjusted p value by multiple hypothesis testing. Thus, the smaller the q value is, the more significant the pathway is; the circle size represented the number of differentially expressed genes annotated with the pathway, the bigger circle size is, the higher number of genes is. The same below.

qPCR Analysis
To confirm RNA-seq data, 26 genes of duck breast and leg muscle obtained at different growth stages were selected for qPCR analysis, respectively. The expression tendency of these genes agreed well with RNA-seq results, which suggested that RNA-seq results were pretty reliable ( Figure 21).

qPCR Analysis
To confirm RNA-seq data, 26 genes of duck breast and leg muscle obtained at different growth stages were selected for qPCR analysis, respectively. The expression tendency of these genes agreed well with RNA-seq results, which suggested that RNA-seq results were pretty reliable ( Figure 21).  Figure 21. qPCR verification of DEGs. "*" was considered significant difference (p < 0.05); "**" was considered extremely significant difference (p > 0.01). Figure 21. qPCR verification of DEGs. "*" was considered significant difference (p < 0.05); "**" was considered extremely significant difference (p > 0.01).

Discussion
Duck growth performance, mainly skeletal muscle growth, provides direct economic benefits to the poultry industry. However, the underlying genetic mechanisms remain unclear. The purpose of this study was to identify candidate genes associated with the growth of duck skeletal muscle and investigate their potential mechanisms. Transcriptome analysis is an efficient and fast tool that has been widely used in animal husbandry. Comparative transcriptome analyses of tissues at different developmental stages provide valuable insights into the question of how regulatory gene networks control specific biological processes [18,19]. In this study, we focus on the skeletal muscle tissue of black Muscovy duck, because it is the prime part of carcass, and the yield of breast and leg meat is one of practical importance to the profitability of production [20]. Besides, muscle fiber numbers increase during embryonic development in poultry, after which myofiber numbers stop to increase, while myofiber volume still increases [8]. Therefore, we studied the gene expression patterns and network pathways at different developmental stages to further understand the molecular mechanism of duck skeletal muscle development. Annotation of the sequence data using the duck genome as a reference revealed expression of 39,401,532 to 66,211,580 genes in the duck skeletal muscle transcriptome, and an average of 49,119,720 high-quality reads per sample were mapped to the Ap genome.

SNP/InDel Analysis
RNA-seq measures not only gene expression, but also structural variations such as fusion transcriptions or mutations. SNP, a type of genome polymorphisms, only involves the variation of a single base. Because SNPs are abundant and have greater stability over generations, genotyping more accurately and easily automates the genotyping processes [21]. Some studies have reported a large number of SNPs associated with animal weight and muscle development traits. These candidate genes could be used as molecular markers in early marker-assisted selection in animal breeding programs [22,23]. InDels are a major class of genomic variation, which are primarily detected from DNA-seq data [24,25]. In order to explore the gene SNPs and InDels related to muscle growth and development, gene levels of the breast and leg muscle tissues were analyzed by high throughput RNA sequencing technology. There were 23,381 to 634,028 SNPs in skeletal muscles at different growth stages, where total numbers of SNPs in the genic region were 20,572 to 574,484, total numbers of SNPs between genes were 2809 to 59,544. The most common change was G/A and C/T, followed by A/G and T/C. The annotations of SNP were mainly distributed in "INTRON" (2,898,440 for breast muscle and 2,103,619 for leg muscle), "SYNONYMOUS_CODING" (1,319,429 for breast muscle and 1,395,604 for leg muscle) and "UTR_3_PRIME" (1,031,194 for breast muscle and 1,050,568 for leg muscle). The annotations of InDels were mainly distributed in "INTRON" (184,929 for breast muscle and 123,569 for leg muscle), "UTR_3_PRIME" (143,626 for breast muscle and 140,766 for leg muscle) and "DOWNSTREAM" (35,847 for breast muscle and 31,447 for leg muscle), indicating that the skeletal muscle SNPs and InDels of black Muscovy ducks at different growth stages were mainly in "INTRON", "SYNONYMOUS_CODING", "UTR_3_PRIME", and "DOWNSTREAM". Inclusion of SNPs and InDels (many of which may be located in genes) in annotated genes will provide a large number of gene centric markers, which will add detailed information for the genetic loci for specific phenotypic traits.

Prediction of AS
AS is a ubiquitous in most eukaryotic genomes, which is a mechanism for organisms to increase their protein pool and regulate physiological and developmental processes/pathways [26][27][28]. Precursor mRNA of AS plays an important role in the regulation of gene expression in higher eukaryotes. Multiple mRNAs can be derived from a single pre mRNA to produce proteins with different functions, which indicates that AS is an important mechanism for regulating life [29]. Isoform expression patterns may provide unique insights into the skeletal muscle transcriptome since it is likely that unique transcript splice variants may play essential roles during development [30]. Due to the lack of detailed full-length cDNA data and high-quality genome annotation, there are few studies about AS in ducks or other birds [31,32]. Yin et al. (2019) found that a total of 199,993 full-length transcripts were obtained from 8 duck tissues by using transcriptome sequencing, and 35,031 AS events were accurately predicted from 3346 genes, which is very useful for the functional research of other birds [33]. In our study, an average of 37,267 AS events were counted from each sample, and the number of AS in each sample were mainly concentrated in TSS and TTS, which indicated that the AS of duck skeletal muscle growth and development mainly existed in the first exon and the last exon.

DEGs Analyzed at All Time Points
The differential expression of growth related genes is considered to be the main cause of genetic variation during duck growth and development, indicating that the regulatory mechanism of growth has changed [13]. In our study, at different time points, the skeletal muscle of black Muscovy duck had 19 to 5377 DEGs. It was worth noting that there were only 19 DEGs in BE34B_vs_BE34L. It is speculated that the genes that regulate the proliferation and differentiation of duck breast and leg muscle were basically the same at this time. It is well known that these processes involved by the DEGs at different developmental stages are essential for maintaining animal muscle growth. Therefore, these genes may play an important role in growth and development of duck skeletal muscle. Further functional studies with these DEGs were warranted to identify key genes influencing muscle growth and development of duck. Besides, some DEGs were found to be closely related to growth and development of skeletal muscle in black Muscovy duck, including transcription factors such as MyoG, FBXO1, MEF2A, and FoxN2, as well as a series of genes related to muscle growth axis, such as FAF1, RGS8, GRB10, SMYD3, and TNNI2. Moreover, some studies have shown that these regulatory transcription factors interacted with each other in regulating animal muscle growth.
Fas-associated factor 1 (FAF1) is involved in diverse bio-chemical processes including cell death, inflammation, cell proliferation, and proteostasis [34]. FAF1 mediates caspase-8 activation via both intrinsic and extrinsic pathways. It suppresses NF-κB activation by interrupting IκB kinase (IKK) complex assembly, and promotes Fas-induced apoptosis [35,36]. FAF1 also inhibits cell cycle by negatively regulating Aurora-A. Moreover, FAF1 regulates the polyubiquitinated protein and valosin containing protein, and inhibits the degradation of ubiquitin dependent proteins [37,38]. Studies have shown that FAF1 antagonized Wnt signal transduction by promoting the degradation of β-Catenin [39]. Wnts regulate cell proliferation and differentiation, and control many biological processes including embryonic development [40,41]. FAF1 regulates Wnt/β-Catenin dependent gene expression in C2C12 myoblasts, including genes involved in osteoblast differentiation [39]. FAF1 is necessary for early embryogenesis. Adham et al. (2008) found that FAF1 gene targeted mice showed embryonic lethality at the two cell stage [42]. Ryu et al. (1999) found Human FAF1 mRNA is abundantly expressed in testis, skeletal muscle, and heart [43]. Fröhlich et al. (1998) identified the avian FAF1 homologue (qFAF) in the pluripotent cells from E0 quail embryos during mesoderm induction in vitro by using mRNA differential display technique, which can be used as the induction gene of fibroblast growth factor (FGF) [44]. Therefore, FAF1 may play an important role in the development of duck skeletal muscle.
Regulator of G protein signaling proteins (RGS) are negative modulators of many G-Protein Coupled Receptor (GPCR) signaling pathways [45]. RGS protein directly binds to the Gα subunit of GTP of activated heterotrimer G protein, which increases the rate of GTP hydrolysis [46]. By this mechanism, RGS proteins rapidly dampen GPCR signal transduction at the level of the active G protein subunits [47]. Regulator of G protein signaling 8 (RGS8) belongs to the R4 subfamily of RGS proteins, and modulates the functioning of G-proteins by activating the intrinsic GTPase activity of the a subunits [48,49]. RGS is involved in many intracellular processes mediated by G protein signal transduction pathway, including cell proliferation, cell differentiation, plasma membrane transport, cell movement, and embryonic development [50]. In our study, RGS8 may regulate G protein by activating GTPase activity of a subunit and participate in development of duck skeletal muscle.
Growth factor receptor-bound protein 10 (GRB10) regulates phosphorylation and activation of the mTORC1 protein, which is a central regulator of cellular metabolism, growth and survival [51,52]. GRB10 is also involved in the regulation of glucose metabolism during fetal and postnatal period [53]. It is also one of the participants in ensuring the metabolic health of fetus to adult. GRB10 binds to Gab1 and participates in the regulation of cell mitosis [54]. Besides, GRB10 is an adaptor protein, which interacts with a number of receptor tyrosine kinases and signaling molecules [55]. GRB10 associates with a variety of growth factors at the cell surface, such as the insulin growth factor receptor (IGFR), and with intracellular protein kinases like Raf1 and MEK1 [56,57]. It is a negative regulator of IGFIR-dependent cell proliferation, and plays a negative regulatory role in the MAPK signaling pathway [58,59]. Notably, IGF and MAPK signaling pathway play an important role in skeletal muscle development. Studies have shown that GRB10 is abundant in brain, fat, muscle and heart of mice [60]. Mutation of GRB10 induces muscle hypertrophy in mice [61], and the embryos and placentas of mice lacking GRB10 were overgrown, such that mutant mice were 30% larger than normal at birth [55]. It was suggested that GRB10 may be very important for skeletal muscle development in duck.
Muscle fibers are composed of myofibrils, and members of the SMYD family play critical roles in myofibril assembly of skeletal and cardiac muscle during development [62]. SET and MYND domain-containing protein 3 (SMYD3) is widely distributed in eukaryotes and participates in epigenetic transcription regulation, development and cell proliferation [63]. SMYD3 is also necessary for regulating skeletal muscle and myocardial development [64]. SMYD3 controls skeletal muscle development and maintenance through transcriptional regulation [65]. Fujii et al. (2011) found that in zebrafish embryos, SMYD3-knockdown led to abnormal expression of myogenic markers including MyoD, which indicated that SMYD3 may play a role in muscle development [64]. In addition, SMYD3 is involved in regulating skeletal muscle atrophy. During the process of dexamethasone-induced skeletal muscle atrophy, the mRNA level of SMYD3 was significantly increased [66]. The differential expression of SMYD3 in duck may be related to greater myogenic potential. Troponin I2 (TNNI2), a muscle growth marker gene [67,68], encodes a subunit of troponin complex [69], which are expressed under muscle type-specific and developmental regulations [70]. Troponin complex is a group of muscle proteins, which is part of the contraction device of rapid skeletal muscle contraction [71]. In the absence of Ca 2+ , TNNI2 prevents muscle contraction by binding actin and tropomyosin [72]. Besides, Yoshimoto et al. (2020) found that expression of TNNI2 was gradually increased as muscle regeneration proceeds, indicating that TNNI2 were excellent indicator to assess myofiber maturity [73]. TNNI2 encodes an isoform of TnI specific to fast-twitch myofibers and troponin I (TnI) mutation with abnormal skeletal muscle structure leads to wing and limb defects in Drosophila [74]. Besides, the expression of slow-twitch TnI is stronger in soleus muscle, while the expression of fast-twitch TnI is stronger in tibial anterior muscle and extensor digitorum longus in neonatal mice [75]. Therefore, TNNI2 may be an interesting candidate gene to explain the phenotypic differences of skeletal muscle development in duck.
The RNA-Seq data were confirmed to be reliable by qPCR. These identified DEGs included many genes significantly related to muscle development.

GO and KEGG Pathway
GO database is a structured standard biological annotation system. It aims to establish a standard system for genes and their products, which is suitable for all species. GO annotation system includes three main branches: Biological Process, Molecular Function, and Cellular Component. In this study, several key DEGs (such as MyoG, FBXO1, MEF2A, FAF1, RGS8, GRB10, SMYD3, and TNNI2) involved in muscle development can be identified in enriched GO terms, which providing a theoretical basis for the skeletal muscle development of black Muscovy duck at different growth stages. The GO analysis of these DEGs showed that the biological processes related to cell growth, muscle development, and cellular activities (such as junction, migration, assembly, differentiation, and proliferation), muscle contraction, as well as glycogen metabolic and biosynthetic processes, were regulated differently at each developmental stage, indicating that the DEGs played an important role in regulating duck skeletal muscle.
In organisms, different genes interact to perform biological functions. Annotation and analysis of pathways are helpful to further understand the functions of genes. KEGG is a database for systematic analysis of gene function and genome information, which can be used as a whole network to study gene and expression information. Duck muscle growth is a complex process influenced by multiple genes and controlled by multiple pathways. In our KEGG analysis, main pathways related to muscle growth were identified, namely, focal adhesion, ECM-receptor interaction, MAPK signaling pathway, neuroactive ligand-receptor interaction, endocytosis, oxidative phosphorylation, ribosome, tight junction, insulin signaling pathway, and regulation of the actin cytoskeleton, of which the focal adhesion, MAPK signaling pathway and regulation of the actin cytoskeleton were the most significantly enriched (Figures S10-S12).
Focal adhesions (FAs) are integrin-containing, multi-protein assemblies crossing the plasma membrane that connect the cellular cytoskeleton to surrounding extracellular matrix. They play a key roles in adhesion and cell signal transduction, and are the main regulators of epithelial homeostasis, tissue response to injury and tumorigenesis [76]. Focal adhesion kinase (FAK) is a multifunctional molecule with the ability to regulate muscle formation, hypertrophy and glucose metabolism [77]. The mitogen-activated protein kinase (MAPK) signaling pathway is a phosphorylation kinase signaling cascade that regulates many cell processes, such as cell division, differentiation, and release of inflammatory mediators [78]. It is well known that MAPK signaling pathway induces protein synthesis and promotes skeletal muscle growth or hypertrophy [79]. The actin cytoskeleton comprises a scaffold of polymeric actin filaments that are assembled and disassembled to organize cell architecture and direct many cell processes. The actin-related proteins control actin filaments reorganization, resulting in significant changes in actin cytoskeleton structure, thus regulating cell processes that affect mitosis, cytokinesis, endocytosis, and cell migration [80,81]. At focal adhesion, the extracellular domain of transmembrane integrin, including α and β subunits, is connected to the extracellular matrix (ECM). The intracellular tail of integrin binds to connexin, and connexin binds to the actin cytoskeleton to perform their related functions [82]. Therefore, the results of this study will allow us to predict the function of new genes and to explore candidate genes that might play a role in muscle growth and development process in duck. Besides, the annotation findings of the current investigation can be used as the selection marker for the genes related to the skeletal muscle growth and might be helpful in better understanding of genetic mechanisms associated with the growth of duck skeletal muscle.

Conclusions
In this study, transcriptome data was generated by RNA-Seq technology, which will help to further understand the molecular sequences and functions of skeletal muscle growth related genes of black Muscovy duck at different stages. There were differences in the expression of genes at different growth stages, including SNPs, InDels, AS, highly expressed genes and pathways. These findings will provide valuable resources for the biological researches of skeletal muscle growth related genes in black Muscovy duck and may also provide clues for understanding the molecular mechanisms in other poultry and mammalian species.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/10/1228/s1, Table S1 The Concentration and RIN value of sample RNA; Table S2 Statistics of sequence comparison between sample sequencing data and Anas platyrhynchos genome; Table S3 The five up-and down-regulated genes among the comparison samples; Table S4 Number

Conflicts of Interest:
The authors declare no conflict of interest.