Transcriptomic characterization of the molecular mechanisms induced by RGMa during skeletal muscle nuclei accretion and hypertrophy

Background The repulsive guidance molecule a (RGMa) is a GPI-anchor axon guidance molecule first found to play important roles during neuronal development. RGMa expression patterns and signaling 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 from our research group 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 profile 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 nuclei accretion induced by RGMa could be due to their capacity to induce the expression of transcripts related to ‘adherens junsctions’ and ‘extracellular-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, especially from sarcolemma and also those associated with striated muscle cell differentiation. Conclusions These results provide comprehensive knowledge of skeletal muscle transcript changes and pathways in response to RGMa. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08396-w.

identified 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 [1][2][3][4] 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 ossification during skeletal development [5], the suppression of endothelial tube formation [6], and inflammatory responses [7,8].
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 affinity to the type I transmembrane neogenin receptor [9,10], which is known as a guidance receptor for migrating neuronal and mesodermal cells [11][12][13]. This domain also harbours a von Willibrand type D structural domain, containing a GDPH autocatalytic site [14]. The RGMa N-terminal 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 [15]. However, RGMa signaling through integrins has not been reported thus far. Notably, N-RGMa presents high affinity to bone morphogenetic proteins (BMP) ligands, making RGMa (and all the members of this family) a modulator of this important signaling pathway [16][17][18][19]. 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 signaling pathway via activation of Smad 1/5/8. RGMa could also integrate neogenin and BMP signaling cascades [5,[20][21][22]. 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) signaling pathway [23].
In previous works, we found RGMa transcripts in the myogenic and satellite cell precursors in the somites during chicken embryonic development [24] and at the sarcolemma and in the sarcoplasm of adult mice muscle cells [25]. RGMa overexpression in C2C12 cells induced the formation of larger myotubes (hypertrophy) with an increased number of myonuclei (nuclei accretion), while its knockdown resulted in the appearance of smaller cells, with a deficient ability to form multinucleated myotubes [25].
Skeletal muscle cell size is known to be determined by the balance between protein and cellular turnover [26][27][28]. Because of cellular turnover, the skeletal muscle cell grows by myonuclei accretion, in a process mediated by cell fusion. The increase of myonuclei into myofibers leads to muscle mass expansion due to the higher rate of transcription given the nuclear turnover [29]. Muscle nuclei accretion is important not only during embryonic development but also during muscle regeneration [29][30][31][32][33][34][35][36]. 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 [27,29]. Although hypertrophy and nuclei accretion are two distinct processes, they frequently occur together [36], and not all signals involved during the proliferation and differentiation of skeletal musculature are known.
Despite having found that RGMa can induce hypertrophy and nuclear accretion in 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. Our hypothesis is that RGMa can modulate the expression of a number of transcripts in skeletal muscle cells, especially those involved with nuclei accretion and striated muscle cell differentiation. 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 first work to show, through RNA-seq analysis, the transcript targets and molecular profile triggered by RGMa during skeletal muscle differentiation and its possible involvement in multiple functions, including cell fusion and hypertrophy.

Overview of the RNA-seq data and differentially expressed transcripts (DETs)
The quality of the generated sequence database was first evaluated to verify the internal consistency and reproducibility of the replicate samples, as well as the disparity among them. The Pearson correlation coefficient (PCC) of the normalized read-counts 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 coefficients above 0.97 among all correlated samples (Fig. 1B). MA-plot analysis revealed that RGMa treatment modulated gene expression in skeletal muscle cells, with very few of them presenting a drastic effect (Fig. 1C).
Twenty-six DETs were exclusively expressed in RGMatreated myoblasts, and 79 DETs had their expression drastically altered by the treatment (Supplementary  Table 1). Differential expression analysis was also performed in gene level. We found that RGMa could modulate the expression of 1,788 genes (DEGs, p < 0.05, Supplementary Table 2). From the 2,195 DETs, 1,091 were also found as DEGs, meaning that 1,104 (~ 50%) were found as differentially expressed only at the transcription level.
The most drastic effects among the DETs were also observed as a heatmap of transcripts with enriched muscle-associated terms (Fig. 1D). 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).
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 (Ap2a1, and Picalm), and with cell cycle progression and apoptosis (Pcgf1, Kifc3 and Cep164) ( Table 1).

