A data-driven analysis of ber type architecture over the entire muscle

Skeletal muscle function is inferred from the spatial arrangement of muscle bers architecture, which corresponds to myober molecular and metabolic features. Myober types can be distinguished by the expression of myosin heavy chain (MyHC) isoforms, representing contraction properties. In most studies, myober typing is determined from a local sampling, typically obtained from the muscle median region. This median region is assumed to represent the entire muscle. However, it remains largely unknown to what extent this local sampling represents the entire muscle. We present here a pipeline to study the architecture of muscle ber type over the entire muscle, from sectioning, staining, imaging to image quantication and data-driven analysis.


Methods
We present here a pipeline to study the architecture of muscle ber type over the entire muscle, from sectioning, staining, imaging to image quanti cation and data-driven analysis.

Results
We reconstructed muscle architecture from consecutive cross-sections stained for laminin and MyHC isoforms. Examining the entire muscle using consecutive cross-sections is extremely laborious, we provide consideration to reduce dataset and yet to cover the entire muscle. Analyses of over 15,000 myo bers, showed spatial variations in myo ber geometric features, myo ber type and the distribution of neuromuscular junctions along the entire muscle.

Conclusions
We suggest that asymmetric spatial distribution of myo ber types, geometric features of myo bers and the neuromuscular junctions along the muscle could affect muscle function. Therefore, instead of a single sampling from a median region, representative regions covering the entire muscle should be investigated in future studies.

