Single-cell spatial transcriptomics reveals a dystrophic trajectory following a developmental bifurcation of myoblast cell fates in facioscapulohumeral muscular dystrophy - PubMed (original) (raw)

Single-cell spatial transcriptomics reveals a dystrophic trajectory following a developmental bifurcation of myoblast cell fates in facioscapulohumeral muscular dystrophy

Lujia Chen et al. Genome Res. 2024.

Abstract

Facioscapulohumeral muscular dystrophy (FSHD) is linked to abnormal derepression of the transcription activator DUX4. This effect is localized to a low percentage of cells, requiring single-cell analysis. However, single-cell/nucleus RNA-seq cannot fully capture the transcriptome of multinucleated large myotubes. To circumvent these issues, we use multiplexed error-robust fluorescent in situ hybridization (MERFISH) spatial transcriptomics that allows profiling of RNA transcripts at a subcellular resolution. We simultaneously examined spatial distributions of 140 genes, including 24 direct DUX4 targets, in in vitro differentiated myotubes and unfused mononuclear cells (MNCs) of control, isogenic D4Z4 contraction mutant and FSHD patient samples, as well as the individual nuclei within them. We find myocyte nuclei segregate into two clusters defined by the expression of DUX4 target genes, which is exclusively found in patient/mutant nuclei, whereas MNCs cluster based on developmental states. Patient/mutant myotubes are found in "FSHD-hi" and "FSHD-lo" states with the former signified by high DUX4 target expression and decreased muscle gene expression. Pseudotime analyses reveal a clear bifurcation of myoblast differentiation into control and FSHD-hi myotube branches, with variable numbers of DUX4 target-expressing nuclei found in multinucleated FSHD-hi myotubes. Gene coexpression modules related to extracellular matrix and stress gene ontologies are significantly altered in patient/mutant myotubes compared with the control. We also identify distinct subpathways within the DUX4 gene network that may differentially contribute to the disease transcriptomic phenotype. Taken together, our MERFISH-based study provides effective gene network profiling of multinucleated cells and identifies FSHD-induced transcriptomic alterations during myoblast differentiation.

PubMed Disclaimer

Figures

Figure 1.

Spatial transcriptomic profiling of differentiated control and FSHD patient/mutant myoblasts by MERFISH. (A) Schematic diagram of the sample preparation process for MERFISH experiments and the genotypes used in the study. (B) Spatial distribution of detected genes by MERFISH. (Column 1) Illustration of the original sample image (cropped). All 140 genes used in the study are illustrated. (Column 2) The magnified example region corresponds to the red square area on the left. (Columns 3–8) Distribution of six example genes in the magnified area. (C) Staining results of the magnified region in B. (Top) DAPI staining of the nuclei in the region; (bottom) cell boundary staining of the region that illustrates the cytoplasm. (D) Automatically segmented mononuclear cells (MNCs) and nuclei inside the magnified region of B. The background is the overlay of cell boundary staining and DAPI staining. (Top) Detected MNCs inside the area. MNCs are circled with red lines. (Bottom) Detected nuclei within the region. Nuclei are circled with red lines. (E) Manually labeled myotube and nonmyotube regions (top), and the corresponding spatial transcriptome illustration with all 140 genes (bottom) in the magnified region in B. (F) Composition of nuclei, MNCs, myotube, and nonmyotube of the different genotypes in this study. (G_–_I) Comparison of MERFISH results with bulk RNA-seq of control myoblasts at day 3 of differentiation. Comparisons between the log-transformed total transcript counts from MERFISH and total transcript counts of all myocytes from bulk RNA-seq (G), intra-myotube transcript counts from MERFISH and total counts from bulk RNA-seq (H), and intra-MNC transcript counts from MERFISH and total counts from bulk RNA-seq (I). Overall MERFISH detections show a high level of correlation with the detection results from bulk RNA-seq and scRNA-seq methods (G) Pearson's r = 0.66, P < 0.001; (H) Pearson's r = 0.78, P < 0.001; (I) Pearson's r = 0.57, P < 0.001; N = 140 genes.

Figure 2.

