A transcriptome atlas of leg muscles from healthy human volunteers reveals molecular and cellular signatures associated with muscle location
- Department of Human Genetics, Leiden University Medical Center, Netherlands
- Division of Cell and Chemical Biology, Leiden University Medical Center, Netherlands
- C.J. Gorter MRI Center, Department of Radiology, Leiden University Medical Center, Netherlands
- Sequencing Analysis Support Core, Leiden University Medical Center, Netherlands
- Orthopedic and Sport Medicine Department, Erasmus MC, University Medical Center Rotterdam, Netherlands
- Orthopedics, Medisch Centrum Haaglanden, Netherlands
- Centre for Molecular and Biomolecular Informatics, Radboud Institute for Molecular Life Sciences, Radboud University Medical Center, Netherlands
- Duchenne Center Netherlands, Netherlands
Figures
Figure 1 with 4 supplements
An overview of biopsies’ location and the study workflow.
(A) A schematic overview of the leg muscles. Arrows point to the muscles that were included in this study. The biopsies, with exception of STM (semitendinosus-middle), were taken from the distal area. (B–D) The study overview includes cryosectioning, RNA-isolation and sequencing (B) data analysis (C) and validations (D). Created with BioRender.
Figure 1—figure supplement 1
Analysis framework.
A flowchart summarizing the analysis framework used to detect molecular signatures characterizing distinct skeletal muscles. Following pre-processing, muscle-specific signatures were identified using three approaches: differential transcript usage analysis (in purple), differential expression analysis (in blue), cell type composition analysis (in yellow), gene co-expression network analysis (in green), and functional enrichment analysis (in red).
Figure 1—figure supplement 2
RNA isolation protocols.
The PCA plot shows the effect of the RNA isolation protocols. The scatter plot of PC1 (x-axis) and PC2 (y-axis) shows GR muscle from five individuals from which RNA was isolated using two RNA isolation protocols.
Figure 1—figure supplement 3
The overview of RNA-seq samples.
An X indicates samples that are not present in the final transcriptome dataset. Fastq files available from the European Genome Archive, Dataset ID: EGAS00001005904.
Figure 1—figure supplement 4
The quality control and batch correction of RNA-seq data.
(A–B) Scatter plots of PC1 (x-axis) and PC2 (y-axis) before (A) and after (B) batch correction. Each dot presents a sample labeled by muscle tissue. Each library preparation batch is shown with a different color. The re-sequenced batch is denoted in black. (C–D) Box plots show the analysis of variance before (C) and after (D) batch correction. Y-axis shows the percentage of variance explained by different factors. The x-axis shows the known biological (muscle and individual, shown in green) and technical (RIN score, concentration, batch, library size, shown in red) factors. The RNA isolation protocol effect is captured in the individual effect.
Figure 2 with 2 supplements
Muscles cluster into three main groups based on cell type composition.
(A) The heatmap shows the mean eigenvalues of genes marking each cell type across all the individuals. Each row shows a muscle, and each column shows a cell type. FAP stands for fibro-adipogenic progenitors. (B) The boxplot shows the eigenvalues for the endothelial cells, fast-twitch, and slow-twitch myofibers per muscle. The boxes reflect the median and interquartile range.
Figure 2—figure supplement 1
Cell type composition pairwise comparison between muscles.
The heatmap shows the differences between each pair of muscles. The statistically significant variation between muscles was tested by ANOVA followed by the post-hoc pairwise comparisons. Each row corresponds to a pairwise comparison, and each column shows a cell type. Color-coded cells show the corresponding t-ratio for the differences in eigenvalue of a cell type in each pairwise comparison. The significant differences (Tukey p-value <0.05) are colored red (significantly higher eigenvalues in muscle M1) or blue (significantly higher eigenvalues in muscle M2). The non-significant differences are colored from pink (relatively higher eigenvalues in M1) to light blue (relatively higher eigenvalues in M1).
Figure 2—figure supplement 2
Cell types’ eigenvalues.
Each boxplot shows the eigenvalues of a cell types across different muscles. The boxes reflect the median and interquartile range.
Figure 3 with 1 supplement
Myofiber type composition is consistent with the expression level of genes marking fast and slow-twitch myofibers.
(A) A representative immunostaining image. The overlay of each myosin heavy chain isoform and laminin are shown separately. Scale bar is 1000 μm. (B) A 3-D scatterplot of the MyHC isoforms MFI. Each dot represents a myofiber. Myofibers in the three largest clusters are denoted with red (Cluster 1), blue (Cluster 2), and green (Cluster 3). The objects with low MFI values for all the isoforms are denoted in yellow (Cluster 4, ~2% of all the dots). In gray are ~4% of myofibers assigned to many small clusters. (C) Density plots show MFI distribution for each MyHC isoform in the three largest clusters. MFI values are scaled (without centering) and transformed. (D–F) Scatterplots show the proportion of the assigned myofibers to each of the largest clusters and the normalized expression of the gene coding the isoform with a relatively higher expression in that specific myofiber cluster. (G) The boxplot shows the proportion of myofibers in the three largest clusters per muscle. Each muscle is depicted with a different color, with G1 muscles in blue, G2 muscles in red and the G3 muscle in gray. The boxes reflect the median and interquartile range.
Figure 3—figure supplement 1
The association of myofiber clusters with MyHC expression, and myofiber composition differences between muscles.
(A) Scatterplots show correlation between the proportion of the assigned myofibers per cluster and MyHC isoform normalized expression levels. MYH2 codes for MyHC2A, MYH7 codes for MyHC1 and MYH1 codes for MyHC2X. (B) The heatmap shows the differences between each pair of muscles. The statistically significant variation between muscles was tested by ANOVA followed by the post-hoc pairwise comparisons. Each row corresponds to a myofiber cluster, and each column shows a pairwise comparison. Color-coded cells show the corresponding t-ratio for the differences in proportions of myofiber in each pairwise comparison. The significant differences (Tukey p-value <0.05) are colored red (significantly higher proportions in muscle M1) or blue (significantly higher proportions in muscle M2). The non-significant differences are colored from pink (relatively higher proportions in M1) to light blue (relatively higher proportions in M1).
Figure 4
Immunostaining confirms higher capillary density in GL compared with STM muscles.
A) Representative images of STM and GL cross-sections immunostained with CD31 (green), ENG (red), and laminin (white). An enlargement of the boxed region is shown on the right: merged and separate channels (CD31 and ENG). Examples of objects recognized as capillaries are depicted by yellow arrows. Examples of objects that were not considered as capillaries due to circularity values are shown by red arrows. Scale bar 1000 μm. (B) The box plot shows the percentage of CD31 positive area in the two muscles. (C) The box plot shows the normalized expression of CD31 gene in the two muscles. (D) The boxplot shows the estimated capillary density in the two muscles. The capillary density was defined as the number of objects (3–51 µm2) with an overlap between CD31 and ENG per unit cross-sectional area of the muscle. The boxes reflect the median and interquartile range (N=19 per muscle). The red dots on the boxes show the mean. **** p-value <1 × 10–6 (linear mixed-model).
Figure 5 with 3 supplements
Gene expression differences between three groups of muscles not driven by cell type composition.
A) Symmetric heatmap plot shows the percentage of DEGs in different pairwise comparisons. Genes with a high Pearson correlation (R>0.5) with the eigenvector of any cell type (cell type related genes) are excluded. (B) Symmetric heatmap plots show the percentage of DETs in pairwise comparisons. Transcripts originating from cell type related genes are excluded. (C) as in (B), but for the number of genes having at least a significant transcripts usage difference. (D) Symmetric heatmap plot shows the number of modules that were not driven by cell type composition and were significantly different in each pairwise comparison. Each row or column in A–D) represents a muscle. (E) The heatmap shows the modules that reflect the intrinsic differences between groups of muscles. Each row represents a muscle, and each column shows a muscle-related module that was not driven by cell type composition. Color-coded cells show the corresponding average of eigenvalues across all individuals (N=20). Modules with an overall higher expression in G1 or G3 are underlined by a blue or gray dashed line, respectively. The red dashed line underlines the modules with an overall higher expression in both G2 and G3. The black asterisks show modules with the largest differences between the three groups of muscles.
Figure 5—figure supplement 1
DEA and WGCNA also clustered muscles in three groups.
(A) Symmetric heatmap shows the proportion of all differentially expressed genes in different pairwise comparisons. Each row or column represents a muscle. (B) As in A, but here for the proportion of all differentially expressed transcripts. (C) As in (A) but here for the number of genes with at least a significant transcript usage difference. (D) as in A, but here for the number of modules that were significantly different in each pairwise comparison.
Figure 5—figure supplement 2
Overview of genes with significant transcript usage difference.
The heatmap shows the number of transcripts with significant usage differences in each pair of muscles. The genes with at least one significant transcript usage difference are listed in rows. The pairwise comparisons with at least one transcript usage difference are shown in columns. Red asterisks highlight genes with a high Pearson correlation (R>0.5) with the eigenvector of any cell type.
Figure 5—figure supplement 3
WGCNA modules.
A) Each row corresponds to a pairwise comparison and each column shows a muscle-related. Color-coded cells show the corresponding t-ratio for the differences in eigenvalue of a module in each pairwise comparison. The significant differences (Tukey p-value <0.05) are colored red (significantly higher eigenvalues in M1) or blue (significantly higher eigenvalues in M2). The insignificant differences are colored from pink (relatively higher eigenvalues in M1) to light blue (relatively higher eigenvalues in M1). (B) The table shows the intersection of the genes in the modules (columns) with genes marking different cell types (rows). Color-coded cells show the corresponding intersection number. The red asterisks show modules driven by cell type composition.
Figure 6
A schematic representation of genes with higher expression in G2 and G3, related to oxidative phosphorylation and metabolic pathways in the mitochondria.
60 (out of the 122) mitochondrial genes with higher expression in G2 and G3 are shown in red. The electron transport chain, lysin and tryptophan catabolism, TCA cycle, and beta-oxidation are shown. The hub genes are underlined and in bold. Created with BioRender.
Figure 7
The expression patterns of HOX genes cluster muscles in the same groups.
A) The graph shows the co-expression subnetwork of HOX genes and genes related to anterior/posterior pattern specification assigned to the M.14 module. Diamonds indicate transcription factors while other genes are indicated by circles. Pink and purple nodes represent the hub genes and non-hub genes, respectively. The genes related to anterior/posterior pattern specification have a black border. The edge thickness reflects the degree of topological overlap. Topological overlap is defined as a similarity measure between each pair of genes in relation to all other genes in the network. High topological overlaps indicate that genes share the same neighbors in the co-expression network. (B) Normalized expression of all HOX genes (scaled by row) represented as a heatmap. The hierarchical clustering was generated using the normalized expression values. Each row represents a gene and each column represents a sample. The side color of columns indicates different muscles. The module in which the gene assigned is given between parentheses. Eleven highlighted HOX genes are assigned into muscle-related modules which showed the largest differences between the groups of muscles (M.14, M.30, and M.32).
Figure 8 with 1 supplement
Distinct expression of HOX genes confirmed by RNAscope.
(A) Representative in situ hybridization images of HOXC10 and HOXA10 in cryosections of STM and GL. The image is a merge image of the channels used for laminin, DAPI, HOXC10 and HOXA10 staining. Scale bar is 100 μm. (B) The boxplots show the average number of foci per myofiber (y-axis) in STM and GL muscles (x-axis). (C) The boxplots show the normalized expression of HOXA10 (top) and HOXC10 (bottom) in STM and GL muscles. The boxes reflect the median and interquartile range (N = 12 per muscle). The red dots on the boxes show the mean. **** p-value < 1 × 10–6 (linear mixed-model).
Figure 8—figure supplement 1
A representative image of negative control probes in the RNAscope experiment in STM muscle cryosections.
Scale bar is 100 μm.
Tables
Table 1
Top enrichment results of muscle-related modules not driven by cell type composition.
| Module | Term | FDR | #Enriched genes | Hub genes |
|---|---|---|---|---|
| Higher expression in G1 (GR, STM, and STD) | ||||
| M.30 (75 genes) | IRE1-mediated unfolded protein response | 4×10–4 | 4 | |
| M.32 (236 genes) | Hormone-mediated signaling pathway | 9×10–3 | 11 | CARM1, WBP2, ZBTB7A |
| M.7 (308 genes) | Negative regulation of nucleobase-containing compound metabolic process | 6×10–4 | 51 | CREBBP, DAB2IP, FOXK2, LARP1, RXRA, THRA |
| Chromatin organization | 1×10–2 | 27 | ARID1B, CREBBP, HUWE1 | |
| M.17 (176 genes) | Chromatin modifying enzymes | 4×10–3 | 11 | HCFC1, SETD1A |
| Higher expression in G2 (RF, VL, and VM) and G3 (GL) | ||||
| M.31 (945 genes) | RNA splicing | 6×10–5 | 54 | SNRNP70 |
| Histone modification | 3×10–3 | 47 | KAT2A | |
| M.33 (538 genes) | Apical junction complex | 2.4×10–2 | 11 | MICALL2 |
| M.13 (300 genes) | Mitochondrion | 4×10–45 | 122 | AIFM1, ATP5F1A, ATP5F1B, CKMT2, COQ9, DLD, DLST, FH, GHITM, HADHA, HADHB, IMMT, MFN2, NDUFS2, PDHA1, PDHB, TRAP1, UQCRC2 |
| M.14 (162 genes) | Anterior/posterior pattern specification | 4×10–3 | 7 | HOXA11 |
| Higher expression in G3 (GL) | ||||
| M.5 (190 genes) | Regulation of lipid metabolic process | 8×10–4 | 15 | ADIPOQ, ADRA2A, CIDEA, LEP, LGALS12, PDE3B, SCD |
| M.25 (188 genes) | Ameboidal-type cell migration | 5×10–3 | 15 | CFL1, PML, TGFB1 |
| Positive regulation of muscle cell differentiation | 9×10–3 | 6 | EHD2, ENG, NIBAN2, TGFB1 | |
| M.11 (299 genes) | Golgi membrane | 2×10–3 | 30 | ASAP2, MAN1A1 |
| Regulation of nervous system development | 8×10–3 | 30 | IQGAP1 | |
Additional files
-
Supplementary file 1
Samples’ metadata.
- https://cdn.elifesciences.org/articles/80500/elife-80500-supp1-v2.xls
-
Supplementary file 2
A combined list of the genes marking different cell types.
- https://cdn.elifesciences.org/articles/80500/elife-80500-supp2-v2.xls
-
Supplementary file 3
The result of differential expression analysis for genes.
- https://cdn.elifesciences.org/articles/80500/elife-80500-supp3-v2.xls
-
Supplementary file 4
The result of differential expression analysis for transcripts.
- https://cdn.elifesciences.org/articles/80500/elife-80500-supp4-v2.xls
-
Supplementary file 5
The result of differential transcript usage analysis.
- https://cdn.elifesciences.org/articles/80500/elife-80500-supp5-v2.xls
-
Supplementary file 6
List of genes – modules.
- https://cdn.elifesciences.org/articles/80500/elife-80500-supp6-v2.xls
-
Supplementary file 7
List of enriched biological processes and molecular functions within modules.
- https://cdn.elifesciences.org/articles/80500/elife-80500-supp7-v2.xls
-
MDAR checklist
- https://cdn.elifesciences.org/articles/80500/elife-80500-mdarchecklist1-v2.docx
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
A transcriptome atlas of leg muscles from healthy human volunteers reveals molecular and cellular signatures associated with muscle location
eLife 12:e80500.
https://doi.org/10.7554/eLife.80500
