Spatially resolved transcriptomics reveals innervation-responsive functional clusters in skeletal 1 muscle.

studies of the mechanisms underlying skeletal muscle homeostasis and innervation.


Graphical Abstract
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022.

36
Skeletal muscle is a complex organ in which several mononucleated cell populations organize within 37 specific anatomic compartments. These domains -nerves, connective, and vasculature -act altogether 38 in coordination with muscle fibers to ensure the seamless execution of each contraction. In response 39 to acute injury specific resident muscle populations, activate and contribute to regeneration by 40 regulating muscle stem cell activity or cooperating with immune cells recruited from the circulatory

55
In particular, most of skeletal muscle mass is composed of long, polynucleated fibers, which, due to 56 the necessary dissociation step, are incompatible with the analysis of cytoplasmic RNAs. Moreover, 57 functional anatomical domains such as NMJs, tendon-to-muscle connections, vessels and connective 58 tissue can have different interactions depending on their relative proximity and locations; cryosection were then placed on Visium Slide capture areas and stained. Each area contains ~5000 86 RNA capture spots characterized by a specific spatial barcode that ensures that transcripts can be 87 mapped to their original histological location. To highlight the different anatomical features of the 88 muscle, we used anti-laminin and anti-collagen-1 antibodies, bungarotoxin (BTX) and Hoechst. This 89 strategy enabled us to visualize fiber boundaries, extracellular matrix scaffolds, NMJs and cellular 90 nuclei, respectively ( Figure 1B and Suppl. Figure 1A and C). After imaging, tissue was permeabilized 91 and polyadenylated RNA captured on the underlying spots. Following manufacturer protocols, we 92 performed reverse transcription, second strand synthesis and amplification. After library preparation 93 and sequencing, reads were processed and spatially resolved using Space Ranger software. The 94 different datasets were then analyzed separately and integrated using Seurat 4 (Hao et al., 2021a).

95
Lastly, data were overlaid onto the acquired image of the tissue sections ( Figure 1B and Suppl. Figure   96 1A and C). Overall, we identified 14,417 features across 2,987 unique capture spots. Unsupervised 97 clustering identified 8 clusters with consistent transcript abundance distribution (Supp. Figure 1E).

98
Tissue projection showed clear spatial segregation for four out of the 8 clusters (cluster 0, 1, 3, 4) 99 ( Figure 1C -Supp. Figure 1B and D). This specific spatial arrangement readily suggested that these 100 clusters could represent the major metabolic domains and general areas of the muscle such as 101 different type of fibers or specific components of the matrix scaffold. Conversely, clusters 2, 5, 6, and 102 7 showed a sparse distribution more in line with small functional elements disperse within the tissue 103 such as NMJ, vessel or nerves. Based on this initial observation we therefore set out to determine the 104 identity of each of these specific cluster based on their gene expression profile (Supplementary Table   105 1) and their overlap with pre annotated anatomical structures in the uninjured muscle.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. Cluster0, was predicted to contain type 2B fibers (fast glycolytic), the portion overlapping with 137 Cluster1 consisted of type 2A fibers (fast oxidative), and type 2X fibers, which are less glycolytic than 138 type 2B fibers ( Figure 2C).

140
Next, we analyzed the molecular signature of the other two clusters that presented a clear spatial 141 distribution, Cluster3 and Cluster4. Upon GO analysis, both Clusters 3 and 4 displayed an enrichment 142 in ECM components (Suppl. Figure 3 Table A). However, while Cluster 3 overlapped considerably with 143 the intramuscular matrix scaffold, Cluster 4 clearly identified the epimysium surrounding the muscle 144 ( Figure 3A and Suppl. Figure 3B). Almost perfect overlap was obtained between the localization of 145 Cluster4 and the localization of collagen 1 surrounding the muscle (validated by Immunofluorescence 146 and Sirius red staining - Figure 3A -inset). As expected, the genes enriched in this cluster are 147 extracellular matrix components such as collagen subunits (Col12a1 and Col1a1) and matrix-148 associated molecules (Fmod and Timp2) ( Figure 3B).

150
Given the complexity of muscle tissue, we wondered whether ST analysis might not only reveal more 151 abundant clusters related to muscle fibers but also identify more restricted morphofunctional zones 152 specific to this tissue. We therefore analyzed the gene expression signature of the sparse clusters 153 (Cluster2 , Cluster 5, Cluster6, Cluster7) and their relative spatial proximity to known morphological 154 structures. Upon visual inspection, Cluster 2 did not present a recognizable distribution pattern nor 155 overlap with clearly detectable anatomical compartments. Furthermore, we could not detect a strong 156 gene expression signature for the cluster. However, it must be considered that each spot contains (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint 2007). Complete denervation of TA and EDL is observed three days after injury (Suppl. Figure 4A Table 2). We then analyzed the specific gene 193 expression profile of each cluster at 3 days and 30 days after denervation ( Figure 5A and Suppl. Table   194 3). Similarly, we identified the upregulation of genes related to atrophy in specific clusters such as the 195 growth arrest and DNA damage-inducible 45α (Gadd45a), the ubiquitin ligase Itch (Nerves) and the 196 homeobox gene Tgif (Blood Vessels) or autophagy related genes Beclin1 (Becn1) and Gabarapl1 in 197 both Nerve and Blood Vessel Clusters ( Figure 5B). Interestingly, we noticed that the induction Atrogin-

204
We confirmed this trend in our observations, via qPCR analysis of laser-microdissected glycolytic (Gly) 205 or oxidative (Ox) areas of the muscle sections (see scheme in Suppl. Figure 4E). MuRF1 and Atrogin-1 206 were more highly expressed in the glycolytic cortex of the denervated TA muscle than in the oxidative 207 core ( Figure 5E). Interestingly, although the reduction in fiber area was rather limited at 3 days of 208 denervation, we observed by separately measuring the fibers in the glycolytic cortex and those in the 209 oxidative core (strategy shown in Suppl. Figure 4E) that the reduction was significantly more