Background
Skeletal muscles facilitate mobilization and stability of the skeleton, which is greatly determined by muscle ber architecture. Muscle ber spatial arrangement is broadly described by their orientation relative to the axis of force generation. In general, three main ber arrangements have been described: 1) bers arrangement parallel to the force-generating axis on both ends of the muscle. 2) bers that are oriented at a single angle relative to the force generating axis, and 3) bers that are oriented at several angles relative to the axis of force generation (Lieber Richard and Ward Samuel, 2011; Lieber and Fridén, 2000; Otten, 1988). It has been recognized that we still far from understand how myo ber properties affect contraction and the speci c movement of muscle groups (Herzog, 2017).
In addition, to muscle ber architecture, contraction properties are also determined by myo ber types.
Fast-and slow-twitch muscle bers can be distinguished by metabolic pathways (Lang et  A simulation of muscle movement has been suggested using architectural models, suggested that even in a muscle with a parallel muscle ber arrangement, as the tibialis anterior (TA), the contraction is not symmetric along the muscle longitudinal axis (Moo et al., 2016;Moo et al., 2020;Otten, 1988;Van Leeuwen and Spoor, 1992). The tibialis anterior (TA) connects the knee to the foot, mobilizing the lower leg (Mathewson et al., 2012). Modelling of myo ber architecture from non-invasive imaging procedures also reported asymmetric myo ber organization along the TA longitudinal (Damon et al., 2011;Sullivan et al., 2019). These reports demonstrate the importance of analysis of myo bers over the entire muscle. However, they do not discriminate between myo ber types.
We investigated myo ber type arrangement over the entire muscle. We present a pipeline to investigate myo ber architecture using imaging and image quanti cation over the entire muscle, including essential steps to reduce dataset size. We show an asymmetric distribution of myo ber geometrical features along the muscle and suggest a pattern that could relate to myo ber-typing and the distribution of neuromuscular junctions.

Mouse and immunohistochemistry procedures
Two TA muscles were excised from two 8-weeks-old male C57BL/6J wild type mice and immediately snap froze in liquid nitrogen in the presence of iso-pentane, as described in (Raz et al., 2019;Riaz et al., 2015). One TA muscle was sectioned in a transverse axis and the other one in a longitudinal axis. Every muscle was sectioned from one end to the other. For each staining, at least three longitudinal sections were used. All 603 transverse sections of 10 µm thick, covering nearly the entire muscle, were stained. The tips at the distal or proximal ends were excluded, as those could not be sectioned. We measured the length of the TA muscle prior to imaging, and it was just over 6 mm ( Fig. 1), indicating that sectioning covered nearly the entire muscle. Sections were made with the CM3050-S cryostat (Leica Germany), and were pasted on "PTFE" printed slides, 12 well -5 mm diameter, (Electron Microscopy Sciences, USA).
Sectioning of the entire muscle was carried out in a single session, single day without removing the muscle from the cryotome. This eliminates differences in the cutting angle (tissue orientation, which will affect myo ber geometry) and batch effect. From the TA muscle of the second mouse, longitudinal sections were generated. Muscle sections were stored at − 20 °C prior to staining.
The immunohistochemical procedure was performed with a staining protocol that is detailed in (Raz et al., 2019). Primary antibodies that were used in this study are: Rabbit-anti-Laminin (1:1000, Abcam, UK).  (Riaz et al., 2016), and α-Bungarotoxin_Alexa uor-488-conjugated (1:1000) was purchased from Thermo Fisher. Slides were mounted with ProLong Gold (Invitrogen, USA). The α-Bungarotoxin antibody was used to stain the neuromuscular junction (NMJ). All cross-sections were simultaneously stained using a single antibody mix, in order to eliminate variations between staining sessions.
The length of the entire muscle or individual muscle bers was measured from longitudinal sections using the microscope scale bar.

Imaging
Imaging was performed with the Pannoramic 250 Flash III slide scanner (3DHISTECH, Hungary). Imaging of all cross-sections was carried out in a single session, eliminating batch effects. Imaging of the α-Bungarotoxin signal was performed with both the slide scanner and a higher resolution with the Leica DM5500. The α-Bungarotoxin signal in NMJs was manually counted from digital magni cations of images. Counting was carried out from three sections, by two researchers. The reported counts are an average of two researchers and three sections.
Image processing, quanti cation and alignment Slide scanner images were stitched, and les were converted to TIFF using CaseViewer (3DHISTECH, Hungary). After conversion, sections were cropped and scaled using custom MatLab (MathWorks Inc.) scripts. Scaling of the images (down-sampling) was crucial to decrease image size and resolution and thereby to increase the processing speed. All sections were exported at a lower resolution scale. Manual investigation of the images shows broken, folded or poorly mounted sections, which were excluded from further analysis (examples are in Fig. S1). This is expected as all wet lab procedures were carried out manually.
The Tiff images were processed and subsequently were segmented and quanti ed using the MuscleJ macro in ImageJ (Mayeuf-Louchart et al., 2018). In brief, myo ber objects were identi ed after segmentation of laminin staining, as reported in (Riaz et al., 2016). An example of tissue outcome is shown in Fig. S2. We applied the following modi cations toe original MuscleJ: the image processing step (prior to segmentation) was changed to match our image resolution, implementation of an automatic loading of the whole image set, and an automatic saving of the segmented mask, circularity, the spatial position (X, Y coordinated) was monitored. The output le included: section number, position, myo ber number, CSA, circularity, spatial position, myHC-2B and MYHC-2A mean uorescence intensity (Table S1). The macro can be obtained upon request. We noticed that this uniform image processing and segmentation settings were not optimal for all sections, but we compromised for a uniform procedure over manual optimization for every section.
We then excluded sections under the following criteria: myo bers were not identi ed in > 70% of the section area, and gaps between consecutive images are not larger than 150 µm. This step was also bene cial to allow proper componential capacity. In total, from the initial 603 Sect. 150 sections remained for further analysis. Representative segmented images are found in Fig. S3). This step was also bene cial to reduce the dataset size otherwise processing took a few days. The remaining sections were then mapped to their physical position along the muscle longitudinal axis. Section position was calculated with the rst section at zero µm and consecutively increasing by 10 µm (sections thickness), including the position of the excluded sections. Median and variance for both cross-sectional area (CSA) and circularity were calculated per section. MyHC mean uorescence intensity (MFI) was obtained after background correction that was separately applied to each image and each uorophore, as described in (Riaz et al., 2016). Additionally, a '2A mask' was generated marking the 2A-positive bers, which were added onto the laminin-segmentation mask. Fibers were considered 2A-positive when their MFI was greater than the MFI mean + one standard deviation of all the bers per section.
The alignment of the images was performed with the '2A masks' in ImageJ using Linear Stack Alignment with the scale-invariant feature transformation (SIFT) algorithm (Lowe, 2004). This method aligns sections consecutively one by one.
Myo ber composition analysis was carried out using R (version 3.5.1). In total, 153930 bers were obtained from the eligible 150 cross-sections and were used in the data-driven analysis. In order to reduce technical differences between images, MFI values were scaled per section. Scaling of myo ber dataset has been demonstrated to be bene cial for myo ber type data-driven analysis (Y et al., 2020). Data-driven analysis of myo bers was carried out after a natural logarithm (ln) transformation of MFI values.
Smoothing lines were calculated and plotted using the smooth.spline function in R, using smoothing parameter = 1. Statistical signi cance was assessed with Pearson's chi-squared test.

