Transcriptomic Characterisation of the Molecular Mechanisms Induced by RGMa During Skeletal Muscle Hyperplasia and Hypertrophy

Background: The repulsive guidance molecule a (RGMa) is a GPI-anchor axon guidance molecule rst found to play important roles during neuronal development. RGMa expression patterns and signalling pathways via Neogenin and/or as BMP coreceptors indicated that this axon guidance molecule could also be working in other processes and diseases, including during myogenesis. Previous works have consistently shown that RGMa is expressed in skeletal muscle cells and that its overexpression induces both nuclei accretion and hypertrophy in muscle cell lineages. However, the cellular components and molecular mechanisms induced by RGMa during the differentiation of skeletal muscle cells are poorly understood. In this work, the global transcription expression prole of RGMa-treated C2C12 myoblasts during the differentiation stage, obtained by RNA-seq, were reported. Results: RGMa treatment could modulate the expression pattern of 2,195 transcripts in C2C12 skeletal muscle, with 943 upregulated and 1,252 downregulated. Among them, RGMa interfered with the expression of several RNA types, including categories related to the regulation of RNA splicing and degradation. The data also suggested that RGMa hyperplasia effects could be due to their capacity to induce the expression of transcripts related to cell-cell adhesion, while RGMa effects on muscle hypertrophy might be due to (i) the activation of the mTOR-Akt independent axis and (ii) the regulation of the expression of transcripts related to atrophy. Finally, RGMa induced the expression of transcripts that encode skeletal muscle structural proteins and members of the signalling pathways associated with GEF and Rho/Rac, common secondary signals of skeletal muscle hypertrophy, and the canonical pathway of the RGMa/Neogenin signalling. Conclusions: These results provide comprehensive knowledge of skeletal muscle transcript changes and pathways in response to RGMa. (Qiagen), using the iTaq Universal Sybr Green Supermix (Bio Rad) and 0.4–0.8 µM of each primer for a nal volume of 10 µl. GAPDH was used as a housekeeping gene. The analysis of differential gene expression was performed using REST 2009 (Relative Expression Software Tool, V.2.0.13) software via randomisation tests (Pair Wise Fixed Reallocation Randomisation Test) [91] with 95% signicance.


