Extraocular muscle stem cells exhibit distinct cellular properties associated with non-muscle molecular signatures

: The muscle stem cell (MuSC) population is recognized as functionally heterogeneous. Cranial muscle stem cells, which originate from head mesoderm, can have greater proliferative capacity in culture and higher regenerative potential in transplantation assays when compared to those in the limb. The existence of such functional differences in phenotypic outputs remain unresolved as a comprehensive understanding of the underlying mechanisms is lacking. We addressed this issue using a combination of clonal analysis, live imaging, and scRNA-seq, identifying critical biological features that distinguish extraocular (EOM) and limb (Tibialis anterior, TA) MuSC populations. Time-lapse studies using a Myogenin tdTomato reporter showed that the increased proliferation capacity of EOM MuSCs is accompanied by a differentiation delay in vitro . Unexpectedly, in vitro activated EOM MuSCs expressed a large array of distinct extracellular matrix (ECM) components, growth factors, and signaling molecules that are typically associated with mesenchymal non-muscle cells. These unique features are regulated by a specific set of transcription factors that constitute a coregulating module. This transcription factor network, which includes Foxc1 as one of the major players, appears to be hardwired to EOM identity as it is present in quiescent adult MuSCs, in the activated counterparts during growth and retained upon passages in vitro. These findings provide insights into how high-performing MuSCs regulate myogenic commitment by active remodeling of their local environment.


INTRODUCTION
Skeletal muscle stem cells are arguably the most diverse in terms of molecular properties compared to other adult stem cells that are distributed throughout the body.Genetic and transcriptomic studies have shown that the MuSC population in any particular anatomical location is heterogeneous.Certain subsets are more prone to self-renewal or differentiation, differ in transplantation efficiency, as well as their stem cell-niche interactions, metabolism, and resistance to stress upon activation (Barruet et al., 2020;Chakkalakal et al., 2014;Cho et al., 2020;Dell'Orso et al., 2019;Dumont et al., 2015;Gayraud-Morel et al., 2012;Hernando-Herraez et al., 2019a;Micheli et al., 2020;Morree et al., 2019;Ono et al., 2012;Rocheteau et al., 2012;Scaramozza et al., 2019;Tierney et al., 2018;Vartanian et al., 2019;Yartseva et al., 2020;Yennek et al., 2014).Despite this diversity, MuSCs share common functions as they are essential for postnatal growth and repair of the skeletal muscle system (Lepper, 2011;Murphy, 2011;Sambasivan, 2011).In response to muscle tissue damage or growth factors in culture media, MuSCs activate, proliferate and undergo cell fate decisions to either self-renew to replenish the MuSC pool (expressing Pax7), or commit to the myogenic program and differentiate into myoblasts and fusioncompetent myocytes that express the commitment and differentiation factors Myod and Myog, respectively (Evano and Tajbakhsh, 2018;Zammit et al., 2006).
An unexpected finding was the discovery that MuSCs in different anatomical locations are programmed with distinct upstream transcription factors (TF) prior to acquiring myogenic identity (Gopalakrishnan et al., 2015;Harel et al., 2009;Kelly et al., 2004;Sambasivan et al., 2009a;Tajbakhsh, 1996).This intriguing observation correlates with reports where subsets of cranial MuSCs are functionally more robust in terms of proliferation, engraftment efficiency, and susceptibility to disease, when compared to those in the limb (Randolph and Pavlath, 2015a).However, the intrinsic and extrinsic regulators conferring robust features to cranial muscles have not been characterized.
EOMs are a subset of head muscles that govern eye movements.EOMs are derived from unsegmented cranial mesoderm and are regulated by distinct transcriptional factors and signalling molecules compared to the somite-derived limb and trunk muscle groups (Grimaldi and Tajbakhsh, 2021;Michailovici et al., 2015;Sambasivan et al., 2011).Notably, when the transcription factor Pax3 is inactivated, no limb muscles develop, yet cranial muscles including the EOMs are present (Sambasivan et al., 2009b;Tajbakhsh et al., 1997).Conversely, 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., 1999a;Zacharias et al., 2011).Additionally, whole muscle transcriptional profiling revealed an extraordinary diversity amongst adult skeletal muscle groups (Kui et al., 2022;Terry et al., 2018).As some muscle subsets, like the EOMs, are preferentially spared in muscular dystrophies and during ageing (Emery, 2002;Formicola et al., 2014;Verma et al., 2017a), it has been proposed that intrinsic properties of MuSCs or the respective myofibers could determine their differentially sensitivity to disease (Randolph and Pavlath, 2015a;Terry et al., 2018).
Similar to MuSCs derived from other craniofacial muscles (Ono et al., 2010a;Terry et al., 2018), EOM MuSCs have a greater proliferative capacity in culture when compared to those in the limb, with this property being conserved in Mdx mice, a model of Duchenne muscular dystrophy (Dmd), and in a preclinical rat model of the disease (Stuelsatz et al., 2015a;Taglietti et al., 2023).Additionally, EOM MuSCs have a higher engraftment potential compared to those in the limb (Stuelsatz et al., 2015a).
Although debated (McLoon and Wirtschafter, 2002;McLoon and Wirtschafter, 2003;Stuelsatz et al., 2015b;Verma et al., 2017b), EOM MuSCs were reported to continuously contribute to myofibres in homeostatic conditions, suggesting a higher level of basal activation of the MuSC population (Keefe et al., 2015;Pawlikowski et al., 2015;Taglietti et al., 2023).In addition, changes in the transcriptomic profile of EOM MuSCs after transplantation into limb muscle (Evano et al., 2020) showed that despite a large transcriptional reprogramming after transplantation, about 10% of EOM specific genes persisted in the heterotypic niche environment.This finding suggests that cell autonomous regulation of MuSC properties predominates to a certain extent in EOM MuSCs (Evano et al., 2020a).Moreover, recent single-cell transcriptomic analysis identified the thyroid hormone signalling pathway as a key factor preventing entry into senescence of EOM MuSCs in Dmd rats (Taglietti et al., 2023).Yet, whether deeply-rooted TF gene regulatory networks exist within MuSCs subsets, and contribute to the maintenance of anatomically-distinct MuSC phenotypes remains largely unexplored.
Here, we used droplet-based scRNA-seq to investigate the transcriptional states that govern the outperformance of MuSC subsets during activation.By integrating genome-wide analysis with cell biology approaches using mouse reporter lines, we identified key regulators that confer distinct mesenchymal-like features to activated EOM MuSCs and expression of extracellular matrix (ECM) components.Computational analyses showed that this pool of EOM cells displays features of a more stem-like state that is actively maintained during in vitro proliferation through a specific set of transcription factors forming a co-regulatory module.

