Single-nucleus RNA-seq and FISH identify coordinated transcriptional activity in mammalian myofibers

Skeletal muscle fibers are large syncytia but it is currently unknown whether gene expression is coordinately regulated in their numerous nuclei. Here we show by snRNA-seq and snATAC-seq that slow, fast, myotendinous and neuromuscular junction myonuclei each have different transcriptional programs, associated with distinct chromatin states and combinations of transcription factors. In adult mice, identified myofiber types predominantly express either a slow or one of the three fast isoforms of Myosin heavy chain (MYH) proteins, while a small number of hybrid fibers can express more than one MYH. By snRNA-seq and FISH, we show that the majority of myonuclei within a myofiber are synchronized, coordinately expressing only one fast Myh isoform with a preferential panel of muscle-specific genes. Importantly, this coordination of expression occurs early during post-natal development and depends on innervation. These findings highlight a previously undefined mechanism of coordination of gene expression in a syncytium.

S keletal muscles constitute 40 to 50% of our body mass and are principally composed of very long myofibers attached to bones via tendons, innervated by motoneurons, and surrounded by vessels and fibroblasts cells. The myofibers themselves are syncytia composed of hundreds of post-mitotic nuclei sharing the same cytoplasm, generated by the fusion of myoblasts during development 1 . However, the protein content of myofibers is not homogenous along their length. The proteins required for the formation of the neuromuscular junction (NMJ) accumulate in the center of the myofiber whereas specific proteins for myotendinous junction specialization (MTJ) accumulates at the periphery of the fiber. Two mechanisms could account for this regionalization: the transport of specific proteins to these sites 2 , and the specialization of transcription within the nuclei in these sites 3,4 . At the NMJ, specific myonuclei express AchR and AchE 5 and at the MTJ-specific myonuclei express Col22a1 6,7 . However, the precise genetic program operating in these specific myonuclei and the transcription factors involved in this specialization during development, compared to the other myonuclei (called body myonuclei), is far from being understood.
In body myonuclei that constitute the vast majority of myonuclei in the fiber, heterogeneity in nuclear protein import 8 and stochastic gene expression 9 have been proposed to occur, suggesting that body myonuclei along the myofiber are not equivalent. Accordingly, gene transcription could occur stochastically (bursts of transcription), resulting in apparent uncoordinated gene expression at a given locus between the nuclei of the same fiber. Furthermore, in other syncytia such as osteoclasts and syncytiotrophoblasts, transcriptional activity is heterogeneous, varying from one nucleus to the other, indicating that their nuclei, although they share the same cytoplasm, are not coordinated 10,11 . Whether the hundreds of body nuclei along the myofiber are transcriptionally active at the same time and express the same set of genes remains debated.
Adult skeletal muscles are composed mainly of slow and fast myofibers and are usually classified by their MYH (myosin heavy chain) expression profile [12][13][14][15] . MYH is the primary determinant of the efficiency of muscle contraction: slow myofibers in mice express Myh7 (MyHC-I), and fast myofibers express Myh2 (MyHC-IIA), Myh1 (MyHC-IIX) or Myh4 (MyHC-IIB). These different Myh isoforms are coded by different genes, and in mammals, fast Myh (fMyh) genes, including Myh2, Myh1, Myh4, embryonic Myh3, neonatal Myh8, and extraocular Myh13 are clustered on a single locus while the slow Myh7 gene lies on another locus 16 . Myofibers can express one, two or three MYH isoforms at the same time at protein and mRNA level [17][18][19][20] . Myofibers expressing more than one MYH isoform are called hybrid fibers in contrast to pure fibers expressing only one isoform. Post-transcriptional mechanisms involving antisense RNA 21 or mir-RNA 22 may participate in the control of the accumulation of specific MYH proteins in both homogenous and hybrid myofibers. Alternatively, hybrid fibers may be constituted of regionalized Myh genes expression or of uniform Myh genes co-expression along the fiber. For the well-studied ß-Globin locus whose organization is reminiscent of the fMyh locus, it was shown that the Locus Control Region (LCR) establishes contacts and activates a single gene at the locus at a given time 23 . It was more recently demonstrated that two adjacent genes can be expressed in the same nucleus at a given time by a shared enhancer 24 . Whether one myonucleus can activate one, two or three fMyh genes at the same time remains to be established.
In this report, we sequence the RNAs and the chromatin accessibility from single nucleus (snRNA-seq and snATAC-seq) isolated from adult skeletal muscles to characterize the transcription profile of known cell populations present in adult skeletal muscles, including Fibro-adipogenic progenitors (FAPs), tenocytes, myogenic stem cells and myonuclei, establishing a blueprint of the transcriptional activity of all its nuclei. We further identify three main populations of myonuclei: NMJ, MTJ, and body myonuclei with specific transcriptional and chromatin accessibility. Our results show that most myonuclei express only one fMyh gene at the same time. By visualizing in vivo in mechanically isolated myofibers the localization of fMyh pre-mRNAs, we validate these results showing a coordination of fMyh expression in the myonuclei of the majority of myofibers of the fast extensor digitorum longus (EDL). We also identify a minority of hybrid myofibers with nuclei expressing several fMyh isoforms at the same time. Those hybrid myofibers are more abundant in slow soleus muscle and during denervation in the fast EDL. Last, we characterize how adult fMyh genes become activated during post-natal development and show the importance of innervation for their coordination.

