Duchenne Muscular Dystrophy human iPSC (DMD-hiPSC) was obtained from DMD patient’s fibroblasts carrying a point mutation in exon 35 (c.4 996C>T; p.Arg1,666X) of the DMD gene that leads to a premature stop codon (Duelen et al., 2022). The DMD isogenic control (DMD-Iso-hiPSC) was generated through CRISPR/Cas9 gene editing from the S. pyogenes system (5’-NGG PAM) as previously described (Ran et al., 2013; Duelen et al., 2022). The healthy control iPSC (HC-hiPSC) line was a gift from Prof. P. Jennings (Medizinische Universität Innsbruck, Austria) to Stem Cell Institute Leuven, KU Leuven and generated by SeV-based reprogramming of male donor fibroblasts (SBAD2). Human iPSC lines were cultured feeder-free on Geltrex LDEV-Free hESC-Qualified Reduced Growth Factor Basement Membrane Matrix and maintained in Essential 8 Flex Basal Medium (Thermo Fisher Scientific) supplemented with Essential 8 Flex Supplement (50x, Thermo Fisher Scientific) and penicillin–streptomycin (0.1%, Thermo Fisher Scientific), at 37°C under normoxic conditions (21% O₂ and 5% CO₂). Colonies were routinely passaged non-enzymatically with 0.5 mM EDTA in Phosphate-Buffered Saline (PBS, Thermo Fisher Scientific). The use of human samples from DMD subjects for experimental purposes and protocols in the present study was approved by the Ethics Committee of the University Hospitals Leuven (S55438 and S65190, respectively).

DMD-, DMD-Iso- and HC-hiPSC lines were differentiated into functional cardiomyocytes (CMs) according to a monolayer-based cardiac differentiation protocol, as previously described (Burridge et al., 2014). Briefly, prior to differentiation, the DMD-, DMD-Iso- and HC-hiPSC lines were suspended into small colonies and subsequently cultured on Matrigel Growth Factor Reduced (GFR) Basement Membrane Matrix layer (Corning) in complete Essential 8 Flex Medium at 37°C under hypoxic conditions (5% O₂ and 5% CO₂) for 3 days, in order to obtain the pre-optimized targeted confluency of 85%. Mesoderm differentiation (day 0) was induced using 6 µM CHIR99021 (Axon Medchem) for 48 h in a chemically defined medium consisting of RPMI 1640 (Thermo Fisher Scientific), 500 μg/ml rice-derived recombinant human albumin and 213 μg/ml L-ascorbic acid 2-phosphate (Sigma-Aldrich). After 24 h of CHIR99021 stimulation, the cells were transferred from hypoxia to normoxia. On day 2 of differentiation, iPSC-derived mesodermal cells were fed with basal medium supplemented with 4 µM IWR-1 (Sigma-Aldrich) for 48 h, to induce cardiac progenitor cell differentiation. From day 4 onwards, medium was refreshed every 2 days with CM Maintenance Medium (RPMI 1640, rice-derived recombinant human albumin and L-ascorbic acid 2-phosphate). Contracting CMs appeared at day 8 or 9 of cardiac differentiation.

A 3% agarose (Invitrogen) gel solution was prepared in PBS. The powder was fully dissolved by heating in microwave oven and the agarose microwells were fabricated in sterile conditions. In brief, the heated agarose solution was added into a custom-made 3D printed micropillar molds (in 24-well plate format). Upon cooling at room temperature for 10 min, the agarose was removed from the molds, thus creating 24 culture inserts each consisting of 137 microwells (diameter × height = 500 × 700 μm). The culture inserts were transferred into a 24-well plate and equilibrated in PBS overnight at 37°C under normoxia conditions (5% O₂ and 5% CO₂).

