ABSTRACT
Skeletal muscle stem cells (MuSCs) are recognised as functionally heterogeneous. Cranial MuSCs are reported to have greater proliferative and regenerative capacity when compared with those in the limb. A comprehensive understanding of the mechanisms underlying this functional heterogeneity is lacking. Here, we have used clonal analysis, live imaging and single cell transcriptomic analysis to identify crucial features that distinguish extraocular muscle (EOM) from limb muscle stem cell populations. A MyogeninntdTom reporter showed that the increased proliferation capacity of EOM MuSCs correlates with deferred differentiation and lower expression of the myogenic commitment gene Myod. Unexpectedly, EOM MuSCs activated in vitro expressed a large array of extracellular matrix components typical of mesenchymal non-muscle cells. Computational analysis underscored a distinct co-regulatory module, which is absent in limb MuSCs, as driver of these features. The EOM transcription factor network, with Foxc1 as key player, appears to be hardwired to EOM identity as it persists during growth, disease and in vitro after several passages. Our findings shed light on how high-performing MuSCs regulate myogenic commitment by remodelling their local environment and adopting properties not generally associated with myogenic cells.
INTRODUCTION
Genetic and transcriptomic studies have shown that the muscle stem cell (MuSC) population in any particular anatomical location is heterogeneous. Some subsets are more prone to self-renewal or differentiation, and differ in their transplantation efficiency, stem cell-niche interactions, metabolism and resistance to stress upon activation (Barruet et al., 2020; Chakkalakal et al., 2014; Dell'Orso et al., 2019; Dumont et al., 2015; Gayraud-Morel et al., 2012; Hernando-Herraez et al., 2019; De Micheli et al., 2020; de Morree et al., 2019; Ono et al., 2012; Rocheteau et al., 2012; Scaramozza et al., 2019; Tierney et al., 2018; Der Vartanian et al., 2019; Yartseva et al., 2020; Yennek et al., 2014). Despite this diversity, MuSCs share common functions that are essential for growth and repair (Lepper, 2011; Murphy et al., 2011; Sambasivan et al., 2011). Myogenic commitment and differentiation involving Myod (Myod1) and Myog, respectively, occur in response to injury or growth factors in culture media after isolation of Pax7-expressing MuSCs (Evano and Tajbakhsh, 2018; Zammit et al., 2004).
Unexpectedly, MuSCs in different anatomical locations were found to be programmed with distinct upstream transcription factors (TFs) before acquiring myogenic identity (Gopalakrishnan et al., 2015; Harel et al., 2009; Kelly et al., 2004; Sambasivan et al., 2009; Tajbakhsh et al., 1997). Extraocular muscles (EOMs) are derived from unsegmented cranial mesoderm and are regulated by distinct transcription factors and signalling molecules compared with the somite-derived limb and trunk muscles (Grimaldi and Tajbakhsh, 2021; Michailovici et al., 2015; Sambasivan et al., 2011). For example, mice lacking the transcription factor Pitx2 do not form EOMs, whereas other cranial and somite derived muscles are unaffected (Diehl et al., 2006; Gage et al., 1999; Zacharias et al., 2011). As some muscle subsets, such as the EOMs, are preferentially spared in muscular dystrophies and during ageing (Emery, 2002; Formicola et al., 2014; Verma et al., 2017), and their MuSCs are functionally more robust in terms of proliferation and engraftment efficiency (Stuelsatz et al., 2015; Hebert et al., 2013), intrinsic properties of MuSCs or myofibers might determine their differential sensitivity to disease (Randolph and Pavlath, 2015; Terry et al., 2018; Kallestad et al., 2011). Notably, single-cell RNA-sequencing (scRNA-seq) analysis has identified the thyroid hormone signalling pathway as a key factor preventing senescence entry of EOM MuSCs in dystrophin deficient (DMD) rats (Taglietti et al., 2023). Moreover, profiling of EOM MuSCs after heterotopic transplantation into limb muscle showed that, despite significant transcriptional changes, ∼10% of EOM-specific genes persisted in the ectopic niche (Evano et al., 2020). This finding suggests that cell-autonomous regulation of MuSC properties predominates to a certain extent in the EOMs. Yet the existence of deeply rooted gene regulatory networks within MuSC subsets and their role in maintaining anatomically distinct phenotypes remains largely unexplored.
Using droplet-based scRNA-seq, we have investigated the transcriptional states that govern the outperformance of MuSC subsets. We identified key regulators that confer distinct mesenchymal-like features to activated EOM MuSCs and extensive expression of extracellular matrix (ECM) components. Activated EOM MuSCs are less prone to differentiation and thus appear to be in a more stem-like state that is maintained in vitro upon passages and in vivo during growth, through a specific set of co-regulated transcription factors.
RESULTS
Temporal heterogeneity in myogenic commitment to differentiation among MuSCs
We isolated EOM and tibialis anterior (TA) MuSCs by FACS using Tg:Pax7-nGFP mice (Sambasivan et al., 2009) and plated them at the same density to follow their proliferation and differentiation dynamics (Fig. 1A,B). Quantification of the number of total and proliferative nuclei [5-Ethynyl-2′-deoxyuridine (EdU) uptake, 2 h pulse] at D3, D4 and D5 showed a 2.9 to 4.4-fold increase in EOM MuSCs compared with TA (Fig. 1C). However, the percentage of total proliferative cells remained unchanged in both conditions at every time point (Fig. 1D), suggesting that the division probability is equal for EOM and TA. Yet, EOM cultures had twice the percentage of PAX7+ cells and fewer fusogenic myoblasts compared with TA at later timepoints (D10, Fig. 1E,F).
Functional differences between EOM and TA MuSCs after activation. (A) Experimental scheme. Adult MuSCs isolated from Tg:Pax7-nGFP mice were plated in growth media (GM) and pulsed with EdU before fixation and immunostaining at day (D) 3, D4, D5 and D10. Media was changed every second day. (B) Immunofluorescence for MF20, PAX7 and EdU detection as in A. Scale bars: 100 µm. (C,D) Bar plots of the total number of nuclei and EdU+ nuclei represented as EOM/TA fold-change (C) and percentage of EdU+ nuclei/total number of nuclei (D) at D3, D4 and D5. N>3 independent experiments with n≥3 mice per experiment. (E) Quantification of the percentage of PAX7+ cells/total number of nuclei at D10. n≥3 mice. (F) Quantification of the Fusion Index % at D10. n=4 mice. (G) Experimental scheme. Live imaging of adult EOM and TA MuSCs from Tg:Pax7-nGFP;MyogntdTom from D3 in vitro. (H) Representative overlaid DIC and red fluorescence channel images at selected time points from experiments outlined in G. Scale bar: 25 μm. (I) Percentage of TOM+ cells over total cell number at each timepoint. n=4 mice. (J) RT-qPCR on in vitro-activated adult EOM and TA MuSCs for Pax7, Myod, Myog and p21 normalised to Rpl13. n=8 mice. (K) Experimental scheme. Adult MuSCs from Tg:Pax7-nGFP;MyogntdTom mice were cultured for 5 days and re-isolated based on GFP and TOM fluorescence for further analysis. (L) RT-qPCR for key myogenic markers of cells isolated as in K. n=4 mice. Data are mean±s.e.m. Two-tailed unpaired Student's t-test; ns, non-significant (P>0.05), *P<0.05, **P<0.01, ***P<0.005, ****P<0.0001. GM, growth media; DM, differentiation media.
Functional differences between EOM and TA MuSCs after activation. (A) Experimental scheme. Adult MuSCs isolated from Tg:Pax7-nGFP mice were plated in growth media (GM) and pulsed with EdU before fixation and immunostaining at day (D) 3, D4, D5 and D10. Media was changed every second day. (B) Immunofluorescence for MF20, PAX7 and EdU detection as in A. Scale bars: 100 µm. (C,D) Bar plots of the total number of nuclei and EdU+ nuclei represented as EOM/TA fold-change (C) and percentage of EdU+ nuclei/total number of nuclei (D) at D3, D4 and D5. N>3 independent experiments with n≥3 mice per experiment. (E) Quantification of the percentage of PAX7+ cells/total number of nuclei at D10. n≥3 mice. (F) Quantification of the Fusion Index % at D10. n=4 mice. (G) Experimental scheme. Live imaging of adult EOM and TA MuSCs from Tg:Pax7-nGFP;MyogntdTom from D3 in vitro. (H) Representative overlaid DIC and red fluorescence channel images at selected time points from experiments outlined in G. Scale bar: 25 μm. (I) Percentage of TOM+ cells over total cell number at each timepoint. n=4 mice. (J) RT-qPCR on in vitro-activated adult EOM and TA MuSCs for Pax7, Myod, Myog and p21 normalised to Rpl13. n=8 mice. (K) Experimental scheme. Adult MuSCs from Tg:Pax7-nGFP;MyogntdTom mice were cultured for 5 days and re-isolated based on GFP and TOM fluorescence for further analysis. (L) RT-qPCR for key myogenic markers of cells isolated as in K. n=4 mice. Data are mean±s.e.m. Two-tailed unpaired Student's t-test; ns, non-significant (P>0.05), *P<0.05, **P<0.01, ***P<0.005, ****P<0.0001. GM, growth media; DM, differentiation media.
We then performed clonal analysis of freshly isolated EOM MuSCs (Fig. S1A) and showed that they displayed a 7-fold higher clonal capacity (mean 1315 cells/clone) compared with those from the TA (186 cells/clone) (Fig. S1B, T1), in agreement with a previous study (Stuelsatz et al., 2015). Surprisingly, the higher clonogenic properties of EOM persisted even when they were pre-amplified for 2 days in vitro, where they yielded 3656 cells/clone compared with 297 cells/clone for TA (12-fold difference; Fig. S1B, T2). Given that EOM and TA MuSCs exhibited a similar proliferative rate in bulk cultures, the clonal data suggest that the greater cellular output of EOM MuSCs is probably due to delayed differentiation. As Myog expression is followed by cell-cycle withdrawal (Andrés and Walsh, 1996; Benavente-Diaz et al., 2021; Guo et al., 1995), a delay in its expression could allow for sustained expansion of EOM MuSCs.
To test this hypothesis, we assessed the differentiation dynamics of EOM and TA MuSCs in vitro by live imaging using the Tg:Pax7-nGFP;MyogntdTom mouse line (Benavente-Diaz et al., 2021; Sambasivan et al., 2009), which allowed simultaneous isolation of MuSCs and monitoring of the onset of Myog expression by a nuclear tdTomato (TOM) reporter. Myogenic cells were tracked continuously from D3 and TOM intensity scored (Fig. 1G). TA myogenic cells initiated reporter gene expression earlier than EOM cells (Fig. 1H). The percentage of TOM+ cells sharply increased in the TA from 5% at 80 h (D3+8 h) to 40% at 96 h (D4) post-plating. Significantly, only ∼7% of EOM cells were TOM+ by 96 h in vitro (Fig. 1I). Western blot (Fig. S1C,D) data were in agreement with the live imaging analysis. For example, expression of the cell cycle inhibitor p21 (Cdkn1a), which regulates cell cycle exit (Guo et al., 1995; Zhang et al., 1999), and expression of phosphorylated p38, which promotes MEF2 transcriptional activity to initiate differentiation (Rugowska et al., 2021), were significantly downregulated in activated EOM MuSCs compared with TA. RT-qPCR showed higher levels of Myog and p21 transcripts in TA compared with EOM MuSCs already at D2 (Fig. 1J). As such, these data suggest that EOM progenitors are less prone to myogenic commitment.
Finally, to test whether EOM cells have a global impairment in myogenic differentiation, we used Tg:Pax7-nGFP;MyogntdTom mice to isolate the upstream (PAX7+/GFP+) and committed (MYOG+/TOM+) fractions at D5 in vitro (Fig. 1K, Fig. S1E). No differences were seen by RT-qPCR among EOM and TA TOM+ cells for Myod, Myog, the fusogenic markers Mymx (myomixer) and Mymk (myomaker) (Sampath et al., 2018), and embryonic myosin heavy chain (Myh3) indicating that differentiation of EOM cells was not generally impeded (Fig. 1L). However, plating EOM and TA TOM+ myoblasts at high density in differentiation media revealed a lower fusion index and higher percentage of EOM myoblasts incorporating EdU on their last division (Fig. S1F-H) (Benavente-Diaz et al., 2021). Therefore, a general delay in the differentiation of EOM MuSCs appears to contribute to the higher number of myogenic cells they give rise to (Fig. S1I).
Distinct signature of EOM MuSCs upon activation
To understand the phenotypic differences between EOM and TA MuSCs, we performed scRNA-seq of in vitro-activated MuSCs using the 10x Chromium platform (Fig. 2A). Unsupervised clustering divided cells into two clusters per sample (Fig. 2B) that were annotated as progenitors (Prog) or differentiating (Diff), based on myogenic marker expression (Fig. 2B,C). Notably, EOM progenitors exhibited reduced expression of Myod mRNA and protein across the entire population at D3 and D5 in vitro, thereby suggesting a lower propensity to differentiate (Fig. 2B, Fig. S2A).
Single cell transcriptome signatures of activated EOM and TA MuSCs. (A) scRNA-seq pipeline for in vitro-activated adult MuSCs. (B) UMAP visualisation of EOM progenitors, EOM differentiating, TA progenitors and TA differentiating clusters (left). Violin plots of the expression of myogenic markers Myf5 and Pax7 (progenitor), Myod (committed) and Myog (differentiating) in each cluster (right). (C) Expression plots of myogenic markers. (D) Heatmap representing the top differentially expressed genes in each cluster and expression levels across all cells. (E,F) Reactome pathway (E) and GO molecular function (F) network analysis on top 100 DEGs of each cluster. Pie charts represent relative contribution of each cluster to this ontology term.
Single cell transcriptome signatures of activated EOM and TA MuSCs. (A) scRNA-seq pipeline for in vitro-activated adult MuSCs. (B) UMAP visualisation of EOM progenitors, EOM differentiating, TA progenitors and TA differentiating clusters (left). Violin plots of the expression of myogenic markers Myf5 and Pax7 (progenitor), Myod (committed) and Myog (differentiating) in each cluster (right). (C) Expression plots of myogenic markers. (D) Heatmap representing the top differentially expressed genes in each cluster and expression levels across all cells. (E,F) Reactome pathway (E) and GO molecular function (F) network analysis on top 100 DEGs of each cluster. Pie charts represent relative contribution of each cluster to this ontology term.
Analysis of differentially expressed genes (DEGs) showed that each cell state had a distinct transcriptional pattern (Fig. 2D, Table S1). As expected, TA progenitors expressed a defined Hox signature (Evano et al., 2020) and inhibitors of differentiation (Id1 and Id2; Jen et al., 1992; Kumar et al., 2009). EOM and TA differentiating cells shared part of their signature (Myog, Mymx and Myl4). However, EOM progenitors displayed markers that were not previously noted in MuSCs, such as Mgp, Bgn, Col1a2 and Acta2 (smooth muscle actin) (Fig. 2D).
We then performed reactome analysis of DEGs (Fig. 2E, Fig. S2B, Table S2). As suggested by the DEG heatmap, TA progenitors shared common pathways with the differentiating clusters. Multiple pathways involved in ECM organisation characterised EOM progenitors, including ‘ECM proteoglycans’, ‘crosslinking of collagen fibrils’ and ‘integrin cell surface interactions’. EOM progenitors were also specifically enriched in Pdgfrb expression, a receptor for platelet-derived growth factor, and with proteins related to insulin-like growth factor binding and integrin signalling (Fig. 2F, Fig. S2C, Table S3). Thus, these analyses uncovered a non-canonical signature in EOM progenitors upon activation.
EOM MuSC molecular identity is partially maintained across cell states
We then compared the signature of EOM and TA myogenic progenitors in vitro with the quiescent counterparts in vivo by scRNA-seq immediately after isolation by FACS (Fig. 3A,B). Although both populations displayed consistent expression of MuSC markers, Myod expression was restricted to TA MuSCs (Fig. 3C,D), suggesting an earlier activation in response to the isolation procedure (van den Brink et al., 2017; Machado et al., 2017; van Velthoven et al., 2017). Differentiation markers such as Myog were absent in both populations (Fig. 3C,D). Visualisation of the top 25 most variably expressed genes highlighted distinct transcriptional programs of EOM and TA quiescent clusters (Fig. 3E, Table S4). Moreover, microarray analysis showed that EOM quiescent MuSCs also stand apart from those in other cranial-mesoderm derived muscles, such as the oesophagus and masseter (Fig. S3A).
Molecular signature of EOM and TA quiescent MuSCs. (A) scRNA-seq pipeline for quiescent adult MuSCs. (B) UMAP visualisation of EOM and TA quiescent clusters. (C,D) Violin plot (C) and expression plot (D) of Myf5, Pax7, Myod and Myog. (E) Heatmap representing the top differentially expressed genes in each cluster and expression levels across all cells. (F,G) Venn diagrams of the overlapping differentially expressed genes between quiescent and global in vitro-activated datasets for EOM (F) and TA (G) MuSCs. (H,I) Expression plots of selected EOM (H) and TA (I) MuSC markers exclusive to quiescence (Unique Q), common to quiescent and activated states (Conserved), and exclusive to activation (Unique Act) from the analysis in F and G.
Molecular signature of EOM and TA quiescent MuSCs. (A) scRNA-seq pipeline for quiescent adult MuSCs. (B) UMAP visualisation of EOM and TA quiescent clusters. (C,D) Violin plot (C) and expression plot (D) of Myf5, Pax7, Myod and Myog. (E) Heatmap representing the top differentially expressed genes in each cluster and expression levels across all cells. (F,G) Venn diagrams of the overlapping differentially expressed genes between quiescent and global in vitro-activated datasets for EOM (F) and TA (G) MuSCs. (H,I) Expression plots of selected EOM (H) and TA (I) MuSC markers exclusive to quiescence (Unique Q), common to quiescent and activated states (Conserved), and exclusive to activation (Unique Act) from the analysis in F and G.
Interestingly, we found several conserved genes between the quiescent and activated scRNA-seq states at each anatomical location, whereas others were unique to the quiescent or activated cell states (Fig. 3F-I). EOM genes unique to the quiescent state included Tshr, encoding the thyroid-stimulating hormone receptor, previously recognised in bulk RNAseq as ‘EOM-resistant’, i.e. retained after engraftment into the limb (Evano et al., 2020), and shown to prevent senescence in EOM MuSCs of DMD rats (Taglietti et al., 2023). Matrix Gla protein (Mgp), a crucial regulator of angiogenesis in multiple organs (Kida and Yamaguchi, 2022), was exclusively upregulated in EOM cells upon activation (Fig. 3H). Among the conserved EOM genes across cell states, we identified Pitx2, a major upstream regulator of EOM development (Gage et al., 1999), Fos, which labels a subset of limb MuSCs with enhanced regenerative capacity (Almada et al., 2021), and Igfbp7, a specific marker of quiescent MuSCs (Fukada et al., 2007) that is upregulated upon exercise (Chen et al., 2020) (Fig. 3F, Fig. S3B,C). The EOM signature also featured Foxc1, which is involved in ocular development (Smith et al., 2000) and is reported to regulate the balance between myogenic and vascular lineages within somites (Lagha et al., 2009; Mayeuf-Louchart et al., 2016) (Fig. S3B). Additionally, the EOM common signature encompassed several ECM components and regulators (e.g. Bgn, Loxl1, Col1a2 and Col6a1), TFs associated with fibrosis and connective tissue development [e.g. Foxp1 (Grimaldi et al., 2022; Shao and Wei, 2018), Egr1 (Havis and Duprez, 2020)] and Acta2 (alpha smooth muscle actin), a marker of smooth muscle, fibroadipogenic progenitors (Joe et al., 2010; Uezumi et al., 2010) and smooth muscle-mesenchymal cells [SMMCs (Giordani et al., 2019)] (Fig. 3F, Fig. S3B,C).
In the TA, Lbx1, Vgll2 and Hox genes were found to be common signatures throughout cell states (Fig. 3G,I; Fig. S3B,C). Lbx1 is a homeobox transcription factor required for the migration of myogenic progenitor cells to the limbs (Brohmann et al., 2000; Gross et al., 2000) and Vgll2 (previously Vito-1) is a key co-factor of the myogenic differentiation program (Günther et al., 2004; Maeda et al., 2002). Lbx1, Vgll2 and several Hox genes had been already identified in limb MuSCs and in whole muscles in the adult (Evano et al., 2020; Honda et al., 2017; Terry et al., 2018). Altogether, our analysis revealed a closer molecular overlap for EOM MuSCs across cell states (quiescence versus in vitro activation), including several TFs and ECM markers that were distinct from TA MuSCs.
Activation of EOM MuSCs involves extensive ECM remodelling
Given the well-established role of ECM synthesis and remodelling on MuSC proliferation and self-renewal (Baghdadi et al., 2018; Rayagiri et al., 2018; Tierney et al., 2016; Yin et al., 2013), we examined these features using the Matrisome database (Matrisome DB), which compiles in silico and experimental data on ECM constituents (Naba et al., 2015). We identified components of the Matrisome DB present in our single-cell dataset and created a global score for ECM component expression (Fig. 4A) (see Materials and Methods). Significantly, the EOM progenitor cluster had the highest matrisome score (Fig. 4A) and expressed the highest number of genes for each matrisome category (Fig. 4B).
EOM MuSC activation is accompanied by ECM remodelling. (A) Violin plots of matrisome scores for each cluster of the in vitro-activated scRNA-seq MuSC dataset. (B) Dot plot visualisation of differentially expressed matrisome components. ECM, extracellular matrix; A-P, affiliated proteins; PG, proteoglycans. (C) Immunofluorescence for fibronectin (FN1), collagen I (COLI), collagen IV (COLIV) and PDGFRβ. Representative of N>2 different experiments. Scale bar: 25 μm. (D) Western blot showing expression of CAV1, CAVIN1, MGP, MMP2, SPARC, IGFBP7 and H3 for normalisation. Cells from n=3 mice pooled per lane. (E,F) Violin plot (E) and UMAP (F) of Pdgfrb expression on the in vitro-activated scRNA-seq dataset. (G) Experimental scheme. MuSCs from adult Tg:Pax7-nGFP mice were activated in vitro, re-isolated based on their PDGFRβ expression and plated for further analysis. (H) Representative images of PDGFRβ-positive and -negative EOM MuSC fractions stained for MYOD, MYOG and EdU as in G. Scale bar: 50 μm. (I-L) Quantification of number of cells/cm2 (I) and percentage of EdU+ (J), MYOD+ (K) and MYOG+ (L) cells. n=3 mice. Data are mean±s.e.m. Two-tailed unpaired Student's t-test; ns, non-significant (P>0.05), **P<0.01, ***P<0.005.
EOM MuSC activation is accompanied by ECM remodelling. (A) Violin plots of matrisome scores for each cluster of the in vitro-activated scRNA-seq MuSC dataset. (B) Dot plot visualisation of differentially expressed matrisome components. ECM, extracellular matrix; A-P, affiliated proteins; PG, proteoglycans. (C) Immunofluorescence for fibronectin (FN1), collagen I (COLI), collagen IV (COLIV) and PDGFRβ. Representative of N>2 different experiments. Scale bar: 25 μm. (D) Western blot showing expression of CAV1, CAVIN1, MGP, MMP2, SPARC, IGFBP7 and H3 for normalisation. Cells from n=3 mice pooled per lane. (E,F) Violin plot (E) and UMAP (F) of Pdgfrb expression on the in vitro-activated scRNA-seq dataset. (G) Experimental scheme. MuSCs from adult Tg:Pax7-nGFP mice were activated in vitro, re-isolated based on their PDGFRβ expression and plated for further analysis. (H) Representative images of PDGFRβ-positive and -negative EOM MuSC fractions stained for MYOD, MYOG and EdU as in G. Scale bar: 50 μm. (I-L) Quantification of number of cells/cm2 (I) and percentage of EdU+ (J), MYOD+ (K) and MYOG+ (L) cells. n=3 mice. Data are mean±s.e.m. Two-tailed unpaired Student's t-test; ns, non-significant (P>0.05), **P<0.01, ***P<0.005.
We then validated some matrisome candidate genes by immunofluorescence or western blot on in vitro-activated EOM and TA MuSCs (Fig. 4C,D). EOM progenitors expressed higher levels of fibronectin (FN1), which promotes MuSC expansion in a cell-autonomous manner (Bentzinger et al., 2013), collagen I (COLI), a major component of the fibrotic ECM (Dulauroy et al., 2012) that has been shown to suppress differentiation of C2C12 cells (Alexakis et al., 2007), and collagen IV (COLIV), which is secreted by MuSCs as well as myoblasts and fibroblasts in culture (Baghdadi et al., 2018; Kühl et al., 1984) (Fig. 4C). In addition, EOM progenitors also displayed differential expression of PDGFRβ (Fig. 4C), a gene identified on the EOM molecular functions (Fig. S2C). PDGFRβ is a tyrosine-kinase receptor commonly expressed by mesenchymal cells and pericytes (Hellström et al., 1999; Levéen et al., 1994; Soriano, 1994). Interestingly, PDGFRβ and Acta2 are markers of smooth muscle-mesenchymal cells (SMMCs), a recently reported subpopulation of Itga7+ Vcam− Pdgfrb+ Acta2+ cells present in adult muscle that exhibits myogenic potential and promotes MuSC engraftment after transplantation (Giordani et al., 2019).
Western blot analysis confirmed that EOM progenitors produce higher levels of caveolin 1 (CAV1) and CAVIN1 (Fig. 4D), which are co-expressed in caveolae and downregulated upon differentiation of rhabdomyosarcoma cells (Faggi et al., 2015). CAV1 is a marker of quiescent and activated mouse MuSCs (Gnocchi et al., 2009), whereas in human, the CAV1+ MuSC subset is associated with ECM organisation, resistance to activation and increased engraftment capacity (Barruet et al., 2020). SPARC, MGP and IGFBP7 were also upregulated in EOM activated MuSCs (Fig. 4D) and their activity appears to be context dependent as they promote or suppress proliferation in different cell types (Ahmad et al., 2017; Artico et al., 2021; Cho et al., 2000; Jing et al., 2019; Li et al., 2020; Melouane et al., 2018). Given the large number of matrisome genes characterising EOM progenitors, we examined the expression of MMP2, a matrix remodelling protein (Gonçalves et al., 2022) whose activation increases the proportion and mobility of MuSCs (Mu et al., 2010). As expected, an enrichment of the active form of MMP2 was observed in EOM progenitors. TA progenitors instead displayed an enrichment of the pro-MMP2 or latent form (Fig. 4D).
In parallel, we assessed to what extent activated EOM MuSCs resemble other cells in skeletal muscle displaying mesenchymal features (Fig. S4A, Table S5). Notably, EOM progenitors displayed a higher score for SMMCs (Giordani et al., 2019), FAPs (Oprescu et al., 2020), myotendinous junction B myonuclei (Kim et al., 2020), Twist2+ cells (Liu et al., 2017), fetal MuSCs (Tierney et al., 2016), developing limb connective tissues (Esteves de Lima et al., 2021) and the skeletal muscle mesenchyme subpopulation identified in human fetal limb (Xi et al., 2020). In contrast, TA progenitors displayed a higher score for myogenic commitment and differentiation (Fig. S4A).
Finally, we focused on PDGFRβ to validate the potential link between mesenchymal EOM features and in vitro expansion capacity. Flow cytometry analysis corroborated the scRNA-seq and immunofluorescence data, where PDGFRβ was enriched in EOM samples (Fig. S4B,C). To distinguish the functional properties of PDGFRβ+ and PDGFRβ− cells, we pulsed them with EdU (Fig. 4G). Notably, PDGFRβ+ cells had a significantly higher proliferative capacity (Fig. 4H-J). Although the percentage of MYOD+ cells was unchanged between the two cell populations, MYOG was more abundant in the PDGFRβ− fraction (Fig. 4K,L). Therefore, the PDGFRβ+ myoblasts are characterised by a higher proliferative potential and decreased differentiation status. Altogether, our analysis shows that EOM progenitors have an unusual transcriptome profile, express a wide range of ECM-related factors and harbour a mesenchymal signature.
EOM transcriptome profile is associated with a unique transcription factor network
To further investigate myoblast heterogeneity, we inferred single cell regulatory networks using pySCENIC (Aibar et al., 2017; Van de Sande et al., 2020), where co-expression patterns and transcription factor binding motifs identify ‘regulons’ (transcription factors and putative targets) (Fig. 5A, Table S6). As expected, regulons associated with myoblast differentiation, such as Myod, Myog and Mef2 family genes, were specifically active in differentiated cells of the EOM and TA (Fig. 5B). Of note, the top five regulons of TA progenitors were also found to be active in EOM progenitors, whereas the top five regulons of EOM progenitors were unique to this cluster (Fig. 5B). EOM progenitors displayed unique regulons involved in connective tissue and ECM remodelling, including Egr1 (Havis and Duprez, 2020) and Creb3l1, a downstream effector of thyroid hormone signalling (García et al., 2017; Sampieri et al., 2019). Top regulons of EOM progenitors are also involved in cell proliferation [Foxc1 (Yang et al., 2017), Sox4 (Moreno, 2019), Fos (Almada et al., 2021), Klf6 (Dionyssiou et al., 2013) and Ebf1 (Györy et al., 2012)], commitment into endothelial and smooth muscle fates [Foxc1 (Han et al., 2017; Mayeuf-Louchart et al., 2016; Whitesell et al., 2019; Yang et al., 2017)] and differentiation into mesenchymal lineages [Ebf1 (El-Magd et al., 2021; Jimenez et al., 2007; Pagani et al., 2021)]. Of note, some EOM regulon TFs were already enriched in quiescent EOM MuSCs, and their expression was upregulated upon activation (Fig. S5A,B).
Distinct gene regulatory networks underlie EOM and TA activation dynamics. (A) Representative top regulon activity for each cluster of the in vitro-activated adult scRNA-seq MuSC dataset overlaid onto UMAP representation. (B) Heatmap of the top five regulons in each cluster with activity level in each cell. Myod_(+) appears twice as it is a regulon of both EOM Diff and TA Diff subpopulations. (C) Transcription factor network highlighting the top regulons of each cluster and common regulons (35 regulons maximum). Proximity of the nodes in the network indicates a higher number of shared edges, highlighting core modules.
Distinct gene regulatory networks underlie EOM and TA activation dynamics. (A) Representative top regulon activity for each cluster of the in vitro-activated adult scRNA-seq MuSC dataset overlaid onto UMAP representation. (B) Heatmap of the top five regulons in each cluster with activity level in each cell. Myod_(+) appears twice as it is a regulon of both EOM Diff and TA Diff subpopulations. (C) Transcription factor network highlighting the top regulons of each cluster and common regulons (35 regulons maximum). Proximity of the nodes in the network indicates a higher number of shared edges, highlighting core modules.
We then built a network restricted to transcription factors (Fig. 5C), where each node (circle) is an active transcription factor and each edge (distance between nodes) is an inferred regulation between two transcription factors (Grimaldi et al., 2022). When placed in a force-directed environment (see Materials and Methods), these nodes aggregated based on the number of shared edges, thereby highlighting associated and co-regulating transcription factor modules. Strikingly, the transcription factors of the most specific regulons of each cluster preferentially organised as tightly related modules (Fig. 5C). In agreement with our previous analyses, known co-regulating transcription factors in differentiated cells (Myod, Myog, Mef2a, Mef2c and Myf6) formed a tight module (Fig. 5C). The TA module was composed of genes required for limb embryonic development (Hox genes and Lbx1) (Gross et al., 2000; Swinehart et al., 2013) (Fig. 5C, Fig. S5B), where Hoxa11 is a determinant of embryonic limb identity (Zakany and Duboule, 2007), and HoxA and HoxC clusters are signatures of adult TA MuSCs (Evano et al., 2020; Yoshioka et al., 2021).
The EOM module included Foxc1, Egr1, Creb3l1, Dmrta2, Sox4, Fos and Egr1, together with Pax7 and Hes1, which support MuSC quiescence (Baghdadi et al., 2018; Mourikis et al., 2012; Olguin and Olwin, 2004; Relaix and Zammit, 2012). However, this network persists in proliferating EOM progenitors, as assessed by the higher levels of PAX7 and CCND1 (cyclin D1), together with EOM TFs (FOXC1, EBF1 and CREB3L1) at the protein level (Fig. S5C,D). Although Fos and Egr1 act as stress signatures after tissue dissociation (Machado et al., 2021), in our dataset, the expression of these genes (i.e. the StressIndex) correlated with Pax7 expression and anti-correlated with Myod (Fig. S5E,F). Regression of the StressIndex or removal of these genes did not alter the general aspect of the data (Fig. S5G), also pointing to a role for these genes in EOM progenitor maintenance.
EOM features are present during the growth phase in vivo and retained upon passages in vitro
Given that activated EOM MuSCs possess an unusual transcriptomic state and delayed myogenic commitment, we examined these features in MuSCs isolated from Tg:Pax7-nGFP;MyogntdTom after passages in vitro (Fig. 6A, Fig. S6A,B). Although the total cell number was reduced with passages for both EOM and TA (Fig. S6A,B), the normalised cellular output was consistently higher for the EOM and correlated with a lower TOM/GFP ratio (Fig. 6B,C). RT-qPCR and protein analysis showed that expression of Pax7 and Hey1, a bHLH transcription factor that is required for MuSC maintenance (Noguchi et al., 2019), EOM specific regulon TFs (Foxc1, Sox4, Ebf1 and Creb3l1) and genes identified by the matrisome or molecular functions (Bgn, Sparc, Igfbp2, Igfbp7 and Pdgfrb) were all retained or even increased in whole EOM cell populations (a mixture of GFP+ and TOM+ cells) after several passages (Fig. 6D-G), and specifically in the GFP+ fraction (Fig. S6C-E). These results show that EOM MuSCs cells retain a cell-autonomous non-canonical signature that is hard-wired even after extended cell culture and that is not present in TA MuSCs.
EOM MuSC properties are retained despite several passages in vitro. (A) Scheme of isolation and passages of EOM and TA adult MuSCs. Cells from Tg:Pax7-nGFP;MyogntdTom were cultured for 3 days (Act) and the entire wells passaged (P1 and P2) every 3 days. (B) Normalised cell number ratios from the experiment in A. The total number of cells per well was counted at day 5 (D5) and upon passages (P), and normalised on the averaged EOM cell number in all wells at each time point. n=3 mice. (C) Ratio of TOM+/GFP+ cells per time point obtained through FACS analysis of Tg:Pax7-nGFP;MyogntdTom EOM and TA adult MuSCs at D5, P1 and P2. n≥3 mice. (D) RT-qPCR at D3 (Act, activated), P1 and P2 on whole-cell populations for key myogenic markers, regulon TFs and matrisome genes identified in EOM progenitors on the in vitro activated scRNA-seq dataset. Cells isolated from adult Tg:Pax7-nGFP mice by FACS. n=4 mice. (E,F) Western blot for TFs (E) and matrisome proteins (F) in the total P2 population from EOM and TA MuSCs isolated from adult Tg:Pax7-nGFP mice. Actin and tubulin were used for normalisation of protein loading. n=3 mice. (G) Densitometric analysis of the western blots in E and F. Data are mean±s.e.m. Two-tailed unpaired Student's t-test; ns, non-significant (P>0.05), *P<0.05, **P<0.01, ***P<0.005.
EOM MuSC properties are retained despite several passages in vitro. (A) Scheme of isolation and passages of EOM and TA adult MuSCs. Cells from Tg:Pax7-nGFP;MyogntdTom were cultured for 3 days (Act) and the entire wells passaged (P1 and P2) every 3 days. (B) Normalised cell number ratios from the experiment in A. The total number of cells per well was counted at day 5 (D5) and upon passages (P), and normalised on the averaged EOM cell number in all wells at each time point. n=3 mice. (C) Ratio of TOM+/GFP+ cells per time point obtained through FACS analysis of Tg:Pax7-nGFP;MyogntdTom EOM and TA adult MuSCs at D5, P1 and P2. n≥3 mice. (D) RT-qPCR at D3 (Act, activated), P1 and P2 on whole-cell populations for key myogenic markers, regulon TFs and matrisome genes identified in EOM progenitors on the in vitro activated scRNA-seq dataset. Cells isolated from adult Tg:Pax7-nGFP mice by FACS. n=4 mice. (E,F) Western blot for TFs (E) and matrisome proteins (F) in the total P2 population from EOM and TA MuSCs isolated from adult Tg:Pax7-nGFP mice. Actin and tubulin were used for normalisation of protein loading. n=3 mice. (G) Densitometric analysis of the western blots in E and F. Data are mean±s.e.m. Two-tailed unpaired Student's t-test; ns, non-significant (P>0.05), *P<0.05, **P<0.01, ***P<0.005.
Next, we assessed whether EOM progenitor features were already present in their activated state in vivo, for example, in fetal development and postnatal stages, where extensive muscle growth and myogenic cell expansion occurs (Gattazzo et al., 2020; Relaix and Zammit, 2012), or whether these features were acquired in adulthood only upon reactivation from the quiescent state in vitro. First, we isolated EOM and TA MuSCs at E18.5 and P21, and showed that here too, EOM cells were less differentiated upon activation in vitro (Fig. S7A). Then we examined activated MuSCs directly in vivo by RT-qPCR. Given that GFP protein persists during myogenic commitment (Sambasivan et al., 2009) and this might introduce a bias in the initial populations, we isolated the GFP+/TOM− fraction from Tg:Pax7-nGFP;MyogntdTom mice at P7-P10 (Fig. 7A, Fig. S7B). Similar to the in vitro scRNA-seq data, RT-qPCR revealed significantly lower levels of Myog and higher transcript levels for EOM-specific regulon TFs (Foxc1, Ebf1, Sox4 and Creb3l1) and matrisome components on the activated EOM progenitors (Fig. 7B). Together with our live-imaging (Fig. 1G-I), FACS (Fig. S1E) and scRNA-seq data (Fig. 2), this analysis suggests that EOM MuSCs repress myogenic commitment and appear to maintain a more ‘stem-like’ state upon activation. This property appears to be conserved from development to adulthood.
EOM activated MuSC features are present during postnatal growth and can be modulated by FOXC1. (A) Scheme of isolation of EOM and TA GFP+/TOM− cells from Tg:Pax7-nGFP;MyogntdTom from postnatal day (P)7-P10 mice. (B) RT-qPCR for key myogenic markers, regulon TFs and matrisome genes identified in EOM progenitors on the in vitro scRNA-seq dataset. n>3 mice. (C,D) Immunostaining for FOXC1, tdTomato (TOM) and dystrophin on cryosections from EOM and TA muscles isolated from P10 (C) and adult (D) Pax7CreERT2; R26tdTom mice. White arrowheads indicate PAX7+FOXC1+ cells. Black arrowheads indicate PAX7+FOXC1− cells. Asterisks indicate myonuclei. Scale bars: 50 µm. (E) Quantification of the experiment in D. n=4 mice. (F) Scheme of lentivirus transduction for FOXC1 gain-of-function (GOF) experiments on adult TA MuSCs. The transduced cells (mCherry+) were re-isolated by FACS or directly analysed in mixed cultures. (G,H) FACS plots of adult Tg:Pax7-nGFP TA MuSCs transduced with control (G) or GOF FOXC1 (H) virus at D5. Both viruses carry an mCherry reporter. (I) Western blot of transduced adult TA MuSCs for mCherry and FOXC1. Vinculin was used for normalisation of protein loading. n≥3 mice. (J) RT-qPCR of transduced adult TA MuSCs for myogenic markers, regulon TFs and key ECM proteins and regulators identified in EOM progenitors. n≥3 mice. (K,K′) Immunostaining of transduced adult TA MuSCs for FOXC1 together with EdU detection and Cell Mask Blue (CM-Blue, cytoplasmic stain). Scale bar: 100 µm. (L) Immunostaining of transduced adult TA MuSCs grown in differentiation media for 48 h and stained for MF20. Scale bar: 100 µm. (M) Quantification of the Fusion Index % for the images in L. n=4 mice. (N-U) Quantification of cellular properties of adult TA MuSCs in pure (N-Q) or mixed (R-U) cultures containing transduced and non-transduced cells. (N,R) Violin plots showing the single cell distribution of FOXC1 fluorescence intensity. (O-Q,S-U) Bar plots of quantification of nucleus area (O,S), total number of nuclei (P,T) and percentage of EdU+ nuclei (Q,U). N-Q, n≥4 mice; R-U, n≥3 mice. Striped bar in T indicates total output in the mixed culture irrespective of FOXC1 expression. (V-X) Quantification of live imaging of adult TA MuSCs in mixed cultures containing transduced and non-transduced cells. Bright-field and red fluorescent channel images were acquired at four different timepoints (+12, +24, +36, +48 h) post-transduction. Quantification of the transduction efficiency (V) at 12 h time point (t) and cell density (W,X) at all timepoints from movie images. n=3 mice. (Y) Proposed model of the role of FOXC1 upon overexpression in TA MuSCs (TA GOF in pure and mixed cultures) and endogenously in EOM MuSCs. Data are mean±s.e.m. Two-tailed unpaired Student's t-test, except in S,U, where P-values were obtained using one-way ANOVA with Tukey's post hoc test; ns, non-significant (P>0.05), *P<0.05, **P<0.01,***P<0.005, ****P<0.0001. GM, growth media; DM, differentiation media.
EOM activated MuSC features are present during postnatal growth and can be modulated by FOXC1. (A) Scheme of isolation of EOM and TA GFP+/TOM− cells from Tg:Pax7-nGFP;MyogntdTom from postnatal day (P)7-P10 mice. (B) RT-qPCR for key myogenic markers, regulon TFs and matrisome genes identified in EOM progenitors on the in vitro scRNA-seq dataset. n>3 mice. (C,D) Immunostaining for FOXC1, tdTomato (TOM) and dystrophin on cryosections from EOM and TA muscles isolated from P10 (C) and adult (D) Pax7CreERT2; R26tdTom mice. White arrowheads indicate PAX7+FOXC1+ cells. Black arrowheads indicate PAX7+FOXC1− cells. Asterisks indicate myonuclei. Scale bars: 50 µm. (E) Quantification of the experiment in D. n=4 mice. (F) Scheme of lentivirus transduction for FOXC1 gain-of-function (GOF) experiments on adult TA MuSCs. The transduced cells (mCherry+) were re-isolated by FACS or directly analysed in mixed cultures. (G,H) FACS plots of adult Tg:Pax7-nGFP TA MuSCs transduced with control (G) or GOF FOXC1 (H) virus at D5. Both viruses carry an mCherry reporter. (I) Western blot of transduced adult TA MuSCs for mCherry and FOXC1. Vinculin was used for normalisation of protein loading. n≥3 mice. (J) RT-qPCR of transduced adult TA MuSCs for myogenic markers, regulon TFs and key ECM proteins and regulators identified in EOM progenitors. n≥3 mice. (K,K′) Immunostaining of transduced adult TA MuSCs for FOXC1 together with EdU detection and Cell Mask Blue (CM-Blue, cytoplasmic stain). Scale bar: 100 µm. (L) Immunostaining of transduced adult TA MuSCs grown in differentiation media for 48 h and stained for MF20. Scale bar: 100 µm. (M) Quantification of the Fusion Index % for the images in L. n=4 mice. (N-U) Quantification of cellular properties of adult TA MuSCs in pure (N-Q) or mixed (R-U) cultures containing transduced and non-transduced cells. (N,R) Violin plots showing the single cell distribution of FOXC1 fluorescence intensity. (O-Q,S-U) Bar plots of quantification of nucleus area (O,S), total number of nuclei (P,T) and percentage of EdU+ nuclei (Q,U). N-Q, n≥4 mice; R-U, n≥3 mice. Striped bar in T indicates total output in the mixed culture irrespective of FOXC1 expression. (V-X) Quantification of live imaging of adult TA MuSCs in mixed cultures containing transduced and non-transduced cells. Bright-field and red fluorescent channel images were acquired at four different timepoints (+12, +24, +36, +48 h) post-transduction. Quantification of the transduction efficiency (V) at 12 h time point (t) and cell density (W,X) at all timepoints from movie images. n=3 mice. (Y) Proposed model of the role of FOXC1 upon overexpression in TA MuSCs (TA GOF in pure and mixed cultures) and endogenously in EOM MuSCs. Data are mean±s.e.m. Two-tailed unpaired Student's t-test, except in S,U, where P-values were obtained using one-way ANOVA with Tukey's post hoc test; ns, non-significant (P>0.05), *P<0.05, **P<0.01,***P<0.005, ****P<0.0001. GM, growth media; DM, differentiation media.
Foxc1 marks the EOM MuSC lineage and plays a role in progenitor cell maintenance
We then focused on Foxc1 for further analysis, as it is one of the top regulons and DEG of the activated and quiescent scRNA-seq dataset (Fig. 5A,B, Fig. S3B,C, Fig. S5B,D), and it was previously identified as a DEG and top regulon of EOM progenitors in the early embryo (E11.5; Grimaldi et al., 2022). Notably, bulk RNA-seq (Terry et al., 2018) also showed higher Foxc1 expression in entire EOMs compared with other adult muscle groups (Fig. S7C). As such, Foxc1 is a good candidate for determining EOM properties throughout the myogenic lineage continuum and across developmental states.
First, we validated FOXC1 protein expression in vivo and in vitro. As expected, immunostaining on tissue sections of E12.5 Myf5Cre;R26tdTom embryos showed expression of FOXC1 in EOM myogenic progenitors but not in limb, back or masseter muscles (Fig. S7D,E). Analysis of Pax7CreERT2;R26tdTom mice showed FOXC1 expression in postnatal (P10) and adult EOM MuSCs (Fig. 7C-E) but not in TA muscles. FOXC1 was also strongly expressed in adult EOM MuSC cultures (Fig. S7F-H). FOXC1 was expressed in about 75% of PAX7+ cells at D5, expression was progressively lost during differentiation and the majority of the FOXC1+ cells incorporated EdU at every time point (Fig. S7I,J). Thus, differences in Foxc1 expression between the EOM and TA appear to arise during development and persist upon activation.
As differences in Myog mRNA between EOM and TA were already evident at D2 after in vitro activation (Fig. 1J), we used this early timepoint to functionally validate a potential role for Foxc1 in progenitor maintenance. Thus, we silenced Foxc1 using siRNA in EOM MuSCs immediately after FACS (Fig. S7K). The total cell number and percentage of PAX7+ cells were not changed 2 days after silencing, despite effective downregulation of FOXC1 at the protein level and transcript level (Fig. S7L-M). Yet, a 2.7-fold increase in Myog expression was detected by RT-qPCR after silencing (Fig. S7M′). Given that the effect of siRNAs was only transient, we transduced EOM MuSCs with lentiviruses expressing different short-hairpin RNAs (shRNA) against Foxc1 to examine later timepoints (Fig. S7N). Immunofluorescence and EdU uptake at D5 showed an efficient depletion of FOXC1 protein and concomitant severe reduction in the total cell number and percentage of EdU+ cells (Fig. S7O-P). Although lower number of myotubes were formed in shRNA FOXC1 TA cells, the fusion index was not significantly altered (Fig. S7Q,R).
To assess whether Foxc1 would confer some EOM features to the TA, we overexpressed this gene in TA MuSC cells (Fig. 7F) using gain-of-function (GOF) FOXC1 lentiviruses carrying mCherry as a reporter, which allowed the re-isolation of the transduced cells (Fig. 7G,H). Notably, there was robust upregulation of FOXC1 (15-30×) at D5 (Fig. 7I,J,K′,N) and of many Foxc1 direct targets, including matrisome components (e.g. Sparc, Pdgfrb and Fbn1) and EOM-specific regulon TFs (e.g. Ebf1, Creb3l1 and Egr1) (Fig. 7J, Table S7). Properties of re-isolated GOF cells (Fig. 7K-Q) included a larger cellular and nuclear area (Fig. 7K,K′,O), and less fusogenic potential (Fig. 7L,M). Unexpectedly, FOXC1 overexpression also resulted in downregulation of Pax7 and Myod, fewer cells/well and reduced EdU uptake (Fig. 7J,P,Q). As these phenotypes might be the result of strong overexpression and cells acquiring a strong mesenchymal character, we thus investigated ‘mixed cultures’, which contain transduced (FOXC1+mCherry+) and non-transduced (FOXC1−mCherry−) TA MuSCs (Fig. 7F,R-X). This strategy allowed us to assess the effect of secreted factors produced by FOXC1+ cells on non-transduced cells. As expected for mixed cultures, FOXC1 fluorescence intensity showed a bimodal distribution, with one population being as the negative control and the other expressing high levels of the protein (Fig. 7R). Interestingly, the total number of nuclei in the mixed GOF wells was higher than controls (Fig. 7T), which contrasts with the result in pure GOF cultures (Fig. 7P). Moreover, although FOXC1+ cells showed higher nuclear area and reduced EdU uptake, as was seen in pure GOF cultures, we observed an augmentation in the percentage of proliferative non-transduced FOXC1− cells in the mixed cultures (compare FOXC1-positive and -negative subpopulations in Fig. 7S,U with Fig. 7O,Q).
Finally, we performed live imaging on FOXC1 GOF mixed cultures (Fig. 7F). When starting with control and GOF cultures with similar transduction efficiencies (88% and 80% mCherry+ cells at 12 h, respectively) (Fig. 7V), we found that the density of mCherry+ cells in GOF wells was reduced to half of that of control wells by 48 h (Fig. 7W, Movies 1 and 2), with a concomitant increase in the density of mCherry− cells (Fig. 7X, Movies 1 and 2). Therefore, our data suggest a role for the FOXC1-induced secretome in cell proliferation of adjacent non-transduced cells (Fig. 7Y). Importantly, overexpression of Creb3l1, another EOM regulon TF, in TA MuSCs (Fig. S7S) did not show alterations in Pax7, Myod, most matrisome genes examined or cellular properties (Fig. S7S-Z). Altogether, our findings support the notion that physiological levels of Foxc1 in EOM MuSCs allow their maintenance in a progenitor-like state, likely through secreted factors.
Transcription dynamics expose EOM and TA disparities in progenitor-state maintenance
We then performed more detailed bioinformatic analysis of the matrisome by correlating the number of matrisome-driving regulons between EOM and TA MuSCs upon activation (Fig. 8A-C, Table S8). EOM MuSCs consistently regulated a higher number of matrisome genes than TA cells (Fig. 8A). The ratio of the number of regulations of matrisome genes between activated EOM and TA MuSCs showed that the difference was maximal when considering the top five regulons, including (as expected) Foxc1, Alx4, Dmrta2, Zmiz1 and Fos (Fig. 8B,C). Interestingly, Foxc1, Dmrta2 and Fos were also active regulons in quiescence, with slight disparities between EOM and TA (Table S9).
RNA velocity reveals distinct population kinetics and potential key regulators of progenitor maintenance. (A) Number of regulatory links between regulons and matrisome genes depending on the number of top regulons in EOM and TA in vitro-activated scRNA-seq MuSC datasets. (B) Ratio of the number of regulations of matrisome genes between EOM and TA adult in vitro-activated MuSCs. Maximal difference for top five first regulons. (C) Heatmap of the top regulons in global in vitro-activated scRNA-seq dataset with activity level in each cell. (D,E) Velocity streams overlaid onto a UMAP representation for EOM (D) and TA (E), along with expression patterns of Myog and Pax7 on in vitro-activated adult scRNA-seq dataset. (F,G) Heatmap of driver gene expression from the velocity streams in D,E. (H) Experimental scheme. MuSCs isolated from adult Tg:Pax7-nGFP mice were cultured for 3 days to obtain EOM- and TA-conditioned media (CM). Hindlimb MuSCs were left untreated (media) or treated with CM from D1 to D4 post-plating and pulsed with EdU 2 h before fixation. (I,J) Quantification of the number of nuclei normalised on initial cell number (I) and percentage of EdU+ nuclei (J) of hindlimb MuSCs treated as in H. FC, fold change. N≥2 independent experiments with n≥2 mice each. Different symbol per experiment. (K) Immunostaining of Myf5nlacZ;PdgfraH2BGFP postnatal day (P) 3 mouse head cryosections for β-galactosidase (β-gal) and GFP. Scale bar: 500 µm. Higher magnifications are in K′-K‴. Scale bars: 100 µm. (L) Quantification of the number of β-gal+ (Myf5nlacZ+) GFP+ cells per mm2 in head muscles. n=3 mice. Data are mean±s.e.m. P-values were obtained using one-way ANOVA with Tukey's post hoc test; ns, non-significant (P>0.05), **P<0.01, ***P<0.005, ****P<0.0001.
RNA velocity reveals distinct population kinetics and potential key regulators of progenitor maintenance. (A) Number of regulatory links between regulons and matrisome genes depending on the number of top regulons in EOM and TA in vitro-activated scRNA-seq MuSC datasets. (B) Ratio of the number of regulations of matrisome genes between EOM and TA adult in vitro-activated MuSCs. Maximal difference for top five first regulons. (C) Heatmap of the top regulons in global in vitro-activated scRNA-seq dataset with activity level in each cell. (D,E) Velocity streams overlaid onto a UMAP representation for EOM (D) and TA (E), along with expression patterns of Myog and Pax7 on in vitro-activated adult scRNA-seq dataset. (F,G) Heatmap of driver gene expression from the velocity streams in D,E. (H) Experimental scheme. MuSCs isolated from adult Tg:Pax7-nGFP mice were cultured for 3 days to obtain EOM- and TA-conditioned media (CM). Hindlimb MuSCs were left untreated (media) or treated with CM from D1 to D4 post-plating and pulsed with EdU 2 h before fixation. (I,J) Quantification of the number of nuclei normalised on initial cell number (I) and percentage of EdU+ nuclei (J) of hindlimb MuSCs treated as in H. FC, fold change. N≥2 independent experiments with n≥2 mice each. Different symbol per experiment. (K) Immunostaining of Myf5nlacZ;PdgfraH2BGFP postnatal day (P) 3 mouse head cryosections for β-galactosidase (β-gal) and GFP. Scale bar: 500 µm. Higher magnifications are in K′-K‴. Scale bars: 100 µm. (L) Quantification of the number of β-gal+ (Myf5nlacZ+) GFP+ cells per mm2 in head muscles. n=3 mice. Data are mean±s.e.m. P-values were obtained using one-way ANOVA with Tukey's post hoc test; ns, non-significant (P>0.05), **P<0.01, ***P<0.005, ****P<0.0001.
Next, we set out to determine whether matrix genes underlie the transition towards progenitors and committed cells during activation using scVelo (Bergen et al., 2020). This method computes local changes in the relative amount of unspliced/spliced variants (La Manno et al., 2018) and identifies candidate ‘driver genes’, i.e. most transcriptionally dynamic and responsible for the inferred velocity (Fig. 8D,E). Two distinct velocity streams stood out in both datasets, towards differentiation (MyogHigh) and towards a progenitor-state (Pax7High). Strikingly, a larger fraction of EOM cells appeared to transition towards a progenitor state compared with TA (Fig. 8D,E). These trajectories did not appear to be specifically correlated with cell cycle phases (Fig. S8A,B), which was shown to influence transcriptomic data in some cases (McDavid et al., 2016). Hence, the velocity streams observed most likely reflect transitions between distinct cell states instead of the cell cycle progression of a homogeneous cell state.
Using scVelo built-in functions, we extracted the top driver genes underlying the velocity towards the progenitor state in both datasets (Fig. 8F,G, Tables S10, S11). Out of the top 100 driver genes, 30 were common to both datasets, including Col5a1, which plays a crucial role in maintenance of quiescence (Baghdadi et al., 2018) (Fig. 8F,G). Interestingly, several matrisome components, including Fn1, other collagens and Igfbp7 (insulin growth factor binding protein 7), were unique to EOM. Of note, IGFBPs bind and regulate insulin-like growth factors (IGFs), thereby repressing myogenic differentiation, but they also have IGF-independent activity (Clemmons, 1997; Engert et al., 1996; Jin et al., 2020). GO molecular functions associated with these driver genes during myogenic progenitor maintenance showed that this transition was characterised by active upregulation of ECM components specifically in EOM, in agreement with our previous results (Fig. S8C).
We then performed conditioned media (CM) exchange experiments to gain functional insights into the role of EOM secreted factors by culturing hindlimb cells in presence of EOM-CM, TA-CM or control media for 3 days (Fig. 8H). In agreement with the FOXC1 GOF experiments (Fig. 7R-Y) EOM-CM induced an approximately 5-fold increase in the total number of nuclei/well compared with control media, whereas no significant differences were observed between TA-CM and control media (Fig. 8I). Furthermore, the percentage of EdU+ cells remained unchanged between EOM-CM and control media (Fig. 8J), suggesting that EOM secreted factors may delay differentiation, thus allowing a greater total output.
To assess whether the mesenchymal features of EOM cells were recapitulated in vivo, we asked whether the bipotent EOM progenitors that give rise to myogenic and connective tissue cells during embryogenesis (Grimaldi et al., 2022) persisted postnatally. Examination of Myf5nlacZ;PdgfraH2BGFP mice, where Myf5nlacZ allows tracing myogenic progenitors and progeny, and GFP labels mesenchymal connective tissue cells, confirmed the presence of βgal+ GFP+ cells in the EOM but not in masseter and tongue muscles in the head (Fig. 8K,L). Surprisingly, analysis of EOM and TA MuSCs and fibroadipogenic cells (FAPs) identified a hybrid population positive for both SCA-1 (FAPS marker; Joe et al., 2010) and GFP, only in EOM samples (8% of total GFP+ population; Fig. S8D-F). We then isolated FAPS, MuSCs and the SCA-1+GFP+ cells from the EOM, and performed RT-qPCR, for Pdfgra, Pdgfrb and Col1a2 (Fig. S8G). Notably, the SCA-1+GFP+ progenitor population showed an intermediate pattern of expression of these matrisome genes between that of MuSCs and FAPS, suggesting a mesenchymal signature that exclusively belongs to EOM MuSCs subsets. Finally, differential gene expression analysis of a publicly available scRNA-seq dataset filtering on our EOM-specific regulon TFs, highlighted Egr1 and Foxc1 as upregulated in dystrophic EOM (Fig. S8H), a condition where EOM MuSCs were described to be proliferative (Taglietti et al., 2023). Altogether, this in vivo data corroborated the results obtained by bioinformatic analysis on our in vitro-activated MuSC scRNA-seq dataset.
DISCUSSION
An unusual feature of skeletal muscle is the reliance on distinct gene regulatory networks (GRNs) in different anatomical locations. Most studies on myogenesis have focused on trunk and limb muscles, and only a handful of transcription factors and signalling pathways have been identified as hallmarks of specific muscle groups. Here, we used multiple approaches and have identified unusual features and specific GRNs that functionally distinguish EOM from TA MuSCs.
EOM MuSCs are more refractory to in vitro differentiation
By monitoring Myog and Myod expression, we showed that EOM MuSCs have a lower propensity to differentiate after activation and persist as a proliferative population. Given the pivotal role of Myod in commitment and differentiation (Vicente-García et al., 2022), the lower levels of Myod transcript and protein in activated EOM MuSCs, even upon upregulation of the MyogntdTom reporter, is intriguing and deserves further investigation. As foetal and early postnatal EOM MuSCs are also refractory to differentiation, we propose that this property might be hardwired by unique GRNs that are retained throughout development and adulthood. Notably, trunk foetal MuSCs were shown to be more resistant to myogenic progression upon in vitro expansion and contribute more efficiently upon transplantation than the adult counterparts (Sakai et al., 2013; Tierney et al., 2016). In addition, limb MuSCs cells isolated at birth displayed prolonged expansion and delayed fusion compared with those in the adult (Gattazzo et al., 2020). Similarly, EOM MuSCs showed a reduced fusion index compared with those in adult limb. Altogether, our data suggest that EOM MuSCs retain features of respective fetal and neonatal myogenic cells.
Major obstacles for cell-based therapies are the large cell numbers required for transplantation and the fact that ex vivo amplification of somite-derived MuSCs leads to a drastic decline in regenerative potential due to commitment to differentiation (Briggs and Morgan, 2013; Ikemoto et al., 2007; Montarras et al., 2005). Although significant progress has been made (Charville et al., 2015; Ishii et al., 2018; L'honoré et al., 2018; Xie et al., 2021), the identification of factors that regulate cell fate decisions in distinct MuSC populations is another resource for advancing knowledge in this context.
EOM progenitors exhibit a mesenchymal genetic signature upon activation
scRNA-seq analyses provided some insights into the transcriptional landscape regulating MuSC quiescence, activation and self-renewal in somite-derived muscles (Dell'Orso et al., 2019; Hernando-Herraez et al., 2019; Machado et al., 2021; De Micheli et al., 2020; Yartseva et al., 2020). For cranial MuSCs, our scRNA-seq analysis of activated EOM and TA MuSCs showed distinct transcriptional profiles that divided myoblasts into two subpopulations: those that resembled progenitors and those differentiating. The EOM progenitor population was characterised by higher Pax7, reminiscent of in vitro reserve cells (Laumonier et al., 2017; Yoshida et al., 1998; Zammit et al., 2004), and ontology analysis of its DEGs showed an enrichment in ECM organisation processes and PDGF signalling.
Transcriptomic analysis of trunk MuSCs highlighted a dynamic profile of Pdgf ligands and receptors during myogenesis (Contreras et al., 2021). Treatment with NOTCH and PDGFRβ ligands (DLL4 and PDGF-BB, respectively) enhanced migration, expression of stem cell markers and perivascular-like features in MuSCs (Gerli et al., 2019) and embryonic myoblasts (Cappellari et al., 2013). Notably also, heterogeneity of activated MuSCs was seen during regeneration, including a transitional Notch2-high state and an ECM-high state that regulates self-renewal (Yartseva et al., 2020). As PDGFRβ and several Notch pathway components (e.g. Notch1, Notch3 and Hey1) are co-expressed in activated EOM progenitors, and at significantly higher levels than in TA muscle, we speculate that crosstalk between these pathways could take place in this subpopulation.
Putative role of ECM secretion and remodelling by activated EOM MuSCs
Self-renewal in the muscle lineage is dependent on cell-autonomous expression, deposition and remodelling of ECM components such as Fn1 (Bentzinger et al., 2013; Lukjanenko et al., 2016; Tierney et al., 2016), Col VI (Urciuolo et al., 2013) and laminin α1 (Rayagiri et al., 2018). EOM progenitors relate more to fetal MuSCs with respect to enrichment in matrisome components, such as Fn1, Fbn1, Vcam and collagens (Tierney et al., 2016). Notably, an extensive and more-complex ECM is observed in EOMs in vivo (McLoon et al., 2018). Therefore, we hypothesise that in vitro-activated EOM MuSCs secrete high amounts of ECM and self-autonomously maintain stemness when removed from their niche. For example, FAPs exert a supportive role in myogenesis by secreted matrix and cytokine components (Joe et al., 2010; Kotsaris et al., 2023; Murphy et al., 2011; Uezumi et al., 2014). Myogenic cells expressing mesenchymal markers were also reported in human embryos and in human iPSCs in vitro (Xi et al., 2020). In the EOM, this plasticity was observed in vivo, where EOM myogenic progenitors transition towards non-myogenic cell fates (Grimaldi et al., 2022). Here, we have extended these observations to the early postnatal period, where considerable proliferation and fusion occurs. We did so by using Myf5nlacZ;PdgfraH2BGFP mice and identified a subpopulation of EOM GFP+ cells co-expressing the mesenchymal stem cell marker SCA-1 in vivo, which is normally used as a negative marker of MuSCs (Liu et al., 2015).
A unique network of transcription factors maintains EOM progenitors
Most of the regulon genes that we identified in EOM progenitors are not typical myogenic TFs. One of the most active regulons is Foxc1, a pro-mitogenic factor in cancer (Yang et al., 2017) and driver of endothelial/smooth muscle fates (Lagha et al., 2009; Mayeuf-Louchart et al., 2016; Whitesell et al., 2019). Another top regulon is Egr1, which promotes expression of many ECM-related genes (Gaut et al., 2016; Milet et al., 2017). Interestingly, Foxc1 and Ebf1 are active transcription factors underlying the transition of myogenic towards non-myogenic cell fates in embryonic EOM (Grimaldi et al., 2022). As Foxc1 is expressed in quiescent and activated EOM MuSCs, this gene might reinforce a progenitor/stem-cell identity throughout the lineage. Interestingly, FOXC1 overexpression in TA MuSCs resulted in upregulation of Egr1 and downregulation of Ebf1 in support of our bioinformatic analysis, indicating that these TFs form a co-regulatory module.
As predicted by our analysis in silico, FOXC1 overexpression in TA MuSCs upregulated several EOM matrisome components and downregulated Myod, in accordance with a previous study (Wright et al., 2021), thereby promoting a less committed state. Given that Pitx2 appears to be a target of Foxc1 in our studies (Table S7) and its overexpression promotes MuSC proliferation and enhances the regenerative potential of dystrophic MuSCs (Vallejo et al., 2018), we expected an increase in nuclear output when it was overexpressed in TA MuSCs. Yet, we observed the opposite phenotype in pure TA GOF cultures, thereby raising several possibilities. First, it has been reported that sustained expression of FOXC1 can induce stem cell quiescence in the skin (Lay et al., 2016; Wang et al., 2016). Alternatively, a strong upregulation of FOXC1 might have induced an irreversible drift towards a mesenchymal phenotype in hindlimb MuSCs, as indicated by the reduced expression of Pax7. Yet, our mixed GOF cultures also suggest that the FOXC1-dependent secretome can promote proliferation of adjacent non-transduced cells. Finally, as Foxc1 is part of a co-regulatory module in the EOM, its singular overexpression might not be sufficient to confer the full EOM phenotype to limb cells. Future studies should assess the role of this TF in EOM progenitor emergence and disease.
Conclusions
Using in-depth bioinformatic analysis, in vitro approaches and analysis of expression patterns in vivo, we propose a model where the outperformance of EOM MuSCs depends on the expression of a tightly associated module of transcription factors regulating a distinct pattern of ECM-remodelling factors, cell receptors and growth factor-binding proteins. These components define the pace at which EOM MuSCs progress through the myogenic lineage and maintain a stem-like population. As such, our study lays the groundwork for elucidating the mechanisms of selective sparing of muscle groups in dystrophic disease by providing information on a unique core GRN within this muscle group.
MATERIALS AND METHODS
Animal care
Animals were handled according to national and European Community guidelines and Institut Pasteur ethics committee approved the protocols (APAFIS#41051-202302204207082). Except when indicated otherwise, males and females of 2-4 months were used. Tg:Pax7-nGFP (Sambasivan et al., 2009), Myf5Cre (Haldar et al., 2007), Myf5nlacZ (Tajbakhsh et al., 1996), MyogntdTom (Benavente-Diaz et al., 2021), PdgfraH2BGFP (Hamilton et al., 2003), Pax7CreERT2 (Mathew et al., 2011) and R26tdTom (Madisen et al., 2009) mouse lines were maintained in a C57Bl/6JRj background. To induce recombination of Pax7CreERT2;R26tdTom mice, a 20 mg/ml stock solution of tamoxifen was prepared in 5% ethanol and 95% sunflower seed oil by thorough resuspension with rocking at 4°C. For adult mice, 2 mg of tamoxifen (Sigma #T5648) were administered by gavage for 5 consecutive days and animals sacrificed 5 days later. To induce recombination of pups, the tamoxifen stock solution was diluted to 15 mg/ml with sunflower seed oil and 0.15 mg were administered daily by subcutaneous injection between P4 and P6 daily. Pups were euthanised at P10 by decapitation and adult mice were euthanised by cervical dislocation.
Muscle stem cell culture, treatment and transfection
MuSC isolation was performed as described previously (Gayraud-Morel et al., 2017) with some modifications indicated in the supplementary Materials and Methods. Cells isolated by FACS were cytospun (3 min at 50 g at room temperature), plated onto Matrigel (Corning, 354248)-coated dishes and cultured in growth media [1:1 mix of Dulbecco's modified Eagle's medium (DMEM) (Gibco, 10569010) and F12 (Gibco, 31765027), 20% FBS (ThermoFisher, 10270), 2% Ultroser (Sartorius, 15950-017), 1% Penicillin/Streptomycin (Gibco, 15140-122)] at 3% O2, 5% CO2, 37°C for the indicated time. Half the volume of medium was changed every second or third day. To assess proliferation, cells were pulsed with 10−6 M of EdU (ThermoFisher, C10640) in cell culture media for 2 h before fixation. To induce myogenic differentiation and fusion, myoblasts were plated at high density (33,000 cells/cm2 in Fig. 7L,M, Fig. S7R,Z or 75,000 cells/cm2 in Fig. S1H) onto Matrigel-coated plates in the growth medium. Once adherent, cells were changed to differentiation medium (DMEM with 5% serum and 1% penicillin and streptomycin).
Conditioned media (CM) experiments were performed by isolating EOM and TA MuSCs by FACS based on GFP fluorescence from Tg:Pax7-nGFP mice, and collecting the media 3 days post-plating. Hindlimb MuSCs were cultured in growth media (control) or treated with CM from day 1 to day 4 post-plating before analysis.
For loss-of-function experiments, freshly isolated satellite cells from Tg:Pax7-nGFP mice were transfected in suspension immediately after FACS with the ON-TARGET plus SMARTpool against FOXC1 (Dharmacon, L-047399-01-0005) or Scramble (Dharmacon, ON-TARGETplus Non-targeting Control siRNA, D0018100205) at 200 nM final concentration using Lipofectamine 3000 (ThermoFisher, L3000001) in Opti-MEM (Fisher Scientific, 11564506) as described by the manufacturer. Briefly, a pre-mix of siRNA/Optimem (1.5 μl/20 μl) and Lipofectamine3000/Optimem (0.3 μl/20 μl) were incubated separately for 5 min at room temperature, mixed at a 1:1 ratio and incubated for a further 15 min at room temperature. Then 2×104 cells in 40 μl of Optimem were incubated with an equal volume of the transfection mix for 2 h in 3% O2 and 5% CO2 at 37°C in Eppendorf tubes in which the caps had been punctured with a needle to allow gas exchange. Two hours after transfection, three volumes of fresh growth medium were added and cells were plated at 10 k/cm2 in Matrigel-coated wells containing growth media. Two days after transfection, wells were processed for immunostaining or RNA collected with TRIzol. Details on immunostaining, RT-qPCR and western blotting are described in the supplementary Materials and Methods. RT-qPCR primers are listed in Table S13.
Lentivirus transduction for gain- and loss-of-function experiments
For gain-of-function (GOF) assays, lentiviral viruses made in our lab [control (pLVX-CAG-P2A-mCherry) and Foxc1-mCherry (pLVX-CAG-Foxc1-P2A-mCherry), CAG promoter] and viruses produced by VectorBuilder (https://en.vectorbuilder.com/) [control (pLVX-EF1A>mCherry), Foxc1-mCherry (pLVX-EF1A>mFoxc1-T2A-mCherry), Creb3l1-mCherry (pLVX-EF1A>mCreb3l1-T2A-mCherry), EF1a promoter] were used. T2A or P2A causes co-translational cleavage of the encoded polypeptide (Table S12). GOF assays of FOXC1 with our own or VectorBuilder-made viruses were used interchangeably as they gave comparable results. For LOF assays, three different Foxc1-shRNA vectors were tested with identical results (U6 promoter for shRNA; hPGK for mCherry reporter).
Freshly isolated satellite cells from Tg:Pax7-nGFP mice were plated in Matrigel-coated wells, cultured overnight and transduced at an MOI (multiplicity of infection) of 100 in 45 µl (for 96-well plates) or 125 µl (for 48-well plates) of MuSC media containing 5 µg/ml of polybrene. After 4 h incubation at 37°C, cells were washed three times with media and cultured for four more days in GM before fixation for immunostaining, protein or RNA collection, or for re-isolation of mCherry-positive cells when needed.
Image analysis
For cells grown on 96 wells, images from the Opera Phenix high-content microscope were quantified using Harmony software and a semi-automated pipeline. Nuclei were detected based on Hoechst signal and mean intensity of fluorescence were automatically quantified on the nuclear region for nuclear markers. The numbers of nuclei positive for a specific nuclear marker were scored on a manually defined threshold. For cells grown on IBIDI eight-well plates, images were quantified with the point counter tool in Fiji (Schindelin et al., 2012) or using TrackMate-StarDist (Ershov et al., 2022).
Immunostaining on ECM markers on cells were scanned using either the Opera Phoenix high-content microscope (Perkin Elmer) or LSM800 microscope (Zeiss) with ZEN software. Imaging of tissue sections was carried out on an LSM800 microscope or a Ti2E Spinning disk (Nikon).
The fusion index was calculated as the fraction of nuclei contained within MF20+ myotubes, which had two or more nuclei, when compared with the number of total nuclei within each image. Figures were assembled in Adobe Photoshop, Illustrator and InDesign (Adobe Systems).
Time lapse microscopy
MuSCs were plated on a microscopy culture chamber (IBIDI, 80826) or 96 wells of TC-Treated Black µCLEAR plates (Greiner Bio-One 655090) and cultured in growth media. The plate was incubated at 37°C in 5% CO2 and 3% O2 (Zeiss, Pecon). A Zeiss Observer.Z1 connected to a Plan-Apochromat 20×/0.8 M27 objective and Hamamatsu Orca Flash 4 camera piloted with ZEN (Zeiss) was used. After activation for 3 days in vitro, EOM and TA cells were imaged every 12 min for ∼48 h. Individual cells were tracked and TOM intensity was scored at selected time points. Tracking of MyogntdTom cells was performed using the Manual Tracking feature of the TrackMate plug-in (Tinevez et al., 2017) in Fiji.
For lentivirus live imaging experiments, detection of total cell numbers and mCherry-positive cells was carried out with a custom-made pipeline. To detect cell contours from the bright field images, we retrained the «livecell» model of CellPose with CellPose 2.0 (Pachitariu and Stringer, 2022). The retraining dataset was created from small crops of the full data with low-density and high-density regions and the different cell phenotypes. To homogenise the data, preprocessing of the bright-field channel was applied to all the data: we enhanced the contrast locally with CLAHE (Contrast Limited Adaptive Histogram Equalization) in Fiji (with a blocksize of 127, default parameters) and normalised the intensity (for each frame). To assess which cells were infected by the virus (mCherry+), the mean mCherry fluorescent intensity was measured in each cell and normalised to the local mean fluorescent intensity in the well (by averaging the fluorescent intensity with a radius of 300 µm). This allowed compensation for differences in illumination, even within the same acquisition, but could slightly underestimate the relative intensity of cells in very dense regions. Finally, the mCherry intensity channel was thresholded after a background removal with the Fiji Triangle threshold method calculated with the histogram of the four time frames together.
The retrained models were then run on the full data (cropped regions or full wells) to detect the cell contours through the Fiji plug-in TrackMate-CellPose (Ershov et al., 2022), which allowed us to directly access measurements of the cells («spots» table in TrackMate).
Data analysis and statistics
Data analysis, statistics and visualisations were performed using Prism (Graphpad Software), or using R (R Core Team, 2014) and the package ggplot2 (Wickham, 2009). Tests are described in the figure legends. Information related to microarray data analysis and scRNA-seq of DMD rat EOM are described in the supplementary Materials and Methods.
scRNA-seq data generation
MuSCs were isolated on BD FACSAria III based on GFP fluorescence and cell viability from Tg:Pax7-nGFP mice. Quiescent MuSCs were manually counted using a hemocytometer and immediately processed for scRNA-seq. For activated samples, MuSCs were cultured in vitro as described above for 4 days. Activated MuSCs were subsequently trypsinised and washed in DMEM/F12 2% FBS. Live cells were re-isolated, manually counted using a hemocytometer and processed for scRNA-seq.
Before scRNA-seq, RNA integrity was assessed using Agilent Bioanalyzer 2100 to validate the isolation protocol (RIN>8 was considered acceptable). 10X Genomics Chromium microfluidic chips were loaded with around 9000 cells and cDNA libraries were generated following the manufacturer's protocol. Concentrations and fragment sizes were determined using Agilent Bioanalyzer and Invitrogen Qubit. cDNA libraries were sequenced using NextSeq 500 and High Output v2.5 (75 cycles) kits. Count matrices were subsequently generated following 10X Genomics Cell Ranger pipeline. After normalisation and quality control, we obtained an average of 5792±1415 cells/condition.
Seurat preprocessing
scRNA-seq datasets were processed using Seurat (https://satijalab.org/seurat/) (Butler et al., 2018). Cells with more than 10% of mitochondrial gene fraction were discarded and 4000-5000 genes were detected on average across all four datasets. Dimensionality reduction and UMAPs were generated by following Seurat workflow. The top 100 DEGs were determined using Seurat ‘FindAllMarkers’ function with default parameters. When processed independently (scvelo), the datasets were first regressed on cell cycle genes, mitochondrial fraction, number of genes and number of UMI after a Seurat dedicated vignette, and doublets were removed using DoubletFinder v3 (McGinnis et al., 2019). A ‘StressIndex’ score was generated for each cell based on the list of stress genes previously reported (Machado et al., 2021) using the ‘AddModule’ Seurat function. Ninety-four out of 98 genes were detected in the combined datasets. UMAPs were generated (1) after StressIndex regression; and (2) after complete removal of the detected stress genes from the gene expression matrix before normalisation. In both cases, the overall aspect of the UMAP did not change significantly (Fig. S5G). Although immeasurable confounding effects of cell stress after isolation cannot be ruled out, we reasoned that our datasets did not show a significant effect of stress with respect to the conclusions of our study.
Matrisome analysis
After subsetting for the features of the Matrisome database (Naba et al., 2015) present in our single-cell dataset, the matrisome score was calculated by assessing the overall expression of its constituents using the ‘AddModuleScore’ function from Seurat (Butler et al., 2018).
RNA velocity and driver genes
Scvelo was used to calculate RNA velocities (Bergen et al., 2020). Unspliced and spliced transcript matrices were generated using velocyto (La Manno et al., 2018) command line function. Seurat-generated filtering, annotations and cell-embeddings (UMAP, tSNE, PCA) were then added to the outputted objects. These datasets were then processed following scvelo online guide and documentation. Velocity was calculated based on the dynamical model [using scv.tl.recover_dynamics(adata) and scv.tl.velocity(adata, mode=’dynamical’)] and differential kinetics calculations were added to the model [using scv.tl.velocity(adata, diff_kinetics=True)]. Specific driver genes were identified by determining the top likelihood genes in the selected cluster. The lists of top 100 drivers for EOM and TA progenitors are given in Tables S10 and S11.
Gene regulatory network inference and transcription factor modules
Gene regulatory networks were inferred using pySCENIC (Aibar et al., 2017; Van de Sande et al., 2020). This algorithm regroups sets of correlated genes into regulons (i.e. a transcription factor and its targets) based on binding motifs and co-expression patterns. The top 35 regulons for each cluster was determined using scanpy ‘scanpy.tl.rank_genes_groups’ function (method=t-test). This function can yield fewer than 35 results depending on the cluster. UMAP and heatmap were generated using regulon AUC matrix (Area Under Curve), which refers to the activity level of each regulon in a given cell. Visualisations were performed using scanpy (Wolf et al., 2018). The outputted list of each regulon and their targets was subsequently used to create a transcription factor network. To do so, only genes that are regulons themselves were kept. This results in a visual representation where each node is an active transcription factor and each edge is an inferred regulation between two transcription factors. When placed in a force-directed environment, these nodes aggregate based on the number of shared edges. This operation greatly reduced the number of genes involved, while highlighting co-regulating transcriptional modules. Visualisation of this network was performed in a force-directed graph using Gephi ‘Force-Atlas2’ algorithm (https://gephi.org/). Of note, a force-directed graph is a type of visualisation technique where nodes are positioned based on the principles of physics that assign forces among the set of edges and the set of nodes. Spring-like attractive forces are used to attract pairs of edges towards each other (connected nodes), while repulsive forces, like those of electrically charged particles, are used to separate all pairs of nodes. In the equilibrium state for this system, the edges tend to have uniform length (because of the spring forces) and nodes that are not connected by an edge tend to be drawn further apart (because of the electrical repulsion).
Gene ontology analysis
Gene ontology analyses were performed on top 100 markers (obtained from Seurat function FindAllMarkers) or on top 100 driver genes (obtained from scvelo), using Cluego (Bindea et al., 2009). ‘GO Molecular Pathway’, and ‘REACTOME pathways’ were used independently to identify common and unique pathways involved in each dataset. In all analyses, an enrichment/depletion two-sided hypergeometric test was performed and P-values were corrected using the Bonferroni step-down method. Only the pathways with P<0.05 were displayed.
Literature scores
The scores for SMMCs (Giordani et al., 2019), FAPs (Oprescu et al., 2020), myotendinous junction B myonuclei (Kim et al., 2020), Twist2+ population (Liu et al., 2017), fetal MuSCs (Tierney et al., 2016), developing limb connective tissues (Esteves de Lima et al., 2021) and the human skeletal muscle mesenchyme (Xi et al., 2020) were calculated by assessing the overall expression of the markers of each population (Table S5) using the ‘AddModuleScore’ function from Seurat (Butler et al., 2018).
Acknowledgements
We gratefully acknowledge the Genomic Platform (Inserm U1016-CNRS UMRS 8104 - Université de Paris Cité) at Institut Cochin, the UtechS Photonic BioImaging (Imagopole), C2RT, Institut Pasteur [supported by the French National Research Agency (France BioImaging; ANR-10–INSB–04; Investments for the Future) and the Center for Translational Science (CRT)-Cytometry] and the Biomarkers Unit of Technology and Service (CB UTechS) at Institut Pasteur for support in conducting this study. We warmly thank Pierre Henri Comere for assistance with cell sorting and Marion Louveaux for establishing the initial image analysis pipeline.
Footnotes
Author contributions
Conceptualization: S.T., G.C.; Methodology: D.D.G., M.B.-D., A.G., G.C.; Software: M.B.-D., A.G., V.L., S.M.; Validation: D.D.G., M.B.-D., M.M., P.T.L., B.E., S.G., G.C.; Formal analysis: D.D.G., M.B.-D., M.M., A.G., P.T.L., B.E., M.K., S.G., V.L., J.-Y.T., G.L., S.M., G.C.; Investigation: D.D.G., M.B.-D., M.M., A.G., M.K., G.C.; Data curation: G.C.; Writing - original draft: D.D.G., M.B.-D., A.G., G.C.; Writing - review & editing: D.D.G., M.M., A.G., P.T.L., B.E., M.K., S.G., V.L., S.M., S.T., G.C.; Visualization: M.B.-D., M.M., A.G., P.T.L., G.C.; Supervision: S.T., G.C.; Project administration: G.C.; Funding acquisition: S.T., G.C.
Funding
We acknowledge funding support from the Institut Pasteur, Agence Nationale de la Recherche (Laboratoire d'Exellence Revive, Investissement d'Avenir, ANR-10-LABX-73 to S.T. and ANR-21-CE13-0005 MUSE to G.C.), Association Française contre les Myopathies (20510 to S.T. and 23201 to G.C.), the European Research Council (Advanced Research Grant 101055234), and theCentre National de la Recherche Scientifique (CNRS). M.B.-D. was supported by a grant from Laboratoire d'Excellence Revive and La Ligue Contre le Cancer. M.K. was supported by a Post-Doctoral Fellowship from Uehara Memorial Foundation.
Data availability
The code that was used to generate the TF network is available at https://github.com/TajbakhshLab/TFnetwork. scRNA-seq datasets are available from the Dryad Digital Repository (Di Girolamo et al., 2024): dryad.b8gtht7k0. Raw scRNA-seq sequencing data have been deposited in the GEO under accession number GSE244964.
References
Competing interests
The authors declare no competing or financial interests.