Transcriptome Sequencing Reveals Pathways Related to Proliferation and Differentiation of Shitou Goose Myoblasts

Simple Summary The Shitou goose, the largest cultivated goose breed in China, has high research value in meat traits. Muscle development is regulated by genes related to myoblast proliferation and differentiation. In this study, the mRNA and lncRNA expression profiles of Shitou goose myoblast in proliferation and differentiation were constructed by Illumina sequencing. A total of 1664 differentially expressed (DE) mRNAs and 244 DE-lncRNAs were identified between the two periods. Functional annotation showed that the DE-mRNAs and DE-lncRNAs were mainly enriched in the Wnt signaling pathway. These results provide new insights into the mechanism of muscle growth and development in large goose breeds. Abstract Chinese Shitou goose is a type of large goose with high meat yield. Understanding the genetic regulation of muscle development in Shitou goose would be beneficial to improve the meat production traits of geese. Muscle development is regulated by genes related to myoblast proliferation and differentiation. In this study, the RNA-seq method was used to construct the mRNA and lncRNA expression profiles of Shitou goose myoblasts and myotubes. A total of 1664 differentially expressed (DE) mRNAs and 244 DE-lncRNAs were identified. The alternative mRNA splicing in proliferation and differentiation stages was also analyzed. Notably, pathways enriched in DE-mRNAs, DE-splicing transcripts, and DE-lncRNAs all point to the Wnt signaling pathway, indicating that the Wnt signaling is a key regulatory pathway of muscle development in Shitou goose. We also constructed the interactive network of DE-lncRNAs and DE-mRNAs and revealed some key genes of lncRNAs regulating the proliferation and differentiation of myoblasts. These results provide new insights for the study of the muscle development of the Shitou goose.


Introduction
In animal husbandry, meat is the main product [1], and most of the selection and breeding of meat-raised animals is focused on improving the percentage of retail cuts of the animals. As the main component of the animal body, skeletal muscle accounts for about 35% of body weight [2], and the development and growth of muscle is the key factor to provide enough meat for humans. Understanding the growth and development of skeletal muscle is important to improve the percentage of retail cuts of livestock. The growth and development of skeletal muscle is an extremely complex process, including the directed differentiation of progenitor cells, the proliferation and differentiation of myoblasts, the fusion of myocytes, and finally the formation of multinucleated muscle fibers with contractile function [3]. There have been many reports on the muscle development of legs and removed the bones, skin, and fascia. The leg muscles were minced with curved scissors in the culture medium and then digested with trypsin at 37 • C for 20-30 min. A DMEM medium containing 20% fetal bovine serum (FBS) (Gibco, Grand Island, USA) was added to terminate the digestion. The mixed medium was centrifuged at 1200 g for 5 min, and then the supernatant was discarded. A DMEM medium with 15% FBS was added to resuspend the cells and then put into a 10 cm diameter petri dish. Serial plating was performed to enrich myoblasts and eliminate fibroblasts. We divided the mixed primary myoblasts into 6 cultures. Three cultures were used to collect proliferating myoblasts (GM), and the other three cultures were induced into differentiation and then used to collect differentiated myotubes (DM). The differentiation of myoblasts was induced by replacing 20% FBS with 2% horse serum (HS) (Biosharp, Hefei, China).

RNA Extraction, cDNA Synthesis, Fluorescence Quantitative PCR
Total RNA was extracted from goose myoblasts using MagZol Reagent (Magen, Guangzhou, China). The Evo M-MLV RT Kit (Agbio, Guangzhou, China) was used for reverse transcription to synthesize cDNA. MonAmp™ ChemoHS qPCR Mix (Monad, Guangzhou, China) was used to perform qRT-PCR. All reactions were set to three replicates. The 2 −∆∆Ct method was used to measure gene expression with GAPDH as the reference gene [20,21].

Library Construction and Illumine Sequencing
The proliferative myoblasts of Shitou geese were cultured to 70% density in the petri dish, then washed with PBS three times, and directly treated with Trizol reagent (Invitrogen, Carlsbad, CA, USA). The differentiated myoblast needs to be cultured with a DMEM of 2% horse serum when the density reaches 80-90%. After three days of differentiation, the samples were collected with Trizol. Total RNA was isolated and purified using Trizol reagent following the manufacturer's procedure. The RNA amount and purity of each sample were quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed by Agilent 2100 with RIN number >7.0. Approximately 5 ug of total RNA was used to deplete ribosomal RNA according to the manuscript of the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, CA, USA). Paired-end sequencing was performed on an Illumina Novaseq™ 6000 (LC Bio, Hangzhou, China) following the vendor's recommended protocol.