After reaching confluency, the DMD-, DMD-Iso- and HC-hiPSC lines were detached using 0.5 mM EDTA at 37°C and resuspended in Essential 8TM medium supplemented with Revitacel^™ Supplement (dilution 1:100, Thermo Fisher Scientific). After cell count, the hiPSCs were resuspended in 1 ml of Essential 8 Flex Basal Medium (Thermo Fisher Scientific) and were plated in agarose inserts at two different cell densities, 5 × 10³ cells/microwell and 1 × 10⁴ cells/microwell respectively. The plates were centrifuged for 10 min at 1,200 rpm to facilitate sedimentation of cells in the microwells. Subsequently, 1 ml of fresh Essential 8 Flex Basal Medium was added to completely cover the microwell area and incubated at 37°C under hypoxic conditions (5% O₂ and 5% CO₂) to promote embryoid bodies (EBs) formation. The medium was refreshed every day for 3 days and cardiac differentiation of the EBs into cardiac organoids (COs) was initiated as described above for the monolayer cardiomyocyte differentiation protocol. On day 5, the COs were transferred from the agarose molds to an ultra-low attachment 6-well plate (Costar, Corning) and dynamic culture was carried out using an orbital shaker at 75 rpm in CM maintenance medium until day 93. The media was changed every 2 days. Contracting COs start to appear from day 8 of the differentiation protocol. The samples were collected on day 10, 14, 28, 56 and 93 for subsequent analyses.

At different time points, the COs were fixed with 4% paraformaldehyde (PFA; Polysciences) for 30 min at room temperature and subsequently embedded in cryogel (Tissue- Tek ^® O.C.T. ^™ Compound). The samples were snap-frozen in liquid nitrogen and stored at −80°C until cryosectioning. The samples were sectioned at the thickness of 6 µm using the HM525 NX Cryostat (Thermo Scientific) and stored at −20°C prior to analysis. For H&E staining, the cryosections were stained in Harris hematoxylin solution (Sigma-Aldrich), counterstained in eosin solution (0.1% erithrosin extra bluish Sigma-Aldrich in 70% ethanol) and mounted with DPX mountant (Sigma) upon dehydration according to routine protocols. For PSR staining, the cryosections were stained for collagen content using the Vitro View^™ Picro-Sirius Red Stain Kit (Cat. No. VB-3017) according to the manufacturer’s instructions (Giarratana et al., 2020). The nuclei were counterstained with Weigert’s Hematoxylin Solution and mounted with DPX mountant (Sigma-Aldrich). Lipid droplets deposition was detected by BODIPY staining. In brief, the BODIPY^™ 493/503 4,4-Difluoro-1,3,5,7,8-Pentamethyl-4-Bora-3a,4a-Diaza-s-Indacene (Invitrogen) powder was dissolved in DMSO at the concentration of 1.3 mg/ml. The cryosections were incubated with the BODIPY solution, diluted 1:2500 in PBS, for 15 min at room temperature and subsequently mounted with Antifade Mounting Medium with DAPI (VECTASHIELD^®). All images were acquired using Axiocam MRm microscope (Zeiss).

For quantitative Reverse Transcription Polymerase Chain Reaction (RT-qPCR) assays, total RNA was isolated through Purelink^® RNA mini kit (Thermo Fisher Scientific) and treated with TurboTM DNA-free kit (Thermo Fisher Scientific) to purify RNA samples. 1 μg RNA was reverse-transcribed using Superscript III Reverse Transcriptase First-Strand Synthesis SuperMix (Thermo Fisher Scientific). Thermal cycler setting: 25°C 10 min, 50°C 30 min, 85°C 5 min, 37°C 20 min incubation with E. Coli RNAse H. A 384-well plate was prepared using Platinum SYBR Green QPCR SuperMix-UDG (Thermo Fisher Scientific) as SYBR Green on 1:5 diluted cDNA. The RT-qPCR was performed by Viia7 384-plate reader (Thermo Fisher Scientific; final primer concentration, 100 nM; final volume, 10 μl; thermal profile, 95°C 15 s, 60°C 60 s, 40×). The oligonucleotide primer sequences are listed in Table 1. Delta Ct (ΔCt) values were calculated by subtracting the Ct values from the genes of interest with the Ct values of the housekeeping genes (GAPDH).

DMD- and DMD-Iso-COs at 14 days were dissociated using Collagenase A (1 U/ml) for 20 min at 37°C followed by 10 min incubation with Accutase^® (Sigma-Aldrich). All flow cytometry procedures were performed according to the manufacturer’s instructions. PBS + 0.1% BSA (Sigma) was used as staining buffer. The cells were stained for the surface markers SIRPA, CD31 and PDGRFa specific for CMs, endothelial cells and fibroblast/adipose progenitors respectively. A total of 50.000 cells was recorded. After single cell gate selection, a cell number total between 15.000 and 20.000 cells was analysed for the used antibodies (CD31, PDGFRα, SIRPA). Fluorescence minus one (FMO) controls and compensations were included for appropriate gating. Samples were analysed using the FACS Canto II HTS (BD Biosciences) and the analysis was performed using FACS Diva Software. Table 2 provides a list of all flow cytometric antibodies used in this study.

