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 there are biochemical and histological assays to study myofiber heterogeneity, efficient methods to analyze the whole transcriptome of individual myofibers are lacking. We have developed single myofiber RNA-Seq (smfRNA-Seq) to analyze the whole transcriptome of individual myofibers by combining single fiber isolation with Switching Mechanisms at 5’ end of RNA Template (SMART) technology. Our method provides high-resolution genome wide expression profiles of single myofibers. Using smfRNA-Seq, we have analyzed the differences in the transcriptome of young and old myofibers to validate the effectiveness of this new method. Using smfRNA-Seq, we performed comparative gene expression analysis between single myofibers from young and old mice. Our data suggests that aging leads to significant changes in the expression of metabolic and structural genes in myofibers. Our data suggests that smfRNA-Seq is a powerful tool to study developmental, disease and age-related dynamics in the composition of skeletal muscle.


INTRODUCTION
Skeletal muscle is composed of a variety of different cell types, including endothelial cells, fibro/adipogenic cells (FAPs), 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 which 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 (MyHC) isoforms (5,7). 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 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 as more than an adequate amount of RNA was recovered, even from a single myofiber, for use with SMART-Seq technology (Table 1).
From the extracted RNA, we successfully 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  (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 ( Figure 3A).
When performing Principal Component Analysis (PCA), single fiber samples cluster together on both axes, but away from the whole muscle on the PCA2 axis ( Figure 3B). 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 ( Figure   3A,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 the whole muscle samples ( Figure 4A-C). The Myh cluster codes for a variety of myosin heavy chain proteins (MyHC), 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. ( Figure   4A). For further confirmation, we also display Ckm, the muscle specific creatine kinase (12) and ACTA1, which codes for skeletal muscle alpha actin (13). These genes have the same pattern as is seen with Myh ( Figure 4B,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 ( Figure 4D). 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 ( Figure 5A,B).
Endothelial cells are also depleted in single fiber samples, since their markers Kdr and PECAM1 are expressed in whole muscle, but not in the single fiber (16,17)( Figure 5C,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)  From this data, we performed a 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 SF, indicating that they originate from non-myogenic cells. Using these criteria, we identified 445 genes that were solely expressed in the WM muscle samples (Table S1).
Some of the more upregulated 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 SF samples, with a q value < 0.01, being expressed at more than 10 RPM and being expressed at least 10 RPM more than in the WM samples. We identified 622 genes that matched these criteria and performed a Gene Ontology (GO) analysis (Table S2) 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 downregulated in old compared to young myofibers ( Figure 7A,B).
Despite not reaching significance, there is also a moderate decrease in Acta1 with a log2fold change of -0.84. Since actins and myosins play a role in myofiber structure, force and contractile properties, we suspect that the downregulation of these genes may play a role in the weakening in contractile strength in aged muscles (28,29). A gene that is significantly upregulated in aged myofibers is Dkk3, which encodes a secreted inhibitor of WNT signaling ( Figure 7C) (30).
This finding is in line with a previous study that demonstrated that increased Dkk3 in aged muscle tissue leads to muscle atrophy through autocrine signaling(30), while also confirming that the myofiber is the source of this secreted factor.
There is also a downregulation in Ndn in the old myofibers ( Figure 7G). Ndn produces Necdin, which is important in muscle regeneration, as it promotes the differentiation of satellite cells, and Necdin knockout mice have a significant impairment in muscle regeneration(31). Necdin was also shown to protect skeletal muscle from tumor-induced muscle wasting, also called cachexia, in part by binding to and inhibiting P53 (32). 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 shows that multiple genes involved in key metabolic processes and maintenance of homeostasis are also deregulated in old compared to young myofibers. For example, we detect a loss in the expression of the long noncoding RNA H19 ( Figure 7E). While the role of H19 during myogenesis has been studied (33), 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 (34,35). We speculate that H19 reduction in old myofibers may therefore be a contributor to metabolic dysfunction, which is a common feature of aging (35).
In old myofibers, we also see an upregulation of Atf3, which is a stress factor that is known to be upregulated in muscle after a disturbance in homeostasis, such as after exercise (36) ( Figure   7D). Interestingly, Atf3 was shown to dampen the expression of pro-inflammatory cytokines after exercise (37). We speculate that its upregulation in aging myofibers may occur in response to an age-related disturbance in homeostasis, such as the stress induced by chronic low-grade inflammation(38). Another gene that is upregulated in aging is Nos1 which produces nitric oxide (NO) and plays a role in skeletal muscle metabolism ( Figure 7F). Notably, Nos1 is known to increase glucose uptake and inhibit mitochondrial respiration(39), which may be a

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

Accession numbers
All gene expression data reported in this study is available through GEO accession number GSE135364 and GSE138591

Commercial kits
The following commercial kits were used in this experiment. RNA extraction buffer is made using 19 µL of the 10X lysis buffer plus 1 µL of RNAse inhibitor from the SMART-Seq HT Kit. 1µL of the previously composed 10X lysis buffer is added to 9 µL of RNAse free water to make 1X lysis buffer

Dissection of the Extensor Digitorum Longus (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 as 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 o C and 5% CO2 for at least one hour.
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 minutes.
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 T0 using 4% paraformaldehyde in PBS for 5 minutes.

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 minutes and then incubated the fiber on ice for 5 minutes 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 #15596018).

cDNA library preparation
cDNA was constructed using the DNA SMART-Seq HT kit and following the manufacturers 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 #A63880) at a 1:1 ratio and quantified using the Quant-iT PicoGreen dsDNA Assay kit (ThermoFisher #P11496) 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 4x volume. The libraries were size selected using AMPure XP beads at 0.85X 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(42), using an index downloaded from the HISAT2 website that jointly indexes the mm10 genome and the ENSEMBL transcriptome definition. (43).

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. Reads were mapped to the mouse genome assembly (mm10) using HISAT2 (42). The number of aligned reads per gene was obtained with HTSeq (44), using gene annotations from GENCODE M23 (45). 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 (46). Log fold changes were calculated, and their associated P-values were corrected by independent hypothesis weighting (47).

Gene Set Enrichment Analysis
To obtain the pathways affected between old and young mouse fiber, we used the implementation of    S1A: List of genes found in the whole muscle but not expressed in myofibers, defined as having an RPM value of at least 10 in the whole muscle and 0 in the single fiber.
Table S1B: GO term analysis of genes not expressed in the myofiber but found in the whole muscle, defined as having an RPM value of at least 10 in the whole muscle and 0 in the single fiber. Used DAVID 6.8 online functional annotation tool with the ENSEMBL_GENE_ID identifier and with a classification stringency of medium. Table S2A: List of genes expressed exclusively in the myofiber, defined as genes with a q value lower than 0.01 between single fiber and whole muscle, more highly expressed in the single fiber with a Log2 fold change of less than -1, more highly expressed than 10 RPM and with a difference in expression between single fiber and whole muscle of at least 10 RPM. Table S2B: Go term analysis of genes exclusively expressed in the single myofiber, defined as genes with a q value lower than 0.01 between single fiber and whole muscle, more highly expressed in the single fiber with a Log2 fold change of less than -1, more highly expressed than 10 RPM and with a difference in expression between single fiber and whole muscle of at least 10 RPM. Used DAVID 6.8 online functional annotation tool with the ENSEMBL_GENE_ID identifier and with a classification stringency of medium. Table S3: Alignment statistics of sequenced reads from young and old single myofibers.             Regulation of APC/C activators between G1/S and early anaphase