Functional differences among MuSCs are accompanied by heterogeneity in myogenic gene expression dynamics
In agreement with previous reports (Stuelsatz et al., 2015b), we show that culture of EOM and Tibialis anterior (TA) MuSCs isolated by FACS using Tg:Pax7-nGFP mice (Sambasivan et al., 2009b), revealed a larger fraction of proliferative cells in the EOM from day 3 (D3) onwards, peaking at D4 and remaining higher even when differentiation and fusion starts to take place at D5 (Figure 1A).To obtain insights into the cellular basis for the outperformance of MuSC subsets upon exit from quiescence, we performed clonal and live imaging analysis in vitro.First, MuSCs from EOM and TA were plated at clonal density in 96-well plates.In parallel, a fraction of the isolated cells was plated as a bulk culture, allowed to activate for 2 days in vitro, re-isolated and then plated at clonal density (Suppl Figure 1A).The total cellular output in both cases was quantified two weeks after seeding.In accordance with a previous study (Stuelsatz et al., 2015a), EOM MuSCs displayed a markedly higher clonal capacity than those from the TA muscle, with a mean of 1315 cells/clone and 186 cells/clone respectively (7-fold difference; Suppl Figure 1B, T1).Surprisingly, the clonogenic properties of EOM MuSCs continued to be elevated after 2 days of activation in vitro where they yielded 3656 cells/clone compared to 297 cells/clone for activated TA MuSCs (12-fold difference; Suppl Figure 1B, T2).
We then monitored their differentiation dynamics in vitro by live imaging using the Tg:Pax7-nGFP;Myog ntdTom mouse lines (Benavente-Diaz et al., 2021;Sambasivan et al., 2009a) to isolate MuSCs by FACS and then monitor the onset of Myog expression using a nuclear localised tdTOMATO fluorescent reporter.After activation for 3 days in vitro, EOM and TA cells were imaged every 12 min for ~48h (Figure 1B).Individual cells were tracked and tdTOMATO intensity was scored at selected time points.EOM and TA MuSCs displayed distinct profiles of lineage progression, where TA myogenic cells initiated reporter gene expression earlier than EOM cells (Figure 1C).The percentage of Myogexpressing cells sharply increased in the TA from 5% at 80h post-plating to 40% at 96h.Significantly, only ~7% of EOM cells were tdTOMATO + by 96h in vitro (Figure 1C).Since Myog expression is followed by cell-cycle withdrawal and terminal differentiation (Andrés and Walsh, 1996;Benavente-Diaz et al., 2021;Guo et al., 1995), a delay in its expression could allow for sustained proliferation of EOM MuSCs.We performed RT-qPCR at 48h post-plating and assessed the expression of key myogenic transcription factors and the cell cycle inhibitor p21, which regulates cell cycle exit during myoblast differentiation (Guo et al., 1995;Zhang et al., 1999).In agreement with the live imaging data, Myog mRNA was significantly upregulated in TA activated MuSCs compared to EOM, concomitant with an increase in p21 mRNA (Figure 1D).No significant differences were found in Pax7 nor Myod expression (Figure 1D).Moreover, EOM progenitors displayed reduced levels of phosphorylated P38 (Suppl Fig1C,D).As p38α/β MAPKs play a central role in promoting MEF2 transcriptional activity and initiating the differentiation program (Rugowska et al., 2021), this data suggests that EOM progenitors are less prone to myogenic commitment.Therefore, our results suggest that a delay in the differentiation of EOM MuSCs might contribute to the higher number of myogenic cells produced compared to TA cells.

Activation of EOM MuSCs include a larger progenitor pool and distinct activation signatures
To gain insights into the potential mechanisms that regulate the phenotypic differences and cell fate decisions between EOM and TA MuSCs, we performed scRNA-seq of in vitro activated MuSCs using the 10x Chromium platform.We first isolated EOM and TA MuSCs by FACS from Tg:Pax7-nGFP mice and cultured them in vitro for 4 days.Cell cultures were subsequently trypsinized, cell-sorted and processed for scRNA-seq (Figure 2A).Unsupervised clustering divided cells into 2 clusters per sample (Figure 2B) that were annotated as progenitors (Prog) or differentiating (Diff) based on the expression of known myogenic makers, such as Pax7/Myf5 and Myod/Myog (Figure 2B, C).
Analysis of differentially expressed genes (DEGs) across these clusters showed that each cell state had a distinct transcriptional pattern (Figure 2D, Suppl Table 1).As expected, TA progenitors expressed a defined Hox-signature (Evano et al., 2020b) and inhibitors of differentiation (Id1, Id2, (Jen et al., 1992;Kumar et al., 2009)).Instead, EOM progenitors displayed markers that had not been previously described in activated MuSC datasets such as Mgp, Bgn, Col1A2 and Acta2 (smooth muscle actin).
EOM and TA differentiating cells shared part of their signature (Myog, Mymx, Myl4).However, TA progenitors shared some markers with the TA differentiating cells (Acta1, Mylpf, Tnnt2), suggesting that Prog and Diff TA clusters are closer at a transcriptional level than their EOM counterparts (Figure 2D).
To get insights into these properties and the prospective functional heterogeneity associated with each transcriptomic signature we performed reactome pathway analysis of the DEGs upon in vitro activation (Figure 2E, Suppl Figure 2A, Suppl Table 2).Multiple pathways involved in ECM organisation were characteristic of EOM progenitors including "ECM proteoglycans", "crosslinking of collagen fibrils" and "integrin cell surface interactions".On the other hand, EOM and TA differentiating cells were both enriched in pathways involved in striated muscle contraction and creatinine metabolism, which are associated with maturation of skeletal muscle.Also, as suggested by the DEG heatmap, TA progenitors shared common pathways with the differentiating clusters.Analysing the molecular function of the DEGs revealed that EOM progenitors were also specifically enriched in Pdgfrb, a receptor for platelet-derived growth factor, and proteins related to insulin like growth factor binding and integrin signaling (Figure 2F, Suppl Figure 2B, Suppl Table 3).Thus, these analyses uncovered a non-canonical signature in the EOM progenitor fraction upon activation.

