Expression profiling by high-throughput sequencing reveals GADD45, SMAD7, EGR-1 and HOXA3 activation in Myostatin (MSTN) and GDF11 treated myoblasts

Abstract

Growth differentiation factor 11 (GDF11) and myostatin (MSTN/GDF8) are closely related members of the transforming growth factor β (TGFβ) superfamily, sharing structural homology. Despite these structural similarities, recent research has shed light on the distinct roles these ligands play within muscle tissue. This study aims to uncover both the differences and similarities in gene expression at the transcriptome level by utilizing RNA sequencing. We conducted experiments involving five distinct groups, each with three biological replicates, using C2C12 cell cultures. The cells were subjected to high-throughput profiling to investigate disparities in gene expression patterns following preconditioning with either GDF11 or MSTN at concentrations of 1 nM and 10 nM, respectively. In addition, control groups were established. Our research revealed concentration-dependent gene expression patterns, with 38 genes showing significant differences when compared to the control groups. Notably, GADD45, SMAD7, EGR-1, and HOXA3 exhibited significant differential expression. We also conducted an over-representation analysis, highlighting the activation of MAPK and JNK signaling pathways, along with GO-terms related to genes that negatively regulate metabolic processes, biosynthesis, and protein phosphorylation. This study unveiled the activation of several genes not previously discussed in existing literature whose full biological implications are yet to be determined in future research.

Keywords:
GDF11; Myostatin; RNA sequencing; myoblasts; gene expression

Introduction

Considering demographic changes in today’s society, preserving muscle function and its regenerative potential are key goals in morbidity and mortality prevention in the aging population (Li et al., 2018; Papadopoulou 2020). In conjunction with the high costs, age-related muscle dysfunction leads to a significant reduction of life quality, progressive weakness, and disability (Janssen et al., 2004; Siparsky et al., 2014; Trombetti et al., 2016). A wide variety of diseases manifesting in muscle dysfunction, with no satisfactory treatment alternatives and no curative therapy options, further increase the clinical urgency and growing interest in muscle aging and regeneration research.

Cytokines of the transforming growth factor β (TGFβ) family contain more than 30 biochemically related but functionally distinct ligands, which are known to play a crucial role in the regulation of developmental patterning, cellular proliferation and differentiation, and the maintenance of tissue homeostasis (Walker et al., 2016; Morikawa et al., 2016). Among others, Growth Differentiation Factor 11 (GDF11) and MSTN are structurally closely related TGFβ family members and share 90% homology in their mature active regions (Walker et al., 2016).

While MSTN is expressed primarily in skeletal muscles, GDF11 is expressed broadly in numerous tissues (McPherron et al., 1997; Nakashima et al., 1999; Gamer et al., 1999) and has been subject to conflicting reports regarding its role in muscle regeneration, aging-process and cancer. Some studies also managed to demonstrate tumor suppressive properties of GDF11 in the Breast and Liver cancer (Bajikar et al., 2017; Gerardo-Ramírez et al., 2019). Additionally, its expression seems to correlate negatively with tumor stage in colorectal cancer and positively with the survival rate in patients with pancreatic cancers (Yokoe et al., 2007; Liu et al., 2018). Prior work also implicates the role GDF11 plays in aging-related phenotypes in the heart, skeletal muscle, and brain (Loffredo et al., 2013; Katsimpardi et al., 2014; Sinha et al., 2014; Egerman et al., 2015; Poggioli et al., 2016; Walker et al., 2016; Kraler et al., 2023) and its expression in epithelial cells of the developing stomach and duodenal epithelium (Harmon et al., 2004).

Despite the high structural similarity between GDF11 and MSTN, mutations in these molecules lead to very different phenotypical manifestations, indicating divergent functionality across numerous tissues. Unlike genetic deficiency of MSTN, which leads to hypermuscularity (McPherron et al., 1997; McPherron and Lee, 1997), homozygous deletion of GDF11 generates defects in axial skeletal patterning and organ development, subsequently leading to perinatal lethality in GDF11-null mice (McPherron et al., 1999).

GDF11 and MSTN are synthesized as precursors and undergo proteolytic processing to produce biologically active mature ligands. Both proteins bind to activin type II receptors (ACVR2A or ACVR2B), and recruit activin receptor-like kinase 4 and/or 5 (ALK4/ ALK5) forming a heteromeric receptor complex to induce downstream intracellular SMAD2/3 mediator pathway signaling via phosphorylation. GDF11 also showed the ability to activate SMAD1/5/9 phosphorylation through the utilization of the ALK1 receptor in different tissues, leading to proliferation and differentiation. Additionally, even though both, GDF11 and MSTN, are able to activate similar receptor types as well as ALK7 receptors, GDF11 has shown to initiate a stronger signal, implicating a higher binding affinity (Walker et al., 2016; Walker et al., 2017).

Despite all similarities in protein sequence or receptor utilization, along with signaling pathways, several studies demonstrated that these two ligands may have different functions across numerous tissues and among others in muscle cells (McPherron et al., 1997; Nakashima et al., 1999; Walker et al., 2016). It has been shown that GDF11 and MSTN are able to initiate some unique activity on muscle, aside from its activation of known (ActRII/ALK/SMAD) signaling. Previous studies suggest that there are noncanonical pathways for activating other non-SMAD proteins, such as p38 MAPK, ERK, and JNK, as well as different mechanisms of modulation of the activity of GDF11 and MSTN through extracellular binding proteins. These extracellular binding proteins, which typically function as antagonists, include follistatin (FST), follistatin-like 3 (FSTL3/FLRG), growth and differentiation factor-associated serum protein 1 (GASP1), GASP2, latent TGF-β binding protein 3 (LTBP3), and decorin (DCN) (Philip et al., 2005; Yang et al., 2006; Wang et al., 2014; Egerman et al., 2015; Biesemann et al., 2015; Walker et al., 2016).

