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

Striated muscle is a highly organized structure composed by well-defined anatomical domains with integrated but distinct assignments. So far, the lack of a direct correlation between tissue architecture and gene expression has limited our understanding of how each unit responds to physio-pathologic contexts. Here, we show how the combined use of spatially resolved transcriptomics and immunofluorescence can bridge this gap by enabling the unbiased identification of such domains and the characterization of their response to external perturbations. Using a spatiotemporal analysis, we followed the changes in the transcriptomics profile of specific domains in muscle in a model of denervation. Furthermore, our approach allowed us to identify the spatial distribution and nerve dependence of atrophic signalling pathway and polyamine metabolism to glycolytic fibers. Indeed, we demonstrate a pronounced alteration of polyamine homeostasis upon denervation. Our dataset will serve as a resource for future studies of the mechanisms underlying skeletal muscle homeostasis and innervation. Graphical Abstract


INTRODUCTION
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.
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 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, expression levels could result in an imbalance in the PA pathway. The decarboxylation of S-238 adenosylmethionine (SAM), catalyzed by Amd1/Amd2, produces an aminopropyl group that acts as a 239 substrate together with Putrescine for the Spermidine synthase and subsequently to form Spermine 240 ( Figure 6A). Thus, reduction of Amd1/Amd2 expression may results in Putrescine accumulation.

241
Indeed, as shown by gas chromatography-mass spectrometry (GC-MS) quantitative analysis, we 242 observed a dramatic change in polyamine ratio 3 days after denervation ( Figure 6F). A mean 25-fold 243 increase of Putrescine was revealed as results of Amd1 limitation. However, although their precursor 244 is significantly increased the amount of spermidine and spermine does not change significantly ( Figure   245 Suppl. 5C). Intriguingly, in support of this observation an induction of Odc1, required to synthesize 246 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

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

320
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.         Methods parameters (normalization.method="SCT",= and reduction = "cca"), and label transfer was performed 478 with the default parameters. GO enrichment analysis was performed with gprofiler2 using as input 479 cluster markers present at least in 50% of each cluster (pct.1>0.5). Pseudobulk data were obatained 480 using "AverageExrepssion" functions and subsequently filtered for expression above 0.1 and lg2FC 481 (3days vs CTR) >-0.3 or < 0.3. KEGG pathway analysis of genes upregulated in 3days vs CTR was 482 performed using Enrichr (Kuleshov et al., 2016). After analysis, for visualization purposes, images were     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).