EOM and TA MuSCs retain distinct molecular signatures after several days in culture
We then set out to determine to what extent the EOM progenitor fraction that was capable of self-renewal in vitro, relates to the quiescent MuSCs present in vivo.Thus, we assessed the profiles of quiescent MuSCs by scRNA-seq of from EOM and TA immediately after isolation by FACS (Figure 3A-B).Despite the homogeneous expression of MuSC markers Myf5 and Pax7 in both populations, Myod expression was restricted to the TA MuSCs (Figure 3C-D), suggesting an earlier MuSC activation in response to the isolation procedure (Brink et al., 2017;Machado et al., 2017;Velthoven et al., 2017).Differentiation markers like Myog were absent in both populations (Figure 3C-D).Visualization of the top 25 most variably expressed genes documented distinct transcriptional programs of EOM and TA quiescent clusters (Figure 3E, Suppl Table 4).Moreover, microarray analysis of adult quiescent Tg:Pax7-nGFP MuSCs revealed that the EOM quiescent signature was also different from that of MuSCs from other cranial-mesoderm derived muscles such as the esophagus and masseter (Suppl Figure 3A).Interestingly, we found several conserved genes between the quiescent and activated sc-RNAseq states at each anatomical location, while others were unique to the quiescent or activated cell states (Figure 3F-I).In the TA, Lbx1, Vgll2 and Hox genes were found as common signature throughout cell states (Figure 3G,I, Suppl Figure 3B, 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., 2000a) and Vgll2 (previously named Vito-1) is a key cofactor of the myogenic differentiation program (Günther et al., 2004;Maeda et al., 2002).All three genes had been already identified in the limb both at the MuSC levels and in entire muscles in the adult (Evano et al., 2020b;Honda et al., 2017;Terry et al., 2018).
Pcdh7, an integral membrane protein belonging to cadherin superfamily that regulates intercellular adhesion (Wang et al., 2020;Yoshida, 2003) is unique to TA quiescent state and the microtubule protein Tubb3 (Duly et al., 2022) is a marker marker of the TA activated state (Figure 3I).EOM genes unique to the quiescent state included Tshr, encoding the thyroid-stimulating hormone receptor, which has been previously identified in bulk RNAseq as "EOM-resistant" kept upon engraftment into the limb (Evano et al., 2020b), and shown to control senescence in MuSCs of DMD rats (Taglietti et al., 2023).Matrix Gla protein (Mgp), which is a critical regulator of angiogenesis in multiple organs (Kida and Yamaguchi, 2022), was identified as exclusively upregulated in EOM cells upon activation (Figure 3H).Amongst the conserved EOM genes across cell states we identified Pitx2, a well-known major upstream regulator of EOM development (Gage et al., 1999b), Fos gene family members (Fos and Fosb), with Fos being recently identified in a subset of limb MuSCs with enhanced regenerative capacity (Almada et al., 2021), and Igfbp7, a specific marker of quiescent MuSCs that is downregulated upon activation (Fukada et al., 2007) but also reported to be upregulated in MuSCs upon exercise (Chen et al., 2020) (Figure 3F, Suppl Figure 3B, C).Additionally, the EOM common signature included several ECM components and regulators (e.g.Bgn, Loxl1, Col1a2, 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)) (Figure 3F,H).Finally, other top markers for EOM MuSCs identified in the quiescent microarray including Foxc1, which is involved in ocular development (Smith et al., 2000) and reported to regulate the balance between myogenic and vascular lineages within somites (Lagha 2009,Mayeuf-Louchart 2016), Eya2, Islr and Igfbp4, were also enriched in the activated EOM sc-RNAseq dataset (Suppl Figure 3B).Altogether, our analysis revealed a closer molecular overlap for EOM MuSCs across cell states, including several TFs and ECM markers that continued to be distinct from TA MuSCs.

EOM MuSC activation involves extensive ECM remodelling and a mesenchymal-like signature
Given the well-established role of ECM synthesis and remodelling on MuSC proliferation and selfrenewal (Baghdadi et al., 2018a;Rayagiri et al., 2018;Tierney et al., 2016a;Yin et al., 2013), we sought to further characterise the EOM and TA signatures using the Matrisome database (Matrisome DB) (Naba et al., 2015), which compiles in silico and experimental data on ECM constituents.We first identified the components of the Matrisome DB present in our single-cell dataset and created a global score for ECM component expression (Figure 4A) (see Material and Methods).Significantly, the EOM progenitor cluster had the highest matrisome score, while no major differences where observed within the other clusters (Figure 4A).Next, we identified the matrisome components that were differentially expressed in each cluster and found the EOM progenitor cluster also expressed the highest number of genes for each matrisome category when compared to the other clusters (Figure 4B).
We then validated some of the matrisome candidate genes at the protein level in activated MuSCs.
EOM and TA MuSCs were isolated by FACS, cultured in vitro for 4 days and protein expression was assessed by immunofluorescence or Western blot (Figure 4C-E).Immunostaining for Fibronectin (FN1) showed a marked difference in fluorescence intensity between EOM and TA cells, confirming that EOM cells express higher levels of this protein (Figure 4D).Given that FN1 secretion by MuSCs promotes MuSC expansion in a cell autonomous manner (Bentzinger et al., 2013), the distinct expression of Fn1 might contribute to maintenance of the progenitor population during EOM MuSC activation.EOM progenitors also had higher levels of Collagen I (COLI) (Figure 4D), a major component of the fibrotic ECM (Dulauroy 2012) shown to suppress differentiation of C2C12 cells and display mutually exclusive expression with Myogenin (Alexakis et al., 2007).Similarly, EOM progenitors expressed higher levels of Collagen IV (COLIV) (Figure 4D), which was reported to be secreted by MuSCs as well as myoblasts and fibroblasts in culture (Baghdadi et al., 2018a;Kühl et al., 1984).In addition, PDGFRβ, a tyrosinekinase receptor commonly expressed by mesenchymal cells and pericytes (Hellström et al., 1999;Levéen et al., 1994;Soriano, 1994), was also differentially expressed in EOM progenitors by immunofluorescence (Figure 4D).Interestingly, Pdgfrβ and Acta2 were described as a markers of Smooth Muscle-Mesenchymal Cells (SMMCs), a novel Itga7+ Vcam-Pdgfrb+ Acta2+ cell subpopulation present in adult muscle that exhibits myogenic potential and promotes MuSC engraftment following transplantation (Giordani et al., 2019).
Western blot analysis confirmed that EOM progenitors produce higher levels of Caveolin1 (CAV1) and Cavin1 (Figure 4E), which are co-expressed in caveolae and shown to be downregulated upon differentiation of rhabdomyosarcoma cells (Faggi et al., 2015).CAV1 was also shown to be a marker of quiescent and Pax7+ activated mouse MuSCs (Gnocchi et al., 2009) whereas in human, the CAV1+ MuSC subpopulation is associated with ECM organization, expression of quiescent markers and display increased engraftment after transplantation (Barruet et al., 2020).SPARC, MGP and IGFBP7 were also found to be upregulated in EOM activated MuSCs by Western blot (Figure 4E).The activity of these proteins seems to be context dependent and they can promote or suppress proliferation in different cell types (Ahmad et al., 2017;Artico et al., 2021;Cho et al., 2000;Jing et al., 2019;Kuronuma et al., 2020;Li et al., 2020a;Melouane et al., 2018;Said et al., 2013).Given the large number of matrisome genes differentially expressed in EOM progenitors, we also examined the expression of MMP2, a matrix remodelling protein (Gonçalves et al., 2022).Western Blot analysis revealed an enrichment in the active form of MMP2 in EOM progenitors whereas the TA progenitors displayed an enrichment of the pro-MMP2 or latent form (Figure 4E).This result is in agreement with the observation that activation of MMP2 in vitro and in vivo increases the proportion and mobility of Pax7+ cells (Mu et al., 2010).
Altogether, this analysis showed that EOM progenitors have an unusual transcriptome profile, express a wide range of ECM-related factors and harbour a mesenchymal signature that may endow them with a higher remodelling capacity of their immediate niche and promotion of an undifferentiated state.