Believed to have very similar functions, the question, of why those two structurally highly homologous proteins have different secretion patterns and effects within various cell types and additionally vary in concentration throughout the lifespan, remains a high priority and subject to ongoing investigations. In this study, we aimed to investigate and expose differences in gene activation patterns, for the first time utilizing RNA sequencing in GDF11 and MSTN-treated C2C12 cell cultures, a subclone from a myoblast cell line, originally isolated by Yaffe and Saxel at the Weizmann Institute of Science in Israel (Yaffe and Saxel, 1977).

Due to the fact that MSTN-null mice survived to adulthood, whereas GDF11-null mice died shortly after birth, a great number of studies used recombinant GDF11 and/ or MSTN proteins in order to reveal the exact function of these ligands. The role of MSTN in satellite cells has been the subject of controversial results, ranging from an inhibitor of C2C12 myoblast proliferation in early studies (Thomas et al., 2000; Taylor et al., 2001), to Rodgers et al. (2014) arguing that recombinant MSTN stimulates C2C12 proliferation. Additionally, the source of the recombinant protein could be an outcome-determining factor on itself (Rodgers et al., 2014).

Sinha et al. (2014) demonstrated that the aged mice treated with recombinant GDF11 protein injections displayed improved metrics in skeletal muscle strength, endurance, muscle regeneration, and even myofibrillar and mitochondrial morphology (Sinha et al., 2014). Further study in fish showed the application of GDF11 recombinant proteins being able to boost antioxidant enzyme activity in muscle, prolonging the lifespan (Zhou et al., 2019). While several studies failed to reproduce those results or even showed conflicting outcomes like impaired muscle regeneration and satellite cell expansion (Egerman et al., 2015; Hinken et al., 2016) our study aimed to provide insights and deeper understanding of the particular function of GDF11 and MSTN in muscle tissue and further outline possible target genes for disease monitoring and interventions.

Material and Methods

RNA sequencing preparation

C2C12 cells, a subclone of the mouse myoblast cell line established by D. Yaffe and O. Saxel (ATCC CRL-1772), biosafety level 1(biosafety classification is based on U.S. Public Health Service Guidelines), were acquired from Thermo Fisher Scientific. Initial experiments were performed at the Harvard University Department of Stem Cell and Regenerative Biology at Wagers Laboratory, Cambridge, MA, US. C2C12 cells are immortalized mouse myoblast cells with very short replication time that can be rapidly differentiated into functional skeletal or cardiac muscle cells (McMahon et al., 1994).

C2C12 cell culturing protocol

C2C12 cells from frozen stock were cultured in a 50 ml flask, and 10 ml special media containing 445 ml of 1 x DMEM - Dulbecco’s Modified Eagle Medium, 50 ml FBS - Fetal Bovine Serum, and 5ml Pen-Strep - Penicillin-Streptomycin (10,000 units penicillin and 10 mg streptomycin/mL). Cultures were split when they reached 80% confluency, approximately every 48 h.

