Gene Co-Expression Network Analysis Provides Novel Insights into Myostatin Regulation at Three Different Mouse Developmental Timepoints

Myostatin (Mstn) knockout mice exhibit large increases in skeletal muscle mass. However, relatively few of the genes that mediate or modify MSTN effects are known. In this study, we performed co-expression network analysis using whole transcriptome microarray data from MSTN-null and wild-type mice to identify genes involved in important biological processes and pathways related to skeletal muscle and adipose development. Genes differentially expressed between wild-type and MSTN-null mice were further analyzed for shared DNA motifs using DREME. Differentially expressed genes were identified at 13.5 d.p.c. during primary myogenesis and at d35 during postnatal muscle development, but not at 17.5 d.p.c. during secondary myogenesis. In total, 283 and 2034 genes were differentially expressed at 13.5 d.p.c. and d35, respectively. Over-represented transcription factor binding sites in differentially expressed genes included SMAD3, SP1, ZFP187, and PLAGL1. The use of regulatory (RIF) and phenotypic (PIF) impact factor and differential hubbing co-expression analyses identified both known and potentially novel regulators of skeletal muscle growth, including Apobec2, Atp2a2, and Mmp13 at d35 and Sox2, Tmsb4x, and Vdac1 at 13.5 d.p.c. Among the genes with the highest PIF scores were many fiber type specifying genes. The use of RIF, PIF, and differential hubbing analyses identified both known and potentially novel regulators of muscle development. These results provide new details of how MSTN may mediate transcriptional regulation as well as insight into novel regulators of MSTN signal transduction that merit further study regarding their physiological roles in muscle and adipose development.

Previous studies have demonstrated that MSTN blockade can functionally improve dystrophic muscle [4] and that MSTN can also play an important role in adipose development [5]. To better understand the molecular mechanisms by which MSTN controls muscle growth and development, microarray data from MSTN-null and wild-type mice were used for gene coexpression network analysis during key timepoints in skeletal muscle development. These timepoints included primary myogenesis at 13.5 days post coitus (d.p.c.), secondary myogenesis at 17.5 d.p.c., and postnatal muscle growth at 35 days after birth (d35) [6][7][8].
In previous studies, the partial correlation and information theory (PCIT) approach and the regulatory impact factor (RIF) algorithm have been used to evaluate the regulatory functions of transcription factors (TFs) on gene expression [9][10][11]. The PCIT algorithm can be used to determine the differential hubbing (DH; [9]) of targets to putative regulatory genes. A regulator that exhibits DH across treatments will have drastically higher connectivity in one treatment compared to the others. The RIF algorithm identifies potential regulators that are "differentially wired" (i.e., differentially correlated) as a function of treatment (Mstn genotype, in the context of this study). Two methods for calculating RIF assign different weights to the target gene's level of expression. The RIF1 score also ranks genes based on their phenotypic impact factor (PIF; a score that combines overall expression and differential expression). RIF1 scores place more emphasis on target transcripts highly differentially wired due to treatment, while RIF2 scores identify regulators whose expression tends to predict changes in expression of target transcripts that are differentially expressed (DE) between treatments. Previous studies evaluating RIF scores demonstrate that genes with relatively stable expression levels in different conditions (i.e., genes not differentially expressed) can still play important roles in biological processes in a coordinated manner. But these studies have limitations: first, these analyses mainly considered DE genes, which were quite limited and may not be representative of all the changes that take place in response to an experimental treatment. Second, only genes with known regulatory functions (i.e., transcription factors) were evaluated in previous studies. Without considering all genes, it is difficult to discover novel regulators of skeletal muscle growth.
In this study, we applied modified versions of the PCIT [12] and RIF algorithms [10,11] to evaluate more than 34,000 annotated genes from a whole microarray dataset as potential regulators of muscle and adipose development in the presence or absence of functional myostatin. We consider the RIF method modified because all genes were tested as possible regulators using RIF scores, whereas previous studies using this methodology considered only TFs as possible regulators. The PCIT method utilized in this anlaysis used the same methodology but was computationally enhanced to complete the analysis in considerably less time than the original algorithm [12]. We calculated RIF scores that contrasted Mstn genotypes at each developmental timepoint to identify potential regulators of Mstn or pathways influenced by MSTN signal transduction. Furthermore, DE genes were compared across Mstn genotypes for differences in the enrichment of TF binding sites, biological processes, and pathways related to MSTN function.

