RNA Sequencing Identifies Upregulated Kyphoscoliosis Peptidase and Phosphatidic Acid Signaling Pathways in Muscle Hypertrophy Generated by Transgenic Expression of Myostatin Propeptide

Myostatin (MSTN), a member of the transforming growth factor-β superfamily, plays a crucial negative role in muscle growth. MSTN mutations or inhibitions can dramatically increase muscle mass in most mammal species. Previously, we generated a transgenic mouse model of muscle hypertrophy via the transgenic expression of the MSTN N-terminal propeptide cDNA under the control of the skeletal muscle-specific MLC1 promoter. Here, we compare the mRNA profiles between transgenic mice and wild-type littermate controls with a high-throughput RNA sequencing method. The results show that 132 genes were significantly differentially expressed between transgenic mice and wild-type control mice; 97 of these genes were up-regulated, and 35 genes were down-regulated in the skeletal muscle. Several genes that had not been reported to be involved in muscle hypertrophy were identified, including up-regulated myosin binding protein H (mybph), and zinc metallopeptidase STE24 (Zmpste24). In addition, kyphoscoliosis peptidase (Ky), which plays a vital role in muscle growth, was also up-regulated in the transgenic mice. Interestingly, a pathway analysis based on grouping the differentially expressed genes uncovered that cardiomyopathy-related pathways and phosphatidic acid (PA) pathways (Dgki, Dgkz, Plcd4) were up-regulated. Increased PA signaling may increase mTOR signaling, resulting in skeletal muscle growth. The findings of the RNA sequencing analysis help to understand the molecular mechanisms of muscle hypertrophy caused by MSTN inhibition.

Blocking the myostatin signaling transduction pathway can induce a highly muscled phenotype [21]. The ability of several myostatin antagonists, such as follistatin, follistatin-related gene, soluble activin receptor IIB and suramin to block MSTN signaling transduction has been tested [22][23][24][25]. Furthermore, some microRNAs can affect the MSTN signaling pathway by targeting the 3' UTR of myostatin [26,27]. We previously generated transgenic mice via the muscle-specific expression of a cDNA sequence (5'-region 886 nucleotides) encoding the MSTN propeptide. The transgene effectively depressed myostatin function, and transgenic mice showed dramatic growth and increases in muscle mass due to muscle fiber hypertrophy [28][29][30]. To further elucidate the mechanisms of muscle hypertrophy caused by MSTN inhibition, we compared the gene expression profiles of the muscle tissue between MSTN propeptide transgenic mice and their littermate wild-type controls. Our results indicate that the Ky/mybph/actins pathway is a novel noncanonical signaling pathway, and phosphatidic acid inhibits mTOR signaling via Dgki, Dgkz and Plcd4 in the hypertrophic skeletal muscle of the MSTN propeptide transgenic mice.

Skeletal Muscle Weights of the Transgenic Mice and Wild-Type Littermate Controls
The mice were sacrificed at four months of age to evaluate the muscle weight and obtain muscle samples. The results were consistent with our previous report; the individual major muscles of the myostatin propeptide transgenic mice were significantly heavier (38% to 95% more muscle) than those of the wild-type littermate controls ( Figure 1). The muscle groups that increased in mass predominantly consisted of fast-twitch muscle fibers, including the gastrocnemius, tibialis anterior, biceps femoris, triceps brachii, longissimus dorsi, and semitendinosus.

Solexa Sequencing Data and Reads Assembly Using Tophat Software
In this study, Solexa sequencing was used to detect the mRNA expression profiles of MSTN propeptide transgenic mice and their littermates. Five RNA libraries were constructed from the total RNA isolated from the gastrocnemius muscle of the transgenic and control mice, which generated 9,710,772 (CN148), 10,371,693 (CN164), 9,991,395 (TN126), 10,258,517 (TN135), and 10,076,991 (TN329) raw reads. After filtering the low quality reads and trimming off the adapted fragments, 8,717,967,9,418,459,9,054,439,9,301,124, and 9,032,266 clean reads were obtained, and the ratios of cleans reads in raw reads were 89.78%, 90.81%, 90.62%, 90.67%, and 89.63% (Table 1). The detailed description of the sequence data is shown in Table 1, and the abundance of mRNA reads mapped to each chromosome is shown in Figure S1. The Tophat aligner software was used to align the sequencing reads to the reference genome download from the ENSEMBEL website ftp://ftp.ensembl.org/pub/release-75/fasta/ mus_musculus/dna/. The mapping rate of each sample exceeded 97% ( Table 2). The statistical alignment data indicated that the sequencing reads were of high quality and the sequencing depth was sufficient for a differential expression analysis between the two groups. Subsequently, HTSeq was selected to control the quality of aligned reads. The quality control result is shown in Figure S2.