Nuclei clusters show differential expression of DUX4 target genes detected by MERFISH. (A) Two-dimensional projections (UMAP) of gene expression profiles of all identified nuclei. (Left) UMAPs overlaid with the corresponding cluster identities identified by scrattch.hicat algorithm; (right) proportion of each cluster in each genotype. (B) Heatmap illustrating DUX4 target gene expression level of the nuclei clusters. The expression profile for each nucleus cluster is the averaged, area-normalized, log-transformed expression across different samples and genotypes. (C) Volcano plot illustrates the differential expression between the two identified nuclei clusters across all samples. About 88% of DUX4 target genes show upregulation in cluster 1 nuclei, whereas one DUX target gene (ZNF596) shows upregulation in cluster 0 nuclei (Log2FC threshold: 0.5; _P_-value threshold: 0.001). (D) Examples of DAPI and cell boundary staining overlaid with nuclei, colored with their corresponding cluster identities. Annotated myotubes are circled in cyan. (E) Total number of intra-myotube and intra-MNC cluster 1 nuclei for each FSHD and DEL5 sample. No significant difference was noted in counts between intra-myotube and intra-MNC cluster 1 nuclei (intra-myotube: 551.0 ± 240.7; intra-MNC: 206.0 ± 66.1; mean ± SE; P = 0.3939, Mann–Whitney U test; N = 6 samples). (F) Proportion of intra-myotube and intra-MNC nuclei to be cluster 1 for each FSHD and DEL5 sample. Intra-myotube nuclei contain a higher percentage of cluster 1 nuclei compared with intra-MNC nuclei for all samples (intra-myotube: 0.3570 ± 0.1130; intra-MNC: 0.0316 ± 0.0098; mean ± SE; P = 0.0087, Mann–Whitney U test; N = 6 samples). (G) UMAP of gene expression profiles of intra-myotube nuclei. UMAPs are overlaid with the clusters identified by scrattch.hicat algorithm. (H) Distribution of intra-myotube nuclei counts for different genotypes. Overall, control myotubes show a significantly higher nuclei count compared with FSHD, whereas DEL5 myotubes do not show significantly different nuclei counts compared with FSHD (370 FSHD myotubes: 15.47 ± 0.60; 213 DEL5 myotubes: 24.60 ± 1.29; 404 control myotubes: 27.42 ± 1.07; mean ± SE; control vs. FSHD: P = 0.0092; DEL5 vs. FSHD: P = 0.0996; linear mixed effects model). (I) Distribution of myotube volume by genotype. Overall, myotubes from control and DEL5 samples show higher volume than those in FSHD (370 FSHD myotubes: 11499 ± 381.9 µm2; 213 DEL5 myotubes: 14652 ± 686 µm2; 404 control myotubes: 15768.0 ± 483.9 µm2; mean ± SE; DEL5 vs. FSHD: P = 0. 0174; control vs. FSHD: P = 2.4655 × 10−6; linear mixed effects model). (J) Comparison of intra-myotube cluster 1 nuclei between FSHD and DEL5 genotype myotubes. Distribution of counts of intra-myotube cluster 1 nuclei per myotube. Only myotubes with cluster 1 nuclei are included. No significant differences are noted between the two genotypes (234 FSHD myotubes: 12.54 ± 0.61; 29 DEL5 myotubes: 11.31 ± 1.38; mean ± SE; FSHD vs. DEL5: P = 0.7145, linear mixed effects model). (K) Proportion of cluster 1 nuclei per myotube. The percentage is calculated as the ratio of cluster 1 nuclei counts to total nuclei counts. Only myotubes with cluster 1 nuclei are included. FSHD myotubes show a significantly higher proportion of cluster 1 intra-myotube nuclei compared with DEL5 (234 FSHD myotubes: 0.7470 ± 0.0150; 29 DEL5 myotubes: 0.5705 ± 0.0535; mean ± SE; FSHD vs. DEL5: P = 0.0011, linear mixed effects model).

Figure 3.