After three PBS washes, the cryosections were permeabilized for 1 h at room temperature using 0.1% Triton X-100 in PBS (Thermo Fisher Scientific). Non-specific antibody binding was blocked by incubation for 30 min with blocking solution containing 5% normal goat serum (NGS, Dako) at room temperature followed by overnight incubation at 4°C with different primary antibodies listed in Table 2. After washing in phosphate-buffered saline (PBS), the samples were incubated with respective secondary antibodies using Alexa Fluor 488 and 555-conjugated secondary antibody (4 μg/ml; Thermo Fisher Scientific). Nuclei were counterstained with Hoechst 33342 (1:1000, Thermo Scientific) for 7 min (Santoni de Sio et al., 2008). The sections were mounted with ProLong ^™ Gold antifade reagent (Invitrogen) and stored in the dark at 4°C until imaging. All images were acquired using Axiocam MRm microscope (Zeiss). Confocal images were acquired using Nikon ECLIPSE Ti Microscope and NIS-Elements AR 4.11 software. Quantifications of immunofluorescent images were performed using ImageJ software tool by calculating the percentage of the signal area/positive nuclei with respect to the total organoid area/total nuclei respectively for each CO.

To assess the contractile properties of DMD-COs and DMD-Iso-COs, 3D cardiac organoids were live-imaged using the Dmi1 Microscope (Leica). The recorded videos were then analysed to determine the CO beating frequency by counting the number of spontaneous contractions per minute. The cardiac organoids growth area was measured at different time points using ImageJ software tool.

For Ca²⁺ imaging experiments, the DMD-, DMD-Iso- and HC-hiPSC monolayers were respectively plated on 35 mm dishes with four Chamber glass bottom. Following 14 days from cardiac induction, the DMD-, DMD-Iso- and HC- CM were incubated with 1 µM Fluo-4 AM solubilized in CM Maintenance Medium. Next, the cells were washed twice with CM Maintenance Medium after which de-esterification was allowed to occur for 45 min at 37°C and 5% CO₂. The Ca²⁺ imaging experiments were performed in pre-warmed (37°C) modified Krebs-Ringer solution (135 mM NaCl, 6.2 mM KCl, 1.2 mM MgCl₂, 12 mM HEPES, pH 7.3, 11.5 mM glucose and 2 mM CaCl₂). Tetracaine was solubilized in the above modified Krebs-Ringer solution at 1 mM final concentration. For the KCl stimulus the modified Krebs-Ringer solution was prepared substituting the NaCl for 140 mM KCl. Imaging was performed using a Nikon eclipse Ti2 inverted fluorescence microscope (Nikon) equipped with excitation filter FF01-378/474/554/635 and dichroic mirror FF01-432/515/595/730 and emission filter 515/30 all from Semrock. Coolled pR-4000 (Coolled) was used for excitation at 470 nm. Acquisition of the fluorescent signal at 520 nM was performed at 10 Hz using a pco.edge 4.2bi sCMOS camera (pCO) (Nakamura et al., 2001). For analysis FIJI software was utilized. In each experiment a region of interest was drawn across spontaneously active cardiomyocytes. The fluorescence intensities were normalized to F0, where the F0 value was obtained after tetracaine administration.

