{'fraction_non_alphanumeric': 0.032082205749696044, 'fraction_numerical': 0.018041338196650585, 'mean_word_length': 4.7648654759213205, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-resolution genome-wide expression analysis of single myofibers using SMART-Seq

Skeletal muscle is a heterogeneous tissue. Individual myofibers that make up muscle tissue exhibit variation in their metabolic and contractile properties. Although biochemical and histological assays are available to study myofiber heterogeneity, efficient methods to analyze the whole transcriptome of individual myofibers are lacking. Here, we report on a single-myofiber RNA-sequencing (smfRNA-Seq) approach to analyze the whole transcriptome of individual myofibers by combining single-fiber isolation with Switching Mechanism at 5′ end of RNA Template (SMART) technology. Using smfRNA-Seq, we first determined the genes that are expressed in the whole muscle, including in nonmyogenic cells. We also analyzed the differences in the transcriptome of myofibers from young and old mice to validate the effectiveness of this new method. Our results suggest that aging leads to significant changes in the expression of metabolic genes, such as Nos1, and structural genes, such as Myl1, in myofibers. We conclude that smfRNA-Seq is a powerful tool to study developmental, disease-related, and age-related changes in the gene expression profile of skeletal muscle.

Skeletal muscle is composed of a variety of different cell types, including endothelial cells, fibro/adipogenic cells (FAPs), 3 adipocytes, mesenchymal cells, and fibroblasts, among others (1)(2)(3)(4). Further heterogeneity of skeletal muscle is also manifested by the diversity in the composition of myofiber types that constitute muscles (5). Skeletal muscle fiber types are often categorized based on their contractile properties, giving two broad categories: fast-twitch muscles and slow-twitch muscles (5,6). These fiber types can be further subcategorized based on metabolic properties and myosin heavy chain isoforms (5,7).
Methods of investigating changes in fiber type in response to different stimuli rely on staining and biochemical analyses of individual fibers, biochemical analyses of the whole muscle, or sequencing of entire muscles (8). Standard bulk RNA-Seq is not suited for the analysis of myofibers at high resolution, because it captures the entirety of the muscle tissue, resulting in the pooling of different fiber types in addition to nonmyogenic cells (8).
Consequently, in such studies, the myofiber-specific gene signature cannot be inferred based on RNA-Seq of whole muscle tissue. The emergence of single-cell technology provides ample opportunity for further investigation into the heterogeneity of single muscle fibers at the transcriptome level. Here, we combine single-myofiber isolation of the extensor digitorum longus (EDL) in Mus musculus with Switching Mechanism at 5Ј end of RNA Template (SMART) technology to analyze the whole transcriptome of individual myofibers ( Fig. 1) (9), in a method called single-myofiber RNA sequencing (smfRNA-Seq).
We describe a robust method to extract RNA from a single myofiber followed by generation of sequencing ready libraries and whole transcriptome analysis. Using this technique, we first determined the genes that are found in the whole muscle that are not produced by the myofiber and are instead produced by nonmyogenic cell types. To demonstrate the effectiveness of this technique, we next went on to analyze the differences in the whole transcriptome between myofibers isolated from young and old mice. smfRNA-Seq proved to be a useful tool in determining gene expression changes that occur in myofibers between different conditions.

Isolation of high quality mRNA from single myofibers
One of the intentions of our novel method is to give researchers the ability to sequence RNA from a single myofiber without the confounding presence of other cell types. A potential source of unwanted signal in the sequencing of myofiber RNA is the presence of muscle stem cells, also known as satellite cells, that are physically associated with the fibers. Using our method, we have found that satellite cells can be almost completely cro ARTICLE removed from the fibers with the addition of trypsin to the collagenase digestion buffer at a final concentration of 0.25%. This process does not damage the fiber itself, with a proper EDL isolation yielding over 200 myofibers per mouse (Fig. S1A). However, this buffer effectively strips the myofibers of their satellite cells as shown by the reduction in the number of PAX7 ϩ cells per fiber (Fig. S1, B-D).
Because muscle fibers are very tough and do not readily break down under normal lysing conditions, extracting the RNA from a single myofiber can prove challenging. With whole muscle, a method to overcome this is by freezing the muscle in liquid nitrogen and grinding it into a powder with mortar and pestle (10). However, this method cannot realistically be done to a single fiber while still collecting all of the RNA. Therefore, we lysed the fiber with lysis buffer in RNase-free water, utilizing osmotic pressure and gentle pipetting to break down the fiber and retrieve the intact RNA. This method proved effective because more than an adequate amount of RNA was recovered, even from a single myofiber, for use with SMART-Seq technology (Table 1).

Single-myofiber RNA-Seq
From the extracted RNA, we successfully generated sequencing-ready cDNA libraries using the DNA SMART-Seq HT kit (Takara Biosciences) in combination with the Nextera XT DNA library preparation kit (Illumina). Final single-myofiber sequencing-ready libraries were of adequate quantity and of ideal size for sequencing, with fragments being of comparable average size as those generated from traditional whole muscle RNA-Seq libraries (Fig. S1, E and F). This implies that our fiber RNA extraction procedure generated high-quality starting material that was compatible with low-input library preparation technologies.

Comparative analysis of whole muscle and single-myofiber RNA sequencing
By evaluating the number of unique reads obtained from smfRNA-Seq samples, we found that sequencing depth of single-myofiber libraries was comparable with whole muscle RNA-Seq ( Fig. 2A). We obtained an average of 24 million unique reads from single-myofiber samples when multiplexing 12 samples per lane on a NextSeq500 and an average of 35 million unique reads when multiplexing 10 samples per lane. This corresponds to an average overall alignment of 82.42%, with an average unique alignment percentage of 64.47% (Table  S3). The overall expression profile of the single myofibers is also very similar to that of the whole muscle, indicating a similar high resolution ( Fig. 2A). When performing principal component analysis (PCA), single-fiber samples cluster together on both axes but away from the whole muscle on the PCA2 axis (Fig. 2B). To compare the efficiency of single-cell technologies with increasing input, we also generated RNA-Seq libraries from five and twenty myofibers. When the number of fibers is increased, samples become less similar to both the whole muscle and the single-myofiber samples, with higher variation between technical replicates (Fig. 2, A and B). This is most likely due to the excessive quantity of sample and could be resolved by scaling up the protocol.
When looking more in depth at individual genes, we see many similarities between the single myofiber and the whole muscle but also crucial differences. Of particular note, we see that muscle-specific genes have similar numbers of reads between the single-fiber and whole-muscle samples (Fig. 3, A-C). The Myh cluster codes for a variety of myosin heavy chain proteins, which are the motor proteins of muscle whose various isoforms are the basis of the different fiber types (7,11). The similar expression between the single fiber and the whole muscle conclusively shows that the RNA sequenced came from a myofiber alone (Fig. 3A). For further confirmation, we also display Ckm, the muscle-specific creatine kinase (12), and Acta1, which codes for skeletal muscle ␣-actin (13). These genes have the same pattern as is seen with Myh (Fig. 3, B and C).
Differences between the single-myofiber and whole-muscle transcriptomes arise from nonmuscle cells. Using smfRNA-Seq, we are capable of completely removing these unwanted cell types to sequence only the myofiber. To demonstrate the removal of cell populations other than the myofiber, we looked at cell-specific genes of a variety of different muscle-resident cell types. We first looked at the expression of the satellite cell marker Pax7 and see that there is no expression of Pax7 in the single-fiber transcriptome despite its presence in whole muscle samples (Fig. 3D). We also analyzed markers for fibroblasts, namely Col1a1 and Thy1 (14,15). As expected for these genes, they are expressed in the whole muscle but not expressed in the single fiber (Fig. 4, A and B). Endothelial cells are also depleted in single fiber samples, because their markers Kdr and Pecam1 are expressed in whole muscle but not in the single fiber (16, 17) ( Fig. 4, C and D). We also analyzed Retn to identify the presence of adipocytes (18), Cd34 for hematopoietic cells (19), Ly6a as a marker for FAPs (2), and Adgre1 for macrophages (20) (Fig. 4, E-J). As expected, in the single-fiber transcriptome there is no expression for any of these genes, all of which were present in the whole muscle. These results clearly demonstrate the purity of smfRNA-Seq in sequencing only the myofiber without confounding cell types.
From these data, we performed a Gene Ontology (GO) term analysis of genes that are expressed solely in the WM samples, defined as having an RPM value of at least 10 in WM and 0 in single-fiber, indicating that they originate from nonmyogenic cells. Using these criteria, we identified 445 genes that were solely expressed in the WM muscle samples (Table  S1). Some of the more up-regulated genes were involved in the formation of the extracellular matrix or immunity (Table  S1).
We also attempted to identify genes that were exclusively produced in the myofiber. We defined these genes as being more highly expressed in the single-fiber samples, with a q value Ͻ 0.01, being expressed at more than 10 RPM and being Single-myofiber RNA-Seq expressed at least 10 RPM more than in the WM samples. We identified 622 genes that matched these criteria and performed a GO analysis (Table S2). From this, we see that some of the top pathways associated with these genes are ribosomal proteins and the respiratory chain. However, this is not an exhaustive list, because of the difficulty in identifying genes expressed solely in myofibers when the WM samples are composed primarily of myofibers. For example, Myh4, which is expressed exclusively in muscle fibers, does not pass the criteria laid down and is not part of the list of genes expressed solely in myofibers.
On the other hand, dystrophin (Dmd), another myofiber-specific gene, does pass our criteria.

Age effect on the transcriptome of single myofibers
To verify the utility of our new method to analyze variation in the myofiber transcriptome under different conditions, we performed smfRNA-Seq on EDL myofibers isolated from young (1 and 3 months old) and aged (19 months old) mice. Using our technique, we see clear differences in the transcriptome of young and old myofibers. In the old myofibers we see the dereg-

Single-myofiber RNA-Seq
ulation of a number of genes, with 181 genes being significantly (p adj Ͻ 0.05) differentially expressed in old compared with young myofibers (Fig. 5, A and B). PCA demonstrates how young myofibers, isolated from 1-month-old mice, cluster together away from the old myofibers on the PC2 axis (Fig. 5E). Furthermore, the myofibers isolated from 3-month-old mice begin to resemble the 19-month-old myofibers while also maintaining similarities with the young 1-month-old myofibers, demonstrating the gradual change in transcriptome of the myofibers as the mouse ages (Fig. 5E). GO term analysis of the genes with the highest variation between young and old myofibers shows that deregulated pathways in aging include the transport

Single-myofiber RNA-Seq
of small molecules such as salts and metal ions and the synthesis of collagens. Additional examples of genes significantly deregulated in aging include Actc1, Myl1, Dkk3, Atf3, H19, Nos1, and Ndn (Fig. 6, A-H). These are genes involved in skeletal muscle structure, metabolism, growth, and maintenance of homeostasis.

Single-myofiber RNA-Seq
Discussion smfRNA-Seq is a powerful new technique that allows highresolution whole transcriptome sequencing of a single myofiber. Here, we have demonstrated that RNA can be successfully isolated from a single myofiber at a suitable concentration and quality to be used with SMART-Seq technology, generating a high-quality cDNA library for sequencing.
What distinguishes smfRNA-Seq from single cell RNA-Seq is that sequencing depth is comparable with large-scale bulk RNA-Seq (Fig. 2). However, as shown by PCA, single-myofiber samples are very similar to whole muscle samples with regards to PCA1 but dissimilar with regards to PCA2. By analyzing individual genes, we can interpret the similarities as being due to muscle-specific genes, whereas the lack of other cell types allows the single fibers to form a distinct cluster away from the whole muscle. Indeed, comparison of single myofiber to wholemuscle samples identified a list of genes expressed only in the whole muscle, which are likely responsible for this distinction (Table S1). This dissociation of the myofiber transcriptome from the transcriptome of other muscle-resident cell types will be of critical importance to tease out the contribution of the myofiber when a gene of interest is expressed in multiple cell types or when a treatment condition has a whole-muscle effect.
All of the fibers analyzed in this study were of the fast type. This is because the fibers were isolated from the EDL, which is composed of over 90% fast fibers (21). We confirmed that they were fast fibers by analyzing the levels of transcript for Tnnt3 and Tnni2, which are expressed in the fast muscle, as well as Tnnt1, which is only expressed in slow muscle fibers (22,23). We see that there is no expression of Tnnt1 and a high expression of the Tnnt3 and Tnni2 genes (Fig. S2, A-C). Because we could identify these fibers as fast-twitch, the method would be appropriate to distinguish differences in the transcriptome between fast and slow fiber types.
One condition that leads to changes in the numbers and in the gene expression profile of various muscle-resident cell types is aging (24,25). Because of the changes occurring in multiple cell types, including immune cells (26,27), FAPs (28), and satellite cells (29,30), whole-muscle RNA-Seq would be inappropriate to distinguish age-induced changes in myofibers alone. Using smfRNA-Seq, we were able to uncover differences in the transcriptome of old and young myofibers, without any added noise from age-related alterations in other cell types. Among the genes deregulated in aging myofibers are genes involved in muscle growth and structure, suggesting that some of these altered genes may be responsible for age-related muscle atrophy, also called sarcopenia. Of note, we see that Actc1 and Myl1 are greatly down-regulated in old compared with young myofibers (Fig. 6, A and B). Despite not reaching significance, there is also a moderate decrease in Acta1 with a log 2 fold change of Ϫ0.84. Because actins and myosins play a role in myofiber structure, force, and contractile properties, we suspect that the down-regulation of these genes may play a role in the weakening in contractile strength in aged muscles (31,32). A gene that is significantly up-regulated in aged myofibers is Dkk3, which encodes a secreted inhibitor of WNT signaling (Fig. 6C) (33). This finding is in line with a previous study demonstrating that increased Dkk3 in aged muscle tissue leads to muscle atrophy through autocrine signaling (33) while also confirming that the myofiber is the source of this secreted factor. There is also a down-regulation in Ndn in the old myofibers (Fig. 6G). Ndn produces Necdin, which is important in muscle regeneration, because it promotes the differentiation of satellite cells, and Necdin knockout mice have a significant impairment in muscle regeneration (34). Necdin was also shown to protect skeletal muscle from tumor-induced muscle wasting, also called cachexia, in part by binding to and inhibiting P53 (35). It would be of interest to see whether Necdin decrease in aging myofibers is also a cause of age-induced muscle atrophy.
In addition to genes involved in muscle structure and growth, our data also show that multiple genes involved in key metabolic processes and maintenance of homeostasis are also deregulated in old compared with young myofibers. For example, we detect a loss in the expression of the long noncoding RNA H19 (Fig. 6E). Although the role of H19 during myogenesis has been studied (36), namely in controlling myoblast proliferation, its role in fully differentiated myotubes is still unclear. However, H19 has also been shown to be involved in glucose metabolism and consequently, insulin sensitivity (37,38). We speculate that H19 reduction in old myofibers may therefore be a contributor to metabolic dysfunction, which is a common feature of aging (38).
In old myofibers, we also see an up-regulation of Atf3, which is a stress factor that is known to be up-regulated in muscle after a disturbance in homeostasis, such as after exercise (39) (Fig.  6D). Interestingly, Atf3 was shown to dampen the expression of pro-inflammatory cytokines after exercise (40). We speculate that its up-regulation in aging myofibers may occur in response to an age-related disturbance in homeostasis, such as the stress induced by chronic low-grade inflammation (41). Another gene that is up-regulated in aging is Nos1, which produces NO and plays a role in skeletal muscle metabolism (Fig. 6F). Notably, Nos1 is known to increase glucose uptake and inhibit mitochondrial respiration (42), which may be a contributor to agerelated metabolic dysfunctions. In addition to these top deregulated genes, GO term analysis revealed that genes involved in collagen synthesis, assembly, and degradation are all perturbed in aging myofibers (Fig. 5D). Because collagens are important components of the muscle extracellular matrix, our data suggest that old myofibers may contribute to changes in extracellular matrix stiffness (43,44).

Single-myofiber RNA-Seq
Altogether, we were able to demonstrate the effectiveness of smfRNA-Seq by sequencing the transcriptome of single myofibers from young and old mice. This method allowed us to conclude that aging does have an effect on myofibers at the level of the transcriptome, without the confounding signal from nonmyogenic cell types. This comparative analysis proves that smfRNA-Seq can thus be used to study other developmental or muscle-wasting disorders, such as muscular dystrophy. Additionally, smfRNA-Seq has the potential to be adapted to a larger scale to perform high-throughput analysis of numerous myofibers, which will be important in the study of myofiber heterogeneity.

Animal care
All procedures that were performed on animals were approved by the McGill University Animal Care Committee.

Accession numbers
All gene expression data reported in this study are available through GEO accession number GSE138591.

Buffers
Myofiber digestion buffer was prepared using 1000 units/ml of collagenase from Clostridium histolyticum (Sigma catalog no. C0130) in unsupplemented DMEM (Invitrogen catalog no. 11995073). Myofiber immunofluorescence blocking buffer is composed of 5% horse serum (Wisent catalog no. 065-250), 2% bovine serum albumin (Sigma catalog no. A8022), 1% Triton X-100 (Sigma catalog no. T9284), in PBS (Wisent catalog no. 311-425-CL). RNA extraction buffer is made using 19 l of the 10ϫ lysis buffer plus 1 l of RNase inhibitor from the SMART-Seq HT kit. 1 l of the previously composed 10ϫ lysis buffer is added to 9 l of RNase-free water to make 1ϫ lysis buffer.

Dissection of the EDL
The EDL was dissected from the hindlimb in the following manner: the skin of the hindlimb was removed by cutting around the ankle with a pair of scissors, and an incision was made along the ventral side of the leg. The epimysium around the tibialis anterior (TA) was removed, and the tendon of the TA was cut at the ankle while making sure to only cut the top tendon because the bottom tendon belongs to the EDL. Using a pair of forceps, the TA was gently peeled off the leg up to the knee and was then cut out as close to the knee as possible with a pair of scissors. To expose the EDL tendon at the knee, the biceps femoris was first removed with a pair of forceps. The EDL was removed by cutting from tendon to tendon with a pair of scissors.

Myofiber isolation
The dissected EDL was placed in a 1.5-ml Eppendorf tube with 800 l of myofiber digestion buffer. Trypsin was added to a final concentration of 0.25% to remove the associated satellite cells. The EDL was incubated at 37°C and 5% CO 2 for at least 1 h.
To disassociate the myofibers, we transferred the EDL to a 6-well plate with 2 ml of unsupplemented DMEM that had previously been coated with 10% horse serum (HS) in DMEM for at least 30 min. The EDL was gently pipetted up and down with a large bore pipette coated in HS until no more myofibers could be retrieved.

Immunofluorescence of myofibers
Briefly, freshly isolated myofibers were fixed at T 0 using 4% paraformaldehyde in PBS for 5 min. They were washed three times with 0.1% Triton X-100 in PBS, then permeabilized with 0.1% Triton X-100 and 0.1 M glycine in PBS for 15 min, and again washed three times with 0.1% Triton X-100 in PBS. Myofibers were blocked for 1 h with blocking buffer, followed by incubation with a Pax7 hybridoma primary antibody (DSHB catalog no. AB_528428) at 1:100 dilution in blocking buffer overnight at 4°C. The next day, the fibers were washed three times with 0.1% Triton X-100 in PBS, before incubating the Alexa Fluor 488 goat anti-mouse IgG1 secondary antibody (Invitrogen catalog no. A21121) at a 1:400 dilution in blocking buffer for 1 h. Fibers were washed three times with 0.1% Triton X-100 in PBS and mounted on a microscope slide with Prolong Gold antifade reagent with 4Ј,6Ј-diamino-2-phenylindole (Invitrogen catalog no. P3695).

Single-myofiber RNA extraction
Using a small bore glass pipette coated with HS, myofibers were transferred to a 6-well plate with 2 ml of PBS to wash the fibers. A single myofiber was transferred by using a coated small bore pipette to a 0.2-ml PCR tube, and the excess PBS was removed using a pipette. Next, we added 10 l of lysis buffer and gently pipetted the myofiber up and down for 3 min and then incubated the fiber on ice for 5 min while periodically vortexing and spinning down the sample. The residual fiber pieces were removed by spinning down the sample and transferring the supernatant to a fresh PCR tube.

Whole-muscle RNA extraction
The whole muscle from the hindlimb of a mouse was dissected, frozen in liquid nitrogen, and ground into a powder using a mortar and pestle. RNA from whole muscle was extracted using TRIzol reagents (Ambion catalog no. 15596018).

cDNA library preparation
cDNA was constructed using the DNA SMART-Seq HT kit and following the manufacturer's recommendations. For a single fiber, we used 12 cycles of PCR amplification on the cDNA. The cDNA was then purified using AMPure XP beads (Beckman Coulter catalog no. A63880) at a 1:1 ratio and quantified using the Quant-iT PicoGreen dsDNA assay kit (ThermoFisher catalog no. P11496).

Single-myofiber RNA-Seq
For the addition of Illumina sequencing adapters, we used 250 pg of cDNA in 1.25 l of water as a starting material and followed the directions provided with the Nextera XT DNA library preparation kit but reduced all quantities by 4ϫ volume. The libraries were size-selected using AMPure XP beads at 0.85ϫ of the sample volume, to remove all fragments below 200 bp.

Sequencing and gene expression analysis from a single fiber and whole muscle
The Illumina NextSeq 500 high-output flow cell was used for sequencing. The sequenced reads were then mapped to the mouse mm10 genome by using HISAT2 (45), using an index downloaded from the HISAT2 website that jointly indexes the mm10 genome and the ENSEMBL transcriptome definition (46).

Differential expression analysis between old and young single myofiber from mouse
RNA-Seq raw reads from young and old single myofiber were included in the analysis. The reads were mapped to the mouse genome assembly (mm10) using HISAT2 (45). The number of aligned reads per gene was obtained with HTSeq (47), using gene annotations from GENCODE M23 (48). Genes with an average read counts smaller than 10 were filtered out. Differentially expressed genes between young and old samples were identified using the R package DESeq2 (49). Log fold changes were calculated, and their associated p values were corrected by independent hypothesis weighting (50).

Gene set enrichment analysis
To obtain the pathways affected between old and young mouse fiber, we used the implementation of Gene Set Enrichment Analysis of the fgsea R package (48). We used the Log fold changes calculated by DESeq2 to create a preranked genelist. The Reactome database was used as a reference for pathways (51).