Results
To develop a methodology reconstructing the 3D-architecture of myo ber type over the entire muscle, we selected the TA of a wild type mouse for its small, which is suitable for such an initial study. The work ow of this methodology is summarized in Fig. 1. One muscle was cross sectioned from the proximal to the distal end, and the other one was sectioned on the longitudinal axis. The entire muscle length was measured to be around 6 mm ( Fig. 2B and 2C), which is in agreement with the TA length that was reported in (Heemskerk et al., 2005) Image processing and quanti cation were carried out on the cross-sections, and data-driven analyses were carried out on data generated from cross-sections. The longitudinal sections were used as reference for coordinates and veri cations of results from the crosssections. A schematic representation of the TA anatomy is shown in Fig. 2A. The muscle ber (named here as myo ber) contours were stained with an anti-laminin antibody, and myo ber types were recognized with an antibody mix to MyHC-2A and MyHC-2B ( Fig. 2C and 2D). Both MyHC isoforms are commonly used for myo ber typing in mouse ( MyHC-2B myo bers encompassed the entire muscle length (Fig. 2C). The length of MyHC-2A myo bers was estimated from three longitudinal sections, (ranged between 160-250 µm), which was used as a criterion for section exclusion.
The wet-lab and imaging steps: cryo-sectioning, immunostaining and imaging were completed within a single batch to avoid technical variations. Sectioning the entire muscle in a single day, without removing the mounted muscle from the machine, ensured that the section angle did not change throughout the entire muscle. Together, this procedure minimizes technical variation that could give rise to potential confounding batch effects. Imaging was carried out in a single batch using a slide scanner, all slides were imaged using the same settings. From all acquired images myo bers were segmented based on the lamina staining using single image processing and image quanti cation procedures. To assess this protocol, the results of lamina segmentation were manually compared to the original images. We noticed that the quality of part of the sections was not su cient for further processing, either due to technical issues (sectioning or folded/broken tissue sections) or bad segmentation. Those sections were excluded under the general guideline for data reduction, allowing a gap of maximum 150 µm between two consecutive sections. This gap size was smaller than the estimated shortest myo ber length (160 µm).
With this criterion we ensured that all myo bers in the muscle tissue will be included in subsequent analyses. In addition, this gap distance was kept constant along the muscle to avoid a bias for over-or under-representative areas.
To investigate changes in myo ber spatial organization along the longitudinal axis we aligned the crosssections using the laminin-based segmentation and MyHC-2A positive bers (Fig. 2E). The MyHC-2A myo bers were de ned as positive when the mean uorescence intensity (MFI) was greater than mean MFI + one standard deviation per section. Reconstruction of muscle ber atchitecture along the longitudinal axis was carried out with the SIFT algorithm. A good alignment was found for consecutive sections over a short distance (AVI. S1). However, the alignment over a larger distance required compromising with reduced dataset for a faster solution. The alignment outcomes that were obtained with gaps between sections were comparable to the alignment without gaps in the same region (AVI. S2, AVI S1, respectively). Since the alignment procedure only accounted for consecutive sections, consecutive alignment of pairs of sections (local) was hampered by the accumulation of errors over sets of sections. Therefore, alignment over the entire muscle was performed in parts. We obtained good alignment results from both proximal and median regions of the muscle, whereas alignments of sections at the distal region were poor (Fig. S1). A good alignment over distance can result from a parallel orientation of muscle bers. A poor alignment over a small distance could indicate bending of the muscle or a notparallel orientation of muscle bers.
The change in myo ber alignment along the longitudinal axis suggests that myo ber properties might alter over the entire muscle. To investigate this, we considered only cross sections with good segmentation (> 85% of the sections were eligible based on segmentation). In total, 153,930 objects from 150 cross-sections were included in the analysis (the full dataset is found in Table S1). From the laminin segmented objects, representing myo ber objects, we extracted cross-sectional area (CSA) and circularity. As the entire muscle was sectioned in a single session, the same cutting angle was kept along the entire muscle, therefore a change in CSA or circularity could not re ect alterations in sectioning orientation but can indicate a biological phenomenon. From each cross section we calculated the median CSA, mean circularity and the variance of both features. A smoothed regression line (curved), which levels off small variations between adjunct cross sections, was applied to assess the change along the entire muscle (Fig. 3). The smoothen regression line enables considering major changes along the longitudinal axis, while disregarding the effect of an individual cross-section or a small group of consecutive sections. The shape of the regression line for the number of myo bers follows the macro shape of the TA muscle: smallest at the proximal end, large at the median part followed by smaller size at the distal end (Fig. 3A).
At the central region we noticed few groups of consecutive sections with distinguished numbers of myo bers (Fig. 3A, 2500-3500 um), since those did not in uence the regression line, they were not considered here. The median CSA plot showed two main regions, in sections at the rst ~ 1/3rd of the muscle the median CSA was smaller than the other part of the muscle. The regression line showed a constant CSA along the other 2/3rd part of the muscle (Fig. 2B). The regression line of the mean circularity showed a similar trend as median CSA along the muscle (Fig. 3C). However, the variance of circularity showed an opposite trend along the muscle (Fig. 3C). At the proximal end the circularity median was low, indicating an elongated orientation of the myo bers, and high variance, indicating heterogeneity. In contrast, at distal and median regions variance in circularity was low and median circularity was high (Fig. 3C). In contrast, the changes in median and variance CSA along the muscle were similar (Fig. S4). Using the intersection between the mean circularity regression line and the circularity variance regression line we de ned the proximal end (Fig. 3C). We then superimposed the CSA and circularity curves in order to assess a spatial relation between these two geometrical features. The intersection between the regression de ned the distal end sub-region (Fig. 3D). Together, the features subregions can be summarised as follows (see also in Fig. 3F): the proximal end (< 1/9 of the muscle) myo bers CSA is small and circularity is low with high variance indicating heterogeneity in myo ber structure. Consequently, in region I (< 1/9 of the muscle) circularity increases but myo ber CSA remains small. Then after, in the median part (~ 5/9th of the muscle length), myo bers have high CSA and circularity values. The small variance in circularity suggests a more homogenous structure of bers. At the distal end, the relations between circularity and CSA were inverse as compared with the preceding region (Fig. 3D). This region is also the one with poor alignment (Fig. S1).
We then investigated if these muscle regions, as de ned by geometrical features of myo bers, can be distinguished by molecular features of muscle biology. We focus on the sarcomeric proteins, MyHC-2B and MyHC-2A, and the distribution of NMJs. The cross-sections were stained with antibody mix to MyHC isoforms, and after imaging the MFI were measured from each myo ber. Main myo ber-type subclasses were determined in the pooled dataset using density plots. The MyHC-2B density plot showed only a single peak, indicating a predominant one myo ber type (Fig. S5A), but two peaks were found in the MyHC-2A density plot (Fig. S5B). This indicated that in TA two myo ber populations are found with high and low MyHC expression, whereas the MyHC-2B myo bers can be described as a homogenous population. For a quantitative assessment of the distribution of the MyHC-2A myo bers along the muscle, we discriminated between the MyHC-2A positive myo bers and the MyHC-2A negative myo bers using density distribution of myHC-2A MFI (Fig. S5B). We then calculated the percentage of 2A-positive myo bers in every section (Fig. 3E). Overall, the regression line of the proportions of the 2A-positive myo bers is consistent with that of the median CSA, the proportions of 2A-positive myo bers were higher at the proximal region, compared with median and distal regions (Fig. 3E). Interestingly, a signi cant difference of MyHC-2A positive myo bers was found between the percentages of the 2A-positive bres in each muscle region (p-value < 2.2*10 − 16 ). The asymmetric distribution of MyHC-2A myo bers along the longitudinal axis suggests that contraction properties are not uniform along the muscle.
Last, we investigated the distribution of NMJs along the muscles and visually assessed a relation to muscle regions. The NMJs were visualized in three longitudinal sections with α-bungarotoxin, which speci cally binds to NMJs (Tse et al., 2014). A representative image is shown in (Fig. S6). The four muscle regions were estimated from the cross-sectional study, using the scale of the muscle longitudinal axis. In this muscle, neuromuscular junctions were absent at the proximal-end region and were rare at the distal end (Fig. 3F). The majority of NMJs (85%) were counted in the median region of the muscle (Fig. 3F).

Discussion
The impact of myo ber type on muscle physiology and function has been mostly interpreted from cross sections at the median part of the muscle. Several reports reinforced the importance to sample the entire muscle in order to understand how myo ber architecture affects muscle function (Moo et al., 2020;Stark and Schilling, 2010 Stark and Schilling, 2010), we also considered the spatial organization of myo ber types over the entire muscle. Our methodological paper, therefore, also reinforces the importance of adequate spatial sampling of whole muscles in future studies.
Modelling muscle ber architecture has been studied using non-invasive methods and ex-vivo imaging procedures. For the mouse TA a parallel myo ber arrangement has been reported, where most bers are parallel to the longitudinal axis at the largest part of the muscle, but close to the tendon at the proximal and distal ends myo bers curve (Heemskerk et al., 2005;Lovering et al., 2013;Moo et al., 2016). This indicates an asymmetric myo ber architecture along the muscle. Ex-vivo measurements of sarcomere force-length demonstrated that the muscle ber architecture affects its biomechanical properties (Moo et al., 2020). Ex-vivo procedures of muscle architecture can in uence the shape of the muscle during the following procedures: 1. pulling when dissecting the muscle out of the mouse. We were careful avoiding pulling tissue, yet this can in uence when comparing between different muscles. 2. As the muscle is not connected to the tendon, the natural muscle-tendon tension cannot be considered. This can affect the muscle tension and length, but how fast this muscle-tendon tension is lost is not well unknown. But a fast tissue freezing, immediately after muscle harvesting, we expect to reduce this loss of tissue tension and length. However, freezing artefacts could also contribute to variations between muscle tissues. 3. Mounting would affect myo ber orientation. Therefore, in our work ow sectioning of the entire muscle was carried out in a single day, without dismounting of the muscle from the cryotome, and the section angle relative to the muscle axis was identical across all sections. Despite those all possible effects, our observations using ex-vivo analysis of myo ber organization over the entire murine TA muscle are in agreement with the in vivo studies. We also found asymmetric myo ber arrangement along the longitudinal axis, and between proximal and distal regions. At the proximal end myo bers were parallel to the longitudinal axis, but at the distal end myo ber orientation was altered, suggesting bending. Ideally, an ex-vivo study of myo ber architecture would be performed in sections of an entire lower leg, in which a muscle is connected to the bone(s) and other tissues exactly as in the full organism. This approach could be investigated in future studies, making use of the work ow we are reporting here.
The distribution of myo ber types along the entire muscle was also previously reported (Wang et al.,

2002)
, but this is the rst study combining immuno uorescence of MyHC isoforms and data-driven analysis. Compared with the ATPase-based analysis of myo ber typing, image analysis of MyHC immuno uorescence can be applied for image quanti cation and more robust downstream analysis of myo ber typing (Mayeuf-Louchart et al., 2018). Furthermore, the four most abundant MyHC isoforms can be investigated in a single study (Riaz et al., 2016;Y et al., 2020). Using ATPase staining, an asymmetric distribution slow-twitch muscle bers was found in rat limb muscles, including the TA (Wang et al., 2002). Here we also report asymmetric myo ber distribution along the murine TA muscle. The MyHC-2B positive myo bers seem to stretch over the entire muscle length, but the MyHC-2A positive myo bers are more abundant in the proximal part of the muscle. The other two MyHC isoforms that are often considered in mouse muscle ber typing are MyHC-1 and MyHC-2X (Rao and Mohanty, 2019; Riaz et al., 2016). We excluded MyHC-1 in the antibody mix, as its expression in TA muscle is limited. We also excluded the MyHC-2X as imaging with the Pannoramic 250 Flash III slide scanner was not optimal for 4 uorophores.
Both MyHC-2B and − 2A isoforms indicate different metabolic processes, their asymmetric distribution may suggest differences in contraction properties along the muscle. We could not conclude whether the MyHC-2B myo bers at the distal end are the curved counterparts of the myo bers seen in the median part of the muscles.
In contrast to most studies that determine myo ber typing from a single region (often median), this study investigates the distribution of myo ber typing along the entire muscle. However, including all cross sections covering the entire muscle is highly laborious and requires an enormous computational power. Here we show that using the shortest myo ber length as a guideline for gaps between consecutive crosssections enabled us to eliminate proportions of the cross-sections without losing spatial information.
Since we show here that myo ber typing changes along the muscle, it potentially implicates that diseaseand age-associated alterations in myo ber type composition could be more affected in speci c muscle regions.
Analytical two-dimensional skeletal muscle models predicted that bers' shape affect intramuscular pressure distribution (van Leeuwen and Spoor, 1996). Other models suggest that pressure distribution along the muscle is not even, the high-pressure region presents half-way along the muscle length, which supports bending of the muscle bres (Van Leeuwen and Spoor, 1992). At the centre of the TA muscle we found that both CSA and circularity have the largest range. This could indicate that the geometry of those myo bers support high pressure. The opposite myo ber geometry was found at the proximal end. The small CSA and circularity could indicate that the proximal end of the TA is a vulnerable area for muscle pressure. The observed asymmetry of geometrical features of myo bers along the entire muscle is consistent with biomechanical models of the TA muscle (Heemskerk et al., 2005;Lovering et al., 2013;Moo et al., 2016). Moo et al., also measured sarcomere force at the median region of a muscle and stressed that their ndings may not be valid to other muscle regions, especially where myo ber lengths in not uniform (Moo et al., 2020). Under the assumption that we aim to describe major changes along the muscle, we applied a smoothed regression line to describe changes in myo ber features along the muscle. In this analysis, we discarded changes in CSA and myo ber numbers were found in only a small number of sections at the median part of the muscle. Combining the regression lines of CSA and circularity we described four muscle regions and suggest a relation with MyHC-2A positive myo bers and the NMJs: the largest muscle part, the median region (region II) is characterized by high circularity and high CSA, less MyHC-2A positive myo bers and higher density of NMJs. The shortest regions are the proximal and distal ends. Both differ in the orientation of myo bers (varied circularity), the percentage of MyHC-2A positive myo bers, and have the lowest proportions of NJMs. Based on circularity, myo bers at the proximal end are more heterogeneous compared to those at the distal end. The second sub-region (region I), could be a transition region between proximal and median regions. It was reported that higher NMJs concentration implicated higher contraction (Lomo, 2003), suggesting that in the TA the signal is mainly found at the median part of the muscle. Yet, our observations were made from a single muscle and conclusions should be cautiously made.

Conclusion
In this methodological paper we describe a pipeline to investigate myo bers features along the entire muscle. Based on geometrical features of myo bers we show asymmetric distribution and a spatial pattern along the muscle and distinguish between four main regions. The distribution of myo ber types and NMJs along the muscle agrees with that of myo ber geometrical features, suggesting an asymmetric myo ber contraction along the muscle, which can play a role in muscle physiology and pathological conditions. Availability of data and materials -All scripts are available upon request. The Myo ber dataset that was used in this study is available in Table S1.

Abbreviations
Competing interests -All authors have no con ict of interests on the submission form.
Funding -the French Muscular Dystrophy Association #26110 to VR.
Authors' contributions -, Wetlab experiments were performed by DB and CSO; imaging, image processing and image analysis were performed by DB and LMV; muscles were collected by MvP; project supervision was carried out by LMV, EvdA and VR; funding was received for VR; All authors read and approved the nal manuscript.  Figure 1 A methodological summary of myo ber data collection and analysis. The work ow of the transverse sections is shown at the left side. The steps from sectioning to imaging were carried out stimulatingly (in brackets) for all sections to avoid technical differences between sections. Exclusion of sections was based on the segmentation step and it was a manual step. The asymmetric distribution of myo bers and expression of markers were veri ed in the longitudinal sections.