Facioscapulohumeral dystrophy weakened sarcomeric contractility is mimicked in induced pluripotent stem cells-derived innervated muscle fibres

Abstract

Background: Facioscapulohumeral dystrophy (FSHD) is a late-onset autosomal dominant form of muscular dystrophy involving specific groups of muscles with variable weakness that precedes inflammatory response, fat infiltration, and muscle atrophy. As there is currently no cure for this disease, understanding and modelling the typical muscle weakness in FSHD remains a major milestone towards deciphering the disease pathogenesis as it will pave the way to therapeutic strategies aimed at correcting the functional muscular defect in patients.

Methods: To gain further insights into the specificity of the muscle alteration in this disease, we derived induced pluripotent stem cells from patients affected with Types 1 and 2 FSHD but also from patients affected with Bosma arhinia and microphthalmia. We differentiated these cells into contractile innervated muscle fibres and analysed their transcriptome by RNA Seq in comparison with cells derived from healthy donors. To uncover biological pathways altered in the disease, we applied MOGAMUN, a multi-objective genetic algorithm that integrates multiplex complex networks of biological interactions (protein-protein interactions, co-expression, and biological pathways) and RNA Seq expression data to identify active modules.

Results: We identified 132 differentially expressed genes that are specific to FSHD cells (false discovery rate < 0.05). In FSHD, the vast majority of active modules retrieved with MOGAMUN converges towards a decreased expression of genes encoding proteins involved in sarcomere organization (P value 2.63e^-12 ), actin cytoskeleton (P value 9.4e^-5 ), myofibril (P value 2.19e^-12 ), actin-myosin sliding, and calcium handling (with P values ranging from 7.9e^-35 to 7.9e^-21 ). Combined with in vivo validations and functional investigations, our data emphasize a reduction in fibre contraction (P value < 0.0001) indicating that the muscle weakness that is typical of FSHD clinical spectrum might be associated with dysfunction of calcium release (P value < 0.0001), actin-myosin interactions, motor activity, mechano-transduction, and dysfunctional sarcomere contractility.

Conclusions: Identification of biomarkers of FSHD muscle remain critical for understanding the process leading to the pathology but also for the definition of readouts to be used for drug design, outcome measures, and monitoring of therapies. The different pathways identified through a system biology approach have been largely overlooked in the disease. Overall, our work opens new perspectives in the definition of biomarkers able to define the muscle alteration but also in the development of novel strategies to improve muscle function as it provides functional parameters for active molecule screening.

Keywords: Facioscapulohumeral dystrophy; Induced pluripotent stem cells; Muscle contraction; Muscle weakening; Pathophysiology; Sarcomere; System biology.

Figure 1

Bosma arhinia and microphthalmia (BAMS) and facioscapulohumeral dystrophy (FSHD) patient's cells progenies express the pathogenic DUX4‐fl transcript and DUX4 target genes. (A) Schematic overview of the samples used and steps of analysis. Induced pluripotent stem cells (hiPSCs) were derived from primary fibroblasts from healthy donors, patients affected with Type 1 FSHD (one FSHD1 patient carrying 7 D4Z4 repeated units (12759) and one patient with mosaicism [25% of the cells carry a 4qA allele with 2 D4Z4 repeated units; 75% of cells carry a 4qA allele with 15 D4Z4 repeated units (17706)], FSHD2 (14586; c.573A>C; p.Q193P), BAMS (BAMS‐1; c.407A>G, p.E136G) (Table S1). Fibroblasts were reprogrammed into human induced pluripotent stem cells (hiPSCs), for FSHD1 mosaic cells, hiPSCs clones carrying the contracted D4Z4 or the healthy D4Z4 allele were isolated separately. Cells were described in Dion et al. All hiPSCs clones were differentiated into functional muscle fibres (MFs). Gene expression analysis was performed by high‐throughput RNA‐sequencing. For each condition, two biological replicates corresponding to two independent differentiation experiments were used. DEGs were selected based on a minimum two‐fold change and statistical significance of FDR, compared with control cells. In complement to classical RNA‐Seq pipeline analysis, we applied MOGAMUN, a recently developed algorithm aimed at revealing active modules in multiplex biological networks. (B) Representative Z stacks of muscle fibres stained for Titin (TTN, red), neurofilament (NF, green), and acetylcholine receptors (AChR) using alpha bungarotoxin coupled with Alexa 555 (White) in control cells or FSHD2 cells at Day 30 post‐differentiation. (C) Schematic representation of the SMCHD1 protein and position of mutations in BAMS (cyan) or FSHD2 patients (red). BAMS‐1 (E136G) carries a missense mutation in the ATPase domain reported as a gain of function. FSHD2 patient #14586 carries a mutation in the ATPase domain (Q193P) reported as a loss of function in the ATPase activity. FSHD2 samples #12051C carrying a synonym mutation at position p.1031 was used as additional control for validation. (D) Venn diagrams for comparison of genes that are differentially expressed in FSHD1, FSHD2 and BAMS MFs compared with controls with a fold‐change −2 < FC > 2 and an FDR < 0.05 in. (E) Analysis of DUX4 expression in the different samples by RT‐qPCR. (F) By comparing our list of DEGs in FSHD1, FSHD2, or BAMS vs. control MFs with DUX4 target genes, only 1 gene, OAS2 (encoding the 2′‐5′‐oligo adenylate synthetase 2), is common between FSHD1, FSHD2, and BAMS MFs, 9 DUX4 target genes are specifically deregulated in FSHD1 MFs, 14 in FSHD2 MFs and 18 in BAMS MFs; 8 are common to FSHD1 and FSHD2. DUX4 target genes that are differentially expressed are listed on the right of the Venn diagram. (G). DUX4 target genes that are differentially expressed in FSHD1, FSHD1 mosaic, and FSHD2 MFs. Red arrows correspond to genes that are up‐regulated and green arrows, to down‐regulated genes. The fold change is indicated for all of them. (H) Validation by RT‐qPCR of DUX4 and selected DUX4 target genes, ZSCAN4 (zinc finger and SCAN domain containing protein 4), MBD3L2 (methyl CpG binding domain protein 3 L2), TRIM43 (tripartite motif containing 43), LEUTX (leucine twenty homeobox), VCAM1 (vascular cell adhesion molecule 1) expression in MFs. Box plots display the results of biological and technical triplicates for each group of samples [controls, FSHD1 short corresponds to the clone containing the contracted D4Z4 allele for the mosaic patient (17706) and FSHD1 long corresponds to its isogenic control; FSHD1 (12759) FSHD2 (14586; 120521C) and BAMS‐1]. Statistical significance was determined using a Kruskal–Wallis statistical test. *P value <0.05, **P value <0.005, ***P value <0.0005, and ****P value <0.00005.