EOM transcriptome profile is associated with a unique transcription factor network
To identify the underlying programs of myoblast heterogeneity between the two muscle groups, we inferred single cell regulatory networks using pySCENIC (Aibar et al., 2017;Sande et al., 2020).SCENIC uses co-expression patterns and transcription factor binding motifs to expose 'regulons" (transcription factors and their putative targets).This pipeline generates a new matrix comprising the activity of each regulon, in each cell (Figure 5A, Suppl Table 6).As expected, regulons associated with myoblast differentiation such as Myod, Myog and Mef2 family were found to be specifically active in differentiated cells of both EOM and TA (Figure 5B).Of note, the top 5 regulons of TA progenitors were also found to be active in EOM progenitors, whereas the top 5 regulons of EOM progenitors were unique to this cluster (Figure 5B).EOM progenitors were characterised by a unique signature of regulons involved in connective tissue/ECM remodelling including Egr1 (Havis and Duprez, 2020b) and Creb3l1, a downstream effector of Thr signaling (García et al., 2017) that plays a critical role during bone development and activates extracellular matrix genes such as Col1a1 and Fn1 (Sampieri et al., 2019).
Interestingly also, for some of these EOM progenitor regulons, the genes themselves were already enriched in quiescent EOM MuSCs and their expression was upregulated upon activation (Suppl Fig

5A, B).
To visualise the potential interactions between these regulons, we built a network restricted to transcription factors.To do so, target genes that were not at the head of regulons themselves were removed (Grimaldi et al., 2022) (see Material and methods) (Figure 5C).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 aggregated based on the number of shared edges.Thus, we were able to highlight associated and coregulating transcription factor modules.Strikingly, the transcription factors of the most specific regulons of each cluster preferentially organised as tightly related modules (Figure 5C).In agreement with our previous analyses, the known co-regulating transcription factors found in differentiated cells (Myod, Myog, Mef2a, Mef2c, Myf6) formed a tight module in this visualisation (Fig 5C, common differentiating category).Interestingly, the TA muscle progenitor module was composed of genes that were reported to be required for limb embryonic development (Hox genes, Lbx1) (Gross et al., 2000b;Pineault and Wellik, 2014;Swinehart et al., 2013).Notably, Hoxa11, is a determinant of limb axial identity during embryonic development (Zakany and Duboule, 2007), and HoxA and HoxC clusters were found as signatures of adult TA MuSCs (Evano et al., 2020b;Yoshioka et al., 2021); Suppl Fig5B).
In this context though, the network seems to maintain a proliferative progenitor state in activated EOM MuSCs cells as assessed by the higher protein levels of PAX7, CCND1 (Cyclin D1), together with other EOM regulon TFs (FOXC1, EBF1 and CREB3L1) (Suppl Fig5C, D).The vicinity of Hes1 to the EOM regulatory module prompted us to explore its direct upstream regulators (one level of ancestry in the network) (Figure 5D).Strikingly, 8 potential activators (including Hes1 itself) were part of the EOM module.By contrast, only 2 upstream regulators were part of the TA regulatory module (Figure 5D).Of note, Fos and Egr1 have also been reported as part of a stress signature following tissue dissociation (Machado et al., 2020).Yet, in our dataset, the expression of these genes (i.e. the StressIndex) was correlated with the expression of Pax7 and anti-correlated with the expression of Myod (Suppl Fig5E,F).
Moreover, when we applied a regression of the StressIndex or removed these genes from the count matrix, the general aspect of the data did not change (Suppl Fig5G), pointing at a role of these genes in the EOM progenitor maintenance independently of the stress response.

EOM features are present during the growth phase in vivo and retained upon passages in vitro
From the analysis above, it appears clear that EOM MuSCs possess an unusual transcriptomic state that may repress myogenic commitment and maintain a more "stem-like" state.We then wondered whether these EOM features are cell-wired and persist upon passage in vitro.Cells were isolated from Tg:Pax7-nGFP;Myog ntdTom mice and analyzed upon in vitro activation or upon one or two passages (Figure 6A, Suppl Fig 6A, B).While the total cell number was reduced with passages for both EOM and TA MuSCs (Suppl Fig 6A , B), the normalized cellular output was consistently higher for the EOM and correlated with a lower tdT/GFP ratio (Figure 6B, C).Real-time qPCR analysis revealed that the expression of a great number of EOM specific regulon TFs (Foxc1, Sox4, Ebf1, Creb3l1) and genes identified by the matrisome or molecular functions (Bgn, Sparc, Igfbp2, Igfbp7, Pdgfrb) were retained or even increased in EOM cultures after several passages compared to activated EOM MuSCs cultured only for 4 days (Figure 6D).These results were confirmed at the protein level by Western blot analysis (Figure 6E, F) where the specific EOM signature described above, that seems to maintain a more proliferative state in activated EOM MuSCs cells, persists even after several passages in vitro.Indeed, EOM MuSCs cells retained higher protein levels of PAX7, EOM regulon TFs (Figure 6E) as well as EOM matrisome or molecular functional genes (Figure 6F) compared to TA MuSCs cells at same passage.These results show that EOM MuSCs cells retain a non-canonical signature that is hard-wired even after extended cell culture.
Next, we asked whether EOM progenitor features are already present in activated fetal development and postnatal stages, where extensive muscle growth and MuSC expansion occurs, or if these features were acquired in adulthood upon reactivation from the quiescent state.To address these possibilities, we first isolated EOM and TA MuSCs from E18.5 and P21 Tg:Pax7-nGFP mice and performed qRT-PCR for Myogenin after 48h in culture (Suppl Figure 7A).As observed for adult MuSCs (Figure 1D), EOM cultures were less differentiated.However, as GFP protein persists during myogenic commitment (Sambasivan et al., 2009a)and this can introduce a bias in the initial populations, we isolated the GFP+/tdT-fractions from Tg:Pax7-nGFP;Myog ntdTom mice at P7 by FACS (Figure 7A, Suppl Figure 7B) and performed RT-qPCR.Similarly to the in vitro sc-RNAseq data, significantly lower levels of Myogenin were detected in EOM GFP+/tdT-cells (Figure 7B).Finally, we confirmed higher transcript levels for the EOM specific regulon TFs (Foxc1, Ebf1, Sox4, Creb3l1) and matrisome or molecular functional components identified by sc-RNAseq of activated MuSCs in vitro (Figure 7B).Altogether, our results indicate that EOM MuSCs repress myogenic commitment and maintain a more "stem-like" state upon activation, and this property is conserved from development to adulthood.

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 sc-RNAseq dataset (Figure 5 A,B, Suppl Figure 3B, Suppl Figure 5B, D) and expressed continuously in the EOM lineage.Notably also, bulk RNAseq from adult muscles (Terry et al., 2018) showed that Foxc1 expression is higher in entire EOMs than in any of the other 10 muscle groups analyzed (Suppl

Fig 7C).
Finally, immunostaining on tissue sections of Pax7 CreERT2 :R26 tdTOM mice displayed colocalization of FOXC1 with tdTOMATO in EOM MuSCs and expression of FOXC1 in EOM myonuclei at P10 and in adulthood but not in the TA (Figure 7C, D).Thus, differences in Foxc1 expression between the EOM and TA seem to arise during development and are kept upon activation, and throughout the myogenic lineage.
As differences in Myogenin expression at the mRNA level are already evident at Day2 upon in vitro activation (Figure 1D), we decided to use this early timepoint to functionally validate a potential role of Foxc1 in the maintenance of the progenitor state.Thus, we silenced Foxc1 via a siRNA approach in EOM activated MuSCs (Figure 7E-G).While the total cell number and percentage of Pax7 cells was not changed upon Foxc1 silencing, immunostaining for FOXC1 showed silencing of 82% of the cells at the protein level (Figure 7F) and RT-qPCR showed a 64% reduction at the transcript level (Figure 7G).
Notably though, a strong correlation between PAX7 and FOXC1 intensity of fluorescence was observed in control cells but this lost upon silencing (Suppl Figure 7D).Moreover, a 1.7 fold increase in Myogenin expression was detected by qPCR following silencing of Foxc1 as early as Day2 (Figure 7G).As Foxc1 direct targets (Suppl Table 7) include matrisome components (e.g.Sparc, Pdgfrb, Fbn1, Nrp1, Figure 4B) and other EOM specific regulon TFs (e.g.Ebf1, Creb3l1, Egr1, Figure 5B), our findings suggest that Foxc1 maintains EOM MuSCs in a progenitor like state through matrix remodeling and cooperation with other TFs.
Different studies have reported that ECM secretion and receptor ligand interactions can direct MuSC behaviour in a cell-autonomous manner (Baghdadi et al., 2018;Tierney et al., 2016).To further investigate the role of Foxc1 targets in EOM MuSCs mesenchymal features and proliferative capacity we decided to focus on Pdgfrβ, a tyrosine-kinase receptor commonly expressed by pericytes (Hellström et al., 1999;Soriano, 1994) that is involved in multiple signalling pathways that regulate cell migration, proliferation and differentiation (Heldin and Westermark, 1999).Since Pdgfrβ is differentially expressed in EOM progenitors (Figure 7H,I), we decided to isolate this subpopulation by FACS from in vitro activated EOM cells using a PDGFRβ antibody.The flow cytometry analysis corroborated the scRNAseq and immunofluorescence data (Figure 4D), where PDGFRβ was enriched in EOM samples (Suppl Figure 7E,F).To distinguish the functional properties of PDGFRβ+ and PDGFRβ-EOM cells, we sorted both subpopulations from EOM cultures, which represent respectively 83 and 17% of total nuclei at D4, and cultured them in growth promoting conditions for additional 24h (Figure 7J-L).Cells were treated with the thymidine analogue 5-Ethynyl-2'-deoxyuridine (EdU) for 2h prior to fixation to monitor cell proliferation.PDGFRβ+ cells had a significantly higher proliferative capacity as measured by the total nuclear output and incorporation of EdU compared to PDGFRβ-cells (Figure 7K,L).Although the percentage of MYOD-expressing cells was not overtly different between the two populations, MYOG was expressed by over 50% of PDGFRβ-cells, whilst being detected in only 7% of PDGFRβ+ cells (Figure 7K,L).Therefore, the PDGFRβ+ myoblast subpopulation is characterised by a higher proliferative potential and decreased differentiation status.