RNA (>10 μg) extracted from DMD-COs and DMD-Iso-COs on day 56 was verified and processed by the Genomics Core (KU Leuven–UZ Leuven). As quality control, the RNA concentration was measured with Nanodrop and quality was checked with Bioanalyzer. The Lexogen QuantSeq 3′ mRNA-Seq library prep kit was used according to the manufacturer's protocol with 500 ng input. After the prep, the libraries were measured with Qubit and put on the Fragment analyser so the libraries could be pooled equimolar to 2 nM. The pool was then quantified with RT-qPCR and a final pool (2 nM) was made for single-read sequencing on the HiSeq4000 (Illumina Inc.). The settings were 51-8-8. The raw sequence files generated (.fastq files) underwent quality control analysis using FastQC v0.11.7 (Andrews, 2010). Adapters were filtered with ea-utils fastq-mcf v1.05 (Aronesty, 2011). Splice-aware alignment was performed with HiSat2 against the human reference genome hg38 using the default parameters. Reads mapping to multiple loci in the reference genome were discarded. Resulting BAM alignment files were handled with Samtools v1.5. (Li et al., 2009). Quantification of reads per gene was performed with HT-seq Count v2.7.14. Count-based differential expression analysis was done with R-based (The R Foundation for Statistical Computing, Vienna, Austria) Bioconductor package DESeq2 (Love et al., 2014). Reported p-values were adjusted for multiple testing with the Benjamini-Hochberg procedure, which controls false discovery rate (FDR). Gene Ontology (GO) and Biological Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were identified using g:Profiler (Raudvere et al., 2019). The GO Biological Process 2018 and KEGG 2016 of each tissue were determined. The significant terms and pathways were selected with the threshold of adjusted p-value < 0.05. Data has been deposited in the NCBI Gene Expression Omnibus (GEO) repository under accession code GSE194297.