MNC clusters with differential expression of myogenesis and muscle function genes identified by MERFISH. (A) Two-dimensional UMAP projection of gene expression profiles of all MNCs. (Left) UMAP overlaid with the corresponding cluster identities; (right) Composition of each cluster in each genotype. (B) Heatmap illustrating DUX4 target gene expression level of the MNC clusters. The expression profile for each MNC cluster is the averaged, normalized, log-transformed expression across different samples and genotypes. (C) Volcano plot visualizes differential expression between the two identified MNC clusters across all samples. Myoblast and myotube marker genes (TTN, MYH3, MYH8, and NEB) show upregulation in cluster B MNCs (Log2FC threshold: 0.5; _P_-value threshold: 0.001). (D) Comparison of expression profiles of the two MNC clusters and myotube populations. (Top) UMAP calculated from the expressions of both MNCs and myotubes from all samples. The myotube population mostly overlaps with cluster 1 MNCs in UMAP space. Each dot represents a single MNC or a myotube. (Bottom) UMAP colored by genotype. Each dot represents a single MNC or myotube. (E) UMAP as in D, but overlaid with the normalized expression level of the myotube target genes, as well as two example DUX4 target genes (CCNA1, ZSCAN4). Cells with higher expression levels of TTN, MYH3, MYH8, and NEB mainly colocalize with cluster B MNCs, further indicating that cluster B MNCs are at a more mature differentiation state. CCNA1 and ZSCAN4 are highly expressed in the area occupied by FSHD and DEL5 myotubes.

Figure 4.

Clustering and pseudotime analysis identify a developmental bifurcation in the myoblast population, associated with a distinct reduction in myotube volume. (A) UMAP plots of MERFISH data from 987 myotubes and 629 nonmyotube regions. The plots are colored according to the clusters determined by the shared nearest neighbor (SNN) algorithm after batch correction using the Harmony algorithm. (B) UMAP from A colored by cell types, including the myotubes and nonmyotube regions of the differentiated control, FSHD and DEL5 cells. (C) Bar plot of the proportion of cell types in each cluster from the UMAP, colored by annotated cell types as in B. The myotubes for FSHD-hi and FSHD-lo are determined by DUX4 target gene expression levels. (D) Average expression profiles of the DUX4 target genes in each cluster. (E) The percentage of FSHD-hi and FSHD-lo cells in FSHD and DEL5 myotubes. (F) UMAP from A colored by designation of FSHD-hi or FSHD-lo versus control myotubes and remaining cells. (G) UMAP from A colored by normalized and scaled expression of example marker genes. (H) Trajectory analysis of batch-corrected MERFISH data using Monocle. The top 18 genes with the highest variability in expression within the myotube and nonmyotube regions are used to construct the pseudotime tree. Myotubes or nonmyotube regions on the tree are colored by pseudotime. (I) The pseudotime tree from H is colored by cell types and displayed individually. (J) Bar plot of the distribution of two types of myotubes, FSHD-hi and FSHD-lo, based on the percentage of cluster 1 nuclei in each myotube. The _x_-axis represents the range of cluster 1 nuclei percentages. The majority of FSHD-lo myotubes have no cluster 1 nuclei, with a few exceptions. To highlight these rare occurrences, the FSHD-lo myotube numbers (1 or 2) are indicated on the bars with red arrows. (K) The number of nuclei per myotube quantified in the FSHD-hi and FSHD-lo groups of FSHD myotubes. (L) The pseudotime heatmap displays the expression level of significantly changed genes (P < 0.01) in the two branches for comparison. Color key from blue to red indicates relative expression levels from low to high.

Figure 5.

Weighted gene coexpression network analysis (WGCNA) of the non-DUX4-targets reveals gene network disruptions in FSHD and DEL5 myotubes. (A) Heatmaps depict the topological overlap matrix (TOM) among the 77 non-DUX4-target genes in the analysis. This matrix indicates similarity between genes by analyzing the extent to which individual genes coexpress with the same subset of genes (including each other). This results in a more robust similarity measure than pure correlation matrices. Lighter colors represent higher overlap, and darker colors represent lower overlap. The three heatmaps are plotted with the same order of genes based on the gene dendrogram of control myotubes. The modules of control are shown along the left side and the top. (B) The bubble plot shows Gene Ontology (GO) enrichment analysis of genes in the modules of control using the online tool g:Profiler. The bubble colors correspond to the module colors in A. The _x_-axis displays log10-adjusted _P_-values used Bonferroni correction and a threshold of 0.05. Bubble size indicates, in each term, the percentage of total genes in the given module. (C) The intramodular _k_-values of the genes in the blue and turquoise modules, identified in control myotubes, are decreased significantly in the FSHD and DEL5 myotubes. (D) The intramodular and intermodular connections of genes change significantly in FSHD and DEL5 myotubes compared with the control. The width of the edges and self-loops is proportional to sum of the weights (topological overlap values) of the edges, thresholded at 0.1. The module size indicates the number of genes within it. Some gene pairs are shown on the right panel as examples.

Figure 6.