Bioinformatics Analysis
RNA-seq data from six geese have been deposited in NCBI Gene Expression Omnibus (GEO) data repository under accession number GSE213147. Cutadapt (https://cutadapt. readthedocs.io/en/stable/, accessed on 3 December 2020) and in-house perl scripts were used to remove reads containing undetermined bases, low-quality bases, or those that contained adaptor contamination. FastQC (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/, accessed on 14 December 2020) was used for quality control. Bowtie2 [22] was used to map reads to the genome of Anser cygnoides (assembly GooseV1.0 (http://asia. ensembl.org/Anser_cygnoides/Info/Index/, accessed on 26 December 2020)). The mapped reads of each sample were assembled using StringTie [23]. Then, all transcriptomes from Shitou goose were merged to reconstruct a comprehensive transcriptome using perlscript and gffcompare. After the final transcriptome was generated, StringTie and Ballgown [24] were used to estimate the expression levels of all transcripts. Differential expression analysis was performed using the edgeR according to the edgeR User's Guide [25]. The differentially expressed mRNAs and lncRNAs were selected with log2 (fold change) > 1 or log2 (fold change) < −1 and with statistical significance (adjusted p-value < 0.05) by R package edgeR. See the supplementary figure for the RNA sequence data processing flow chart (Supplementary Figure S1).
The lncRNA prediction software CPC (coding potential calculator) [26] and CNCI (coding non-coding index) [27] were used to predict the lncRNA transcripts larger than bp. All transcripts with CPC score < −1 and CNCI score < 0 were removed. The codes for Gene differential splicing analysis can be referred to Supplementary File S1.
The alternative splicing transcripts in RNA-seq data were analyzed by rMATS [28]. Differential splicing genes were selected with log2 (fold change) > 1 and abs(lncLeveldiff) > 0.1. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of the enriched genes was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/tools.jsp/, accessed on 22 March 2021) and R Studio. Genes expressed in our samples were used as the background. The R Studio code can refer to Supplementary File S2.