The PPI network of differentially upregulated genes in DMD-COs was constructed by feeding a list of gene symbols and their log₂fold changes into the NetworkAnalyst platform (http://www.networkanalyst.ca/) using the IMEx interactome database with Steiner Forest Network (SFN) reduction algorithm. Subsequently, the gene-miRNA interactions (Rotini et al., 2018) for the selected KEGG pathways were constructed based on the miRTarBase (v8.0) database, and the network was reduced using the SFN algorithm. The degree of each node was calculated based on its number of connections to other nodes. In the network, the area of an individual node indicates the degree, and the color represents the expression. The identified top five miRNAs were mapped out in the KEGG pathways to show their interactions with the genes of a particular pathway.

Data were statistically analysed using GraphPad Prism. All data were reported as mean ± standard deviation (SD). Differences between groups were examined for statistical significance using ANOVA and two-way ANOVA. Significance of the differences was indicated as follows: *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001.

Following the 2D monolayer differentiation protocol, we generated cardiomyocytes (CMs) from both DMD patient-derived hiPSC (DMD-CMs), the isogenic control hiPSC (DMD-Iso-CMs) and healthy control hiPSC (HC-CMs) monolayer cultures. Figure 1A shows representative images of 2D culture morphology of DMD-CMs and DMD-Iso-CMs on day 22. These cells started to develop contractile phenotype around day 8 and were morphologically similar in both conditions. On day 22, RT-qPCR analysis showed that the DMD-CMs expressed significantly lower dystrophin (DMD) than the isogenic controls and HC-CMs (Figure 1B), confirming the restoration of DMD expression in the isogenic controls, as described in Duelen et al. (2022). The DMD-CMs also expressed significantly lower sarcomeric α-actinin (ACTN2), the pacemaker gene HCN4, and the troponin-related genes (TNNI1, TNNC1 and TNNI3) but not the myosin light (MYL7, MYL2) or heavy chain (MYH6, MYH7) genes, than the isogenic controls (Figure 1C). Next, we established a cell physiological analysis of 14-day differentiated DMD-CMs and DMD-Iso-CMs. A hallmark of functional CMs is their ability to generate cytosolic Ca²⁺ signals that are driven by ryanodine receptors (RyRs), intracellular Ca²⁺ -release channels residing at the sarcoplasmic reticulum of CMs. Therefore, cytosolic Ca²⁺ imaging was performed in single-cell CMs loaded with Fluo-4. In the presence of extracellular Ca²⁺ (1.5 mM CaCl₂), spontaneous Ca²⁺ oscillations were observed in both DMD-CMs and DMD-Iso-CMs that could be blocked by tetracaine, an inhibitor of RyR channels. However, spontaneous Ca²⁺ oscillations appeared to display a lower frequency and amplitudes in DMD-CMs compared to isogenic and healthy controls, indicating a defect in physiological Ca²⁺ signalling in dystrophic CMs that is corrected in the isogenic controls. Moreover, DMD-CMs displayed a lower Ca²⁺ response to KCl, which provokes membrane depolarization, compared to DMD-Iso-CMs (Figures 1D,E). These findings validated the dystrophic properties of DMD-CMs and their defects in physiological Ca²⁺ signalling, whereby both deficiencies could be reverted in DMD-Iso-CMs generated in this study.

We adapted the cardiomyocyte monolayer differentiation protocol to generate COs by direct differentiation of the embryoid bodies (EBs) (Figure 2A). By using the agarose microwell culture inserts, we could promote self-aggregation of the DMD-hiPSCs and DMD-Iso-hiPSCs into EBs at cell seeding number of 5,000 and 10,000 cells per microwell (Figure 2B). This allowed us to generate 137 EBs per insert per well of 24-well plate. On day 5 of cardiomyocyte differentiation, the resulting DMD-COs and DMD-Iso-COs were transferred to 6-well plate on orbital shaker for dynamic culture in the CM maintenance medium (Figure 2C). Contractile cardiomyocyte protrusions (Figure 2C, arrow) and self-organized cellular structures (Figure 2C, arrowheads) at the organoid periphery, both with specific spatial distribution of NKX2.5 and αACTN positivity, could be observed. The non-translucent organoid structure (#) was negative for both NKX2.5 and αACTN. Immunofluorescence staining showed abundant DYS localization in DMD-Iso-COs, which was undetectable in the DMD-COs (Figure 2D). The representative images of CO over 28 days of dynamic culture showed morphological changes and variations in organoid size (Figure 2E). Quantification of the organoid surface area over 28 days of dynamic culture showed no significant differences on the organoid size between the two cell seeding numbers within each cell line, however, the size of DMD-Iso-COs was significantly smaller than DMD-COs on day 14 and 28, respectively (Figure 2F). The DMD-COs displayed contraction on day 8 (± 19 per minute) which decreased over time and stopped contraction between day 14 and 18 (Figure 2G). The DMD-Iso-COs displayed contraction on day 12 (± 18 per minute) which persisted till day 28. FACS analysis was performed on DMD- and DMD-Iso-COs after 14 days to determine SIRPA (APC), CD31 (PE/Cy7) and PDGFRα (PE) positive cells (Supplementary Figures S1A,B). The analysis disclosed that both DMD- and DMD-Iso-COs accounted for 40/50% of SIRPA⁺ cells, while DMD-Iso-COs showed the double amount of the CD31⁺ cell (1.4% ± 0.2) population and the half of the PDGFRα⁺ cells (1.7% ± 0.2) compared to DMD-COs.

We performed immunofluorescence staining for α-sarcoglycan (SCGA), sarcomeric α-actinin (αACTN), and NKX2.5 on day 10, 14 28, 56 and 93 in order to assess cardiac differentiation and contractile protein development within the organoids. The results showed abundant SCGA expression in DMD-COs on day 10, which became low on day 14 and undetectable from day 28 onwards (Figures 3A,C and Supplementary Figure S2A). Conversely, the SCGA expression in DMD-Iso-COs persisted till day 93. As reported in Supplementary Figure S3A also the β-, γ- and δ-sarcoglycans started to disappear from day 14 following the α-sarcoglycan trend. A transient expression of the early cardiac differentiation marker NKX2.5 was observed up to day 28 in both DMD-COs and DMD-Iso-COs, which became undetectable on day 56 and 93 (Figure 3C). Additionally, abundant αACTN, a cardiac contractile protein, was observed in both DMD-COs and DMD-Iso-COs on early time points, which remained detectable on day 93 (despite at lower expression level) in both conditions (Figures 3B,C and Supplementary Figure S2B). There was no distinguishable difference in the SCGA, αACTN and NKX2.5 expression between organoids generated from the two cell seeding numbers within a cell line. These results demonstrated a progressive loss of SCGA protein expression in DMD-COs (a member of the dystrophin associated complex, DAC) as compared to the isogenic controls. The sarcomeric pattern is shown at higher resolution and magnification in Supplementary Figure S3B which reports representative confocal images of αACTN and SCGA. In order to better characterize the cardiomyocyte population within the COs we performed the RT-qPCR for atrial/ventricular gene markers and TNNT2 at all the time points (Supplementary Figures S4A,B). Both atrial and ventricular markers showed a peak at the initial stage followed by a decline at later stages. Intriguingly IRX4, a ventricle-specific transcription factor, appeared to be always expressed in both samples until day 28 and declined at late stages as the majority as the other atrial and ventricle markers. However, the two major ventricular myosin heavy-chain genes MYL2 and MYH7 were sustained longer during the differentiation compared to the atrial ones (Supplementary Figure S4A). Likewise, TNNT2 expression reached the peak at day 14 and then decreased without showing striking differences between the two cell lines (Supplementary Figure S4B). Additionally, RT-qPCR analysis showed a significant upregulation of genes related to cardiac contractility in DMD-COs as compared to DMD-Iso-COs, in particularly from day 56 onwards (Figure 3D). These included ACTN1, IRX4, MYBPC3, MYL2, MYOM1, TNNC2 and TPM1.

We examined cell proliferation and apoptotic conditions within the DMD-COs and DMD-Iso-COs by immunostaining of the proliferation marker Ki67 and apoptotic marker cleaved caspase 3 (CCASP3). The results showed low Ki67 staining in DMD-COs but relatively higher signal in DMD-Iso-COs on day 10, while the signal became comparable on day 28 and 93 (Figures 4A,E and Supplementary Figure S2C). Though initially CCASP3 signal area was slightly more in DMD-Iso-COs, from day 28 it became significantly higher in DMD-COs (Figure 4E). However, no significant difference in Ki67 and CCASP3 staining was observed at both cell seeding densities for both COs conditions (data not shown). These data suggest that the DMD-COs was lacking an initial proliferative capacity at early time point while at later time points, they were more apoptotic. We then assessed the metabolic activity within the COs by immunostaining of the glycolytic marker phosphoglycerate kinase 1 (PGK1). The results showed high and comparable PGK1 staining in both COs conditions at all time points (Figures 4B,E and Supplementary Figure S2D), which was independent from the cell seeding densities (data not shown). These results suggest the glycolytic condition of immature COs in both CO conditions. The cellular stress was assessed by immunostaining of two known endoplasmic reticulum (ER) stress markers ARCN1 and GORASP2. Interestingly, we detected relatively higher level of ARCN1 in DMD-COs than DMD-Iso-COs at all the time points (Figures 4C,E and Supplementary Figure S2E), whereas GORASP2 increased progressively over the 28 days in DMD-COs, independently from the cell seeding densities (data not shown) and at higher level than that in DMD-Iso-COs at all the time points (Figures 4D,E and Supplementary Figure S2F). This finding indicated a high level of ER stress occurred in DMD-COs. To assess oxidative stress due to high NOX4 expression, we performed immunofluorescence staining on the organoid cryosections for NOX4. As shown in Supplementary Figure S5A, high NOX4 protein level was detected in the DMD-COs on day 14 which was barely present in DMD-Iso-COs.

We performed histological examination to assess any cytoarchitecture changes and DMD-related pathological progression within the COs over 93 days. The DMD-COs displayed normal cardiomyocyte-like structures similar to that of DMD-Iso-COs on day 10, which deteriorated on day 14 (indicated as “#”) and developed fibrotic-like structures (indicated as “f”) at later time points (Figure 5A; H&E staining on day 56, and Figure 5B; Picro-Sirius red staining for collagen deposition on day 93). These findings were corroborated by a significant upregulation of gene markers associated with fibrosis COL1A2, COL3A1 and FN1 in DMD-COs on day 56 and 93 as compared to DMD-Iso-COs (Figure 5C). Additionally, H&E staining also revealed adipose tissue formation in DMD-COs on day 28 (Figure 5D), confirmed by the detection of lipid droplets via BODIPY staining and immunolabelled PDGFRα⁺ cells (an adipocyte marker) in DMD-COs on day 28 and 56 (Figure 5E). Interestingly, GDF10 protein (an adipogenesis inhibitor) was also detected near the PDGFRα⁺ cells in DMD-COs (Figure 5F). These findings suggest that DMD-COs displayed an initial normal cardiac phenotype which deteriorated progressively and exhibited pro-fibrotic and adipogenic phenotypes upon long-term culture, resembling pathologic events associated with DMD cardiomyopathy.

Principal component analysis of the RNA transcriptomic data showed a distinct separation between DMD-COs and DMD-Iso-COs clusters (PC1: 90%) with low intra-condition variance (PC2: 5%) (Figure 6A). Based on the enhanced volcano plot, out of 22,371 gene variables, 1,518 and 554 genes were differentially upregulated in DMD-COs and DMD-Iso-COs, respectively (Cut-off: log₂ fold change = 1.5; −Log₁₀ p = 10^–16) (Figure 6B). Among the top 30 most differentially upregulated genes in DMD-COs (Figure 6C), MGP, MYL1, COL1A2, HAPLN1 and OGN were the five most significant upregulated genes (Figure 6B and Table 3). In addition, we validated the upregulation of the majority of the genes in DMD-COs showed in Figure 6B by RT-qPCR in DMD-, DMD-Iso- and HC-COs at day 56 as additional control (Supplementary Figure S6A). Based on gProfiler analysis, gene ontologies that were significantly enriched for molecular function in extracellular matrix regulation (i.e. collagen and glycosaminoglycan; GO:MF), cardiac tissue structure formation (i.e. external encapsulating structure such as sarcolemma; GO:MM), and cardiovascular development (GO: BP) could be identified in DMD-COs (Figure 6D). Additionally, KEGG pathways associated with protein digestion and absorption, dilated and hypertrophic cardiomyopathy, ECM-receptor interaction, and cGMP-PKG signalling pathway (known to positively modulates cardiac contractility, hypertrophy and protects against apoptosis (Takimoto, 2012) were significantly enriched in DMD-COs (Table 4 and Figure 6E (i)). These findings were corroborated by the analysis on human phenotype ontology, whereby ontology related to abnormal cardiovascular system physiology, including abnormal left ventricular function, abnormal endocardium morphology, atrial arrhythmia and fibrillation, supraventricular arrhythmia, myopathy and cardiac arrest, as well as abnormal adipose tissue morphology and lipodystrophy were significantly enriched in DMD-COs as compared to DMD-Iso-COs (Figure 6E (ii)). Moreover, the gProfiler analysis also identified three top miRNA regulators for the differentially upregulated genes in DMD-COs, namely hsa-mir-335-5p, hsa-mir-29a-3p and hsa-mir-29b-3p. Altogether, the RNA sequencing data validated the histological observations described above on cardiomyocyte deterioration, adipogenesis and fibrosis at the transcriptomic level.

PPI analysis of the differentially upregulated genes in DMD-COs revealed a gene network consisted of 2,289 nodes and 2,288 edges. According to the degree level (d), the top five hub nodes were HNF4A (d = 257), UBC (d = 108), UBD (d = 66), APP (d = 38) and EGR1 (d = 31) (Figure 7A). By exploring the miRNA database (i.e. mirTarBase v8.0), the top three miRNA regulators of this gene network were hsa-mir-335-5p, hsa-mir-124-3p, and hsa-mir-26b-5p. Together with the hsa-mir-29b-3p and hsa-mir-29a-3p identified by gProfiler2, we mapped out these miRNAs on the gene-miRNA regulatory networks for the selected KEGG pathways relevant to the DMD-COs phenotypes: 1) Hypertrophy cardiomyopathy, 2) Dilated cardiomyopathy, 3) Arrhythmogenic right ventricular cardiomyopathy (ARVC), 4) PPAR signalling pathway (for adipogenesis), and 5) PI3K-Akt signalling pathway (for cardiac fibrosis (Table 5) (Qin et al., 2021)). The results showed that hypertrophy and dilated cardiomyopathy networks shared the same gene set (50 nodes, 49 edges), miRNA interactions (Figure 7B) and 16 genes similarity with the ARVC network (Figure 7C). Except hsa-mir-124-3p, the other four top miRNAs were mapped in these three networks, respectively. The PPAR signalling gene-miRNA network consisted of 33 nodes and 43 edges (Figure 7D). In addition to hsa-mir-26b-5p and hsa-mir-355-5p, the hsa-mir-124-3p was mapped in the network and found interacts with the gene ACSL5 and ACADL. The has-mir-29b-3p, hsa-mir-26b-5p and hsa-mir-355-5p were the main miRNA regulators in the PI3K-Akt signalling (147 nodes, 146 edges), which interact with one of the two hub genes CCND2 (Figure 7E).