Figure 2

Fi Differentiation of hiPSCs from FSHD1, FSHD1 mosaic, FSHD2, and Bosma arhinia and microphthalmia (BAMS) into functional muscle fibres revealed defects in sarcomere‐related functions in cells affected by FSHD. Overrepresentation test analyses were performed using enrichGO from the R package clusterProfiler (v3.10.15). We identified biological processes (BP) with a false discovery rate (FDR) <0.05. (A). Biological pathways (BPs) corresponding to enrichment analysis of DEGs in FSHD1 vs. control MFs filtered on −2 < FC > 2 and FDR < 0.05. Bar plot in the left represents the percentage of genes that are deregulated and associated with a GO‐term shown in the right column. Light grey bars in the right represent the enrichment score (Log10 of FDR) for each GO‐term. (B) BP corresponding to enrichment analysis of DEGs in FSHD1 mosaic vs. its isogenic control MFs filtered on −2 < FC > 2 and FDR < 0.05. Bar plot in the left represents the percentage of genes that are deregulated and associated with a GO‐term shown in the right column. Light grey bars in the right represent the enrichment score (Log10 of FDR) for each GO‐term. (C) BP corresponding to enrichment analysis of DEGs in FSHD2 vs. control MFs filtered on −2 < FC > 2 and FDR < 0.05. Bar plot in the left represents the percentage of genes that are deregulated and associated with a GO‐term shown in the right column. Light grey bars in the right represent the enrichment score (Log10 of FDR) for each GO‐term. (D) BP corresponding to enrichment analysis of DEGs in BAMS vs. control MFs filtered on −2 < FC > 2 and FDR < 0.05. Bar plot in the left represents the percentage of genes that are deregulated and associated with a GO‐term shown in the right column. Light grey bars in the right represent the enrichment score (Log10 of false discovery rate) for each GO‐term. (E) Heatmap for genes associated with muscle GO terms with dendrogram branches associating the different groups of samples based on transcripts per million (TPM) values. Distance was determined by Manhattan and Clustering, Ward.D2. (F) Venn diagram for FSHD MFs. FSHD1 and FSHD2 DEGs were obtained by comparisons with control MFs, FSHD1 mosaic (contracted D4Z4 allele only). By comparison to the clone with the normal D4Z4 allele (15 RUs), a list of 132 DEGs common to the three conditions was obtained. (G) Cellular components analysis for the list of 132 DEGs in common between FSHD1, FSHD1 mosaic, and FSHD2 MFs. Genes with the higher FC are related to muscle contraction. Twenty genes belong to ‘supramolecular fibres’ cellular component (P value, 2.10e⁻⁰⁶), 17 genes to ‘myofibril’ (P value, 2.19e⁻¹²), 16 genes to ‘sarcomere’ (P value, 2.63e⁻¹²) and 12 to ‘actin cytoskeleton’ (P value, 9.4e⁻⁰⁵).