Background
Repulsive guidance molecule a (RGMa) comprises the rst repulsive glycoprotein member identi ed in the family of repulsive guidance molecules [1]. It was originally identi ed as a repulsive clue in the orientation of axonal growth in the central and peripheral nervous system and as an important target for neuronal survival [2][3][4][5] However, RGMa action domains were found to go beyond the processes related to neurogenesis and could be extended to different processes, including the induction of endochondral ossi cation during skeletal development [6], the suppression of endothelial tube formation [7], and in ammatory responses [8,9].
These diverse functions can be performed by RGMa because it can signal through different receptors and work as a modular protein. The RGMa C-terminal domain (C-RGMa) harbours a GPI-anchor and presents a nity to the type I transmembrane neogenin receptor [10,11], which is known as a guidance receptor for migrating neuronal and mesodermal cells [12][13][14]. This domain also harbours a von Willibrand type D structural domain, containing a GDPH autocatalytic site [15]. The RGMa Nterminal domain (N-RGMa) harbours a signal peptide, an additional neogenin-binding site, and an RGD motif, that is known to be important in cell-cell adhesion processes mediated by integrins [16]. However, RGMa signalling through integrins has not been reported thus far. Notably, N-RGMa presents high a nity to bone morphogenetic proteins (BMP) ligands, making RGMa (and all the members of this family) a modulator of this important signalling pathway [17][18][19][20]. N-RGMa shares the same binding site on the BMP ligand with the ectodomain of the BMP type I receptor A (BMP-R1A), meaning that RGM can induce the BMP canonical signalling pathway via activation of Smad 1/5/8. RGMa could also integrate neogenin and BMP signalling cascades [6, [21][22][23]. Finally, RGMa was recently found to promote astrogliosis and glial scar formation in a rat model of middle cerebral artery occlusion/reperfusion by forming a complex with ALK5 and Smad2/3, which are the main members of the transforming growth factor β1 (TGFβ1) signalling pathway [24].
In previous works, we found RGMa transcripts in the myogenic and satellite cell precursors in the somites during chicken embryonic development [25] and at the sarcolemma and in the sarcoplasm of adult mice muscle cells [26]. RGMa overexpression in C2C12 cells induced the formation of larger myotubes (hypertrophy) with an increased number of myonuclei (hyperplasia), while its knockdown resulted in the appearance of smaller cells, with a de cient ability to form multinucleated myotubes [26].
Skeletal muscle cell size is known to be determined by the balance between protein and cellular turnover [27][28][29]. Because of cellular turnover, the skeletal muscle cell grows by myonuclei accretion (hyperplasia), in a process mediated by cell fusion. The increase of myonuclei into myo bers leads to muscle mass expansion due to the higher rate of transcription given the nuclear turnover [30]. Muscle hyperplasia is important not only during embryonic development but also during muscle regeneration [30][31][32][33][34][35][36][37]. In contrast, because of protein turnover, the skeletal muscle cell grows by upregulating protein synthesis pathways, consequently increasing the level of protein within the muscle tissue [28,30]. Although hypertrophy and hyperplasia are two distinct processes, they frequently occur together [37], and not all signals involved during the proliferation and differentiation of skeletal musculature are known.
Despite having found that RGMa can induce hypertrophy and hyperplasia of skeletal muscle cells cultivated in vitro, the molecular mechanisms that are induced by this axon guidance molecule in these particular cells have not been investigated thus far. In this work, C2C12 cells were treated with RGMa recombinant protein to investigate the molecular mechanisms that are modulated by this axon guidance molecule during myogenic differentiation. This was the rst work to show, through RNAseq analysis, the transcript targets and molecular pro le triggered by RGMa during skeletal muscle differentiation and its possible involvement in multiple functions, including cell fusion and hypertrophy.

Results
Overview of the RNA-seq data and differentially expressed transcripts (DETs) The quality of the generated sequence database was rst evaluated to verify the internal consistency and reproducibility of the replicate samples, as well as the disparity among them. The Pearson correlation coe cient (PCC) of the normalised readcounts revealed a perfect positive linear correlation between all RGMa-treated samples and an extremely strong correlation among the control ones (Fig. 1A). The analysis also revealed a subtle difference between treated and control samples, as there was a positive linear correlation showing Pearson r coe cients above 0.97 among all correlated samples (Fig. 1A). MA-plot analysis revealed that RGMa treatment modulated gene expression in skeletal muscle cells, with very few of them presenting a drastic effect (Fig. 1B).
The expression of 23,855 transcripts could be detected after normalisation, and 2,195 were found as differentially expressed transcripts (DETs, p < 0.05, Fig. 1B, grey dots), with 943 upregulated and 1,252 downregulated by RGMa treatment compared to the control (Fig. 1B, blue dots with Log2(FC) > 0 and Log2(FC) < 0, respectively). Twenty-six DETs were exclusively expressed in RGMa-treated myoblasts, and 79 DETs had their expression drastically altered by the treatment (Supplementary Table 1). The most drastic effects among the DETs were also observed as a heatmap of transcripts with enriched muscle-associated terms (Fig. 1C). The heatmap also allowed the observation that the expression of the majority of the transcripts did not change considerably between the control and treated samples. The 20 most upregulated and 20 most downregulated transcripts by RGMa treatment in C2C12 cells are shown in Table 1.
Notably, another isoform of the Pou2F1 transcription factor (Pou2F1-205, ENSMUST00000111427.9) was found to be one of the most downregulated genes by RGMa treatment (Table 1). Note. The twenty most highly downregulated (Log2(Fold Change) < 0) and twenty most highly upregulated (Log2(Fold Change) > 0) Differentially Expression Transcripts (DET -with a false discovery rate (FDR) < 0,05) modulated in C2C12 cells treated with RGMa during differentiation. Access number and transcript name identi ed in the Ensembl database; log2(FD) < 0 corresponds to the fold change of the downregulation and log2(FD) > 0, of the upregulation of each transcript after RGMa treatment.
Among the others, the downregulation of genes from the same categories included those associated with regulators of muscle mass and structure, such as Tsc1 and Tnpo3, with the formation of clathrin-coated vesicles (Kifc3, Ap2a1, and Picalm), and with cell cycle progression (PCGF1) (  2). Among the 943 upregulated DETs, 786 (83.3%) were protein coding, 115 (12.2%) were processed transcripts, 36 (3.8%) were NMDs, 5 (0.5%) were pseudogenes, and 1 (0.1%) was a TEC (Fig. 2). Overall, this data revealed that most of the RNA biotypes that were modulated by RGMa treatment were ORF-containing RNAs, while the remaining were composed of RNAs mainly associated with the regulation of gene expression, including the NMD category, which was composed of transcripts containing a premature stop codon, and processed transcripts, a category composed of lncRNA, ncRNA, antisense, and intron-retained RNAs.
GO pathway enrichment analysis of the non-coding RNA found as DETs The non-coding RNA found as DETs that were upregulated by RGMa treatment were mostly involved with the 'regulation of RNA splicing,' 'stress bre,' 'myoblast fusion,' and 'integrin biding' (Fig. 3A), while 'regulation of RNA transport' and 'peptide biosynthetic process' were enriched among the downregulated non-coding DETs (Fig. 3B).
GO pathway enrichment analysis of the protein coding RNA found as DETs DETs were characterised based on the Gene Ontology (GO) terms to identify the pathways that were enriched among the upand downregulated transcripts. The enriched GO terms for the protein coding upregulated DETs were mostly related to the following biological processes: 'morphogenesis,' 'metabolism,' and 'developmental regulation of muscle cell' (Fig. 4A). Related to the cellular components, RGMa treatment could induce the upregulation of transcripts associated with 'cytoskeleton,' 'cell projection,' 'endomembrane system,' 'adherens junction,' 'nucleus,' and 'nucleoplasm' (Fig. 4B); and related to molecular function, transcripts were grouped as 'nucleic acid binding,' 'transcription factor binding,' and 'regulation of GTPase' and 'Ras GTPase activity' (Fig. 4C).
A different pattern was found with the classi cation of the protein coding downregulated DETs; terms related to 'metabolism' and 'tissue survival' were the most downregulated after RGMa treatment. 'Purine nucleoside triphosphate metabolic process,' 'peptide biosynthetic process,' 'translation,' and 'positive regulation of apoptotic signalling pathway' were the most enriched terms of biological processes (Fig. 5A). Cellular components were mostly associated with 'mitochondria protein complex,' 'actin cytoskeleton,' 'cytosolic ribosome,' 'proteasome complex,' and 'vesicle coat' (Fig. 5B), while those found to be associated with molecular function were grouped in 'activity of nucleoside-triphosphatase,' 'ATPase,' and 'positive regulation of catalysis' categories ( Fig. 5C).

Discussion
Although originally identi ed as a guidance clue for axonal growth, RGMa has been identi ed as playing roles in a number of different biological processes, including during myogenesis. RGMa transcripts could be found in chicken somites at the origin site of the muscle and satellite cell precursors [25]. In adult muscle, RGMa was found regions of the sarcolemma and sarcoplasm, with an expression pattern similar to sarcomeric proteins [26]. Initial functional studies revealed that RGMa can induce myonuclear accretion and hypertrophy of myotubes, suggesting that this axon guidance molecule might be involved with the mechanisms that modulate skeletal muscle cell size [26]. However, the molecular mechanisms induced by RGMa during these important muscle phenotypes have not been clari ed thus far. RGMa exerts its canonical effects through the type-I transmembrane neogenin receptor [7,8,10,38], but it can also work as a bone morphogenetic protein (BMP) co-receptor, as it shares the same binding site in BMP-R1A with BMP ligands [16]. Notably, both signalling pathways seem to be active in skeletal muscle cells, inducing similar phenotypes in controlling the cell size, but these effects were never investigated in the context of having RGMa as a possible ligand. Using RGMa recombinant proteins in C2C12 cells, we could not clearly elucidate if RGMa effects were induced via neogenin and/or BMP signalling pathways, possibly because these receptors do not have RGMa as an exclusive ligand [39]. For this reason, in this work, the transcriptome of C2C12 cells was sequenced after being treated with RGMa recombinant protein during the late differentiation stage to detect the transcripts that had their expression modulated by this axon guidance molecule during the induction of hypertrophy and hyperplasia of skeletal muscle cells.
A database composed of 23,856 transcripts expressed during C2C12 differentiation was generated. Sequenced biological triplicates from treated and control groups were found to be homogeneous, conferring internal consistency and reproducibility of replicate samples. Three biological replicates are considered a su cient for a reliable quantitative inferential analysis [40]. From these expressed transcripts, 2,195 were modulated by RGMa treatment, with 943 upregulated and 1,252 downregulated.
From this database, it was noted that RGMa was able to modulate the expression of ve RNA biotypes. The most frequent RNA biotype modulated by RGMa was 'protein coding,' which included ORF containing transcripts. However, a signi cant portion of DETs were included in categories involved with the regulation of gene expression, including in the 'nonsense mediated decay' (composed of transcripts with a premature stop codon) and 'processed transcript' (composed of 'retained intron RNA', 'antisense' and 'ncRNA') biotypes. According to Wong et al. (2013), a number of transcripts must be destroyed to permit developmental transitions during differentiation [41]. Therefore, this data suggests that RGMa treatment induced the regulation of the genes that were being expressed during the differentiation stages using these particular molecular mechanisms, allowing the adaptation of these cells to reach terminal differentiation.
Additionally, our analysis also revealed that RGMa could differentially induce the expression of alternatively spliced transcripts. Pou2F1 (Pou Class 2 Homeobox 1, also known as Oct-1) isoforms were found to be the most upregulated DETs (Pou2f1-208, ENSMUST00000160260.9), as well as the most downregulated DETs (Pou2f1-205, ENSMUST00000111427.9) by RGMa treatment in skeletal muscle cells. Pou2F1 is a member of the Pou transcription factor family and is associated with a plethora of processes, including the activation of some snRNA, histone H2B, immunoglobulins, and other housekeeping genes [42], the regulation of the circadian clock [43], and glycolytic metabolism [44]. In skeletal muscle cells, this transcription factor was associated with the activation of pro-in ammatory immune response in patients with myalgia [45] and with MyHC IIB expression, when associated with MEF2 and the serum response factor (SRF) [46][47][48]. Pou2F1 was also identi ed on a slow skeletal muscle troponin I promoter in Gaoyou duck skeletal muscle [49]. In addition to Pou2F1, multiple isoforms for myoferlin (Myof), myosin heavy chain 10 (Myh10), myosin IXB (Myo9b), titin (Ttn), tensin 2 (Tns2), supervillain (Svil), and chromodomain helicase DNA binding protein 2 (Chd2) were also found as modulated by the RGMa treatment in C2C12 cells; these genes are of wide importance for development, differentiation, and maintenance of skeletal muscle cells.
The protein coding DETs were analysed to determine how RGMa induces hypertrophic and hyperplasic effects on skeletal muscle cells.
Our transcriptome database showed that RGMa induced the expression of the mammalian target of rapamycin (mTOR) transcript, which is a common factor from different pathways that culminate with skeletal muscle hypertrophy [28, [50][51][52][53][54]. RGMa could speci cally induce mTOR transcript isoform 202 (ENSMUST00000103221.10), suggesting a new mechanism for this isoform in these cells. The effect of RGMa on mTOR expression was also con rmed by qPCR. mTOR exerts its effects as part of two complexes, termed mTORC1 and mTORC2. Increased mTORC1 activity can positively regulate muscle protein synthesis via S6K1 and also inhibit its negative regulation when working via 4EBP1 [28,55]. TSC1, in a complex with TSC2, is responsible for the negative regulation of mTORC1 signalling, inhibiting the nutrient-mediated or growth factor-stimulated phosphorylation of S6K1 and 4EBP1 [54]. Furthermore, TSC1-204 (ENSMUST00000113870.3) was also highly downregulated by RGMa treatment in skeletal muscle cells. The inhibition of TSC1/2 protein synthesis resulted in rapid activation of mTORC1 signalling independent of Akt [54,56]. The hypertrophic effects observed by RGMa treatment could then be a result of the inhibition of the TSC1 transcript and of the induction of mTOR expression, which are both crucial for muscle growth. Additionally, although the TSC1/2 complex is not physically associated with mTORC1, it is required for mTORC2 activation and consequently, for Akt phosphorylation, in a manner that is independent of its GTPase-activating protein activity toward Rheb [57]. Thus, the inhibition of TSC1 by RGMa suggests that RGMa simultaneously works to prevent mTORC2 activation. The fact that TSC1 inhibition contributes to mTORC1 activation independently of Akt, as well as to mTORC2 inhibition, resulting in the loss of Akt stimulation [56], might explain why Akt was not induced by RGMa in skeletal muscle cells. Our results suggest that mTOR upregulation in response to RGMa is independent of Akt phosphorylation.
Other factors associated with the mTORC pathway were also dysregulated by the RGMa treatment and could contribute to including this axon guidance molecule in an alternative muscle hypertrophic pathway. For example, RGMa could induce the expression of the phospholipase D1 (Pld1-202, ENSMUST00000120834.8) transcript, which was found to be an activator of mTORC1 [51,58].
RGMa could also induce the upregulation of members of the Myocyte Enhancer Factor 2 (Mef2) family, speci cally Mef2a-204 (ENSMUSG00000030557.17) and Mef2d-204 (ENSMUSG00000001419.17) isoforms. Mef2 transcription factors activate many muscle-speci c growth factor-induced genes and regulate muscle cell differentiation and muscle embryonic development [59][60][61]. Mef2 can also act as a nodal point for remodelling programs in metabolic gene expression, bre-type switching, and skeletal muscle regeneration [59,60,62]. Mef2a upregulation can also contribute to terminal differentiation and myoblast fusion, which is also consistent with the present GO  Myof, for example, is a member of the Ferlin protein family, highly expressed in myoblasts during the pre-fusion phase of differentiation and in myo bers, especially during regeneration after injury [67-69]. It is associated with fusion events and intracellular tra cking in muscle, including myoblast fusion, vesicle tra c, membrane repair, and endocytic recycling [70].
Myh9 and Myh10 are equally fundamental for the positive regulation of cell-cell adhesion and myoblast fusion [71]. Myh9 is known as non-muscle myosin heavy chain IIa (NMMHC-IIA), while Myh10 is the non-muscle myosin heavy chain IIb [72]. These myosins are expressed in most cell types, working as motor proteins in a variety of processes requiring contractile force, such as cytokinesis, cell migration, polarisation and adhesion, maintenance of cell shape, and signal transduction [73][74][75]. In skeletal muscle cells, non-muscle myosins drive myoblasts to align and fuse to form multinucleated myotubes [71,76]. The knockdown of these myosins inhibit the change of the myoblast shape, interfering with cell-cell adhesion and fusion [71].
Dab2 plays an important role as a modulator of cell-cell interactions, as it is a clathrin adaptor and can mediate integrin signalling [77].

Conclusion
The current work allowed us to unravel some molecular mechanisms that were altered in skeletal muscle cells after treatment with RGMa, especially those associated with muscle hyperplasia and hypertrophy. Our analysis suggested that RGMa induced cell hypertrophy via (i) upregulation of hypertrophic markers, (ii) downregulation of inhibitors of hypertrophic pathways, (iii) downregulation of transcripts related to the positive regulation of muscle atrophy, and (iv) upregulation of transcripts that negatively regulate atrophy. At the same time, transcripts associated with known hyperplasic pathways were also found to be modulated by RGMa, mainly those related to cell-cell adhesion pathways.
Our results provide comprehensive knowledge of skeletal muscle transcriptional changes and pathways in response to RGMa treatment.

Cell culture and differentiation
The lineage of immortalised mouse myoblasts C2C12 (ATCC® CRL1772™) was cultured at 37°C and 5% CO 2 in growth medium (GM), composed of DMEM (Dulbecco's Modi ed Eagle's Medium) with high glucose and L-glutamine (Gibco), supplemented with 10% foetal bovine serum (FBS, Gibco) and 1% penicillin, streptomycin, and amphotericin B solution (Gibco). Myogenic differentiation was induced in differentiation medium (DM), composed of DMEM, supplemented with 2% horse serum (Gibco) and 1% penicillin, streptomycin, and amphotericin B. For growth or differentiation conditions, the medium was replaced every 2 days.
RGMa recombinant protein treatment C2C12 cells were seeded at 2x10 4 cells per well in 24-well plates and cultivated in GM at 37°C and 5% CO 2 . After reaching 90-100% con uency, cells were induced to differentiate in DM for 72 h. DM was then replaced with fasting medium (FM), composed of DMEM supplemented with 0.2% FBS, and cells were incubated at the same conditions for 3 h. Subsequently, C2C12 were treated with 50 ng/ml mouse RGMa recombinant protein (R&D Systems) in FM and incubated for an additional 48 h, as previously described [39]. The recombinant protein was omitted in the control samples.
Total RNA isolation and cDNA synthesis Cells were harvested in TriReagent (Sigma Aldrich) as pools of three wells in triplicate. Total RNA isolation was performed according to the manufacturer's instructions. Sample integrity, purity, and concentration were evaluated by electrophoresis in 1% agarose gel and in NanoDrop® ND-1000 UV/Vis Spectrophotometer, respectively.
The quality of the total RNA was also evaluated in a Bioanalyzer (Agilent) before being submitted to sequencing. Values for RNA integrity number (RIN) ranging from 8 to 10 were considered suitable for RNA-seq.

RNA-seq library preparation and next-generation sequencing (NGS)
For cDNA library construction, 2 µg of total RNA were treated with 1U of DNaseI ampli cation grade (Invitrogen) and puri ed according to the Illumina protocol (http://grcf.jhmi.edu/hts/protocols/mRNA-Seq_SamplePrep_1004898_D.pdf), using magnetic microspheres for messenger RNA separation. The puri ed mRNA was fragmented in Illumina buffer. Superscript III (Invitrogen) and oligo(dT) were used for reverse transcription of the rst cDNA strand. The second strand was synthesised using the enzymes RNase H and DNA Polymerase I (Illumina). Molecule ends were treated with T4 DNA Polymerase and Klenow DNA Polymerase (Illumina), making them blunt. The 3' end of the synthetised cDNA was phosphorylated with T4 PNK (Illumina) and adenylated with Klenow exo (Illumina). Adaptors were bound to cDNA ends, and the samples were puri ed and selected by size of 200 bp ± 25 bp after fractioning in agarose gel electrophoresis (QIAquick Gel Extraction Kit, QIAGEN). Puri ed cDNA was quanti ed by RT-qPCR using adaptor-speci c oligonucleotides (Illumina).
Sequencing was performed using Illumina HiSeq 2500, according to the manufacturer's recommendations, and using the paired-end reads protocol. Each sample was sequenced until it reached around 34 million reads/library.

Mapping RNA-seq data
Transcript quanti cation analysis was performed based on Salmon (version 0.13.1), an open-source and freely-licensed software (available at https://github.com/COMBINE-lab/Salmon [87]). Raw reads were used as an input to quantify transcripts in mapping-based mode. The current version of the mouse transcriptome (available at https://www.gencodegenes.org/mouse/release_M20.html) was used as a reference.
Statistical RNA-seq Statistical analysis was performed using the DESeq2 package of R Bioconductor [88]. An adjusted p-value with a false discovery rate (FDR) correction (Benjamini and Hochberg, 1995) of 5% was calculated and used to control false-positive signi cance in transcript expression variation. Log2(fold change) > 0 and log2(fold change) < 0 were selected as the threshold to show an increase or decrease in transcript expression of treated groups relative to the control group.
Transcript expression pattern and RNA-seq quality analysis Transcript-speci c normalisation was performed to remove disparities in the base means correlations and to eliminate the noise of transcripts with low expression.
Normalised transcripts were plotted in MA form using the DESeq2 package to generate a scatter plot of log2 fold changes < 0 and > 0 versus the mean of normalised counts of transcripts, considering DE those with FDR < 0.05. The correlation of each sample and the clustering of the treated and control groups was performed by calculating the PCC of normalised read-counts.
Functional RNA-seq analysis The GO enrichment analysis of DETs was performed using the plugin ClueGO v. 2.5.4 [89] for Cytoscape v3.7.1 [90]. The node colours represent the functional groups, and the node size represents the term enrichment signi cance. Only the most signi cant term in the group was labelled. The edges connecting the nodes were based on the Kappa statistic (Kappa Score Threshold of 0.4). The right-sided hypergeometric test was used to identify overrepresented GO terms, and the Benjamini-Hochberg method was used for the correction of the p-values (p < 0.001). The functionally grouped network of enriched categories for DETs was annotated for GO terms (cellular component, biological process, molecular function, and immune system process) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) terms. The Ensembl Transcript ID of the DETs was used as input for ClueGO analysis.
The heatmap graph was obtained using the D3Heatmap package (https://www.rdocumentation.org/packages/d3heatmap/versions/0.6.1.2), using ID ensemble transcripts as an input (of clue go output for muscle associated terms) and the correlated base mean expression.
Before primer design, the multiline interface (http://multalin.toulouse.inra.fr/multalin/) was used for the alignment of different isoforms of the same gene. The non-consensus sequences were selected to obtain an amplicon of up to 250 bp. Primer 3.0 software was used for primer design. Manual primers were designed for small speci c strings.
qPCR was performed in the Rotor-Gene RT-qPCR system (Qiagen), using the iTaq Universal Sybr Green Supermix (Bio Rad) and 0.4-0.8 µM of each primer for a nal volume of 10 µl. GAPDH was used as a housekeeping gene.    Functional analysis of the non-protein coding RNA differentially regulated by RGMa treatment. For this analysis, we considered upregulated DETs that do not encode proteins. Pie analysis of the GO enrichment, showing the most frequent terms, including cellular component, biological process, molecular function, and immune system process, and KEGG GO terms that were (A) upregulated and (B) downregulated. The right-sided hypergeometric test was used in statistical inference, and the Benjamini-Hochberg method was applied for a p-value correlation (p < 0.05). The analysis was conducted using the plugin ClueGO (v.2.5.4) for Cytoscape (v3.7.1).

Figure 4
Functional analysis of the protein coding RNA upregulated by RGMa. For this analysis, we considered the DETs that encode proteins that were found to be upregulated (FC > 1) by the treatment with RGMa, compared to the control. (A-C) Pie chart analysis of the three GO categories used to classify the upregulated protein coding transcripts. The right-sided hypergeometric test was used in statistical inference, and the Benjamini-Hochberg method was applied for a p-value correlation (p < 0.0001). The analysis was conducted using the plugin ClueGO (v.2.5.4) for Cytoscape (v3.7.1).

Figure 5
Functional analysis of the protein coding RNA downregulated by RGMa. For this analysis, we considered the DETs that encode proteins and were found to be downregulated (FC < 1) in the RGMa treated group, compared to the control one. (A-C) Pie chart analysis of the three GO categories for downregulated DETs. (D) Functionally grouped network of enriched categories for expressed transcripts, annotated for 'biological process,' 'cellular component,' and 'molecular function' GO terms. The rightsided hypergeometric test was used in statistical inference, and the Benjamini-Hochberg method was applied for a p-value correlation (p<0.001). The analysis was conducted using the plugin ClueGO (v.2.5.4) for Cytoscape (v3.7.1).

Figure 6
Muscle-related enriched terms from the functional analysis of all DETs in response to rcRGMa. Functionally grouped network of enriched categories for all types of DET annotated as 'cellular component,' 'biological process,' and 'molecular function' GO terms. For this analysis, we considered the DETs that were upregulated (A) and downregulated (B) in the RGMa-treated group, compared to the control group. Only the most signi cant term in the group was labelled. GO terms are represented as nodes, and the node size represents the term enrichment signi cance. The edges connecting the nodes are based on the Kappa statistic (Kappa Score Threshold of 0.4), which measures the overlap of shared genes between terms. The right-sided hypergeometric test was used in statistical inference, and the Benjamini-Hochberg method was applied for a p-value correlation (p < 0.05). The analysis was conducted using the plugin ClueGO (v.2.5.4) for Cytoscape (v3.7.1).

Figure 7
Validation of the RNA-seq expression pro les by qPCR. A subset of twelve DETs that were upregulated and downregulated by RGMa treatment during muscle differentiation were used to validate the obtained RNA-seq expression data. Transcripts were selected by their expression and their known association with muscle hyperplasic or hypertrophic phenotypes. Expression patterns indicate agreement between the two methods and *, signi cance of p-adj < 0.05.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.