Differential correlation between DUX4 target genes in FSHD and DEL5 myotubes reveals network hubs showing differential association with non-DUX4 targets. (A) Pearson correlation analysis of DUX4 target genes in FSHD and DEL5 myotubes. Correlation between genes is determined by normalized and batch-corrected expression values of three DUX4 target genes versus all 20 DUX4 target genes. The top panel shows data for FSHD myotubes, and the bottom panel shows data for DEL5 myotubes. Some specific results are highlighted, with red indicating a very strong positive relationship, green indicating a very weak relationship, and gray indicating a _P_-value > 0.05. (B) The expression values of 20 DUX4 target genes (enclosed in the block box) and 14 non-DUX4-target genes in FSHD and DEL5 myotubes are analyzed by WGCNA. Red edges indicate gene pairs that show strong coexpression in both cell lines, with weights greater than 0.4 and mean values exceeding 0.5. Dark gray and light gray edges represent positive correlations between genes that do not meet the threshold of the red edges. Dark gray edges indicate that the correlation values between FSHD and DEL5 are not significantly different. Blue edges indicate negative correlations between genes. The thickness of each edge represents the average absolute weight of the corresponding gene pair in both cell lines. For the correlations between DUX4 target genes, only weights greater than 0.2 are displayed, except for DUXA, ZNF296, and PRAMEF19. These three DUX4 target genes have a lower threshold of 0.06 and are presented with their largest mean weight. For the correlations between DUX4 target genes and non-DUX4-target genes, the network represents positive correlations with weights above 0.06 or negative correlations above 0.025 in both cell lines. The correlations between non-DUX4-target genes are not shown. (C) Three subfigures extracted from B. The LEUTX-CCNA1-KHDC1L groups, ZSCAN4-KDM4E-RFPL4B groups, and six genes related to muscle function or development are classified into three groups, respectively, with the width of the edges proportional to the sum of the weights of the edges between two individual genes in B. To highlight the differences in edge weights, the scale used in these subfigures differs from that of B. The red self-loop is not scaled in the same manner as the other edges. In the top left panel, the two DUX4 target gene groups and other DUX4 target genes connected to them are displayed. SMCHD1 and DUXA are also included in this panel. The top right panel presents seven non-DUX4-target genes along with the DUX4 target genes/groups they are connected to by dark edges. The bottom left panel displays all three gene groups and shows the DUX4 target genes that are negatively correlated with the muscle genes group. The connections between DUX4 target individual genes and the groups and in top right and bottom left panels are represented by very light edges. For all three panels, the edges between individual DUX4 target genes are not shown, except for DUXA in the top left panel.

References

1. Athwal VS, Pritchett J, Martin K, Llewellyn J, Scott J, Harvey E, Zaitoun AM, Mullan AF, Zeef LAH, Friedman SL, et al. 2018. SOX9 regulated matrix proteins are increased in patients serum and correlate with severity of liver fibrosis. Sci Rep 8: 17905. 10.1038/s41598-018-36037-4 -DOI -PMC -PubMed
1. Banerji CRS, Panamarova M, Hebaishi H, White RB, Relaix F, Severini S, Zammit PS. 2017. PAX7 target genes are globally repressed in facioscapulohumeral muscular dystrophy skeletal muscle. Nat Commun 8: 2152. 10.1038/s41467-017-01200-4 -DOI -PMC -PubMed
1. Benedetti R, Romeo MA, Arena A, Gilardini Montani MS, Di Renzo L, D'Orazi G, Cirone M. 2022. ATF6 prevents DNA damage and cell death in colon cancer cells undergoing ER stress. Cell Death Discov 8: 295. 10.1038/s41420-022-01085-3 -DOI -PMC -PubMed
1. Blondel VD, Guillaume J-L, Lambiotte R, Lefebvre E. 2008. Fast unfolding of communities in large networks. J. Stat Mech Theory Exp. 2008: P10008. 10.1088/1742-5468/2008/10/P10008 -DOI
1. Bolcato-Bellemin AL, Lefebvre O, Arnold C, Sorokin L, Miner JH, Kedinger M, Simon-Assmann P. 2003. Laminin α5 chain is required for intestinal smooth muscle development. Dev Biol 260: 376–390. 10.1016/S0012-1606(03)00254-9 -DOI -PubMed

Single-cell spatial transcriptomics reveals a dystrophic trajectory following a developmental bifurcation of myoblast cell fates in facioscapulohumeral muscular dystrophy - PubMed (original) (raw)