Ethics Statement
This study was carried out in strict accordance with the recommendations in the Animal Guide of the United States Department of Agriculture. The protocol was approved by the Institutional Animal Care and Use committee of Iowa State University (Protocol Number: 11-05-6013-M).

Animals and sample collection
Animals were reared, sacrificed by CO 2 asphyxiation (embryonic timepoints) or cervical dislocation (d35), and tissues were collected from whole embryos (13.5 d.p.c.), hindquarters (17.5 d.p.c.), or pectoralis muscle (d35). The pectoralis muscle was chosen at day 35 because it is the muscle with the greatest change in muscle mass due to Mstn deletion and had a high level of myostatin expression as described in the original paper characterizing myostatin function in skeletal muscle [2]. After collection, muscles were flash-frozen in liquid nitrogen as previously described [13].

RNA isolation and target preparation
The frozen tissues from whole embryo or hindquarter only were immersed in RLT lysis buffer (Qiagen, Valencia, CA) and homogenized for 45 seconds, then centrifuged at 3600g for 10 min to pellet cell debris. Total RNA was isolated using the RNase Midi Kit (Qiagen) with on-column DNase digestion according to manufacturer's protocol. Isolation from pectoralis muscle was conducted as described [13]. The concentration of isolated RNA was evaluated based on absorbance at 260 nm and the quality was determined by agarose gel electrophoresis. For microarray analysis, 10 μg of RNA from each of three individual animals was pooled to give five pools representing 15 indviduals for each genotype/timepoint combination. A final concentration of 1 μg/ml was yielded from the pooled RNA, the quality of which was reconfirmed using the 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). Microarray target preparation, hybridization, and scanning were conducted as described [13].

Microarray processing
Target RNA was hybridized to both chip A and chip B of the GeneChip Mouse Expression Set 430 (MOE430; Affymetrix, Santa Clara, CA), which includes over 45,000 probe sets and 39,000 transcripts and variants, containing over 34,000 well-annotated mouse genes. Robust multiarray analysis (RMA) was performed to normalize the raw data [14] by using R/affy package [15]. All of the microarray data is available at NCBI GEO under accession number GSE63154.

Statistical analysis
The normalized gene expression data were used for further analysis. To detect differentially expressed genes in each genotype, expression data with five replicates were analyzed as a completely randomized design by using R/limma package [16]. Because the tissues are different across timepoints, the analysis was performed separately for each timepoint. The model used for analysis was, where y ij is the normalized gene expression, μ is the overall mean, G i is the ith genotype, and ε ij is the error for the jth replicate within the ith genotype. DE genes were identified at a false discovery rate (FDR) <0.05 using the q-value package from Bioconductor [17].

RIF analysis and differential hubbing
Using all of the Affy ProbesetIDs after excluding the control IDs and housekeeping genes, we applied a modified partial correlation and information theory approach (PCIT) [10,12] to compute the co-expression correlation between each gene pair in both MSTN-null and wildtype conditions. Co-expression analyses were done only within timepoint since tissue types were not identical across different developmental stages. As described in Reverter et al. [11], we calculated the regulatory impact factors (RIF), which compute the differential wiring from the difference in co-expression correlation of each pair of genes. A key difference in this analysis is that we used all genes for the RIF analysis to try to identify novel regulators of muscle growth in response to lack of functional MSTN. Differential hubbing calculations were conducted using a Perl script that captured only the significant partial correlations identified by the PCIT algorithm. Correlates were filtered to include only those genes with |r| 0.90 to ensure that only strong correlations were captured, given the small sample size and large amount of data. The concept of RIF, PIF, and differential hubbing can be challenging to understand in the context of biological processes. Hudson et al. [9] provide a nice description of how these co-expression methods can be used to understand the genetic regulation of muscle mass.