Figure 3

Active modules analysis in muscle fibre (MF) derived from FSHD1 hiPSCs reveals defects in sarcomeric protein network. (A–H) Representative active modules sampled out of 23 nodes for FSHD1 datasets using the MOGAMUN algorithm. Integration of protein–protein interactions (blue lines), biological pathways (orange lines), and co‐expression data (yellow lines) with our lists of differentially expressed genes (DEGs) enabled the identification of active modules for each category of samples. Up‐regulated nodes are coloured in red, and down‐regulated ones are in green. The intensity of the colour reflects the fold change. Thickness of the dark line around each rectangle reflects the level of significance [false discovery rate (FDR) <0.05 and 1 < FC > 1]. Each active module contains between 15 and 16 genes. Genes corresponding to each nodes were analysed using g:Profiler to define corresponding molecular function and P value. (A) Extracellular matrix (ECM) organization, P value 7.5e⁻⁹. (B) Extracellular matrix organization, P value 1.1e⁻⁹. (C) Actin–myosin filament sliding, P value 9.8e⁻²⁴. (D) Actin–myosin filament sliding, P value 7.9e⁻³⁵. (E) Transmembrane receptor tyrosine kinase signalling, P value 3.3e⁻¹⁸. (F) Actin–myosin filament sliding, P value 7e⁻⁹. (G) Phosphatidyl inositol‐mediated signalling, P value 2.13e⁻¹⁴. (H) Phosphatidyl inositol‐mediated signalling, P value 1. 3e⁻¹⁴.

Figure 4

Network‐based analysis of hiPSCs derived into functional muscle fibres highlights active modules that are specific to each pathology. Representative active modules sampled from the accumulated Pareto front of 30 runs for FSHD2 datasets using the MOGAMUN algorithm. The colour of the nodes represents the fold change, where green and red nodes correspond to down‐regulated and up‐regulated genes, respectively. Nodes with bold black border correspond to genes significantly differentially expressed (FDR < 0.05 and −1 > FC > 1). The colour of the edges corresponds to the layer of the multiplex where the interaction comes from. Blue for protein–protein interactions, orange for biological pathways, and yellow for co‐expression. Each active module contains between 15 and 18 genes. Genes corresponding to each nodes were analysed using g:Profiler to define corresponding molecular function and P value. (A–D) Selection of representative active modules for FSHD2 cells. (A) Actin–myosin filament sliding, P value 2.3e⁻²⁸. (B) Actin–myosin filament sliding, P value 7.5e⁻²¹. (C) Positive regulation of macromolecules synthesis, P value 4.2e⁻¹¹. (D) Actin–myosin filament sliding, P value 1.04e⁻¹⁵. (E–H) Selection of representative active modules for BAMS cells. (E) Actin–myosin filament sliding, P value 2.3e⁻²⁸. (F) Extracellular matrix organization, P value 8.5e⁻¹⁵. (G) Transmembrane receptor tyrosine kinase signalling, P value 5.1e⁻¹⁴; ERBB2 signalling, P value 1.85e⁻¹³. (H) Transmembrane receptor tyrosine kinase signalling, P value 1.85e⁻¹⁴; ERBB2 signalling, P value 8.5e⁻¹³.

Figure 5

Sarcomeric dysfunction and calcium release in vivo and innervated muscle fibres derived from patients affected with FSHD. (A) RT‐qPCR was performed in technical triplicates for each biopsy sample. For each group, two different biopsies from two different individuals were used (Table S1). Expression was normalized to three housekeeping genes (GAPDH, HPRT, and PPIA). Statistical significance was determined by Kruskal–Wallis statistical test. *P value <0.05, **P value<0.005, ***P value <0.0005, and ****P value<0.00005. (B) Muscle contractions were recorded for 30 s with images taken every 200 ms in live innervated muscle fibres (Movie S1). (C) Scattergram displays the length of contraction determined as the shifting over time of randomly distributed point (n = 1200 data points per condition) along independent muscle fibres (y‐axis) analysed using the Imaris software in the different conditions. Data were analysed using a Brown–Forsythe statistical test with a Games‐Howell correction, ****P value <0.0001. (D) Intracellular calcium signalling was measured in hiPSC‐ derived live myofibres using an Axio observed microscope at day 45 post‐differentiation. Different regions were tracked to record the fluorescent intensity at a 10× magnification during 30 s and an interval of 200.0 ms corresponding to 151 frames. Regions of interest (ROIs) corresponding to muscle fibres were selected and intensity of fluorescence was acquired using the Imaris software. (E) Scattergram displays the amplitude of calcium uptake and release in 50 fibres randomly selected in different fields for each sample (y‐axis). Data were analysed using a Kruskal–Wallis non‐parametric test (ns: non‐significant; ****P value <0.0001).