Activation of 2/3 SMAD and/or p-SMAD (sc-6032, Lot#G2214 and sc-11769, Lot# H2014, Antibodies, Santa Cruz Biotechnology) pathways was confirmed using the Protein Simple Wes system.

mRNA extraction

All experimental groups had three biological replicates (A-C), with a total number of 15 samples. C2C12 cells were incubated in the growth media until approximately 70% confluency was reached, then washed with PBS and incubated in a starvation media for 3 hr. The relatively short starvation allowed to avoid reduced cell survival and increased apoptosis as a reaction to prolonged starvation. It also provided enough time for the adaptation processes in cell metabolism and the induction of cell cycle synchronization.

Starvations media contained 5 mL Pen-Strep = Penicillin-Streptomycin (10,000 units penicillin and 10 mg streptomycin/mL), 495 mL DMEM, and 200 uL FBS (0.02%).

Recombinant pure proteins GDF11, and MSTN (PEPROTECH, USA, REF #120-11 and #120-00) were prepared using following buffer: 10 mM HCL, 0.1% protease free BSA, and PBS. Subsequently, starvation media was changed, and the cultures were divided into five experimental groups and co-incubated with recombinant proteins GDF11, and MSTN at concentrations of 1 nM and 10 nM, as well as “buffer only” vehicle group, for the duration of one hour. Based on previous studies the concentration of both proteins varied across different tissues and during the lifespan. Additionally, there is a known difference in receptor affinity for both proteins. To account for both factors in our experiments, we utilized two different concentrations for each protein.

Plates were twice washed in ice cold PBS, then total RNA was extracted from all 20 experimental samples using TRIzol Reagent Protocol (Thermo Fisher: Pub. No. MAN0001271). The concentration of RNA was confirmed using NanoDrop and Qubit. BioAnalzyer was used to assess the quality of the probes. RNA Integrity Number (RIN) values ranged from 9.9 - 10.

Library preparation

Sequencing was performed at Bauer Core Facility, Harvard University, Cambridge, USA. Libraries were prepared using a SciClone G3 NGSx workstation (Perkin Elmer) using the Kapa mRNA HyperPrep kit (Roche Sequencing). Polyadenylated mRNAs were captured using oligo-dT-conjugated magnetic beads (Kapa mRNA HyperPrep kit, Roche Sequencing) from 500 ng of total RNA on a Perkin Elmer SciClone G3 NGSx automated workstation. Poly-adenylated mRNA samples were immediately fragmented to 200-300 bp using heat and magnesium. First-strand synthesis was completed using random priming, followed by second-strand synthesis and A-tailing. A dUTP was incorporated into the second strand to allow strand-specific sequencing of the library. Libraries were enriched and indexed using nine cycles of amplification (Kapa mRNA HyperPrep kit, Roche Sequencing) with PCR primers, which included dual 8 bp index sequences to allow for multiplexing (IDT for Illumina unique dual 8 bp indexes). Excess PCR reagents were removed through magnetic bead-based cleanup using KAPA Pure magnetic beads on a Sciclone G3 NGSx workstation (Perkin Elmer). The resulting libraries were assessed using a 4200 TapeStation (Agilent Technologies) and quantified by qPCR (Roche Sequencing). Libraries were pooled and sequenced using paired-end, 75 bp reads. On average, 27.4 M reads were obtained per library (SD: 3.3 M).

RNA sequencing data analysis

Sequence reads were processed with fastp (v0.20.1) to remove sequences of sequencing adapters and low quality (Phred quality score below 15) sequences from the 3’-end of the sequence reads (Chen et al., 2018). Hereafter, reads were aligned to the mouse reference assembly (GRCm39.104) using STAR (v2.7.9a) (Dobin et al., 2013). Differential expression was assessed with DESeq2 (Love et al., 2014). A gene was considered to be significantly differentially expressed if the corresponding false discovery rate (FDR) did not exceed a value of 0.05 and the absolute value of the logarithmic fold change (log2FC) was 1 or higher. The detection of GO-terms over-represented in sets of differentially expressed genes was performed using GOrilla (database version: Mar 6, 2021) (Eden et al., 2009) in combination with the Gene Ontology database (Ashburner et al., 2000).

Results

High-throughput gene expression profiling was conducted on C2C12 cell cultures, comprising five distinct experimental groups, each with three biological replicates. Within each group, cells were pre-treated with either GDF11 or MSTN at concentrations of 1 nM and 10 nM. Control groups treated with a vehicle were also included in the study.

We conducted separate assessments of genes that exhibit upregulation or downregulation when compared to the control group. Although we did not anticipate or investigate linear correlations, we observe clear distinctions among the treatment groups (^{Figure S1}). Thirty-five genes were significantly differentially expressed in the comparison of GDF11 10 nM and the control. 28 of these were higher expressed in GDF11 10 nM, whereas seven were higher expressed in the control. In the comparison of GDF11 1 nM, MSTN 10 nM and MSTN 1 nM 26, 12, and 2 genes were differentially expressed respectively, whereas the genes differentially expressed in the comparisons including MSTN were wholly contained in the sets of genes of the comparisons including GDF11 (Figure 1). In addition we performed direct comparisons between all treatment groups. In the comparison of GDF11 1 nM and MSTN 1 nM, three genes (Smad7, Wnt9a and Hbegf) are significantly differentially expressed. In the comparison GDF11 10 nM and MSTN 1 nM, there are six differentially expressed genes (Smad7, Wnt9a, Hbegf, Bhlhe40, Dusp2 and CCn2). No differentially expressed genes could be detected in the other direct comparisons.

Figure 1 -
Intersections of genes that are significantly differentially expressed in comparison to the control. The two respectively 12 genes that are differentially expressed between MSTN and the control are also differentially expressed when GDF11 is compared to the control. Contrary to the general trend, there appear to be three genes that are differentially expressed solely between GDF11 1 nM and the control. However, DUSP6 and TNFSF15 both satisfy the FDR criterion and miss the foldchange criterion (|log2FC| >= 1) in GDF11 10 nM compared to control only narrowly with log2 foldchanges of -0.9 and 0.97.

While at first, it might seem that there are 12 or 13 genes that play a role solely in GDF11 10 nM and GDF11 1 nM, respectively, a closer look reveals that gene expression instead undergoes a gradual change from control to MSTN 1 nM, MSTN 10 nM, GDF11 1 nM, and then GDF11 10 nM (Figure 2). The comparatively small changes between GDF11 1 nM and GDF11 10 nM could indicate that, in the case of GDF11, a concentration of 1 nM already produces an effect that can be increased at most slightly at higher concentrations. Apart from the concentration dependence described above, no expression differences suggesting different regulatory mechanisms of GDF11 and MSTN were observed.

Figure 2 -
Row Z-Scores of the normalized expression of all 38 genes that were significantly differentially expressed in at least one comparison. The expression of most genes increases (n=24) or decreases (n=4) continuously from the Control over MSTN 1 nM, MSTN 10 nM and GDF11 1 nM to GDF11 10 nM. The corresponding gene symbols are printed in bold. Additionally, genes which were exclusively differentially expressed in both of the two GDF10 groups are marked with an asterisk.

Over-representation analysis (ORA) was performed based on the comparison of GDF11 10 nM and the control. Genes corresponding with the Gene Ontology (GO) -terms associated with MAPK activation (GO:1900745, GO:0000185, GO:0032874, GO:0046330, GO:0043410) and JNK (GO:0046330) were significantly enriched and displayed the smallest FDR. Furthermore, the analysis showed GO-terms correlated to negative regulation of metabolic processes (e.g., GO:0010605, GO:0009892), biosynthesis activity (e.g., GO:0009890), and protein phosphorylation (e.g., GO:0001933).

The sequence reads have been submitted to the European Nucleotide Archive (ENA) as part of the study PRJEB57932. For accession numbers please refer to ^{Table S1}.

Discussion

To our knowledge, this is the first study to report gene expression pattern comparison between GDF11 and MSTN in C2C12 cell culture using RNA sequencing (RNA-seq) technology. In previous microarray analysis of hSkMDCs treated with GDF11 or MSTN performed by Egerman et al. (2015) the human primary muscle cells were exposed to a 24-hour stimulation with 300 ng/ml of MSTN, GDF11, or buffer alone (serving as a negative control). Subsequent analysis of the log fold change in gene expression demonstrated that GDF11 and MSTN induced very similar, insignificantly different, genes expression patterns (Egerman et al., 2015). However, RNA-seq provides counts of aligned sequence reads, resulting in a very broad dynamic range, improving the possibility of detecting rare transcripts and conclusively increasing its sensitivity and accuracy. Additionally, RNA-seq does not use probes or primers, therefore, the data suffer from much lower biases allowing more robust results. Mentioned above also results in RNA-seq outrunning microarrays for their high degree of reproducibility (Wang et al., 2009; Wilhelm and Landry, 2009; Zhao et al., 2014; Liu et al., 2015; Rao et al., 2019). Rather than restricting the comparison of the gene activation pattern of GDF11 and MSTN directly, we took the approach of specifically looking into the expression patterns in comparison to the vehicle controls, which allowed us to further eliminate noise, increasing the sensitivity of our experiment.

These are two structurally very similar and closely related ligands that naturally circulate in the blood in their inactive form. Recent studies demonstrate that the molar concentration of circulating MSTN seemed to be around 500 times higher than that of GDF11 (Rodgers and Eldridge, 2015). Furthermore, being ~90% identical in their mature domains, there is only ~52% identity in their pro-domains (Walker et al., 2017), leading to the conclusion that their activity depends on structural changes, previously described regulatory mechanisms and receptor utilization in target tissue (Zimmers et al., 2002; Walker et al., 2016). Walker et al. (2017) described MSTN as having lower binding affinity to all type 1 receptors, conclusively GDF11 inducing a greater SMAD3-dependent signal (Walker et al., 2017). Although the biological consequences of the described differences remain to be determined, our study clearly underlined differences in binding affinity, indirectly confirming these findings. As visualized in Fig. 2, the gradual increase in signal correlates with the respective increase in ligand concentration, with MSTN requiring higher concentration to induce comparable strength in downstream activation.

Upon utilization of the string network analysis, we could demonstrate that there are described relationships between 38 differentially expressed genes in GDF11/ MSTN vs. controls, with a subset of genes clustered in groups. Even though our study did not specifically aim to explain the function of the 38 genes found to be significantly differentially expressed, there is a rationale for deriving possible functions based on clustered groups and already published data. Here we present a subset of these genes, that might offer a deeper understanding of functioning of these two ligands. In our belief, supported by previous studies, these selected genes represent the most promising targets for in-depth analysis, which could reveal unknown functions and interactions of these two ligands in future investigations.

Gadd45 α/ β/ γ network

GADD45 (growth arrest and DNA-damage inducible) is a gene family of three highly homologous small nuclear proteins: GADD45α, GADD45β, and GADD45γ, that have been linked to regulation of many cellular functions including DNA repair, cell cycle control, genotoxic stress, as well as senescence and apoptosis. Due to the pro-apoptotic activities, GADD45 is believed to play an essential role in oncogenesis, functioning as a tumor suppressor. All three members of the family are found in muscle tissue with GADD45β (also known as Myd118) being the least abundant of the three (Tamura et al., 2012).

MSTN is a potent negative regulator of skeletal muscle growth and its inactivation, as previously described, induces skeletal muscle hypertrophy, while its overexpression leads to muscle atrophy. The role of GDF11 in muscle has been subject to conflicting reports, and its function in muscle tissue is still to be determined. Interestingly our results show equal dose-dependent activation of the GADD45 network with both ligands, pointing towards a similar function. Bullard et al. (2016) and Ebert et al. (2012) both described the role of GADD45α in muscle atrophy either by activating stress-inducible, pro-atrophy transcription factor ATF4 or the interaction with MAP3K4, a mitogen-activated protein kinase complex (Ebert et al., 2012; Bullard et al., 2016). Over-representation of GO-terms linked to negative regulation of metabolic processes (e.g., GO:0010605, GO:0009892), biosynthesis activity (e.g., GO:0009890), as well as protein phosphorylation (e.g., GO:0001933), further supports the theory of negative regulation of skeletal muscles through GADD45 network (^{Figure S2}). To our knowledge, there is no previous reporting on GADD45 activation through GDF11/ MSTN, which could be a missing link between these two ligands and their activity in muscles and an exciting study subject for future investigations.

SMAD7

Forbes et al. (2006) described Smad7 as being involved in the auto-regulation of the MSTN promotor by the negative feedback loop (Forbes et al., 2006). The auto-regulation triggered by MSTN appears to be dose-dependent, as shown by our results, demonstrating a correlation between expression averages and the MSTN concentration used (Forbes et al., 2006). Furthermore, according to a more recent study, the loss of Smad7 leads to various changes in muscle physiology, resulting in decreased muscle mass, impaired force generation, delayed recovery from injury, and many more (Cohen et al., 2015). Conclusively, a distracting negative loop leads to adverse effects of MSTN overexpression in muscle tissue.

Moreover, according to our results, GDF11 triggered an even stronger SMAD7 response leading to the conclusion that this same protein may be involved in regulating GDF11 expression. Previous studies explored different regulatory mechanisms of GDF11/ MSTN activity, including Growth and differentiation factor (GDF) -associated serum protein-1 (GASP-1) and GASP-2, an inhibitor of mature MSTN and GDF11 (Lee and Lee, 2013; Walker et al., 2016). However, SMAD7 has not been previously reported in conjunction with GDF11 and muscle tissue, warranting further investigation.

EGR1

Early growth response transcription factor (EGR1) has been subject to reports linked to its negative role in cardiovascular disease progression (Khachigian, 2021), profibrotic response, and physiologic and pathological connective tissue remodeling (Bhattacharyya et al., 2011). Given its manifold role in various tissue, our study discovered the activation of this transcript being weaker with GDF11 compared to MSTN.

MSTN has been shown to directly regulate fibrosis in muscle tissue by stimulating the proliferation of muscle fibroblasts through activation of Smad, p38 MAPK, and Akt pathways (Li et al., 2008). Inhibition of MSTN was demonstrated to reverse muscle fibrosis by reducing fibroblasts’ resistance to the apoptosis (Li et al., 2012). The utilization of GO-term analysis using the ORA, demonstrated positive activation of the p38MAPK and MAP3K cascade, indirectly confirming the results of Li et al. (2008) and pointing towards fibrosis regulation through EGR-1 triggered p38MAPK signaling (^{Figure S2}). Having both physiologic and pathologic functions, EGR-1 should be further studied towards its role in MSTN-provided muscle tissue fibrosis and its activation through GDF11.

HOXA3

Gradual inhibition of the Hoxa3 genes was noted in both GDF11 and MSTN-stimulated C1C12 cells in our study. It also seemed to be inversely correlated with concentrations in both ligands. Turner et al. (2020) did great work describing epigenetic changes in skeletal muscle tissue by investigating the methylation status of various subsets of genes in their recent study. Interestingly, Hoxa3 showed an inverse relationship between gene expression and methylation, demonstrating significantly downregulated expression through hypermethylation in aged cells. Further investigations by this group have also determined that increased levels of endurance exercises were able to reverse methylation and increase gene expression of Hoxa3, practically reversing age-related epigenetic changes in muscle cells (Turner et al., 2020). Another study concluded that GDF11 could help to prevent cardiomyocyte pyroptosis (inflammation-dependent type of programmed cell death) via GDF11/Smad2/3/HOXA3/NLRP3 (Nucleotide-binding oligomerization domain-like receptor pyrin domain-containing 3) signaling pathway, by leading to overexpression of Hoxa3 with subsequent inhibition of NLRP3 (Li et al., 2020). Our study could not identify an increase in HOXA3 gene expression upon exposure to GDF11/ MSTN. In summary, there is little known about GDF11/ MSTN triggered HOXA3 activity in skeletal muscle tissue, therefore, it is our belief that our results should prompt further investigation that might uncover the additional mechanism of action of these two ligands.

One significant limitation in our experimental design may stem from the duration of the experiment. Specifically, the treatment with GDF proteins was conducted for the period of one hour, with sequential harvesting for further analysis. However, even though some biological processes, particularly those related to differentiation, could need extended time frames, the main goal our experiment was the genes expression patterns analysis. Satisfying time for these processes was assured by confirmation of SMAD activation utilizing Western Blot.

Another limitation of our study that is to be mentioned is the absence of validation with qPCR. Past experiences with microarrays suggested the necessity of such validation of genome-scale expressions using RT-qPCR analysis in order to mitigate issues with reproducibility and bias (Zhang et al., 2006; Balázsi and Oltvai, 2007). Confirmation of RNA sequencing results with quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR) is a common practice in molecular biology and genomics research, and it serves a valuable purpose in many cases. However, there are situations where confirmation with qRT-PCR may not be necessary or may not provide substantial additional benefits. High-throughput RNA sequencing platforms have become increasingly accurate and reliable. These platforms are designed to generate precise gene expression data, and many researchers have confidence in the quality of their results, especially when using well-established protocols and bioinformatics pipelines. A number of studies have specifically addressed the correlation between results obtained with RNA seq and qPCR, showing similar results as those published by Everaert et al. (2017). Furthermore, as Coenye T. outlined in his recent article the RNA-seq methods and data analysis approaches are robust and therefore do not always require confirmation through other approaches or methods such as RT qPCR (Coenye, 2021).

Even though our study demonstrated similar direct signaling of GDF11 and MSTN, and confirmed the results of the microarray analysis performed by Egerman et al. (2015), it is still possible that GDF11 and MSTN might have some unique divergent activity on myoblasts and muscle tissue. As several studies have already stressed, the activity of both molecules in different tissues is not only determined through the activity of the mature ligand itself but also through the cleaving activity of proteases that cleave off the pro domains and additionally by modulation and inhibition by different extracellular binding proteins, as well as the availability of specific binding domains (Kondás et al., 2008; Walker et al., 2016). Described mechanisms could, technically speaking, change the outcome of the experiment in vivo since using active ligands and cell culture does not allow any further modification steps and changes that might happen in tissue and be crucial for these ligands’ function.

Furthermore, the actual source of the ligands used was previously already described as a possible important factor that could alter experimental results (Rodgers et al., 2014), implicating that there is an actual rationale in determining expression patterns utilizing ligands from various sites and comparing the downstream activation. Assuming that the source of the ligands could play a vital role gives a high probability of conflicting results if performing experiments with ligands from different vendors.

Conclusions

In conclusion, in part contracting studies highlight the need for further investigation to fully delineated and uncover the functions of GDF11 and MSTN in muscle tissue. Our study did not reveal any statistical difference in gene expression patterns between these two ligands. ORA confirmed the activation of the non-canonical pathway by the enrichment of GO-terms for MAPK (GO:0043410) and JNK activation (GO:0046330). Furthermore, we were able to show similar activation patterns and outline novel relationships between genes not yet described in connection with GDF11/ MSTN in muscle tissue, pointing out important questions for future investigations.

Acknowledgements

This work was made possible by the Wagers Laboratory. We thank Jill Goldstein, Jordan Lewandowski, Kathleen Messemer, Bryan Peacker and Naftali Horwitz for scientific discussions and technical assistance. We would also like to thank the Bauer Core Facility at Harvard University and specifically Emma White and Claire Bailey Hartmann for technical assistance. Furthermore, we thank the German Academic Scholarship Foundation for the financial support. We also acknowledge financial support from the Open Access Publication Fund of UKE - Universitätsklinikum Hamburg-Eppendorf.

References

• Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT et al (2000) Gene ontology: Tool for the unification of biology. The Gene Ontology Consortium . Nat Genet 25:25-29.
• Bajikar SS, Wang C-C, Borten MA, Pereira EJ, Atkins KA and Janes KA (2017) Tumor-suppressor inactivation of GDF11 occurs by precursor sequestration in triple-negative breast cancer. Dev Cell 43:418-435.
• Balázsi G and Oltvai ZN (2007) A pitfall in series of microarrays: The position of probes affects the cross-correlation of gene expression profiles. Methods Mol Biol 377:153-162.
• Bhattacharyya S, Wu M, Fang F, Tourtellotte W, Feghali-Bostwick C and Varga J (2011) Early growth response transcription factors: Key mediators of fibrosis and novel targets for anti-fibrotic therapy. Matrix Biol 30:235-242.
• Biesemann N, Mendler L, Kostin S, Wietelmann A, Borchardt T and Braun T (2015) Myostatin induces interstitial fibrosis in the heart via TAK1 and p38. Cell Tissue Res 361:779-787.
• Bullard SA, Seo S, Schilling B, Dyle MC, Dierdorff JM, Ebert SM, DeLau AD, Gibson BW and Adams CM (2016) Gadd45a protein promotes skeletal muscle atrophy by forming a complex with the protein kinase MEKK4. J Biol Chem 291:17496-17509.
• Chen S, Zhou Y, Chen Y and Gu J (2018) Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34:i884-i890.
• Coenye T (2021) Do results obtained with RNA-sequencing require independent verification? Biofilm 3:100043.
• Cohen TV, Kollias HD, Liu N, Ward CW and Wagner KR (2015) Genetic disruption of Smad7 impairs skeletal muscle growth and regeneration. J Physiol 593:2479-2497.
• Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M and Gingeras TR (2013) STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29:15-21.
• Ebert SM, Dyle MC, Kunkel SD, Bullard SA, Bongers KS, Fox DK, Dierdorff JM, Foster ED and Adams CM (2012) Stress-induced skeletal muscle Gadd45a expression reprograms myonuclei and causes muscle atrophy. J Biol Chem 287:27290-27301.
• Eden E, Navon R, Steinfeld I, Lipson D and Yakhini Z (2009) GOrilla: A tool for discovery and visualization of enriched GO terms in ranked gene lists. BMC Bioinformatics 10:48.
• Egerman MA, Cadena SM, Gilbert JA, Meyer A, Nelson HN, Swalley SE, Mallozzi C, Jacobi C, Jennings LL, Clay I et al (2015) GDF11 increases with age and inhibits skeletal muscle regeneration. Cell Metab 22:164-174.
• Everaert C, Luypaert M, Maag JL V, Cheng QX, Dinger ME, Hellemans J and Mestdagh P (2017) Benchmarking of RNA-sequencing analysis workflows using whole-transcriptome RT-qPCR expression data. Sci Rep 7:1559.
• Forbes D, Jackman M, Bishop A, Thomas M, Kambadur R and Sharma M (2006) Myostatin auto-regulates its expression by feedback loop through Smad7 dependent mechanism. J Cell Physiol 206:264-272.
• Gamer LW, Wolfman NM, Celeste AJ, Hattersley G, Hewick R and Rosen V (1999) A novel BMP expressed in developing mouse limb, spinal cord, and tail bud is a potent mesoderm inducer in Xenopus embryos. Dev Biol 208:222-232.
• Gerardo-Ramírez M, Lazzarini-Lechuga R, Hernández-Rizo S, Jiménez-Salazar JE, Simoni-Nieves A, García-Ruiz C, Fernández-Checa JC, Marquardt JU, Coulouarn C, Gutiérrez-Ruiz MC et al (2019) GDF11 exhibits tumor suppressive properties in hepatocellular carcinoma cells by restricting clonal expansion and invasion. Biochim Biophys Acta Mol Basis Dis 1865:1540-1554.
• Harmon EB, Apelqvist AA, Smart NG, Gu X, Osborne DH and Kim SK (2004) GDF11 modulates NGN3+ islet progenitor cell number and promotes beta-cell differentiation in pancreas development. Development 131:6163-6174.
• Hinken AC, Powers JM, Luo G, Holt JA, Billin AN and Russell AJ (2016) Lack of evidence for GDF11 as a rejuvenator of aged skeletal muscle satellite cells. Aging Cell 15:582-584.
• Janssen I, Shepard DS, Katzmarzyk PT and Roubenoff R (2004) The healthcare costs of sarcopenia in the United States. J Am Geriatr Soc 52:80-85.
• Katsimpardi L, Litterman NK, Schein PA, Miller CM, Loffredo FS, Wojtkiewicz GR, Chen JW, Lee RT, Wagers AJ and Rubin LL (2014) Vascular and neurogenic rejuvenation of the aging mouse brain by young systemic factors. Science 344:630-634.
• Khachigian LM (2021) Early Growth Response-1, an integrative sensor in cardiovascular and inflammatory disease. J Am Heart Assoc 10:e023539.
• Kondás K, Szláma G, Trexler M and Patthy L (2008) Both WFIKKN1 and WFIKKN2 have high affinity for growth and differentiation factors 8 and 11. J Biol Chem 283:23677-23684.
• Kraler S, Balbi C, Vdovenko D, Lapikova-Bryhinska T, Camici GG, Liberale L, Bonetti N, Canestro CD, Burger F, Roth A et al (2023) Circulating GDF11 exacerbates myocardial injury in mice and associates with increased infarct size in humans. Cardiovasc Res 119:2729-2742.
• Lee Y-S and Lee S-J (2013) Regulation of GDF-11 and myostatin activity by GASP-1 and GASP-2. Proc Natl Acad Sci U S A 110:E3713-E3722.
• Li R, Xia J, Zhang XI, Gathirua-Mwangi WG, Guo J, Li Y, McKenzie S and Song Y (2018) Associations of muscle mass and strength with all-cause mortality among US older adults. Med Sci Sports Exerc 50:458-467.
• Li Z, Xu H, Liu X, Hong Y, Lou H, Liu H, Bai X, Wang L, Li X, Monayo SM et al (2020) GDF11 inhibits cardiomyocyte pyroptosis and exerts cardioprotection in acute myocardial infarction mice by upregulation of transcription factor HOXA3. Cell Death Dis 11:917.
• Li ZB, Kollias HD and Wagner KR (2008) Myostatin directly regulates skeletal muscle fibrosis. J Biol Chem 283:19371-19378.
• Li ZB, Zhang J and Wagner KR (2012) Inhibition of myostatin reverses muscle fibrosis through apoptosis. J Cell Sci 125:3957-3965.
• Liu Y, Morley M, Brandimarto J, Hannenhalli S, Hu Y, Ashley EA, Tang WHW, Moravec CS, Margulies KB, Cappola TP et al (2015) RNA-Seq identifies novel myocardial gene expression signatures of heart failure. Genomics 105:83-89.
• Liu Y, Shao L, Chen K, Wang Z, Wang J, Jing W and Hu M (2018) GDF11 restrains tumor growth by promoting apoptosis in pancreatic cancer. Onco Targets Ther 11:8371-8379.
• Loffredo FS, Steinhauser ML, Jay SM, Gannon J, Pancoast JR, Yalamanchi P, Sinha M, Dall’Osso C, Khong D, Shadrach JL et al (2013) Growth differentiation factor 11 is a circulating factor that reverses age-related cardiac hypertrophy. Cell 153:828-839.
• Love MI, Huber W and Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15:550.
• McMahon DK, Anderson PA, Nassar R, Bunting JB, Saba Z, Oakeley AE and Malouf NN (1994) C2C12 cells: Biophysical, biochemical, and immunocytochemical properties. Am J Physiol 266:C1795-802.
• McPherron AC and Lee SJ (1997) Double muscling in cattle due to mutations in the myostatin gene. Proc Natl Acad Sci U S A 94:12457-12461.
• McPherron AC, Lawler AM and Lee SJ (1997) Regulation of skeletal muscle mass in mice by a new TGF-beta superfamily member. Nature 387:83-90.
• McPherron AC, Lawler AM and Lee SJ (1999) Regulation of anterior/posterior patterning of the axial skeleton by growth/differentiation factor 11. Nat Genet 22:260-264.
• Morikawa M, Derynck R and Miyazono K (2016) TGF-β and the TGF-β family: Context-dependent roles in cell and tissue physiology. Cold Spring Harb Perspect Biol 8:a021873.
• Nakashima M, Toyono T, Akamine A and Joyner A (1999) Expression of growth/differentiation factor 11, a new member of the BMP/TGFbeta superfamily during mouse embryogenesis. Mech Dev 80:185-189.
• Papadopoulou SK (2020) Sarcopenia: A contemporary health problem among older adult populations. Nutrients 12:1293.
• Philip B, Lu Z and Gao Y (2005) Regulation of GDF-8 signaling by the p38 MAPK. Cell Signal 17:365-375.
• Poggioli T, Vujic A, Yang P, Macias-Trevino C, Uygur A, Loffredo FS, Pancoast JR, Cho M, Goldstein J, Tandias RM et al (2016) Circulating growth differentiation factor 11/8 levels decline with age. Circ Res 118:29-37.
• Rao MS, Van Vleet TR, Ciurlionis R, Buck WR, Mittelstadt SW, Blomme EAG and Liguori MJ (2019) Comparison of RNA-Seq and microarray gene expression platforms for the toxicogenomic evaluation of liver from short-term rat toxicity studies. Front Genet 9:636.
• Rodgers BD and Eldridge JA (2015) Reduced circulating GDF11 is unlikely responsible for age-dependent changes in mouse heart, muscle, and brain. Endocrinology 156:3885-3888.
• Rodgers BD, Wiedeback BD, Hoversten KE, Jackson MF, Walker RG and Thompson TB (2014) Myostatin stimulates, not inihibits, C2C12 myoblast proliferation. Endocrinology 155:670-675.
• Sinha M, Jang YC, Oh J, Khong D, Wu EY, Manohar R, Miller C, Regalado SG, Loffredo FS, Pancoast JR et al (2014) Restoring systemic GDF11 levels reverses age-related dysfunction in mouse skeletal muscle. Science 344:649-652.
• Siparsky PN, Kirkendall DT and Garrett Jr WE (2014) Muscle changes in aging: Understanding sarcopenia. Sports Health 6:36-40.
• Tamura RE, de Vasconcellos JF, Sarkar D, Libermann TA, Fisher PB and Zerbini LF (2012) GADD45 proteins: central players in tumorigenesis. Curr Mol Med 12:634-651.
• Taylor WE, Bhasin S, Artaza J, Byhower F, Azam M, Willard DHJ, Kull FCJ and Gonzalez-Cadavid N (2001) Myostatin inhibits cell proliferation and protein synthesis in C2C12 muscle cells. Am J Physiol Endocrinol Metab 280:E221-8.
• Thomas M, Langley B, Berry C, Sharma M, Kirk S, Bass J and Kambadur R (2000) Myostatin, a negative regulator of muscle growth, functions by inhibiting myoblast proliferation. J Biol Chem 275:40235-40243.
• Trombetti A, Reid KF, Hars M, Herrmann FR, Pasha E, Phillips EM and Fielding RA (2016) Age-associated declines in muscle mass, strength, power, and physical performance: Impact on fear of falling and quality of life. Osteoporos Int 27:463-471.
• Turner DC, Gorski PP, Maasar MF, Seaborne RA, Baumert P, Brown AD, Kitchen MO, Erskine RM, Dos-Remedios I, Voisin S et al (2020) DNA methylation across the genome in aged human skeletal muscle tissue and muscle-derived cells: The role of HOX genes and physical activity. Sci Rep 10:15360.
• Walker RG, Poggioli T, Katsimpardi L, Buchanan SM, Oh J, Wattrus S, Heidecker B, Fong YW, Rubin LL, Ganz P et al (2016) Biochemistry and biology of GDF11 and myostatin: Similarities, differences, and questions for future investigation. Circ Res 118:1125-1142.
• Walker RG, Czepnik M, Goebel EJ, McCoy JC, Vujic A, Cho M, Oh J, Aykul S, Walton KL, Schang G et al (2017) Structural basis for potency differences between GDF8 and GDF11. BMC Biol 15:19.
• Wang C, Gong B, Bushel PR, Thierry-Mieg J, Thierry-Mieg D, Xu J, Fang H, Hong H, Shen J, Su Z et al (2014) The concordance between RNA-seq and microarray data depends on chemical treatment and transcript abundance. Nat Biotechnol 32:926-932.
• Wang Z, Gerstein M and Snyder M (2009) RNA-Seq: A revolutionary tool for transcriptomics. Nat Rev Genet 10:57-63.
• Wilhelm BT and Landry J-R (2009) RNA-Seq-quantitative measurement of expression through massively parallel RNA-sequencing. Methods 48:249-257.
• Yaffe D and Saxel O (1977) Serial passaging and differentiation of myogenic cells isolated from dystrophic mouse muscle. Nature 270:725-727.
• Yang W, Chen Y, Zhang Y, Wang X, Yang N and Zhu D (2006) Extracellular signal-regulated kinase 1/2 mitogen-activated protein kinase pathway is involved in myostatin-regulated differentiation repression. Cancer Res 66:1320-1326.
• Yokoe T, Ohmachi T, Inoue H, Mimori K, Tanaka F, Kusunoki M and Mori M (2007) Clinical significance of growth differentiation factor 11 in colorectal cancer. Int J Oncol 31:1097-1101.
• Zhang L, Yoder SJ and Enkemann SA (2006) Identical probes on different high-density oligonucleotide microarrays can produce different measurements of gene expression. BMC Genomics 7:153.
• Zhao S, Fung-Leung W-P, Bittner A, Ngo K and Liu X (2014) Comparison of RNA-Seq and microarray in transcriptome profiling of activated T cells. PLoS One 9:e78644.
• Zhou Y, Ni S, Song L, Wang X, Zhang Y and Zhang S (2019) Late-onset administration of GDF11 extends life span and delays development of age-related markers in the annual fish Nothobranchius guentheri. Biogerontology 20:225-239.
• Zimmers TA, Davies MV, Koniaris LG, Haynes P, Esquela AF, Tomkinson KN, McPherron AC, Wolfman NM and Lee S-J (2002) Induction of cachexia in mice by systemically administered myostatin. Science 296:1486-1488.

• Data Availability

  Sequence data was submitted to the European Nucleotide Archive (ENA) and made publicly available under the study PRJEB57932 on the 15th of February 2024. ^{Table S1} in the Supplementary material lists all accession numbers separately.

Supplementary material

The following online material is available for this article:

Figure S1 - Boxplots.

Figure S2 - Vulcano plots.

Table S1 - Accession Numbers.

Edited by

Associate Editor:

Emmanuel Dias Neto

Data availability

Sequence data was submitted to the European Nucleotide Archive (ENA) and made publicly available under the study PRJEB57932 on the 15th of February 2024. ^{Table S1} in the Supplementary material lists all accession numbers separately.

Publication Dates

• Publication in this collection
  15 July 2024
• Date of issue
  2024

History

• Received
  01 Nov 2023
• Accepted
  08 May 2024