RNA categories among DET
Given the reliability of the transcriptome data, we next classified all 2,195 DETs by RNA biotypes to determine which were the main RNA categories influenced by RGMa treatment. Among the 1,252 DETs that were found to be downregulated, 917 (73.2%) were protein coding, 246 (19.6%) were processed transcripts, 67 (5.35%) were NMDs, and 22 (1.75%) were pseudogenes (Fig. 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 fibre, ' 'myoblast fusion, ' and 'integrin binding' (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 characterized based on the Gene Ontology (GO) terms to identify the pathways that were enriched among the up-and 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 classification 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 signaling pathway' were the most enriched terms of biological processes (Fig. 5A). Cellular components were mostly associated with 'mitochondria protein complex, ' 'actin cytoskeleton, ' 'cytosolic ribosome, ' Table 1 The most drastically altered transcripts in RGMa-treated C2C12 myoblast, during myogenic differentiation 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 identified 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 '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 identified as a guidance clue for axonal growth, RGMa has been identified 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 [24]. In adult muscle, RGMa was found in regions of the sarcolemma and sarcoplasm, with an expression pattern similar to sarcomeric proteins [25]. 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 [25].
However, the molecular mechanisms induced by RGMa during these important muscle phenotypes have not been clarified thus far. RGMa exerts its canonical effects through the type-I transmembrane neogenin receptor [6,7,9,37], 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 [15]. Notably, both signaling 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 signaling pathways, possibly because these receptors do not have RGMa as an exclusive ligand [38]. 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 differentiation 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 technical replicates are considered a sufficient for a reliable quantitative inferential analysis [39]. 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 five RNA biotypes. The most frequent RNA biotype modulated by RGMa was 'protein coding, ' which included ORF containing transcripts. However, a significant 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 [40]. 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 DET (Pou2f1-208, ENS-MUST00000160260.9), as well as the most downregulated DET (Pou2f1-205, ENSMUST00000111427.9) by RGMa treatment in skeletal muscle cells. Although the specific functions of each of these isoforms have not been described thus far, it is known that Pou2F1 is an ubiquitously expressed 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 [41], the regulation of the circadian clock [42], and glycolytic metabolism [43]. In skeletal muscle cells, this transcription factor was associated with the activation of pro-inflammatory immune response in patients with myalgia [44] and with MyHC IIB expression, when associated with MEF2 and the serum response factor (SRF) [45][46][47]. Pou2F1 was also identified on a slow Fig. 6 Muclei accretion and muscle-related enriched terms from the functional analysis of all DETs in response to RGMa. GO enrichment and the network analysis of DETs was performed using the software ClueGO. Terms were selected for network analysis related to nuclei accretion A and to muscle differentiation and structure B. The right-sided hypergeometric test was used in statistical inference, and the Benjamini-Hochberg method was applied for a p-value correlation (p < 0.001). The network was designed using the ForceAtlas2 algorithm and node size represents network centrality which was calculated using Eigenvector Centrality algorithm Fig. 7 Validation of the RNA-seq expression profiles 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 *, significance of p-adj < 0.05 skeletal muscle troponin I promoter in Gaoyou duck skeletal muscle [48]. 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.

RGMa treatment modulated the expression of muscle hypertrophic markers
The protein coding DETs were analyzed to determine how RGMa induces hypertrophic and nuclear accretion 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 [27,[49][50][51][52][53]. RGMa could specifically 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 confirmed 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 [27,54]. TSC1, in a complex with TSC2, is responsible for the negative regulation of mTORC1 signaling, inhibiting the nutrientmediated or growth factor-stimulated phosphorylation of S6K1 and 4EBP1 [53]. Furthermore, TSC1-204 (ENS-MUST00000113870.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 signaling independent of Akt [53,55]. 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 [56]. 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 [55], 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 [50,57].
RGMa could also induce the upregulation of members of the Myocyte Enhancer Factor 2 (Mef2) family, specifically Mef2a-204 (ENSMUSG00000030557.17) and Mef2d-204 (ENSMUSG00000001419.17) isoforms. Mef2 transcription factors activate many muscle-specific growth factor-induced genes and regulate muscle cell differentiation and muscle embryonic development [58][59][60]. Mef2 can also act as a nodal point for remodeling programs in metabolic gene expression, fiber-type switching, and skeletal muscle regeneration [58,59,61]. Mef2a upregulation can also contribute to terminal differentiation and myoblast fusion, which is also consistent with the present GO term analysis and with the RGMa muscle phenotype [25,38]. Mef2a, Mef2c, and Mef2d deleted in combination in satellite cells abolished skeletal muscle regeneration after cardiotoxin injury [59].
Our RNA-seq database suggested other hypertrophic mechanisms that could be regulated by RGMa treatment, including the upregulation of Sirtuin 1 (Sirt1), which is known to regulate protein degradation via FoxO inhibition [62]; the upregulation of Nos1, which interacts with Sirt1 [63]; or the downregulation of genes that promote muscle protein degradation, such as the activating transcription factor 4 (Atf4) [64,65].

RGMa treatment also modulated the expression of genes associated with nuclei accretion
We have also searched for genes associated with myonuclear accretion that were modulated by RGMa treatment in C2C12 cells. Among these, cadherin2 (Cdh2, Myof, for example, is a member of the Ferlin protein family, highly expressed in myoblasts during the prefusion phase of differentiation and in myofibers, especially during regeneration after injury [66][67][68]. It is associated with fusion events and intracellular trafficking in muscle, including myoblast fusion, vesicle traffic, membrane repair, and endocytic recycling [69]. Myh9 and Myh10 are equally fundamental for the positive regulation of cell-cell adhesion and myoblast fusion [70]. Myh9 is known as non-muscle myosin heavy chain IIa (NMMHC-IIA), while Myh10 is the non-muscle myosin heavy chain IIb [71]. 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 [72][73][74]. In skeletal muscle cells, non-muscle myosins drive myoblasts to align and fuse to form multinucleated myotubes [70,75]. The knockdown of these myosins inhibit the change of the myoblast shape, interfering with cell-cell adhesion and fusion [70].
Dab2 plays an important role as a modulator of cellcell interactions, as it is a clathrin adaptor and can mediate integrin signaling [76]. In the musculature, Dab2 was detected during early myogenic differentiation [77,78]. Shang et al. (2020) showed that Dab2 expression is upregulated in C2C12 myoblast during the differentiation in myotubes, and its knockdown resulted in reduced myoblast fusion and fewer myotubes. Besides, Dab2 overexpression could enhance the myotube formation and also restore the myotube differentiation capacity of its knockdown [79].
The calcium voltage-gated channel subunit alpha1 S (Cacna1s-202, ENSMUST00000112068.10) encodes one of the five subunits of the L-type voltage-dependent calcium channel in skeletal muscle cells. In the musculature, calcium is generally related to muscle contraction and muscle relaxation [80][81][82]. However, the regulation of calcium influx into muscle cells plays a critical role in muscle differentiation [82,83]. Intracellular calcium is able to regulate transcription factors necessary for myotube fusion [83,84], while its reduction inhibits myoblast differentiation [85]. The upregulation of Cacna1s in response to RGMa treatment suggests an association with the regulation of intracellular calcium, which is important for the myoblast fusion process and myotube contraction.

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 nuclei accretion 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 myoblast fusion 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 immortalized mouse myoblasts C2C12 (ATCC ® CRL1772 ™ ) was cultured at 37 °C and 5% CO 2 in growth medium (GM), composed of DMEM (Dulbecco's Modified Eagle's Medium) with high glucose and L-glutamine (Gibco), supplemented with 10% fetal 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 2 × 10 4 cells per well in 24-well plates and cultivated in GM at 37 °C and 5% CO 2 . After reaching 90-100% confluency, cells were induced to differentiate in DM for 72 h (Fig. 1A). 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 [38]. 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 amplification grade (Invitrogen) and purified according to the TruSeq Stranded mRNA Sample Prep LS Protocol of Illumina (http:// grcf. jhmi. edu/ hts/ proto cols/ mRNA-Seq_ Sampl ePrep_ 10048 98_D. pdf), using magnetic microspheres for messenger RNA separation. The purified mRNA was fragmented in Illumina buffer. Superscript III (Invitrogen) and oligo(dT) were used for reverse transcription of the first cDNA strand. The second strand was synthesized 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 synthetized cDNA was phosphorylated with T4 PNK (Illumina) and adenylated with Klenow exo (Illumina). Adaptors were bound to cDNA ends, and the samples were purified and selected by size of 200 bp ± 25 bp after fractioning in agarose gel electrophoresis (QIAquick Gel Extraction Kit, QIAGEN). Purified cDNA was quantified by RT-qPCR using adaptor-specific oligonucleotides (Illumina).
Sequencing was performed using HiScanSQ (Illumina),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 quantification analysis was performed based on Salmon (version 0.13.1), an open-source and freely-licensed software (available at https:// github. com/ COMBI NE-lab/ Salmon [86]). 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. genco degen es. org/ mouse/ relea se_ M20. html) was used as a reference, which includes all RNA categories used to classify the transcripts obtained in this work.

Statistical RNA-seq
Statistical analysis was performed using the DESeq2 package of R Bioconductor [87]. 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 significance 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-specific 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 and the network analysis of DETs was performed using the software ClueGO v.2.5.4 [88] and Gephi v 0.9.2 (https:// ojs. aaai. org/ index. php/ ICWSM/ artic le/ view/ 13937). The right-sided hypergeometric test was used to identify overrepresented GO terms and the BenjaminiHochberg method was used for the correction of the p-values (p < 0.001). The Ensembl Transcript ID of the DETs was used as input for ClueGO analysis. Terms were selected for network analysis by related to nuclei accretion (Fig. 6A) and related to muscle differentiation and structure (Fig. 6B). The network was designed using the ForceAtlas2 algorithm and node size represents network centrality which was calculated using Eigenvector Centrality algorithm.
The heatmap graph was obtained using the D3Heatmap package (https:// www. rdocu menta tion. org/ packa ges/ d3hea tmap/ versi ons/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.

Primer design and qPCR
qPCR was performed for the specific upregulated isoforms of 204 and 203,208 and 214, Parva-203, Pou2f1-208 and for the downregulated isoforms of Cep164-170, Kifc3-202 and Pcgf1-201, that were selected due to their importance for muscle phenotypes, as well as by their relevance between the more enriched terms.
The multiline interface (http:// multa lin. toulo use. inra. fr/ multa lin/) was used for the alignment of genes with some specific isoforms up and others downregulated by RGMa. The non-consensus sequences among them were selected to avoid undesired isoforms and the consensus ones were used to obtain an amplicon of up to 250 bp for the chosen isoforms. Primer 3.0 software was used for primer design. Manual primers were designed for small specific strings.
cDNA was synthesized using 1 μg of total RNA following the recommendations of the RevertAid ™ H Minus First Strand cDNA Synthesis kit (Fermentas). 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 final 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) [89] with 95% significance.