Transcription dynamics expose EOM and TA disparities in progenitor-state maintenance
We first attempted to correlate the number of regulons responsible for most of the matrisome gene expression differences between EOM and TA upon activation (Figure 8A-C, Suppl Table 8).We observed that activated EOM MuSCs consistently regulated a higher number of matrisome genes than TA cells independently of the number of regulons considered (Figure 8A).However, the ratio of the number of regulations of matrisome genes between EOM and TA activated MuSCs showed that the difference was maximal when considering the top 5 regulons, Alx4, Dmrta2, Foxc1, Zmiz1 and Fos (Figure 8B, C).Notably, Dmrta2, Foxc1 and Fos were also active regulons in quiescence, with slight disparities between EOM and TA (Figure 8D, Suppl Table 9).
Second, we set out to determine whether matrix genes underlie the transition towards progenitors and committed cells during activation.To do so, we inferred RNA velocity using scVelo (Bergen et al., 2020).This method computes local changes in the relative amount of unspliced and spliced variants, which depend on the rates of transcription, degradation and splicing (Manno et al., 2018a).ScVelo identifies candidate "driver genes", i.e. the most transcriptionally dynamic genes in a given cluster that are responsible for most of the inferred RNA velocity, and thus underlie its transitional trajectory.We applied the scVelo pipeline to EOM and TA datasets independently, and examined the directional trajectories underlying commitment and progenitor maintenance (Figure 8E,F).Two distinct velocity streams stood out in both datasets, towards differentiation (Myog High ) and towards a progenitor-state (Pax7 High ).
Strikingly, a larger fraction of EOM cells appeared to transition towards a progenitor state compared with the TA (Figure 8E,F).Conversely, most cells in the TA appeared to be directed towards differentiation, consistent with an overlapping transcriptomic profile with differentiating cells, shown as top DEGs (Figure 2D) and regulatory networks (Figure 5C).These trajectories did not appear to be specifically correlated with cell cycle phases (Suppl Figure 8 A,B), which was shown to influence transcriptomic data in some cases (McDavid et al., 2016).Hence, the velocity streams observed are most likely to 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 (Figure 8G-I, Suppl Table 10, 11).Out of the top 100 driver genes, 30 were common to both datasets, including ColV, which plays a critical role in maintenance of quiescence (Baghdadi et al., 2018a) (Figure 8G,H).We assessed the GO molecular functions associated with these driver genes to identify the unique and shared pathways employed by both datasets during myogenic progenitor maintenance (Figure 8J).In agreement with our previous results, this transition was characterised by the active upregulation of ECM components specifically in the EOM.
In summary, our results establish a new regulatory network permitting an enhanced expansion capacity and delayed commitment to differentiation of EOM MuSCs upon activation.This regulatory module is hardwired to EOM identity and relies on the active maintenance of progenitor MuSC features.

DISCUSSION
One of the fundamental, yet unanswered questions in stem cell and regenerative biology regards the mechanisms that allow an orchestrated balance between proliferation, differentiation and self-renewal.
An unusual feature of skeletal muscle stem cells is their reliance on distinct GRNs in different anatomical locations.How these pathways confer the unique functional properties of each stem cell population in cranial and limb/limb muscles remains poorly understood.Indeed, most studies on adult and developmental 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.In this study we performed scRNA-seq accompanied by experimental validations to characterise the GRN underlying the functional heterogeneity of MuSCs derived from different muscles, taking as EOM and TA as archetype muscles for high and low performing stem cell populations.