Annotation and pathway analysis
Annotation of ProbesetIDs is based on the annotation file from the Affymetrix website (http:// www.affymetrix.com), while the IDs without annotation were annotated using Biomart (Ensembl, http://www.ensembl.org/biomart/martview/c56c1fc10cb311e0b5fe1c05dbd3f217) and DAVID [18]. Differentially expressed genes and genes co-correlated to top RIF were used to identify over-represented pathways across different developmental timepoints in web-based DAVID, while the significant regulators of these gene lists were analyzed in Sub-Network Enrichment Analysis (SNEA) via Pathway Studio [19].
Discovering significantly over-represented motifs using DREME Among those genes up-and downregulated between genotypes as well as co-expressed with top RIF regulators, we applied a motif discovery method for non-redundant, statistically significant, and discriminative motifs using Discriminative Regular Expression Motif Elicitation (DREME) [20] implemented in iPlant (http://www.iplantcollaborative.org). We used 1,500 bp upstream flanking sequence in FASTA format. Generally, we ran the analysis using the default settings. First, the sequences of genes of interest from up-and downregulated groups were used for positive input sequences, while the negative sequences were the sequences of well-annotated genes from both MOE430A and MOE430B chips after excluding those genes of interest. Second, we chose 3 as the minimum width of core motif and 8 as the maximum. Third, 0.05 was chosen as the motif e-value cutoff, where the e-value is the motif statistic p-value times the number of candidate motifs tested, and the counts created after removing sites that matched motifs found previously were used to calculate the p-value.
To evaluate the discovered motifs, we used TOMTOM [21] implemented in MEME Suite (http://meme.nbcr.net/meme/intro.html) to compare the motifs with the database of known transcription factor motifs in JASPAR (JASPAR_CORE_2009.meme; 476 motifs) [22] and UniProbe (uniprobe_mouse.meme; 386 motifs) [23] using default settings with the significance threshold 0.05. Table 1 and Table 2 present the top 20 DE genes ranked by significance at 13.5 d.p.c. and d35, respectively. The DE genes were identified using R/limma (S1 Table) and a significance threshold of FDR < 0.05. At 13.5 d.p.c., 283 genes were DE, 17 of which were downregulated and 266 of which were upregulated in WT. Surprisingly, we did not identify any differentially expressed genes (FDR < 0.05) at 17.5 d.p.c. At d35, 2034 genes were DE, 928 of which were downregulated and 1106 of which were upregulated in wild-type mice.

Differentially expressed genes
DAVID annotation and Pathway Studio analysis of DE genes at 13.5 d.p. c. and d35 Pathway analysis results from DAVID and Pathway Studio 9 are listed in S2 and S3 Tables. The 283 genes DE at 13.5 d.p.c. were significantly enriched for 46 Gene Ontology (GO) terms [24] and several PANTHER pathway terms and SP_PIR_KEYWORDS [18]. Most of the 13.5 d.p.c. DE genes were involved in the regulation of transcription, RNA processing, negative regulation of cellular biosynthetic process and transcription, nucleotide binding, transcription repressor activity, and WNT signaling pathway. For the 2034 genes DE at d35, 58 GO terms, 19 SP_PIR_KEYWORDS, and several KEGG and PANTHER pathway terms were over represented. The DE genes at d35 were largely involved in muscle organ/tissue development, fatty acid metabolism, and mitochondria. Pathway Studio results at d35 identified 100 significantly over-represented regulator-target relationships (p < 0.05); PPARG, LEP, calmodulin, and GH1 were linked to the largest numbers of differentially expressed genes. MEF2A, MEF2C, creatine kinase, and UCP2 were also significantly enriched in the regulator-target relationships.
Over-represented DNA motifs in differentially expressed genes

Regulatory impact factor analysis
Tables 3-6 list the top ten most extreme RIF scores, ranked by Z-score, at 13.5 d.p.c. and d35. All RIF scores are provided in S4 . Several genes were identified with top-ten RIF scores at both timepoints (i.e., Beta-s, GM10420, March-2). In addition, Tnni2 was identified in the top ten at d35 for both RIF1 and RIF2. Target genes highly connected to Tnni2 (Pearson correlation coefficient based on PCIT, |r| 0.9) were further investigated using ontology enrichment analysis to identify potential pathways involved in muscle development. Genes highly correlated with Tnni2 in wild-type mice were enriched in the KEGG pathways viral myocarditis, immune/cancer, and proteolysis (adj.P, Benjamini < 0.1; S5 Table). However, no KEGG terms were significantly enriched in Tnni2 correlates in MSTN-null mice.

Differential hubbing
Differential hubbing scores are listed in S6 Table. Fig. 4 illustrates extreme DH observed for Ctnna1, Bmp8b, and Myl3 in wild-type vs. MSTN-null mice for highly correlated target genes (|r| 0.9). For instance, there were over 2,000 genes highly co-expressed (correlated) with Ctnna1 in MSTN-null mice, but only one gene was highly co-expressed with Ctnna1 in wild type at 13.5 d.p.c. Interestingly, Ctnna1 was not differentially expressed (q = 0.56). Bmp8b was also found to be differentially hubbed between wild-type and MSTN-null mice at 13.5 d.p.c.; over 3,000 genes were highly co-expressed with Bmp8b in wild-type mice, while less than 1,000 genes were highly co-expressed with Bmp8b in MSTN-null mice. Fig. 5 contrasts the differences in expression patterns of genes highly correlated to Myl3 across genotypes at d35. Expression data used in Fig. 5 are listed in S7 Table. The 5,500 genes highly correlated to and co-expressed with Myl3 in wild-type mice were almost always downregulated, whereas the 1,189 co-expressed genes in MSTN-null mice were almost exclusively upregulated at d35. In MSTN-null mice, the Myl3 differentially hubbed correlates were enriched for the KEGG pathway term mmu04360: "Axon guidance". In contrast, those genes highly correlated with Myl3 in wild-type mice were significantly enriched for the KEGG terms actin cytoskeleton, T/B cell receptor signaling, viral myocarditis, Huntington's disease, cancer pathways, and ubiquitinmediated proteolysis (Benjamini p < 0.10; Table 7).

Discussion
Interpretation of genes differentially expressed during primary myogenesis and postnatal skeletal muscle growth Microarray analysis of tissues from MSTN-null and wild-type mice identified DE genes at both 13.5 d.p.c. and d35, whereas no DE genes (q < 0.05) were identified at 17.5 d.p.c. Ideally, comparisons would have been made using the same tissue sample from each of the timepoints. However, obtaining single muscles from 13.5 d.p.c. and 17.5 d.p.c. embryos was not feasible Co-Expression Network Analysis of Myostatin Regulation at 3 Timepoints due to their small size at these stages of development. An important consequence is that no comparisons can be made across timepoints in this study since the tissue types are not comparable. It is also important to note that results from 13.5 and 17.5 d.p.c. may reflect gene expression changes in tissues other than muscle, even though myostatin is predominantly expressed in muscle. We would also like to note that this study is a meta-analysis of a previous study and that, although there is significant overlap, the number of DE genes identified at d35 in this study differs from the results of the previous study in our lab [13]. This is because different statistical methods were used in these two analyses. Despite the difference in the methodology used to analyze differential expression in the current study, all of the DE genes validated by real-time qPCR in our previous study were also significantly DE in this study (q<0.05; S1 Table) [13].
Ontology-based enrichment analysis identified biological processes related to myostatin-mediated changes in skeletal muscle growth and metabolism.
DAVID ontology enrichment analyses identified significantly over-represented regulators relevant to both muscle and adipose development. At 13.5 d.p.c., over-representation of the WNT KEGG pathway was identified for genes DE by genotype. This result is consistent with a previous study from our lab that also identified a major role for WNT signal transduction in MSTN-mediated skeletal muscle growth at the d35 postnatal timepoint [13]. The DE genes at d35 were enriched for the DAVID terms muscle organ/tissue development, fatty acid MicroRNA-regulated gene networks were over represented at 13.5 d.p.c.
Pathway Studio identified potential regulators within DE gene lists based on text mining of previously published results. Several microRNAs were identified as potential regulators from the sub-network analysis at 13.5 d.p.c., including Mir29b1, Let7g, Mir139, and Mir122. Since microRNAs are involved in many early developmental processes, it is possible that they may be involved in muscle and adipose development in pathways upstream or downstream of MSTN signaling.
Fatty acid metabolism, myogenic, and muscle fiber type pathways were over represented at d35 Pathway enrichment with Pathway Studio at d35 after birth identified a number of significantly enriched pathways (p<0.05) relevant to muscle and adipose development. For example, several myogenic factors' sub-networks were significantly enriched, including networks regulated by MEF2A, MEF2C, MYH1, and TNNT2 (S3 Table). MEF2C was previously found to be involved in the formation of myofibers [25]. In addition, several sub-networks involved in fatty acid metabolism were also significantly enriched at d35, including networks regulated by ADIPOQ, CEBPA, FABP4, LDL, LEP, PGC1A, PPARG, SREBF1, and VLDL. PPARG, LEP, calmodulin, and GH1 were linked to the largest numbers of differentially expressed genes at d35. UCP2, a mitochondrial membrane transporter involved in energy balance regulation, was significantly enriched as well (p = 0.00573; S3 Table). Notably, WNT4 pathway members were over represented at d35 (p = 0.0227), which is consistent with previously experimentally confirmed results from this study [13]. The transcription factor GATA5, which is an important regulator in muscle development [26], was identified by Pathway Studio because its targets, including the muscle fiber type development genes Myl2, Tnnc1, Myh7, Myl3, and Myh11, were significantly over represented (p = 0.00133). These results indicate that genes DE due to MSTN loss of function at 13.5 d.p.c. and d35 include several important regulators of muscle and adipose development. All Pathway Studio results are presented in S3 Table. Known and novel transcriptional regulatory motifs were differentially enriched in the presence or absence of functional myostatin at 13.5 d.p. c. and d35 Previously, SMAD transcription factors have been implicated in MSTN signal trasnduction [27]; however, not all promoters of MSTN-regulated genes contain SMAD binding sites. Using DREME, several over-represented motifs were identified in DE genes at both 13.5 d.p.c. and d35 ( Figs. 1 and 2, Figures A and B in S1 File). Zinc-finger protein family (ZFP) binding site was the most significant motif identified at 13.5 d.p.c. Binding sites for PLAGL1, a known biomarker of disease [28], were over represented in genes upregulated in WT at 13.5 d.p.c. Neither Plagl1 (q = 0.72) nor Zfp187 (q = 0.54) was DE at 13.5 d.p.c. However, Plagl1 itself was significantly 0.67-fold downregulated in WT at d35 (q = 0.01). In contrast, ZFP187 TF binding sites were over represented in genes downregulated in WT at 13.5 d.p.c., while the gene itself was significantly 1.21-fold upregulated in WT at d35 (q = 0.03). It is interesting that PLAGL1 binding sites were identified as over represented in genes that are more highly expressed in the WT genotype, since the similar gene Plag1 was recently identified in a haplotype that had undergone a selective sweep in double-muscled Belgian Blue cattle [29].
At d35, SMAD3 binding site was found to be the most over-represented motif in DE genes upregulated in MSTN-null mice at d35, indicating its importance in MSTN signal transduction during postnatal growth. SMAD3 has proved to be necessary and sufficient for MSTN signal transduction during muscle development [30]. In addition, SP1 was identified as an over-represented motif in DE genes upregulated in WT at d35. Sp1 was a DE gene and 1.29-fold upregulated in wild-type mice at d35. The SP1 transcription factor was reported to play an important role in the transcriptional activity of genes [31] and induction of lipogenesis via FOXO1 [32]. Further, SP1 is known to be important in skeletal muscle proliferation through transcriptional regulation of fibroblast growth factor receptor 1 (FGFR1) [33,34], and differentiation by regulation of ERK5/Sp1/Klf [35].
An important unanswered question addressed in this study is how MSTN signals to all of the genes that it regulates. Besides SMAD3 and SP1, additional over-represented motifs were identified in DE genes between MSTN-null and wild-type mice. These results are important because they suggest novel transcriptional regulators that MSTN may signal through to directly or indirectly alter transcriptional regulation of genes within its signal transduction pathway. However, these sequence motifs have no known transcription factors that bind to them. It is possible that currently unknown sequence motifs could regulate MSTN-reactive genes at these sites. These results provide novel insights into possible mechanisms of transcriptional regulation in the context of MSTN signal transduction and indicate that there is still much to learn about MSTN signaling.

Co-expression methods identified genes known to interact with myostatin as well as novel putative pathway components
Application of the RIF and PIF algorithms [11] to all genes in this dataset (over 34,000 annotated genes) allowed genes to be ranked as regulators, not just transcription factors. This allowed novel regulatory genes to be identified in this study. All PIF scores are provided in S8 Table. Genes with the 20 highest PIF scores at 13.5 d.p.c. (Table 8) were shown by DAVID analysis to be enriched for several RNA and protein modification terms, including protein acetylation, methylation, and RNA-binding (Benjamini p<0.05; S9 Table). Several cell cycle and cancer-related genes were also over represented. The highest PIF value genes from d35 (Table 9) appear to be related to fiber type changes between MSTN-null and wild-type mice (e.g., Tpm1, Ckm, Mylpf, Actn3, Myh4, Tnnt3, Actc1, Myh1, and Tnni2). A DAVID enrichment analysis of the top 20 PIF genes confirms that these genes are involved in a multitude of muscle-specific structural and metabolic functions consistent with fiber type changes that occur due to MSTN inactivation (S9 Table).
Several known regulators of skeletal muscle growth were identified using the RIF method. Potentially novel regulators that fit closely with known MSTN functions were also identified. The following passages present several novel possibilities regarding regulation of MSTN signal transduction.
At 13.5 d.p.c., voltage-dependent anion channel 1 (Vdac1) and thymosin B4 (Tmsb4x) were identified among the top ten RIF2 scores. The Vdac1 gene encodes a mitochondrial porin that appears to regulate mitochondrial function by regulating the metabolites and ions that can enter the mitochondria [36,37]. A recent study reports that treatment of HeLa cells with MSTN leads to increases in VDAC1-mediated apoptosis via BAX [38]. Interestingly, Vdac1 knockout mice also have impaired exercise capacity and reduced glucose tolerance [36]. The Tmsb4x gene has been proposed to be involved in wound healing and muscle regeneration. Studies in C2C12 and satellite cell-derived primary myocytes and myoblasts indicate that Tmsb4x is more highly expressed at muscle wound sites and that it acts as a myoblast chemoattractant during muscle regeneration [39]. In the mdx mouse model of muscular dystrophy, chronic administration of TMSB4X was also found to increase the number of regenerating skeletal muscle fibers [40].
The Sox2 gene was identified as a top-ten ranked RIF1 gene at 13.5 d.p.c. Several studies have documented the role of SOX2 in SMAD-mediated TGFβ signal transduction during embryonic development [41]. In addition, Sox2 sustained expression is promoted by TGFβ signal transduction to promote embryonic stem cell self renewal [42]. Therefore, it is plausible that altered MSTN signal transduction, which acts through TGFβ, may impact SOX2 signal transduction and processes associated with stem cell renewal. Since whole embryos were used for 13.5 d.p.c., it is not possible to determine which tissues are actually affected.
At d35, Apobec2 and Mmp13 were identified as potential regulators based on high RIF2 and RIF1 scores, respectively. The APOBEC2 protein is an RNA-editing enzyme that is expressed almost exclusively in striated muscle and appears to be important in fiber type determination. Mouse knockouts of Apobec2 tend to have a greater ratio of slow to fast twitch muscle, 15-20% reduction in body mass, and evidence of mild myopathy in aged mice [43]. Vonica   differentiation of muscle, and inhibition of TGFβ signal transduction [44]. Based on studies in zebrafish, it has been suggested that APOBEC2 may mediate its effects in part by altering methylation patterns [44,45]. The role of MMP13 in skeletal muscle has been suggested to center around myoblast migration and differentiation during muscle regeneration. Pharmacological blockade of MMP13 in C2C12 cells leads to reduced myoblast migration and differentiation, whereas overexpression leads to improved myoblast migration. Expression of Mmp13 appears to be increased during muscle repair as well as following myoblast fusion during myotube formation. Although its exact role is unknown, it has been suggested to alter extracellular matrix remodeling during muscle regeneration [46]. The identification of Atp2a2 (a.k.a. Serca2a) as the top-ranked RIF2 gene at d35 is particularly noteworthy, because it was previously identified as a potential modifier of MSTN in a large F2 mouse study [47]. The ATP2A2 calcium pump affects sarcoplasmic reticulum function and muscle fiber type. Transgenic Atp2a2 mice exhibit cardiac hypertrophy [47][48][49][50]. Furthermore, pharmacological inhibition of SMAD2 in cardiomyocytes indicates that TGFβ/SMAD2 signaling regulates ATP2A2 function [51]. Therefore, it is very plausible that MSTN may regulate ATP2A2 in skeletal muscle to mediate changes in muscle fiber type and size.
Differentially hubbed genes that may act as regulators or key mediators in muscle and adipose development in the context of MSTN are listed in S1 Table. Ctnna1 was identified as differentially hubbed at 13.5 d.p.c. despite not being DE, which indicates that MSTN may affect some regulators both directly and indirectly so as to cause or influence the differential Co-Expression Network Analysis of Myostatin Regulation at 3 Timepoints expression of target genes. Similarly, Bmp8b was identified to have over 2,000 genes differentially hubbed in wild-type versus MSTN-null mice (Fig. 4). A recent report demonstrated that BMP8B directly regulates thermogenesis in brown adipose tissues [52], which indicates that MSTN may alter adipose metabolism through BMP8B. In addition, recent research identified a mutation in Myl3 that can cause hypertrophic cardiomyopathy in infants [53]. In our study, over 4,000 genes were differentially hubbed with Myl3 in wild-type compared to MSTN-null mice. Interestingly, significant enrichment of biological pathways involved in muscle development and disease occurred only in wild-type mice ( Table 7, adj.P, Benjamini < 0.1), which is consistent with previous findings that MSTN plays an important role in phenotypes such as muscle hypertrophy or disease caused by muscle dysfunction [54,55].

Conclusions
Differentially expressed genes were identified between wild-type and MSTN-null mice during primary myogenesis and postnatal muscle growth, but not during secondary myogenesis. In addition, co-expression analyses identified several known and novel regulators or pathway effectors of myostatin. Regulators with known roles in muscle growth, TGFβ signaling, and mitochondrial function were identified by the RIF method. Another important finding was that genes not DE or a priori considered regulators were identified as important in differentiating wild-type and MSTN-null mice. These potential regulators would not have been identified if analyses had been limited to considering only transcription factors as regulators. Furthermore, motif enrichment analyses identified several transcription factor binding sites in genes differentially expressed due to Mstn genotype. This study demonstrates that all genes must be considered as potential regulators when trying to identify critical control points of pathway circuitry.
Supporting Information S1 File. Supporting figures. Figure A in S1 File. Over-represented motifs in genes upregulated in WT at d35. Figure B in S1 File. Over-represented motifs in genes upregulated in MSTN-null at d35. (PDF) S1