Differentially Expressed Genes of Skeletal Muscle between Transgenic Mice and Their Littermate Controls
The cufflinks software was used to assemble the individual transcripts that had been aligned to the genome and quantify the expression level of each transfrag in the sample. After cufflinks treated the reads, cuffmerge was used to parsimoniously merge the assembled transfrags. The annotation reference file of Mus musculus integrated into the merged assembly was downloaded from ENSEMBEL (ftp://ftp.ensembl.org/pub/release-75/gtf/mus_musculus/) [31]. Cuffdiff was used to calculate the differential expression and analyze the significances of observed changes between two groups. A cluster analysis was performed to elaborate the expression patterns of genes in MSTN propeptide transgenic mice and their controls ( Figure 2A). The results showed that approximately 70% genes were up-regulated in the transgenic mice. CummeRbund volcano plots ( Figure 2B) were available to intuitively observe the genes that were differentially expressed between the two groups. Finally, the expression levels of 132 genes were identified to be significantly different between the transgenic mice and control mice; among these genes, 97 were up-regulated and 35 were down-regulated (FDR < 0.05). Of the 132 genes, 117 were annotated in the database. The differentially expressed genes with more than two-fold changes and an FPKM value more than 5 are listed in Table 3. Several genes have not been reported to be involved in skeletal muscle hypertrophy, including myosin binding protein H (mybph), beta, gamma and alpha actin (Actb, actc1, actg1), protein tyrosine phosphatase receptor type c (Ptprc), zinc metallopeptidase STE24 (zmpste24), and genes related to phosphatidic acid (PA) pathways (Dgki, Dgkz and Plcd4). Moreover, the gene kyphoscoliosis peptidase (Ky), which plays an important role in muscle growth, was also up-regulated in the skeletal muscle of the transgenic mice. The bar plots of the differentially expressed genes are shown in Figure S3, and all differentially expressed genes identified by RNA-seq are summarized in Table S1. Each row represents a gene ID and the column represents a sample; (B) The volcano plots reveal genes that differ significantly between two groups. The black dots represent the expression levels of genes that do not differ and the red dots represent the expression level differences.

Function Annotation of the Differentially Expressed Genes
To gain insight into the function of the differentially expressed genes in the hypertrophic skeletal muscle of the MSTN propeptide transgenic mice, the DAVID (Database for Annotation, Visualization and Integrated Discovery) website (http://david.abcc.ncifcrf.gov/) was used to identify the gene function with the following parameters: Count = 2 and EASE = 0.01. The gene ontology functional classification is shown in Figure 3. The functions of genes in biological processes, cellular components and molecular processes were annotated based on the GO categories. Most genes were enriched in two types of cellular components: the extracellular region and extracellular space. The genes enriched in biological processes were mainly involved in biological adhesion and cell adhesion. The genes enriched in molecular function were mainly associated with calcium ion binding.

Validation of the Differentially Expressed Genes by RT-qPCR
To verify the reliability of the sequencing data, five genes were randomly selected from Table 3 for RT-qPCR analyses ( Figure 4). As expected, the gene expression patterns validated by RT-qPCR strongly correlated with the sequencing results (Pearson's r correlation coefficient is 0.97591, Figure 5). Both RT-qPCR and solexa sequencing indicated that MSTN propeptide, Cbr2, Tnfaip2 and Mybph were up-regulated, while Mss51 was down-regulated in MSTN propeptide transgenic mice. Among the differentially expressed genes, the MSTN propeptide mRNA level increased by 12.78-fold according to solexa sequencing ( Figure S4) and 43-fold according to the RT-qPCR analyses ( Figure 4A). We also detected the expression level of endogenous MSTN; the endogenous myostatin expression did not significantly differ (p > 0.05) between the transgenic mice and their littermate wild-type controls ( Figure 4B).

Signaling Pathway Analysis
To identify pathways that are involved in the regulation of muscle hypertrophy in the MSTN propeptide transgenic mice, we employed the DAVID website to analyze the differentially expressed genes with the following parameters: Count = 2 and EASE = 0.1. The pathways and related genes are listed in Table 4. This analysis identified the hypertrophic cardiomyopathy (HCM) and dilated cardiomyopathy signaling pathways, both of which regulate muscle growth and development. In addition, pathways associated with the maintenance of muscle structure and myocytes were identified, including the ECM-receptor interaction and focal adhesion pathways. Notably, the gene expression of three key enzymes, Dgki, Dgkz and Plcd4, of the phosphatidic acid (PA) signaling system was significantly enhanced in the skeletal muscle of the transgenic mice compared with their littermate controls ( Figure 6). The hematopoietic cell lineage genes Itga4, cd3d and cd5 were also identified.

Discussion
High-throughput sequencing technologies are currently emerging as a powerful tool to initially elaborate the genetic mechanisms responsible for a phenotype. Previous studies have detected genes related to the hypertrophic muscle observed in myostatin propeptide transgenic mice by employing microarray analysis [30]. The results identified 52 unique genes that are differentially expressed, and myostatin propeptide promotes adult muscle growth and maintenance by up-regulating the expression levels of myogenic regulatory factors and extracellular matrix components and down-regulating protein degradation-related genes as well as mitochondrial ATP synthesis-related genes [30].
To further illustrate the molecular mechanism and identify potential regulatory genes associated with the muscle hypertrophy caused by myostatin inhibition, we employed RNA sequencing technologies to study differences in the expression patterns between MSTN propeptide transgenic mice and wild-type littermate control mice at four months of age; these mice were younger than those used in the microarray study. In the current study, 132 differentially expressed genes were identified in MSTN propeptide transgenic mice compared with control mice based on the RNA sequencing technology and Cufflinks analysis method. Among these 132 genes, most were up-regulated in transgenic mice, and only 35 genes were down-regulated. Functional annotation showed that these differentially expressed genes mainly correlated with muscle development, lipid metabolism, energy metabolism and the immune system. Intriguingly, both of the abovementioned studies revealed that the function of some differentially expressed genes involves the extracellular region. Additionally, pathway analysis found four pathways in both studies: focal adhesion, ECM-receptor interaction, the regulation of the actin cytoskeleton, and cardiac muscle contraction.
In the current study, we found that the expression levels of cytoskeletal and muscle growth and development-related genes were changed in the transgenic mice. Strikingly, Myh6 was down-regulated in the transgenic mice compared with the wild-type controls. A previous study showed that protein synthesis and muscle hypertrophy were enhanced when the expression of myh6 was repressed [37][38][39], suggesting that myh6 is a negative regulator of muscle mass. This finding is consistent with our results, which showed that Myh6 expression was decreased in muscle hypertrophy mice. Moreover, a previous study reported that the mitochondria number, the ratio of mitochondrial DNA to nuclear DNA, and the level of oxidative enzymes used to synthetize ATP all decreased in myostatin-deficient muscle [40]. Other studies showed that Myh6 had higher ATPase activity compared with other myosins [38,41]. Thus, the energy metabolism of myostatin propeptide transgenic mice may decrease due to a decrease in the expression of Myh6, which reduces ATP consumption. Unlike Myh6, the muscle growth-related genes Ky, Actb, Mybph, Actc1, RMST and Actg1 were up-regulated in our study. Among these genes, Ky (kyphoscoliosis peptidase) plays a vital role in muscle growth; the absence of Ky protein leads to muscular dystrophy [42,43]. A Ky mutation results in degenerative myopathy and loss of muscle hypertrophy responses [42]. Ky protein can interact with myosin-binding protein and filamin C [44]. Moreover, Myosin-binding protein H (Mybph) reportedly interacts with Rho kinase 2 (ROCK2) to affect the interaction of actin-myosin [45][46][47]. Additionally, Actb, Actc1 and Actg1 are three types of actins that participate in the maintenance of the cytoskeleton [48]. A functional analysis of these cytoskeletal and muscle development-related genes suggested that myostatin propeptide can enhance muscle mass via the Ky/Mybph/Actins pathway.
Lipids play a critical role in signal transduction; they often act as a secondary messenger to activate downstream pathways [49]. We found that the overexpression of myostatin propeptide activated the phosphatidic acid (PA) signaling pathway, which generates signaling lipids and participates in the cell cycle and cell growth [50,51]. In the current study, we found that Plcd4 and Dgkz, two key enzymes of the PA signaling pathway, were up-regulated in myostatin propeptide transgenic mice. This result suggested a connection between the myostatin signaling and PA signaling pathways. Plcd4 is a member of the delta class of phospholipase C enzymes. Phospholipase C enzymes can catalyze hydrolyzing phosphatidylinositol 4,5-bisphosphate (PIP2) to two intracellular secondary messengers, inositol 1,4,5-trisphosphate (IP3) and diacylglycerol (DAG) [32]. Previous studies reported that DAG could be phosphorylated by Dgkz, and the phosphorylated DAG could be digested to phosphatidic acid (PA). PA can reportedly bind to mTOR to induce muscle fiber hypertrophy [52][53][54][55][56]. Other studies reported that myostatin mediates myoblast differentiation and myotube hypertrophy by inhibiting Akt/mTOR/p70S6 protein signaling [18,20,57]. Overall, we may infer that myostatin can negatively regulate the PA signing pathway to inhibit the mTOR signaling pathway (Figure 4). In addition to activating the mTOR signaling pathway, PA also activates p21-activated kinase 1 (PAK1) to initiate the release of RhoGDI from Rac1, and RhoGD1 subsequently changes the actin dynamics [58].
In biology, redox reactions can frequently store and release biological energy by electron transfer processes. In our study, we found that heme-binding protein Mss51, which is related to redox reactions, was down-regulated in myostatin propeptide transgenic mice. In yeast, Mss51 can regulate the biogenesis of cytochrome c oxidase (COX), the terminal electron transport chain (ETC) oxidase, thereby controlling energy production [59,60]. Previous studies showed that cytochrome c oxidase subunit VIc (Cox6c) expression is decreased in myostatin propeptide transgenic mice [30]. Although we did not detect differential Cox6c expression in our study, the expression of its upstream gene, Mss51, was decreased. Thus, we can infer that myostatin may influence ATP synthesis via Mss51. In addition, the gene Cbr2, which participates in neutralizing excessive reactive oxygen species (ROS), is also regulated by myostatin propeptide [61][62][63]. Myostatin induces the production of ROS, which decreases muscle mass by regulating antioxidant gene expression and mitochondrial biogenesis [64][65][66]. In combination with our results, these findings suggest that myostatin leads to muscle dystrophy by negatively regulating Cbr2 to increase ROS.

Muscle Tissue Sample Collections
Myostatin propeptide cDNA under the control of rat myosin light chain 1 (MLC1) with its enhancer, SV40 poly A signal sequence, was constructed and served as the transgene construct, and transgenic mice were generated by standard microinjection techniques as previously described [28][29][30]. Male mice (hemizygous genotype for the transgene) from the high-expressing line were mated with B6SJL F1 wild-type females. The mice were housed in cages with a 12-h light/dark cycle, and the room temperature was maintained at 22 °C. The mice were weaned at four weeks of age and given free access to a standard diet (10% kcal fat, ME 3.85 kcal/g). All procedures and animal care were in accordance with the institution guidelines and approved by the Institutional Animal Care and Use Committee of the University of Hawaii, Honolulu, HI, USA. Male mice were sacrificed at four months of age after 8 h of fasting to sample and dissect the muscle tissue. Gastrocnemius muscle samples from both legs were immediately dissected from the carcass, quickly frozen in liquid nitrogen, and stored at −80 °C. Three transgenic male mice and two wild-type littermate male control samples were used for this study.

RNA Isolation and RNA Sequencing
Total RNA was extracted from freshly frozen muscle tissues using the TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. A magnetic bead homogenizer was used to homogenize the tissue and TRIzol Reagent. The quality and concentration of RNA was detected by NanoDrop ND2000 (Thermo Fisher Scientific, Waltham, MA, USA) spectrophotometry and gel electrophoresis. Solexa sequencing was used to detect the differentially expressed genes in the gastrocnemius of the transgenic and control mice. Five cDNA libraries were constructed using the TruSeq Stranded Total RNA LT Sample Prep Kit (Illumina, Santiago, CA, USA). The library construction and solexa sequencing were performed by Genergy biological technology Limited company (Shanghai, China).

Quality Control for Raw Sequencing Data
The raw sequencing data were in a FASTQ file, and the data with clean reads were obtained by trimming the adapter contaminants and filtering the low-quality reads. A quality control tool, HTSeq (https://pypi.python.org/pypi/HTSeq), which can provide a quick impression of high throughput data, was utilized in our study. The quality control results are presented in the supplementary information.

Bioinformatics Analysis
The software Bowtie2 (http://bowtie-bio.sourceforge.net/index.shtml/) was used to merge the chromosome genomes of Mus musculus downloaded from the ENSEMBL database (ftp://ftp.ensembl.org/pub/release-75/fasta/mus_musculus/dna/) and convert them to the FM index, which can store massive genomic data and rapidly be searched. Tophat (http://ccb.jhu.edu/software/ tophat/index.shtml), which uses Bowtie as an engine to align the sequencing reads to reference genome, was used to align the clean reads of each individual to the large genomes. Subsequently, the transcriptome was assembled with the Cufflinks software (http://cufflinks.cbcb.umd.edu/), and these assemblies were merged using Cuffmerge in the presence of the reference genome annotation (ftp://ftp.ensembl.org/pub/release-75/gtf/mus_musculus/). Cuffmerge provided a uniform basis to calculate the expression level of transcripts and genes at different conditions. Cuffdiff, a separate program contained in Cufflink, was used to calculate the expression level of each individual and test the statistical significance between transgenic and control mice. Cuffdiff used FPKM (fragments per kilobase of transcript per million mapped fragments) method to calculate the expression level. In addition, CummeRbund was used to visualize and integrate all results provided by the Cuffdiff analysis. DAVID Bioinformatics Resources were used for gene annotation and pathway analysis.

Differential Expression Gene Validation by Real Time RT-PCR
Quantitative real-time PCR was performed to validate the differentially expressed genes identified by the Solexa sequencing data. Five differentially expressed genes were randomly selected for quantitative real-time PCR analysis. The Primer Premier 5.0 (PREMIER Biosoft International, Palo Alto, CA, USA) software was used to design the quantitative primers based on the cDNA sequence. RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) was used for RNA reverse transcription following the manufacturer's instructions, and DNase I and RNase-free (Thermo Fisher Scientific, Waltham, MA, USA) were used to treat the RNA sample to avoid genomic DNA contamination. The qPCR reaction was performed on a LightCycler 480 Real-Time PCR instrument (Roche, Penzberg, Germany) as follows: single cycle of denaturation at 95 °C for 5 min, 40 cycles of denaturation at 95 °C for 15 s, annealing at 60 °C for 15 s and extension at 72 °C for 15 s. SYBR Green (BIO-RAD, Hercules, CA USA) was used for the qPCR reaction, and all reactions were performed in triplicate. The 2 −ΔΔCt method was used to calculate the relative expression levels between transgenic and control mice, and all expression levels were normalized to the expression of 18S RNA. Significant differences between groups were assessed with a two-tailed t-test. The difference between groups was considered significant when p < 0.05.

Conclusions
In this study, gene expression profiles associated with hypertrophic skeletal muscle in myostatin propeptide transgenic mice were investigated by RNA sequencing. A total of 132 genes were identified and used for function annotation and gene enrichment analysis. The novel up-regulated genes were related to hypertrophic cardiomyopathy, focal adhesion, and PA signaling pathway, and ROS removal processes, suggesting that myostatin directly affects these processes during the regulation of skeletal muscle growth and mass. Two new signaling pathways, the Ky/Mybph/Actins and PA signaling pathways, were uncovered, which may account for myostatin propeptide-induced muscle hypertrophy. Our results provide information to help further understand how myostatin inhibition causes skeletal muscle hypertrophy.