EOM MuSCs are more refractory to in vitro differentiation
Although subsets of cranial derived MuSCs were previously described to be more proliferative and to have a higher engraftment capacity than trunk MuSCs (Ono et al., 2010b;Randolph and Pavlath, 2015b;Stuelsatz et al., 2015a), the mechanisms responsible for these phenotypic differences remained elusive.
As these features depend on the balance between proliferation and differentiation, one possibility is that different upstream regulators define the pace at which MuSCs progress in the lineage in diverse anatomical positions.
Here, by monitoring Myog expression, we showed that EOM MuSCs have a lower propensity to differentiate following activation and thus persist as a proliferative population.As foetal and early postnatal EOM MuSCs are also refractory to differentiation, we propose that this property might be hardwired to some extent by unique GRN that are retained throughout development and adulthood.
Notably, foetal MuSCs were shown to be more resistant to myogenic progression upon in vitro expansion and contribute more efficiently upon transplantion than the adult counterparts (Sakai 2013,Tierney 2016).In addition, a recent study showed that limb MuSCs cells isolated at birth displayed prolonged expansion rate and a delay in differentiation and fusion compared with those isolated at later stages of postnatal growth (P7 and P15) and adulthood (Gattazzo et al., 2020).Therefore, is tempting to speculate that EOM MuSCs retain features of the respective foetal and neonatal precursor cells.
EOM MuSCs are strong candidates for cell-based therapies given their robust engraftment capacity in vivo (Stuelsatz et al., 2015b).So far, major obstacles in the clinic are the large cell numbers required for transplantation and the fact that ex vivo amplification of trunk MuSCs leads to a drastic decline in regenerative potential due to commitment to differentiation (Briggs and Morgan, 2013;Ikemoto et al., 2007).Substantial advances have been made by modifying the MuSC culture conditions (Charville et al., 2015;L'honoré et al., 2018), or by the use of teratoma-derived skeletal myogenic progenitors (Xie et al., 2021).As such, the identification of factors that regulate cell fate decisions in distinct MuSC populations serves as a resource for advancing knowledge in the context of regenerative medicine.

EOM progenitors exhibit a mesenchymal genetic signature upon activation
Recent scRNA-seq analyses have 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., 2019b;Machado et al., 2021;Micheli et al., 2020;Yartseva et al., 2020).To characterise these processes in cranial MuSCs, we performed a comparative scRNA-seq analysis of activated EOM and TA MuSCs.Distinct transcriptional profiles divided the myoblast pool into two subpopulations: those that resembled progenitors and those that were differentiating.The EOM progenitor population was characterised by a higher Pax7 levels, and thus reminiscent of in vitro reserve cells (Laumonier et al., 2017;Yoshida et al., 1998;Zammit et al., 2004) and ontology analysis of its DEGs revealed a surprising enrichment in ECM organisation processes and Pdgf signalling.
By surveying bulk RNA sequencing and single-cell profiling data of trunk skeletal muscle stem cells, a recent study highlighted a dynamic profile of Pdgf ligands and receptors during myogenesis (Contreras et al., 2021).In addition, treatment with NOTCH and PDGFRβ ligands (DLL4 and PDGF-BB respectively) was shown to enhance migration, expression of stem cell markers and perivascular-like features in MuSCs (Gerli et al., 2019) and embryonic myoblasts (Cappellari et al., 2013).Given that PDGFRβ and several Notch pathway components are co-expressed in EOM progenitors, and at higher levels than the TA, we speculate that cross-talk between these pathways could take place in this subpopulation.
In addition to Pdgfrβ, which is a well-known pericyte marker (Lindahl et al., 1997), EOM progenitor cells express Acta2 which is involved in the contractile apparatus of smooth muscle and also used as a SMMCs (Pdgfrb+ Itga7+ Vcam-Pdgfrb+ Acta2+) and FAPs (Giordani et al., 2019;Joe et al., 2010;Uezumi et al., 2010) marker.However, while the EOM progenitor fraction expressed Fbn1 and Loxl1, which are markers of a new subtype of FAPs identified by scRNA-seq of whole muscle (Rubenstein et al., 2020), they were negative for canonical FAPs and pericyte markers such as Pdgfrα (Joe et al., 2010;Uezumi et al., 2010) and Cspg4 (Ng2) (Birbrair et al., 2013), respectively.In addition, the EOM progenitor fraction also related more to fetal MuSCs than the TA counterparts in terms of enrichment in matrisome components such as Fn1, Fbn1, Vcam and collagens (Tierney et al., 2016b).It is well established that cell autonomous ECM deposition and remodelling are key instructive cues governing cell fate decisions in MuSCs (Baghdadi et al., 2018b;Lukjanenko et al., 2016;Mashinchian et al., 2018;Tierney et al., 2016b).Thus, we hypothesize that EOM progenitors might exploit these features to preserve self-renewal potential and mimic their in vivo stem cell niche even after in vitro culture.In the EOM, this plasticity has already been observed in vivo, where EOM myogenic progenitors transition towards non-myogenic cell fates in areas deprived of neural crest cells during mouse embryogenesis (Grimaldi et al., 2022).Additionally, a resident subpopulation expressing reduced levels of canonical myogenic markers but increased levels of mesenchymal markers has been detected in embryonic and foetal human limb muscle (Xi et al., 2020b).This plasticity seems to be present also in vitro, where directed myogenic differentiation of hiPSCs generates heterogeneous cell types including both myogenic and non-myogenic cells (Xi et al., 2020b).

A unique network of transcription factors maintains EOM progenitors
Analysis of single cell regulatory networks (SCENIC) revealed a unique set of regulons whose activity was restricted to EOM progenitors.Of note, the majority of the identified head regulon genes in EOM progenitors are not typical myogenic TFs.One of the most active regulons in this cluster was Foxc1, which was described as a pro-mitogenic factor in the cancer field (Yang et al., 2017) and as driver of endothelial/smooth muscle fate alone (Whitesell et al., 2019b) or in combination with Foxc2 (Mayeuf-Louchart et al., 2016).Another top regulon was Egr1, which was shown to promote the expression of many ECM-related genes in different tissues (Gaut et al., 2016;Guerquin et al., 2013;Milet et al., 2017;Wu et al., 2017).Interestingly, Foxc1 and Ebf1 were found as active transcription factors underlying the transition of myogenic towards non-myogenic cell fates in the EOM anlage during embryonic development (Grimaldi et al., 2022).Both of these transcription factors were found in EOM progenitors in this study, suggesting that they might be continuously active from development to adulthood specifically in the EOM.Interestingly, expression of Foxc1 by hair follicle stem cells activates quiescent gene networks upon activation (Wang et al., 2016).As Foxc1 is expressed both by quiescent and activated EOM MuSCs, this gene might reinforce a progenitor/stem-cell identity throughout the lineage.
Notably, Pitx2 appears as a target of Foxc1 in our analysis (Suppl Table 7) and its overexpression promotes MuSC proliferation in vitro and enhances the regenerative potential of dystrophic MuSCs in vivo (Vallejo et al., 2018)Other non-canonical myogenic TF identified in EOM progenitors are the Kruppel-like factor family members Klf4, Klf6, and Klf9 with Klf6 identified as top regulons.This gene family participates in the development and homeostasis of numerous tissues (McConnell and Yang, 2010), and Klf4 is pioneer factor well-known for its reprogramming capabilities (Schmidt and Plath, 2012) When assessing the interactions between transcription factors, the EOM Progenitor regulon TFs clustered together in a regulatory module including Pax7 and the Notch target Hes1, supporting our hypothesis that these cells are in a more stem/progenitor state.This module also included the transcription factor Fos/AP-1, whose expression in a subset of freshly isolated trunk MuSCs is related to enhanced regenerative properties such as rapid cell cycle entry, higher clonal capacity and more efficient engraftment (Almada et al., 2021).Further, by means of an inducible FOS expression system, it was reported that elevated FOS activity in recently activated muscle progenitor cells perturbs cellular differentiation by altering the 3D chromosome organization near critical pro-myogenic genes (Barutcu et al., 2022).Altogether, it is tempting to speculate that a common regulatory module allows the establishment and maintenance of a MuSC subpopulation in the EOMs that displays mesenchymal features, partially overlaps with the quiescent state, and is refractory to differentiation.
To investigate the maintenance of this stem state during EOM activation we used RNA velocity to infer the different directional trajectories and their driver genes.The proportion of cells transitioning to a more progenitor state was higher in the EOM when compared to TA, and this progenitor stream was driven by a specific gene set.Amongst these genes there was an enrichment of ECM-related factors, underscoring its active role in the maintenance of a less differentiated state.Interestingly, we identified Igfbp7 (Insulin growth factor binding protein 7), as unique EOM driver gene, which has been found to be a specific marker of quiescent MuSCs (Fukada et al., 2007).By binding IGFs (insulin growth factors), IGFBPs regulate the action of IGFs and repress myogenic differentiation (Damon et al., 1997;JONES and CLEMMONS, 1995;Li et al., 2020b).IGFBP7 has low affinity for IGF and higher affinity for Insulin and Activin (Jin 2020).IGFBP7 also binds Collagen IV and Heparan sulphate on the cell surface (Jin 2020).This suggests a unique mode of repression of myogenesis through the regulation of IGF intake in EOM progenitors.Taken together, our results highlight an active maintenance of a progenitor pool that partially overlaps with the quiescent state during EOM MuSC activation, possibly through the steady induction of inhibitory differentiation signalling cues and promotion of ECM remodelling.