227
It is noteworthy that, concomitantly with a reduction in overall expression in the glycolytic cluster, we 228 observed a partial loss of the spatial restriction of these genes ( Figure 6B and D). Intriguingly, at 30 229 days after the induction of nerve damage-when innervation is restored-Amd1, Amd2 and Smox 230 expression in the muscle is still decreased compared to CTR levels ( Figure 6C). This could suggest that 231 more time or more extensive maturation of NMJs is needed to re-establish the correct expression of 232 these markers.

233
We validated this observation through laser microdissection; indeed, while the expression of Amd1,

234
Amd2 and Smox was high in the glycolytic portion of the TA in control tissue, after denervation, due 235 to their reduced expression in glycolytic fibers, the relative levels of these genes were similar between (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. Putrescin from Ornithine, is upregulated in denervated muscle (Suppl. Figure 5D).

248
Collectively, these data show how ST could be used to identify variation in the spatial patterning of (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

310
Here, we only focused on TA and EDL muscle groups alone which are almost exclusively composed of 311 fast type fibers. In the future it would be interesting to use TS to analyze atrophy in different muscle 312 groups with a more divers fiber composition to assess whether the same fiber types activate the same 313 transcriptomic program in different muscles and to highlight differences with slow twitch muscle 314 fibers (e.g. predominant in Soleus muscle). In addition, exploring the regeneration process at later 315 time after denervation would be useful to more accurately define the molecular events that occur   333 nonetheless these approaches often require scRNA-seq or bulk RNA-seq data as initial references.

334
Eventually, the advent of more dense spot arrays will be overcome this technical bottleneck. In 335 muscle, the limit in resolution is partially compensated by the large size of fibers that make less likely 336 the colocalization of different fibers within a single spot. Despite the current limitations we foresee 337 that the use of ST technology will be instrumental to understand the molecular dynamics occurring at (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint in terms of tissue localization, will contribute to elucidate the response of specific histological 340 functional units to external perturbations such as traumatic injury or chronic inflammatory conditions.

342
In summary, our study provides a detailed molecular profiling that will serve as a steppingstone to 343 investigate the mechanism underpinning muscle specific response to denervation at morpho-344 functional level -as in the case of differential fiber sensitivity to atrophy -.Furthermore, by providing 345 a spatially resolved transcriptomic reference of the unperturbed muscle will allow the identification 346 of new cellular or structural interactions in response to any perturbation of the homeostatic balance.

361
The authors declare no competing interests.

363
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

371
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.  (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

391
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.   (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

414
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

427
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.  (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

495
Total RNA was extracted from total muscle or dissected specimens (look above) using Qiagen RNeasy 496 mini kits following the manufacturer's protocol. Total RNA was quantified with a NanoDrop ONE c 497 spectrophotometer (Thermo Scientific). First-strand cDNA was synthesized from the total RNA using 498 a PrimeScript TM RT Reagent Kit with gDNA Eraser (Takara) following the manufacturer's protocols. The 499 generated cDNA was used as a template in real-time PCR with Fast Q-PCR Master Mix (SMOBIO), which 500 was run on a QuantStudio 7 Flex (Applied Biosystems, Thermo Scientific) machine with three-step (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint standard 1,6-diaminohexane and adjusted to pH≥12 with 0.5 mL of 5 M NaOH. The samples were then 529 subjected to sequential N-ethoxycarbonylation and N-pentafluoropropionylation. For DM2 samples, 530 biopsies were also resuspended in 0.2 M HClO4 and processed as described above. GC-MS analyses 531 were performed with an Agilent 7890B gas chromatograph coupled to a 5977B quadrupole mass 532 selective detector (Agilent Technologies, Palo Alto, CA). Chromatographic separations were carried 533 out with an Agilent HP-5ms fused-silica capillary column. Mass spectrometric analysis was performed 534 simultaneously in TIC (mass range scan from m/z 50 to 800 at a rate of 0.42 scans s-1) and SIM mode 535 (put, m/z 405; spd, m/z 580, N1-acetyl-spm, m/z 637; spm, m/z 709).   (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted April 1, 2022. ; https://doi.org/10.1101/2022.03.31.486563 doi: bioRxiv preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.   Suppl. Figure 4: A) Representative immunostaining for Synaptophysin (green), bungarotoxin (red) and Caveolin-3 (cyan), the first image is an overlay among synaptophysin and bungarotoxin. Scale 0.5 mm B) Relative expression of Fbox32 Trim63 during denervation in TA and EDL muscle (n = 3, values represent mean SD, ***P < 0.001, ****P < 0.0001; by 1-way ANOVA Tukey's multiplecomparison test). C) Pseudo-bulk KEGG pathway enrichment analysis of whole Crush 3 days sample. D) Violin plot showing the expression of Fbox32 Trim63 in the different ST timepoints. E) Muscle section compartmentalization strategy for microdissection. Immunostaining for Laminin (red). Scale 1 mm Areas highlighted in red have been used to calculate CSA from oxidative and glycolytic parts of the muscle Suppl. Figure 5: A) Relative expression levels of Amd2 scale 2mm (white dotted lines highlights EDL muscle) B) Serial section stained with Hematoxylin and Eosin illustrating muscle localization C) GC-MS quantitative analysis of Spermidine and Spermine levels in control (CTR) and denervated (Crush 3 days) muscles (n = 4, values represent mean SD by 2-tailed Student's tstudent test). D) Relative expression of Odc1 during denervation in TA (n = 3, values represent mean SD,*P < 0.05 by 2-tailed Student's t-student test).