Construction of DE-lncRNAs-DEGs Interaction Network
We searched for all of the DEGs 100 kb upstream and downstream of differentially expressed lncRNA by the UCSC Genome Browser database [29]. These genes that are genomically adjacent and coexpressed in the expression pattern are likely to be the cistarget genes of the lncRNA. The lncRNA-gene Interaction was constructed by Cytoscape software (https://cytoscape.org/, accessed on 14 June 2021). Functional analysis of the cis-target genes was conducted by BLAST2GO [30]. Differences at p < 0.05 were regarded as statistically significant.

Statistical Analysis of qPCR Results
For the qPCR experiment, each sample was analyzed based on results that were repeated in triplicate and analyzed by Graphpad 8.0 (GraphPad Software, San Diego, CA, USA). The statistical significance of differences between the two groups was determined by a standard Student's t-test. In all cases, differences at p < 0.05 were regarded as statistically significant (Significance is coded as * for p < 0.05 and ** for p < 0.01).

Isolation and Transcriptome Sequencing of Shitou Goose Myoblasts
To identify the genes involved in the development of Shitou goose muscle during the proliferation and differentiation of myoblasts, we used the Illumina HiSeq platform to sequence the transcriptome of cells in the proliferation (GM) and differentiation (DM) stages ( Figure 1A). The box diagram ( Figure 1B) showed that the sequencing quality of the six sequencing samples is relatively consistent. The percentage of unique Mapped reads of the samples is greater than 70%, the Q30 value is greater than 98.4%, and the GC base content is between 47-48% (Table 1). There were 1664 differentially expressed genes (DEGs) between the proliferation and differentiation stages of Shitou goose myoblasts, including 685 up-regulated DEGs and 845 down-regulated DEGs ( Figure 1C). The top ten differentially expressed genes are listed in Table 2.

Functional Analysis of DEGs between Shitou Goose Myoblasts and Myotubes
We then performed GO and KEGG analysis on the DEGs to identify important genes and pathways during goose muscle development. GO functional enrichment analysis found that down-regulated DEGs were significantly enriched in muscle development, such as positive regulation of cell proliferation, positive regulation of cytosolic calcium,

Functional Analysis of DEGs between Shitou Goose Myoblasts and Myotubes
We then performed GO and KEGG analysis on the DEGs to identify important genes and pathways during goose muscle development. GO functional enrichment analysis found that down-regulated DEGs were significantly enriched in muscle development, such as positive regulation of cell proliferation, positive regulation of cytosolic calcium, cell adhesion molecule binding, and actin binding (Figure 2A). The up-regulated DEGs were significantly enriched in cellular calcium ion homeostasis, skeletal muscle tissue development, and calcium ion transmembrane transport ( Figure 2B). KEGG pathway enrichment analysis found that the DEGs were significantly enriched in the pathways Animals 2022, 12, 2956 6 of 14 related to cell growth and metabolism. The down-regulated DEGs enriched pathways include the PI3K-Akt signaling pathway, the JAK-STAT signaling pathway, and the ECMreceptor interaction ( Figure 2C). The up-regulated DEGs enriched pathways include the calcium signaling pathway, cardiac muscle contraction, and hypertrophic cardiomyopathy ( Figure 2D). calcium signaling pathway, cardiac muscle contraction, and hypertrophic cardiomyopathy ( Figure 2D).
Next, Gene set enrichment analysis (GSEA) of the RNA-seq data demonstrated that the differentiation of goose myoblast led to the positive enrichment of muscle cell differentiation, muscle organ development, and cell adhesion ( Figure 3A). Genes related to cell cycle and DNA duplication were enriched in proliferative cells ( Figure 3B). We randomly selected some DEGs for qPCR verification. Muscle differentiation-related genes such as MYH15 and MyoM are highly expressed in the differentiation stage. Cell cycle-related genes, such as CDKN1A and CCNB1, are significant differences in expression between the two stages ( Figure 3C). These results are consistent with the RNA sequencing results (Figure 3D).  Next, Gene set enrichment analysis (GSEA) of the RNA-seq data demonstrated that the differentiation of goose myoblast led to the positive enrichment of muscle cell differentiation, muscle organ development, and cell adhesion ( Figure 3A). Genes related to cell cycle and DNA duplication were enriched in proliferative cells ( Figure 3B). We randomly selected some DEGs for qPCR verification. Muscle differentiation-related genes such as MYH15 and MyoM are highly expressed in the differentiation stage. Cell cycle-related genes, such as CDKN1A and CCNB1, are significant differences in expression between the two stages ( Figure 3C). These results are consistent with the RNA sequencing results ( Figure 3D).  The data in C are mean ± S.E.M. with three replicates per group. One-sample t-test was used to assess the difference between the two groups. * p < 0.05; ** p < 0.01.

Differential Splicing Analysis of mRNA During Goose Myoblast Differentiation
When DNA is transcribed into RNA, it is first transcribed into Pre-mRNA by RNA polymerase. Pre-mRNA is spliced into mature mRNA by the spliceosome, and the mature RNA can be further translated into protein [31,32]. Different types of cells can produce with three replicates per group. One-sample t-test was used to assess the difference between the two groups. * p < 0.05; ** p < 0.01.

Differential Splicing Analysis of mRNA during Goose Myoblast Differentiation
When DNA is transcribed into RNA, it is first transcribed into Pre-mRNA by RNA polymerase. Pre-mRNA is spliced into mature mRNA by the spliceosome, and the mature RNA can be further translated into protein [31,32]. Different types of cells can produce different splicing variants through variable splicing [33,34]. Therefore, we next analyzed the differentially splicing genes (DSGs) between goose myoblast and myotube (Table 3). A total of 854 DSGs were found between the two groups (Supplementary File S3), and the main splicing type is skipped exon (SE), which accounts for 67.70% ( Figure 4A). GO and KEGG analyses were performed on the 854 DSGs. A KEGG pathway analysis found that DSGs were mainly enriched in the Wnt signaling pathway and fatty acid biosynthesis ( Figure 4B). A GO enrichment analysis found that the DSGs were enriched in skeletal muscle tissue development, sequence-specific DNA binding, and regulation of alternative mRNA splicing, via spliceosome ( Figure 4C). We also analyzed genes that were both differentially expressed and differentially spliced. A total of 43 genes were significantly differentially expressed and differentially spliced between the two groups ( Figure 4D, Table S1). The functional analysis found that the 43 genes were significantly enriched in an extracellular matrix and actin binding (Table 4).

Differential Expression Analysis of LncRNA during Goose Myoblasts Differentiation
The development of myogenesis is regulated by epigenetic factors, including DNA methylation, RNA methylation, histone modification, and noncoding RNA [35][36][37]. There have been many reports on the research of lncRNA in poultry [38,39]. In this study, we examined the differential expression analysis of lncRNA in the proliferation stage and differentiation stage of Shitou goose myoblasts. A total of 244 DE-lncRNAs were found between the two groups ( Figure 5A). A KEGG enrichment analysis found that DE-lncRNAs were enriched in the cGMP-PKG signaling pathway, the calcium signaling pathway, the Wnt signaling pathway, and the cAMP signaling pathway ( Figure 5B). These pathways are related to muscle development [40,41]. A GO analysis of DE-lncRNA found that the functions of lncRNAs were mainly enriched in epigenetic regulation, such as transcriptional regulation and transcriptional regulation of an RNA polymerase promoter ( Figure 5C). Next, we constructed a putative DE-lncRNA-DEGs interactive network involved in goose myogenesis. Among them, 316 lncRNAs were found to target 150 cis-DEGs, with a total of 631 pairs of lncRNA-mRNA connections (Figure 6, Supplementary  Figure S2). The expression relationship of 378 pairs of lncRNA-mRNA is positive, and the expression relationship of 253 pairs of lncRNA-mRNA is negative.

Differential Expression Analysis of LncRNA during Goose Myoblasts Differentiation
The development of myogenesis is regulated by epigenetic factors, including DNA methylation, RNA methylation, histone modification, and noncoding RNA [35][36][37]. There have been many reports on the research of lncRNA in poultry [38,39]. In this study, we examined the differential expression analysis of lncRNA in the proliferation stage and differentiation stage of Shitou goose myoblasts. A total of 244 DE-lncRNAs were found between the two groups ( Figure 5A). A KEGG enrichment analysis found that DE-lncRNAs were enriched in the cGMP-PKG signaling pathway, the calcium signaling pathway, the Wnt signaling pathway, and the cAMP signaling pathway ( Figure 5B). These pathways are related to muscle development [40,41]. A GO analysis of DE-lncRNA found that the functions of lncRNAs were mainly enriched in epigenetic regulation, such as transcriptional regulation and transcriptional regulation of an RNA polymerase promoter ( Figure 5C). Next, we constructed a putative DE-lncRNA-DEGs interactive network involved in goose myogenesis. Among them, 316 lncRNAs were found to target 150 cis-DEGs, with a total of 631 pairs of lncRNA-mRNA connections ( Figure 6, Supplementary Figure S2). The expression relationship of 378 pairs of lncRNA-mRNA is positive, and the expression relationship of 253 pairs of lncRNA-mRNA is negative.

Discussion
In this study, we identified the genes and pathways related to the proliferation and differentiation of Shitou goose myoblasts through transcriptomics and spliceomics analysis. Since the proliferation and differentiation of myoblasts affect muscle growth and hypertrophy [42], the determination of related genes and pathways has important reference value for subsequent meat goose breeding. We also found that there was different splicing of multiple genes during myoblast proliferation and differentiation, indicating that one gene may play different roles in the regulation of muscle development.
Animal meat production mainly depends on the increase in the number of muscle fibers and muscle fiber hypertrophy. Previous studies have shown that the number of muscle fibers does not increase after birth and acquired muscle growth mainly depends on muscle fiber hypertrophy [43]. As a large goose species, the Shitou goose has faster muscle development before hatching, and it has more and larger muscle fibers than the small goose species [8]. Similarly, a chicken with high weight selection grows faster in the later stage of embryonic development than a chicken with low weight selection [44]. These results suggest that more and faster muscle fiber development before hatching is critical for the postnatal growth of the Shitou goose. The sequencing samples in this study are myoblasts isolated from 15 embryonic Shitou goose embryos. The myoblasts obtained in this period have active proliferation and differentiation potential [45], better reflecting the proliferation of myoblasts in the embryonic stage of Shitou goose and showing the process of differentiation and fusion of myoblasts to form polynuclear muscle fibers.
In recent years, it has become more and more popular to use bioinformatics to analyze changes in mRNA expression profiles of growth and development in non-model animals. As an excellent meat goose variety, it is very important to study the mRNA expression profile for the regulation of Shitou goose muscle development. In birds, muscle growth and development are regulated by a variety of signaling pathways, including the JAK-STAT signaling pathway [46], the MAPK signaling pathway [47,48], the PI3k

Discussion
In this study, we identified the genes and pathways related to the proliferation and differentiation of Shitou goose myoblasts through transcriptomics and spliceomics analysis. Since the proliferation and differentiation of myoblasts affect muscle growth and hypertrophy [42], the determination of related genes and pathways has important reference value for subsequent meat goose breeding. We also found that there was different splicing of multiple genes during myoblast proliferation and differentiation, indicating that one gene may play different roles in the regulation of muscle development.
Animal meat production mainly depends on the increase in the number of muscle fibers and muscle fiber hypertrophy. Previous studies have shown that the number of muscle fibers does not increase after birth and acquired muscle growth mainly depends on muscle fiber hypertrophy [43]. As a large goose species, the Shitou goose has faster muscle development before hatching, and it has more and larger muscle fibers than the small goose species [8]. Similarly, a chicken with high weight selection grows faster in the later stage of embryonic development than a chicken with low weight selection [44]. These results suggest that more and faster muscle fiber development before hatching is critical for the postnatal growth of the Shitou goose. The sequencing samples in this study are myoblasts isolated from 15 embryonic Shitou goose embryos. The myoblasts obtained in this period have active proliferation and differentiation potential [45], better reflecting the proliferation of myoblasts in the embryonic stage of Shitou goose and showing the process of differentiation and fusion of myoblasts to form polynuclear muscle fibers.
In recent years, it has become more and more popular to use bioinformatics to analyze changes in mRNA expression profiles of growth and development in non-model animals.
As an excellent meat goose variety, it is very important to study the mRNA expression profile for the regulation of Shitou goose muscle development. In birds, muscle growth and development are regulated by a variety of signaling pathways, including the JAK-STAT signaling pathway [46], the MAPK signaling pathway [47,48], the PI3k pathway [49], the focal adhesion pathway [50], the TGF-β signaling pathway [51] and so on. In this study, we found that the JAK-STAT signaling pathway and PI3k pathway were significantly enriched in DEGs during goose proliferation (down-DEGs). Many studies on these signaling pathways related to muscle have been reported. For example, the JAK-STAT signaling pathway is related to satellite cell proliferation and muscle atrophy [52]. In the mouse model, the JAK-STAT signaling pathway is related to Critical Illness Myopathy (CIM), and Chronic JAK/STAT activation promoters loss of muscle mass and function [53]. The P13K-Akt signaling pathway at the GM stage is significantly enriched (Figure 2C), indicating that it plays an important role in the proliferation of Shitou goose myoblasts and skeletal muscle development. The PI3K Akt signaling pathway is highly conserved in animals [54] and plays an important role in normal cell activities related to cell growth, proliferation, movement, survival, metabolism, and apoptosis [55]. The PI3K-Akt signaling pathway is also associated with muscle-related diseases and muscle development [56], indicating that this pathway may play an important regulatory role in Shitou goose muscle development.
This study not only analyzed the DEGs based on transcriptome sequencing data but also analyzed the DSGs during goose myoblast differentiation. Splicing depends on the dynamic assembly of ribonucleoprotein complexes, including spliceosomes and a large number of helper proteins. The spliceosomes influence the fate of RNA and cell developmental processes by selecting alternative splicing sites [9]. At the single-cell level, the dynamic splicing of RNA is the key to the bioinformatics evaluation of cell fate [57]. This suggests that differential splicing may be an important factor for the smooth progress of cell differentiation during the muscle development of the Shitou goose. Some pieces of literature have pointed out that some signaling pathways can indirectly affect cell fate by affecting splicing factors, such as the MAPK signaling pathway [58] and the Wnt signaling pathway [59]. The data of this study also showed that the DEGs and DSGs during Shitou goose myoblasts differentiation were significantly enriched in the Wnt signaling pathway, indicating that this pathway may regulate the proliferation and differentiation of myoblasts by affecting gene signaling transduction and mRNA splicing.

Conclusions
In conclusion, we constructed the mRNA and lncRNA expression profiles of Shitou geese during myoblast differentiation. A total of 1664 DEGs, 854 DSGs, and 244 DE-lncRNAs were identified between myoblasts and myotubes. By analyzing the biological functions of DEGs, DSGs, and DE-lncRNAs, a set of signaling pathways, such as the Wnt signaling pathway, was found to be highly related to the muscle development of Shitou geese. The results of this study can provide a basis for a follow-up study on the mechanism of goose muscle growth and development.  Table S1: 43 differentially expressed genes with simultaneous differential splicing; Supplementary File S1: Codes for Gene differential splicing analysis; Supplementary File S2: codes for R Studio; Supplementary File S3: Differentially splicing genes in the main splicing type.

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets supporting the conclusions of this article are available in GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE213147, public on 15 September 2022. The data analysis pipeline (Linux code and R scripts) is provided as a zip file in the supplemental materials.