Conclusion
Using in vitro and in silico approaches, we propose a model where the outperformance of EOM MuSCs depends on the expression of a tightly associated module of transcription factors regulating a characteristic and 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 maintenance of 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 gene regulatory network within this muscle group.In addition, our study could serve as basis for the optimization of in vitro protocols that would exploit cranial MuSC properties in the context of cell-based muscle therapies.

Animal care
Animals were handled according to national and European Community guidelines and an ethics committee of the Institut Pasteur (CETEA) in France approved protocols.Except when indicated otherwise, males and females of 2-4 months were used.Tg:Pax7-nGFP (Sambasivan et al., 2009b), Myog ntdTom (Benavente-Diaz et al., 2021), Pax7 CreERT2 (Mathew et al., 2011)and R26 tdTom (Madisen et al., 2009)mouse lines were maintained in a C57Bl/6JRj background.To induce recombination of Pax7 CreERT2 :R26 tdTom 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 during 5 consecutive days and animals sacrificed 5 days later.To induce recombination of pups, the Tamoxifen stock solution was diluted to 15mg/ml with sunflower seed oil and 0.15 mg were administered daily by subcutaneous injection between P4 and P6 daily.Pups were sacrificed at P10 by decapitation and adult mice by cervical dislocation.

MuSC isolation by FACS
Muscles were dissected and minced in ice-cold DMEM.Samples were then incubated in DMEM (Fisher Scientific, 11594446), 0.08% Collagenase D (Sigma, 11088882001), 0.2% Trypsin (ThermoFisher, 15090) and 10 µg/mL of DNAse I (Sigma, 11284932) for 25 minutes at 37°C under gentle agitation for 5 rounds of digestion.After each round, samples were allowed to precipitate for 5 minutes, the supernatant was collected on 4mL of Foetal Bovine Serum (FBS) on ice and fresh digestion buffer was added to the remaining muscle pellet.The collected supernatants were centrifuged for 15 minutes at 550 g at 4°C, resuspended in DMEM 2% FBS and filtered through a 40µm strainer (Corning, 352235) before cell sorting.
For myoblast isolation based on PDGFRβ expression, dishes were washed twice with PBS prior to trypsinization.Cells were subsequently incubated with α-PDGFRβ antibody in HBSS (Fisher Scientific, 11560616) media with 2% FBS for 30 minutes on ice.After two washes with HBSS 2% FBS, cells were incubated with streptavidin-coupled BV786 and incubated on ice for 25 minutes.After one wash with HBSS 2% FBS, the cell suspension filtered through a 40µm strainer (Corning, 352235) before cell sorting.Cells were isolated based on size, granulosity and GFP or BV786 fluorescence on using an Aria III (BD Biosciences) machine.

Muscle stem cell culture, treatment and transfection
Cells isolated by FACS were plated onto Matrigel® (Corning, 354248) coated dishes and cultured in growth media at 3% O2, 5% CO2, 37°C for the indicated time.To assess proliferation, cells were pulsed with 10 -6 M of EdU (ThermoFisher, C10640) in cell culture media for 2h prior to fixation.
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 (Gibco) as described by the manufacturer.Briefly, a pre-mix of siRNA/Optimem (1.5ul /20ul) and Lipofectamine3000/Optimem (0.3ul / 20ul) were incubated separately 5 min RT, mixed at 1:1 ratio and incubated 15 min more at RT. 2.10 4 cells in 40ul of Optimem were incubated with an equal volume of the transfection mix for 2h at 3% O2, 5% CO2, 37°C in Eppendorf tubes whose caps had 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 10k/cm2 in Matrigel coated wells containing growth media.Two days upon transfection wells were processed for immunostaining or RNA collected with Tryzol as above.

RT-qPCR
RNA from in vitro myoblast cultures was extracted using a Trizol-based kit (Zymo Research, R2061) and reverse transcribed using SuperScriptIII (Invitrogen, 18080093).RT-qPCR to assess for mRNA relative expression performed with SYBR green master mix (Roche, 04913914001) in an Applied biosciences machine.Data analysis was performed using the 2-ΔΔCT method (Livak and Schmittgen, 2001) and mRNA expression was normalized with Rpl13.

Immunofluorescence
For immunostaining of ECM markers, cells were fixed in warm (37°C) 3% paraformaldehyde (PFA, Electron Microscopy Sciences) 2% sucrose.Upon incubation in the fixative for 15min at RT, cells were washed three times with PBS and incubated with NH4Cl 100mM in PBS to quench aldehyde groups.
Cells were washed three times in PBS, permeabilized in 0.5% Triton X-100 (Merck, T8787) for 5 min RT, washed again three times in PBS and blocked with 10% goat serum (GIBCO).Cells were incubated with the indicated primary antibodies overnight at 4°C in PBS with 2% goat serum.Upon PBS washes, cells were incubated with secondary antibodies and Hoechst (ThermoFisher, H1399) for 45 min at RT. EdU was detected after antibody staining if needed according to manufacturer instructions (ThermoFisher, C10640).
Following immunostaining, 96-well plates (DUTSCHER DOMINIQUE, 655090) were imaged on an Opera Phenix high-content microscope (Perkin Elmer) using a 20x or 40x objective.Quantification of acquired images was performed using Harmony (Perkin Elmer) or Columbus software, using an automated pipeline.First, nuclei were detected based on Hoechst signal.Number of nuclei was assessed and mean and maximum intensity were automatically quantified on the nuclear region for nuclear markers.For acquisition of ECM markers either the Opera Phenix high-content microscope (Perkin Elmer) or Zeiss LSM800 microscope with ZEN software were used.