Results
SnRNA-seq identifies several cell populations in adult skeletal muscles. To characterize the transcriptional profile of myonuclei within myofibers, we developed single nucleus RNA-seq (snRNA-seq) experiments using a droplet-based platform on purified nuclei from different skeletal muscles of adult mice ( Fig. 1a and Supplementary Fig. 1a-d). By using Seurat software 25 and Uniform Manifold Approximation and Projection (UMAP) clustering we were able to detect a total of 19 different cell populations clusters present in the different muscles ( Fig. 1b and Supplementary Fig. 2b, c). The origin of the nuclei was identified by specific markers (Fig. 1b Fig. 2d). To look at the heterogeneity of the myonuclei, we have then reclustered the myonuclei expressing Titin (Ttn) gene to remove all the non-fiber nuclei (Fig. 1c). Several populations of myonuclei where identified, among which MTJ myonuclei (expressing Col22a1), NMJ myonuclei (expressing Ache), a yet unknown myonuclei subpopulation (expressing Myh9, Flnc, Runx1), and body myonuclei accounting for 94% of all myonuclei (Figs. 1c and 2a). By comparing gene expression in body nuclei and specialized nuclei, we characterized numerous known (Ache and Etv5 for NMJ 5 ; Coll22a1 for MTJ 6 and unknown genes specific to these specialized myonuclei (Etv4 for NMJ, Maml2 for MTJ) (Fig. 2b). Unidentified, NMJ and MTJ myonuclei also expressed Myh genes (Fig. 2c). Slow myonuclei expressing Myh7 clustered independently from the fast Myh1/2/4 myonuclei. The fast myonuclei formed a large cluster that could be subdivided into several subpopulations, all expressing sarcomeric genes like Ttn and various Myh isoforms (Figs. 1c and 2a). Myh4 + myonuclei represented 51% of this fast myonuclei population. We also compared the expression of genes in slow and fast nuclei expressing the different types of Myh (Fig. 1d). Fast Myh2 and Myh1 nuclei showed a similar transcription program, contrasting with slow Myh7 and very fast Myh4 nuclei each showing a distinct transcriptional program (Fig. 1d). Interestingly, three main categories of Myh4 nuclei were identified, identifying transcriptional heterogeneity among these nuclei ( Supplementary Fig. 3a- d   Myh7  Tnnc1  Myl2  Tnnt1  Tnni1  Atp2a2  Myh7b  Tpm3  Mhrt  Lgi1  Gm28653  5730419F03Rik  Ablim1  Astn2  Rp1  Tsga10  Myl3  Gm122  Inpp4b Car3 Enox2 Adamtsl1 Ttn SnATAC-seq identifies heterogenous transcriptional programs in myonuclei. Recent studies profiled the genome-wide chromatin accessibility of adult skeletal muscles 26 . However, this global approach did not discriminate the chromatin status of each cell population present in skeletal muscles. Thus, we set up a protocol to perform single nucleus ATAC-seq (snATAC-seq) with purified nuclei from adult soleus and quadriceps using a droplet-based platform. Briefly, nuclei were transposed using Tn5 transposase and then encapsulated into gel beads emulsion using 10X protocol. After quality assessments and filtering (Supplementary Fig. 4a-c), 132 966 peaks were identified from 6037 sequenced nuclei. To identify the different populations of nuclei, we classified nuclei from ATAC-seq based on snRNA-seq experiments from Fig. 1b: we used methods for cross-modality integration of snRNA-seq and snATAC-seq data and label transfer from Signac software ( Fig. 3a and Supplementary Fig. 4d, e). This classification was robust, as shown by the chromatin accessibility of the Myh genes present only in myonuclei ( Fig. 3b and Supplementary Fig. 5a, b). Using chromVar we characterized the transcription factors (TF) variations in motif accessibility between the different types of myonuclei (Fig. 3c). Knowing that this analysis frequently does not distinguish between related TFs of the same family (usually sharing similar motifs) 27 , we could characterize that the NFAT and SOX binding motifs accessibility are enriched in slow Myh7 myonuclei, ESRRβ and NR4A2 motifs in fast Myh2 and Myh1 myonuclei and SIX and NRL/MAF motifs in fast Myh4 myonuclei (Fig. 3c, d and Supplementary Fig. 3b). These distinct combinations of TF could participate in the genetic code controlling the expression of specific fast-subtypes and slow sarcomeric genes and of genes associated with myofiber specialization (Fig. 1d). Equivalently to Fig. 1d presenting genes differently expressed in slow and fast myonuclei (from snRNA-seq data), we observed specific chromatin accessibility between slow and fast myonuclei in genes differently expressed in slow and fast myonuclei. For exemple, the gene Pvalb expressed only in fast Myh4 myonuclei (Fig. 1d) presented a strong chromatin accessibility only in fast Myh4 myonuclei ( Supplementary Fig. 5d). At the opposite a strong chromatin accessibility in the Atp2a2 was observed in Myh7 nuclei only where this gene is expressed ( Fig. 1d and Supplementary Fig. 5e). These results showed with a high-resolution how slow and fast-subtypes genetic programs are associated with different sets of TF controlling the expression and the chromatin accessibility of a set of co-expressed genes. Interestingly, a specific set of TF motifs is enriched in MTJ myonuclei as compared to body myonuclei, suggesting their participation in regionalized expression of specific genes such as Col22a1 ( Fig. 3e and Supplementary Fig. 5c).
Most myonuclei express a single fMyh gene. Fast Myh genes are located side by side on the same locus and it is presently unknown whether they can be transcribed simultaneously (Fig. 4a). Using snRNA-seq, we found that only rare mixed nuclei co-expressed simultaneously Myh7 and Myh2 or Myh1 and Myh4 (Fig. 4b). Mixed nuclei co-expressing Myh2 and Myh1 were more numerous. Other combinations of coexpression were not detected ( Supplementary Fig. 6a). This showed that the vast majority of adult myonuclei analyzed were pure and expressed only one Myh gene. To visualize the expression of fMyh genes, we performed fluorescent in situ hybridization (RNAscope) experiments using intronic probes to detect fMyh pre-mRNA in myofibers from EDL and soleus muscles ( Fig. 4c-e). Myonuclei expressed fMyh genes in a mono or bi-allelic manner (Fig. 4c). Expression of two distinct fMyh genes from identical or distinct alleles was detected in some nuclei (Fig. 4d). Expression of two fMyh genes by the same allele (single bicolor dot) indicated that two Myh genes were transcribed simultaneously on the same locus ( Fig. 4e-g). In some myonuclei, no fMyh pre-mRNA was detected. Quantification of RNAscope experiments is provided in (Fig. 4f, g): in EDL, the majority (>80%) of nuclei expressed only one Myh isoform and <5% co-expressed two isoforms. The number of nuclei co-expressing two isoforms of Myh was highest in the soleus, where~25% co-expressed Myh1 and Myh2 (Fig. 4f). The co-expression was in most cases from the same allele (Fig. 4g), suggesting that Myh2 and Myh1 can be activated simultaneously from one allele. The results of RNAscope experiments thus confirmed snRNA-seq results, showing that almost all myonuclei were pure, expressing only one Myh isoform at a time, except in the soleus. Then we tested whether different populations of myonuclei express either metabolic, calcium handling or sarcomeric genes, or if all myonuclei co-express simultaneously these different categories of genes. SnRNA-seq data showed that most Myh4 myonuclei expressed the calcium handling gene Pvalb and the glycolytic gene Pfkfb3, that a majority of fMyh myonuclei expressed Tbc1d1 involved in muscle glucose uptake, and most Myh7 nuclei expressed the sarcoplasmic reticulum Calcium transporting 2 gene Atp2a2 ( Fig. 5a-d). These results were confirmed by RNAscope visualizing the localization of the pre-mRNA of AldolaseA, associated with glycolytic metabolism in myofiber 28 , and of Idh2, associated with oxidative metabolism 29 . Accordingly, more than 70% of Myh4 expressing nuclei also expressed the glycolytic enzyme AldolaseA pre-mRNA, while the TCA cycle isocitrate dehydrogenase Idh2 pre-mRNA was undetectable. At the opposite most Myh2 expressing nuclei expressed Idh2 pre-mRNA ( Fig. 5e-h), while AldolaseA pre-mRNA was not detected. Thus, myonuclei co-expressed coordinated specialized sarcomeric, calcium handling and metabolic genes. Altogether these results showed that fMyh gene expression, and more generally muscle gene expression within the myofiber syncytium is not random but precisely controlled, each myonucleus showing a precise pattern of co-expressed genes.
Transcription is coordinated in myonuclei along the myofiber. We next evaluated the coordination of fMyh genes expression in distinct myonuclei of the same fiber. In EDL, coordinated Myh expression was detected in all myonuclei in a majority of fibers ( Fig. 6a, upper panel). In contrast, a large number of uncoordinated fibers was detected in soleus ( Fig. 6b). Different nuclei within the same soleus fiber were detected expressing different fMyh isoforms, while some nuclei co-expressed two isoforms. Quantification of these experiments is provided in (Fig. 6c, d). We never observed hybrid coordinated myofibers: meaning that we did not observe that all nuclei of a given fiber express the same combination of 2 or 3 Myh genes. In quadriceps and EDL, most myofibers (90%) showed Myh gene expression coordination between their myonuclei expressing only Myh4 (Fig. 6c, d). Myh4 mRNA accumulated homogenously in the coordinated fibers expressing Myh4 pre-mRNA ( Supplementary Fig. 7a). Interestingly, half of the fibers from the soleus were hybrid (Fig. 6c, d). In hybrid Myh4 and Myh1 fibers, Myh4 mRNA accumulated less efficiently, and did not spread far from the nuclei expressing Myh4 pre-mRNA ( Supplementary Fig. 7a), demonstrating that mRNA diffusion was confined within myonuclear domains 8 . The transcriptional program of Myh4-positive myonuclei may differ in hybrid myofibers from pure myofibers, explaining the heterogeneity of Myh4 myonuclei observed by snRNA-seq (Supplementary Fig. 3b, c). Myh4 pre-mRNA was detected also in NMJ myonuclei, as identified by AchE mRNA expression (Supplementary Fig. 7b).
Innervation controls myonuclei coordination along the myofiber. To understand the basis behind pure and mixed myonuclei and their coordination along the myofiber, we tested whether innervation, known to modulate adult muscle fiber type 14,15,30 , was involved in the coordination process. At the protein level, 1 week after denervation, we observed a preferential atrophy of fast MYH4 and MYH1-positive fibers and 4 weeks after denervation a fast to slow transition in all muscles ( Supplementary Fig. 8a). In agreement, 1 week after denervation, the number of nuclei negative for fast Myh genes was increased (Figs. 4f and 6a). The percentage of uncoordinated myofibers after 1-week denervation in the EDL was increased more than five-fold but did not significantly differ in soleus (Fig. 6d). These observations showed that innervation is absolutely required to coordinate fast Myh genes expression in the nuclei within myofibers in EDL (  [32][33][34] . Knowing that in adult myofibers the homeoprotein SIX1 controls the expression of many fast/glycolytic muscle genes among which fMyh genes 35,36 and that SIX binding motifs accessibility are enriched in Myh4 myonuclei (ref. 26 and Fig. 3c), we looked at the nuclear accumulation of this protein in fast myonuclei along EDL myofibers. As shown in Fig. 6e we observed robust SIX1 accumulation in all myonuclei from an extremity to the other of the fiber. To further test the involvement of SIX1 in fMyh gene coordination along the myofiber we tested the consequences of its absence in mutant myofibers from muscle-specific Six1 conditional knock-out (Fig. 6f, g). In absence of Six1, EDL myofibers appeared slower and presented an increased number of non-coordinated hybrid fibers (Fig. 6f, g), showing that SIX transcription complexes participate to coordinate the expression of fMyh genes in adult myonuclei.
Coordination of myonuclei during development. We next addressed how fMyh expression coordination appears during development. A sequential transition from embryonic to neonatal and to adult myosin heavy chain isoforms has already been characterized in skeletal muscles during mouse post-natal development 37 . We further analyzed the expression pattern of mRNAs or pre-mRNAs of all sarcomeric Myh (embryonic, neonatal, slow and adult fast) on isolated forelimbs myofibers at different stages. At E15.5, when myofibers were already innervated, embryonic primary myofibers expressed Myh3 and Myh8 mRNA ( Supplementary Fig. 9a). At E18.5, secondary fetal myofibers (small fibers) displayed a higher accumulation of Myh3 than of Myh8. Myh3 mRNAs accumulated strongly in MTJ regions of myofibers at all stages, (Supplementary Fig. 9a, b) and after birth also at the center of myofibers, where NMJs are suspected to be localized ( Fig. 7a and Supplementary  Fig. 9a). These results suggested the specific transient activation of Myh3 in NMJ and MTJ myonuclei, or its transient activation in newly accreted myonuclei in these regions of growing myofibers as proposed previously 38 . In contrast, Myh8 mRNA was not regionalized along the myofiber at any stage. Adult fMyh pre-mRNAs were first detected in a few nuclei at E15.5 and their number increased drastically after E18.5 (Fig. 7b, c). At these different stages, most of the myonuclei expressed only one adult fMyh gene (Fig. 7d) without obvious regionalization along the myofiber. We observed mono-allelic expression before E18.5, and bi-allelic expression after P2 (higher transcription rates) (Fig. 7e). At P2 and P5, more than 90% of the nuclei were coordinated along the myofibers and expressed one fMyh gene (Fig. 7f). This percentage was similar to that of adult EDL, suggesting that fMyh gene coordination is established early when the expression is initiated. Myh2 and Myh4 mRNAs accumulated strongly from P5 but some areas of the fiber did not present mRNAs accumulation (Supplementary Fig. 10a). Slow Myh7 mRNAs accumulated in all primary fibers at the MTJ level at E15.5 (Fig. 7g). After E18.5, the expression of Myh7 was restricted to slow fibers, which homogeneously accumulated this mRNA specifically (Supplementary Fig. 10b). These experiments showed that the majority of myonuclei started to express adult fMyh genes at P2, a single gene being activated by nucleus, and that myonuclei in the growing myofiber were coordinated from their onset.

Discussion
By combining single-nucleus RNA-seq/ATAC-seq and FISH experiments, we established an atlas of the gene expression pattern and chromatin landscape of the different cell populations of adult skeletal muscles, including mononucleated cells and multinucleated myofibers. From our experiments, we estimated that myonuclei constitute~68% of the nuclei present in adult skeletal muscles of 2 months old mice at rest. Recent single-cell RNA-seq studies of skeletal muscle 39,40 have characterized the gene expression pattern of mononucleated cells of adult skeletal muscle, but failed to identify the transcription signature of multinucleated myofibers. The composition and gene expression of cell types identified by these reports is in agreement with our data, providing the same global signature of skeletal muscle at rest, identifying distinct populations of FAPs, tenocytes, immune cells, pericytes, adipocytes, smooth muscle cells, and myogenic stem cells. Expanding these studies, our results showed that myofibers are composed of three different types of myonuclei localized differently and expressing distinct sets of genes (Fig. 2), namely the MTJ, NMJ, and body myonuclei. The full transcriptome of these nuclei was to this day unknown, and we characterized numerous genes expressed specifically in these specialized myonuclei. We also characterized an unidentified population of Titin + myonuclei expressing Runx1, Dysf, Ehd4, Flnc, and Myh9. As Filamin C is involved in early stages of myofibrillar remodeling and repair, this myonuclei population, which had never been described could correspond to myonuclei in area of myofiber damage to support fast repair 41,42 . Our results demonstrated that most myonuclei express a single fMyh gene, and that in the majority of myofibers myonuclei coordinately express a specific set of genes, including this Myh gene (Fig. 8). Our RNAscope experiments also identified a feature of transcriptional control in myonuclei of adult myofibers. In fast skeletal muscles, the transcriptional activity of nuclei within each myofiber is finely coordinated. We showed that this coordination is dependent upon innervation and is established early during development.
Several reports have shown that apparently homogenous populations of cells present an important variability of gene expression. This variability can be the consequence of stochastic gene expression, different stages of cell cycle, or distinct subcategories of cells in the studied population 43 . In syncytia like osteoclasts and syncytiotrophoblasts, gene expression can be heterogeneous in nuclei that share the same cytoplasm 10,11 . We found that the transcription of adult fMyh genes in myonuclei along given myofibers was finely controlled during development, allowing the coordinated expression in their nuclei of a single fMyh gene. This transcriptional coordination could be a consequence of various processes. In Drosophila, one founder myonucleus was shown to reprogram newly fused myonuclei, inducing coordination of expression between nuclei of the same fiber 44 . Alternatively, distinct populations of preprogrammed  Fig. 1a. Each dot corresponds to one myonucleus and the axis correspond to the gene expression level. Atp2a2 is expressed only in most Myh7 and in few Myh2 nuclei. b Same as a with Tbc1d1 expression. Tbc1d1 is expressed only in fast and not in slow nuclei. c Same as a with Pvalb expression. Pvalb is expressed only in Myh1 and Myh4 nuclei. d Same as a with Pfkfb3 expression. Pfkfb3 is expressed preferentialy in Myh4 nuclei. e RNAscope on isolated fibers from EDL showing the localization of Myh4 (red) and AldoA (green) pre-mRNAs. Most myonuclei expressed at the same time Myh4 and AldoA (white arrowhead). However, some nuclei expressed Myh4 without AldoA (red arrowhead), and some nuclei are negative for both genes. f Percentage of the nuclei co-expressing Myh4 and AldoA or expressing only Myh4 or AldoA and negative nuclei in EDL fibers. g RNAscope on isolated fibers from soleus showing the localization of Myh2 (red) and Idh2 (green) pre-mRNAs. Like in e most myonuclei expressed at the same time Myh2 and Idh2 (white arrowhead). However, some nuclei expressed Myh2 without Idh2 (red arrowhead) and other nuclei are negative for both genes. h Percentage of the nuclei co-expressing Myh2 and Idh2 or expressing only Myh2 or Idh2 and negative nuclei in soleus fibers. For e.g., scale bar: 100 μm. For f, h the graphs show data pooled from three animals and fifteen fibers per animal. Numerical data are presented as mean ± s.e.m. *P < 0.05, **P < 0.01, ***P < 0.001. Source data for f, h are provided in the Source Data file.
NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18789-8 ARTICLE myogenic stem cells with distinct potential [45][46][47][48][49] could fuse (homotypic fusion 50 ), activating a single Myh gene along the myofiber. Moreover, extracellular signals from surrounding cells may also intervene 51 . The de-coordination of myonuclei after denervation observed in our study supports the importance of innervation in the coordination process during development. A better understanding of the matching between motoneuron subtypes and myofiber subtypes that operates during development 52 may help understand how neuromuscular specialization contributes to the onset of coordination.
Adult myofibers are plastic and can adapt their contractile and metabolic properties depending on external stimuli by activating These changes occur in a sequential order recapitulating the organization of the genes at the fast Myh locus (Myh4-Myh1-Myh2) 15 . In this report, we showed that most myonuclei expressed only one Myh isoform at a given time (Fig. 4). This robust and exclusive expression is reminiscent of the control of the ß-Globin genes by a locus control region 23 . The regulation of the fMyh locus is not fully understood and could be controlled like the ß-Globin genes by a locus control region, such as a super enhancer 53,54 , activating the expression of associated embryonic, fetal and adult Globin genes in a temporal order. Transient transfection of adult mouse TA and soleus with GFP reporters under the control of 1 kb upstream of each Myh2, Myh1 and Myh4 promoters recapitulate their specific expression in the corresponding MYH2, MYH1 and MYH4 myofibers, respectively 55,56 . This suggests that each fMyh promoter possesses its own fiber specific cis-regulatory elements. Accordingly, in Myh4 + nuclei, we observed in snATAC-seq experiments a higher chromatin opening in Myh4 than in Myh1 and Myh2 genes. In Myh4 + nuclei, we failed to detect chromatin opening at the Myh1 and Myh2 promoter regions, suggesting the absence of TF binding. In soleus myofibers with a high proportion of Myh1 and Myh2 co-expression (where Myh1 and Myh2 pre-mRNAs are detected in many myonuclei), we mainly observed Myh1 and Myh2 pre-mRNAs as two dots in the myonuclei, suggesting that each allele transcribed a single gene, and that these two adjacent fMyh genes could not be transcribed simultaneously. Whether this exclusive expression reflects the competition for a shared enhancer at the locus, as already illustrated for the ß-Globin genes 23 , is an interesting possibility. The coordinated expression of a single isoform of Myh could be due to a similar combination of transcription activators and repressors concentration in these myonuclei, affecting all the genes that they express, or only those with high expression, excluding stochastic expression as previously proposed 9 . In support of this, snRNA-seq experiments showed that specific genes are co-expressed with Myh4, such as genes coding for other sarcomeric proteins, proteins associated with metabolism or associated with calcium handling. The snRNA-seq results showed harmonious non-random gene expression pattern in each nucleus of the myofiber, with some variability between muscles with distinct subprograms among Myh4 + myonuclei. In the mouse, denervation induces muscle atrophy, associated with a fast to slow transition of Myh at the protein level 32 . These changes occur also at the pre-mRNA level, as shown in the present study. We observed that transitions after denervation were sequential, Myh4 + fibers expressing Myh1 but not Myh2, and some nuclei switching their Myh expression profile earlier than others. Thirty days after denervation, myofibers remained uncoordinated showing that the increased number of hybrid fibers observed in EDL 7 days after denervation was not due to a transitional state but that innervation was required to maintain the coordination of fMyh expression. The atrophy of myofibers induced by denervation involves the activation of protein degradation pathways and the repression of protein synthesis 57 . We observed an increased number of fMyh-negative nuclei as well as a decreased bi-allelic expression of fMyh after muscle denervation, suggesting that downregulation of fMyh expression may participate as well in this atrophy.
The loss of coordination in myonuclei of adult myofibers after denervation observed in the present study is consistent with the role of motoneurons in controlling adult muscle physiology 14,15,30,32,58,59 . Calcium is a diffusible signal distributed throughout myofibers and a potential candidate for the control of specific muscle gene coordination. Resting Ca 2+ concentration differs in fast and slow fibers and is modulated by slow or fast motoneuron firing 60-62 . Ca 2+ transients following slow motoneuron stimulation activate specific signaling pathways 14,63 among which Calcineurin phosphatase that controls the nucleo-cytoplasmic localization of NFATc1 in slow myofibers important for the expression of slow muscle genes 64,65 . MYOD activity as well is modulated by motoneuron activity and phosphorylation of its T115 by slow motoneuron activation impairs its DNA binding activity in the soleus 66 . In contrast, fast-twitch motoneurons, characterized by phasic highfrequency firing do not activate Calcineurin phosphatase, suggesting that other signaling pathways relay motoneurons firing information in fast myofibers. Furthermore, fast-subtypes myofibers are innervated by fast-subtypes motoneurons, including fast fatigable (FF, innervating MYH4 myofibers), fast intermediate (FI, innervating MYH1 myofibers) and fast resistant (FR, innervating MYH2 myofibers) [67][68][69] . How these distinct classes of fast motoneurons control the fast-muscle subtypes phenotype and myonuclei coordination is presently unknown 32,59,67,70 . Whether Ca 2+ concentrations are similar along the myofiber and how denervation may perturb this basal intracellular Ca 2+ concentrations are potential important issues to understand myonuclei coordination. Several proteins controlling intracellular Ca 2+ concentration in myofibers (Calsequestrin 1 et 2, Serca1 and 2, Parvalbumin, Calmodulin, Phospholamban, Myozenin 1 and 2, Sarcalumenin, IP3R) are differentially expressed in fast and slow myofibers, and induce a variety of local Ca 2+ concentration gradients that may modulate the activity of different transcription factors [71][72][73][74] . Here, we showed that the transcription factor SIX1 accumulates in all myonuclei of EDL myofibers and participates in the coordination of fMyh expression. The mechanisms controlling SIX1 nuclear accumulation specifically in fast myofibers remain to be established. We observed by snRNA-seq an increased expression of the SIX cofactor EYA4 in Myh4 + myonuclei. Whether the SIX1-EYA4 transcription complex activity may be controlled by motoneuron firing frequency and participates in Fig. 6 Coordination of fast Myh expression along myofibers. a RNAscope on isolated fibers from EDL showing the localization of Myh4 (red), Myh1 (green), and Myh2 (yellow) pre-mRNAs. Up: EDL myofiber showing coordination of Myh4 expression between nuclei (pure myofiber). Down: EDL myofiber after 1 week of denervation. After denervation, myonuclei are no more coordinated in the myofiber and express either Myh4 or Myh1 (green arrowhead) or co-express both isoforms (hybrid myofiber). In most myonuclei, fMyh pre-mRNAs are no more detected (white arrowhead). b Same as a with Myh4 (yellow), Myh1 (green), and Myh2 (red) pre-mRNAs in a soleus hybrid myofiber. The white arrowhead shows a mixed Myh1, Myh2 myonucleus, whereas the red arrowhead shows a pure Myh2 myonucleus. c Percentage of the different types of fibers observed in quadriceps, EDL and soleus, and 1 week after denervation. d Percentage of hybrid (de-coordinated) fibers in control EDL and soleus and after 1 week of denervation. e SIX1 accumulates in all myonuclei of adult EDL fibers as identified by immunohistochemistry with SIX1-specific antibodies. f Percentage of the different types of EDL fibers observed in WT and Six1 flox/flox ;HSA-CRE mice by RNAscope experiments. g Percentage of hybrid EDL fibers (de-coordinated) in WT and Six1 flox/flox ;HSA-CRE mice. For a, b (up) scale bar: 100 μm. For b (down) scale bar: 20 μm. For c, d, n = 3 and 20 fibers per animal. For f, g the graphs show data pooled from three animals and more than fifty fibers in total. Numerical data are presented as mean ± s.e.m. *P < 0.05, **P < 0.01, ***P < 0.001. Source data for c, d, f, g are provided in the Source Data file.
myonuclei coordination is an interesting possibility that may explain its role in the control of the adult fast phenotype 75 .
The genome-wide chromatin landscape of adult muscles had not been evaluated previously at the myofiber level, due to the lack of adapted experimental protocols 26 . We show here that snATAC-seq is a powerful tool to overcome this limitation and that it allows to detect all the cell populations identified by snRNA-seq in adult skeletal muscles, including myofibers. We identified chromatin accessibility domains coherent with the defined patterns of gene expression, specific of fast-subtypes or slow myonuclei, of MTJ or body myonuclei, as well as of other cell types in skeletal muscle. This analysis also allowed to identify enriched TF-binding sites in ATAC-peaks specific for each myonuclei population. We identified TF-binding motifs known to participate in muscle genes regulation in slow and fast myofibers, among which the NFATc 65 , MEF2 76  nuclear receptors of the ERR-ß/γ 79 and the NR4A1 80 families. The differential TF ATAC-peaks enrichment observed between fast-subtypes myonuclei sheds some light on the basis of their specialization. We identified TF DNA binding motifs, which had never been implicated in muscle genes regulation, such as those recognized by the b-ZIP family of NRL/MAF transcription factors 81 . Accordingly, we identified c-Maf mRNAs enriched in the Myh4 + myonuclei in our snRNA-seq, suggesting that this TF may contribute to the fine-tuning of muscle gene expression in specialized fast myonuclei, in synergy with SIX homeoproteins. Serria et al showed earlier that c-Maf expression is upregulated in myogenic cells during their differentiation in culture, but the role of c-Maf in adult muscle physiology is presently unknown 82 .
Last, in addition to defining the exact transcriptional landscape within skeletal muscle, the snRNA-seq and snATAC-seq methods presented here open new perspectives to characterize in the future the crosstalk between myofibers and their environment, such as with associated myogenic stem cells in pathophysiological conditions, to date poorly investigated due to technical limits.

Methods
Animals. Animal experimentations were carried out in strict accordance with the European STE 123 and the French national charter on the Ethics of Animal Experimentation. Protocols were approved by the Ethical Committee of Animal Experiments of the Institut Cochin, CNRS UMR 8104, INSERM U1016, and by the Ministere de l′education nationale, de l′enseignement et de la recherche, no. APAFIS#15699-2018021516569195. We used 6-to 8-week-old C57black6N mouse female for most of our experiments. Mice were anesthetized with intraperitoneal injections of ketamine and xylazine and with subcutaneous buprecare injections before denervation was performed by sectioning of the sciatic nerve in one leg.
FACS sorting of nuclei from adult skeletal muscle for single-nucleus RNA-seq. For the mix condition 6 tibialis, 6 gastrocnemius, 6 soleus, 6 plantaris, and 6 EDL were pulled together. For the soleus condition, 24 soleus, for the tibialis condition 6 tibialis and for the quadriceps condition 6 quadriceps were pulled together. Muscles were dissected and pulled in cold phosphate-buffered saline (PBS) with 0.2 U/μl RNase inhibitor. PBS was removed and muscles were minced in 1 ml cold lysis buffer (10 mM Tris-HCl pH7.5, 10 mM NaCl, 3 mM MgCl 2 , and 0.1% NonidetTM P40 in Nuclease-Free Water) with scissors. After 2 min, 4 ml of cold lysis buffer were added and muscles were lysed for 3 min at +4°C. Nine milliliter of cold wash buffer (PBS, BSA 2% and 0.2 U/μl RNase inhibitor from Roche) was added and the lysate was dounced with 10 strokes of loose pestle avoiding too much pressure and air bubbles. The homogenate was filtered with 70 and 40 μm cell strainers. Nuclei were pelleted by centrifugation for 5 min at 500 × g at +4°C. Nuclei were washed to remove ambient RNAs in 1 ml of cold wash buffer, transferred in a 1.5 ml Eppendorf tube and centrifuged 5 min at 500 × g at +4°C. The pellet was resuspended in 250 μl of wash buffer and stained during 15 min in the dark at +4°C with DAPI (10 μg/ml). Then the nuclei were washed with 1 ml of wash buffer, centrifuged 5 min at 500 × g at +4°C, resuspended in 300 μl of wash buffer and filtered with 30 μm cell strainers. Nuclei were then FACS sorted to exclude debris with a BD FACSAria III and the BD FACSDIVA software. More than 200,000 nuclei were collected in a tube containing 200 μl of wash buffer.
Single-nucleus RNA-seq from skeletal muscle. After FACS sorting, nuclei were centrifuged 5 min at 500 × g at +4°C, resuspended in 15 μl of wash buffer and counted with a hemocytometer. The concentration of nuclei was then adjusted to 1000 nuclei/μl with wash buffer. We loaded around 4000 nuclei per condition into the 10x Chromium Chip. We used the Single-Cell 3′ Reagent Kit v3 according to the manufacturer's protocol. GEM-Reverse Transcription (GEM-RT) was performed in a thermal cycler: 53°C for 45 min, 85°C for 5 min. Post GEM-RT Cleanup using DynaBeads MyOne Silane Beads was followed by cDNA amplification (98°C for 3 min, cycled 12 × 98°C for 15 s, 67°C for 20 s, 72°C for 1 min). After a cleanup with SPRIselect Reagent Kit and fragment size estimation with High-SensitivityTM HS DNA kit runned on 2100 Bioanalyzer (Agilent), the libraries were constructed by performing the following steps: fragmentation, endrepair, A-tailing, SPRIselect cleanup, adaptor ligation, SPRIselect cleanup, sample index PCR, and SPRIselect size selection. The fragment size estimation of the resulting libraries was assessed with High-SensitivityTM HS DNA kit runned on 2100 Bioanalyzer (Agilent) and quantified using the QubitTM double stranded DNA (dsDNA) High-Sensitivity HS assay (ThermoFisher Scientific). Libraries were then sequenced by pair with a HighOutput flowcel using an Illumina Nextseq 500 a b  with the following mode: 26base-pairs (bp) (10x Index + UMI), 8 bp (i7 Index), and 57 bp (Read 2).
Single-nucleus RNA-seq analysis. A minimum of 50,000 reads per nucleus were sequenced and analyzed with Cell Ranger Single-Cell Software Suite 3.0.2 by 10x Genomics. Raw base call files from the Nextseq 500 were demultiplexed with the cellranger mkfastq pipeline into library-specific FASTQ files. The FASTQ files for each library were then processed independently with the cellranger count pipeline with parameter set to default. The single-nucleus RNA-seq assay captures unspliced pre-mRNA, as well as mature mRNA. To include these intronic reads in our analysis, we created a custom "pre-mRNA" reference package, listing each gene transcript locus as an exon with cellranger mkgtf and cellranger mkref according to 10x genomics instruction (https://support.10xgenomics.com/single-cell-geneexpression/software/pipelines/latest/advanced/references). We used a reference genome built against mouse mm10, Sequence: GRCm38 Ensembl 93. We used STAR21 to align cDNA reads against our reference genome. Once aligned, barcodes associated with these reads-cell identifiers and Unique Molecular Identifiers (UMIs), underwent filtering and correction. Reads associated with retained barcodes were quantified and used to build a transcript count table. The subsequent visualizations, clustering and differential expression tests were performed in R (v 3.4.3) using Seurat (v3.0.2) 25 . Quality control on aligned and counted reads was performed keeping cells with >200 and <2500 nFeature RNA and <5% mitochondrial genes. We get 6477 nuclei from the mix of tibialis, EDL, gastrocnemius, plantaris and soleus, 1335 nuclei from quadriceps, 4838 nuclei from tibialis, 2517 nuclei from soleus. We detected the expression of approximately one thousand genes per nucleus in each of these muscles ( Supplementary Fig. 1).
Single-nucleus ATAC-seq from skeletal muscle. We used the 10x genomic nuclei Isolation for Single-Cell ATAC Sequencing protocol (CG000169|Rev B) with the following modifications: 12 quadriceps and 12 soleus were dissected and pulled in cold PBS. PBS was removed and muscles were minced 2 min in 1 ml of cold ATAC-lysis buffer (10 mM Tris-HCl pH7.4, 10 mM NaCl, 3 mM MgCl 2 , 1% BSA and 0.1% Tween-20 in Nuclease-Free Water). Six milliliter of cold ATAC-lysis buffer were added and muscles were lysed on ice. After 3 min the lysate was dounced with 10 strokes of loose pestle avoiding too much pressure and air bubbles. After douncing, 8 ml of wash buffer were added and the homogenate was filtered with 70, 40, and 20 μm cell strainers. Nuclei were pelleted by centrifugation for 5 min at 500 × g at +4°C. Next, we used the Chromium Single-Cell ATAC kit according to the manufacturer's protocol. Nuclei were resuspended in nuclei buffer from the kit, transposed 1 h at 37°C. We loaded around 6000 nuclei into the 10X Chromium Chip. GEM incubation and amplification were performed in a thermal cycler: 72°C for 5 min, 98°C for 30 s, and 12 repeated cycles of 98°C for 10 s, 59°C for 30 s, and 72°C for 1 min. Post GEM Cleanup using DynaBeads MyOne Silane Beads was followed by library construction (98°C for 45 s, cycled 12 × 98°C for 20 s, 67°C for 30 s, 72°C for 1 min). The libraries were constructed by adding sample index PCR, and SPRIselect size selection. The fragment size estimation of the resulting libraries was assessed with High-SensitivityTM HS DNA kit runned on 2100 Bioanalyzer (Agilent) and quantified using the QubitTM dsDNA High-Sensitivity HS assay (ThermoFisher Scientific). Libraries were then sequenced by pair with a HighOutput flowcel using an Illumina Nextseq 500.
FISH with amplification (RNAscope) on isolated fibers. RNAscope® Multiplex Fluorescent Assay V2 was used to visualize fMyh pre-mRNAs and mRNAs, and AldolaseA and Idh2 pre-mRNAs. Twenty different pairs of probes against the first intron of each fMyh, AldolaseA and Idh2 pre-mRNAs were designed by ACDbio. Muscles were dissected and immediately fixed in 4% PFA at +4°C for 30 min.
After fixation muscles were washed three times in PBS for 5 min. Myofibers were dissociated mechanically with small tweezers and fixed onto Superfrost plus slides (Thermo Fischer) coated with Cell-Tak (Corning) by dehydration at +55°C during 5 min. Slides were then proceeded according to the manufacturer's protocol: ethanol dehydration, 10 min of H 2 O 2 treatment and 30 min of protease IV treatment. After hybridization and revelation, the fibers were mounted under a glass coverslip with Prolong Gold Antifade Mountant (Thermofischer). Myofibers were imaged with a Leica DMI6000 confocal microscope composed by an Okogawa CSU-X1M1 spinning disk and a CoolSnap HQ2 Photometrics camera. Images were analyzed with Fiji Cell counter program.
Immunohistochemistry. For immunostaining against MYH4, MYH2, MYH7, and Laminin, adult legs without fixation and without skin were embedded with Tis-suTEK OCT, and directly frozen in cold isopentane cooled in liquid nitrogen.
Muscles were conserved at −80°C and cut with Leica cryostat 3050s with a thickness of 10 μm. Cryostat sections were washed three times 5 min with PBS and then incubated with blocking solution (PBS and 10% goat serum) 30 min at room temperature. Sections were incubated overnight with primary antibody solution at +4°C, then washed three times for 5 min with PBS and incubated with secondary antibody solution 1 h at room temperature. Sections were further washed three times for 5 min and mounted with Mowiol solution and a glass coverslip. Images were collected with an Olympus BX63F microscope and a Hamamatsu ORCA-Flash 4.0 camera. Images were analyzed with ImageJ program. The references of the antibodies used are listed in Supplementary Table 1.
Statistical analysis. The graphs represent mean values ± SEM. Significant differences between mean values were evaluated using two-way ANOVA with Graphpad 6 software and student t test for Fig. 6g.