Image analysis
Cell tracking was performed using the Manual Tracking feature of the TrackMate plug-in (Tinevez et al., 2017) in Fiji (Schindelin et al., 2012).Fiji and Harmony (Perkin Elmer) were used for image analysis.
Figures were assembled in Adobe Photoshop, Illustrator and InDesign (Adobe Systems).

Data analysis and statistics
Data analysis, statistics and visualisations were performed using Prism (Graphpad Software) or using R (Team, 2014) and the package ggplot2 (Wickham, 2009).For comparison between two groups, two tailed unpaired Student's t-tests were performed to calculate p values and to determine statistically significant differences (see Figure Legends).

scRNAseq data generation
MuSCs were isolated on BD FACSAria™ III based on GFP fluorescence and cell viability from Tg:Pax7-nGFP mice (Sambasivan et al., 2009b).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 four days.Activated MuSCs were subsequently trypsinized and washed in DMEM/F12 2% FBS.Live cells were re-sorted, manually counted using a hemocytometer and processed for scRNA-seq.
Prior to scRNAseq, 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 and cDNA libraries were generated following 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.
Following normalisation and quality control, we obtained an average of 5792 ± 1415 cells/condition.
Cells with more than 10% of mitochondrial gene fraction were discarded.4000-5000 genes were detected on average across all 4 datasets.Dimensionality reduction and UMAPs were generated 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, number of UMI following 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 2020) using the "AddModule" Seurat function.94 out of 98 genes were detected in the combined datasets.
UMAPs were generated after 1. StressIndex regression, and 2. after complete removal of the detected stress genes from the gene expression matrix before normalization.In both cases, the overall aspect of the UMAP did not change significantly (Figure S5).Although immeasurable confounding effects of cell stress following 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., 2018b).

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 (Manno et al., 2018b) 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 Suppl Tables 10 and 11.

Gene regulatory network inference and transcription factor modules
Gene regulatory networks were inferred using pySCENIC (Aibar et al., 2017;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).Note that this function can yield less 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.Visualizations 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 2 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.Visualization of this network was performed in a force-directed graph using Gephi "Force-Atlas2" algorithm (https://gephi.org/).

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 a p-value lower than 0.05 were displayed.

Microarray data and analysis
EOM, TA, Esophagus, Masseter and Diaphragm muscles from three adult TgPax7nGFP mice were digested and processed for FACs as described before.GFP+ cells were collected and RNA extracted using the RNeasy Micro Kit (Qiagen).After validation of the RNA quality with Bioanalyzer 2100 (using Agilent RNA6000 pico chip kit), 450 pg of total RNA is reverse transcribed following the Ovation Pico WTA System V2 (Tecan).Briefly, the resulting double strand cDNA is used for amplification based on sc-RNAseq.AG and MBD run the bioinformatic analysis with help from SM and VL.AG, MBD and GC wrote the manuscript with input from DDG and ST.GC and ST provided funding.Supplementary Tables Table 1.Top 100 marker genes of EOM and TA of in vitro activated MuSCs    F-G.Venn diagrams of the overlapping differentially expressed genes between quiescent and global in vitro activated datasets for EOM (F) and TA (G).

FIGURES AND
H. Expression plots of selected EOM markers exclusive to the quiescence state (Unique Q), common to quiescent and activated states (conserved) and exclusive to activation (Unique Act) from the analysis in F.

Figure 1 .
Figure 1.EOM and TA MuSCs exhibit functional differences following activation A. MuSCs were isolated by FACS based on GFP fluorescence from Tg:Pax7-nGFP mice, plated and pulsed with EdU before fixation and immunostaining at Day (D) 3, D4 and D5.B. MuSCs from EOM and TA were cultured for 72h and imaged for 48h.Representative overlayed DIC and red fluorescence channel images at selected time points showing a temporal increase of tdTOMATO fluorescence.Scale bar 25 μm.C. Quantification of selected timepoints in (B).Graph represents mean percentage of tdTOMATO + cells over total number of cells at each timepoint for each muscle.Percentage of tdTOMATO + cells was compared between

Figure 2 .
Figure 2. Single cell transcriptome signatures of activated EOM and TA MuSCs

Figure 3 .
Figure 3. Molecular signature of EOM and TA quiescent MuSCs A. Experimental scheme of sc-RNAseq pipeline for quiescent EOM and TA MuSCs.B. UMAP visualization of EOM and TA quiescent clusters.C-D.Violin plot (C) and expression plot (D) of Myf5, Pax7, Myod and Myog.E. Heatmap representing top differentially expressed genes in EOM and TA quiescent clusters and their expression levels across all cells.

Figure 5 .
Figure 5. Distinct gene regulatory networks underlie EOM and TA activation dynamics A. Top 4 regulon activity for each cluster of in vitro activated sc-RNAseq dataset overlaid onto UMAP representation.B. Heatmap of Top 7 regulons in each cluster of the activated sc-RNAseq dataset with activity level in each cell C. Transcription factor network highlighting top regulons of each cluster as well as common regulons (35 regulons maximum) of activated sc-RNAseq dataset.Proximity of the nodes in the network indicate a higher number of shared edges, highlighting core modules.D. Filtered network showing Hes1 direct regulators (nodes with 1 level ancestry).

Figure 6 .
Figure 6.EOM properties are retained despite several passages in vitro.A. Scheme of isolation and passages of EOM and TA MuSCs.Cells were isolated by FACS based on GFP fluorescence from Tg:Pax7-nGFP;Myog ntdTom and cultured for 3 days (Act).The entire well was passaged (P1, P2) every 3 days.Following each passage, a fraction of cells was analysed by FACS to assess GFP and tdTom ratios.B. EOM/TA normalized cell number ratios.Total cell numbers per well were counted upon trypsinization at day 3 (D3), days 5 (D5) and at passages P1 and P2.EOM/TA normalized cell number ratio was generated for each time point analysed (n=3 mice).For normalization, the ratio between cell numbers counted per well and averaged EOM cell number in all wells at each time point was calculated.C. FACS analysis of Tg:Pax7-nGFP;Myog ntdTom EOM and TA MuSCs upon 5 days in culture (D5), one or two passages (P1, P2) (n>3 mice).Graph represents ratio of TdT + /GFP + cells for each time point.Two-tailed unpaired Student's t-test; *p-value<0.05,** p-value<0.01.D. qRT-PCR at D3(Amp), P1, P2 on whole populations present in dish for key myogenic markers, regulon TFs and Matrisome/Molecular function genes identified in EOM progenitors (n=4 mice).E. Western blot of total P2 population of EOM and TA MuSCs (n=3 mice; one mouse per lane).Note higher expression of Pax7 and EOM progenitor specific TF regulons (FOXC1, CREB3L1, EBF1) identified in activated EOM MuSCs by sc-RNAseq.ACTIN was used for normalization of protein loading.

Figure 8 .
Figure 8. 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 sc-RNAseq datasets.