Comprehensive analysis of m6A regulators characterized by the immune microenvironment in Duchenne muscular dystrophy
Journal of Translational Medicine volume 21, Article number: 459 (2023)
Duchenne muscular dystrophy (DMD) is an X-linked, incurable, degenerative neuromuscular disease that is exacerbated by secondary inflammation. N6-methyladenosine (m6A), the most common base modification of RNA, has pleiotropic immunomodulatory effects in many diseases. However, the role of m6A modification in the immune microenvironment of DMD remains elusive.
Our study retrospectively analyzed the expression data of 56 muscle tissues from DMD patients and 26 from non-muscular dystrophy individuals. Based on single sample gene set enrichment analysis, immune cells infiltration was identified and the result was validated by flow cytometry analysis and immunohistochemical staining. Then, we described the features of genetic variation in 26 m6A regulators and explored their relationship with the immune mircoenvironment of DMD patients through a series of bioinformatical analysis. At last, we determined subtypes of DMD patients by unsupervised clustering analysis and characterized the molecular and immune characteristics in different subgroups.
DMD patients have a sophisticated immune microenvironment that is significantly different from non-DMD controls. Numerous m6A regulators were aberrantly expressed in the muscle tissues of DMD and inversely related to most muscle-infiltrating immune cell types and immune response-related signaling pathways. A diagnostic model involving seven m6A regulators was established using LASSO. Furthermore, we determined three m6A modification patterns (cluster A/B/C) with distinct immune microenvironmental characteristics.
In summary, our study demonstrated that m6A regulators are intimately linked to the immune microenvironment of muscle tissues in DMD. These findings may facilitate a better understanding of the immunomodulatory mechanisms in DMD and provide novel strategies for the treatment.
Duchenne muscular dystrophy (DMD) is an X-linked, incurable, degenerative neuromuscular disease caused by mutations in the DMD gene coding for dystrophin protein. The absence of dystrophin compromises the integrity of the sarcolemma and leads to uncontrolled inflammation, which is followed by extensive degeneration of the muscle fibers . Currently, there is no cure for DMD, although numerous therapeutic strategies have been developed to improve survival. Glucocorticoids remain the standard of therapy, but their use is limited by the occurrence of side effects such as Cushing’s syndrome. Several promising therapeutic strategies aimed at the restoration of dystrophin production, including gene therapy and stem cell therapy, have been hampered by the few benefited population and the hosts' immune response [2,3,4,5]. Therapies designed to ameliorate inflammation in the muscle microenvironment represent a feasible therapeutic avenue to both prevent muscle deterioration and enhance the tolerability of emerging approaches . Therefore, further characterization of the muscle microenvironment and extensive exploration of the immunomodulatory mechanisms is indispensable to develop effective therapies.
Compelling evidence suggests the crosstalk between the immune system and DMD [7, 8]. A previous study has identified that aberrant signaling pathways regulate immune processes leading to the degenerative process of DMD . Enhanced expression of inflammatory genes and increased infiltration of activated immune cells are evident early in the progress of DMD [10, 11]. Since many unknown factors could influence the immune status, the regulatory mechanisms responsible for immunity are not fully elucidated. N6-methyladenosine (m6A) modification is the most prevalent internal transcript modification of RNA in eukaryotes, which is dynamically mediated by specific m6A regulatory enzymes, including “methyltransferases” (mainly METTL3 and METTL14), “reading proteins” and “demethylases” (ALKBH5 and FTO) . m6A modification is widely involved in various physiological and pathological processes [13, 14]. Emerging evidence indicates that aberrant expression and mutation in the m6A regulators were related to abnormal processes, including metabolism abnormality, dysregulated cell cycle and proliferation, etc. [15, 16]. Recently, several studies have demonstrated that m6A regulators have a close relationship with immunological regulation . For instance, the deletion of m6A reader YTHDF2 enhances the activation of NF-κB and MAPK signaling pathways to upregulate the expression of osteoclast-associated gene and immunity processes . Besides, the m6A writer METTL3 facilitates M1 macrophage to M2 macrophage polarization by STAT1 methylation . To our knowledge, however, few studies have explored the relationship between m6A modification and the immune microenvironment in DMD until now.
In this work, we studied the characteristics of the immune environment in the muscle tissues of DMD based on the next-generation sequencing data. Flow cytometry (FCM) of the muscle tissues in mdx mice (a mouse model of DMD) and immunohistochemistry (IHC) in DMD patients were introduced to validate the infiltration of dominant immune cells. Then, we performed a systematic assessment of the DMD m6A modification pattern and revealed the close relationship between m6A regulators and the immune microenvironment in the muscle tissues of DMD. In summary, our findings uncovered the potential role of m6A modification in the immune microenvironment of DMD and may provide new potential therapeutic avenues for this disease.
Materials and methods
Mdx and C57BL/6 mice were purchased from the Nanjing Biomedical Research Institute of Nanjing University (Nanjing, China). All experiments were conducted based on protocols and approved by the Second Hospital of Hebei Medical University Animal Care and Use Committee (approval number: 2022-AE283). The sample sizes of mice in the experiment were established according to previous experience and the analyses were terminated when the differences between each group were considered statistically significant .
Human tissues were collected from patients with suspected muscle disease admitted to the Second Hospital of Hebei Medical University. All of the patients signed written informed consents to allow the collection of muscle samples and agreed to use these samples/cells for research purpose. The diagnosis of DMD was confirmed by genetics. Negative muscle samples included patients referred for muscle discomfort who had normal histology, histo-enzymology, and immunohistochemistry at the time of muscle biopsy assessment (Additional file 7: Table S1). The analyses were stopped after analyzing 4 patients’ biopsies because a clear statistical difference between DMD patients and non-DMD controls (n = 4) was observed.
Microarray datasets collection and data process
Microarray datasets were retrieved from the Gene Expression Omnibus (GEO) database, including GSE109178, GSE6011, GSE38417, and GSE1004 [21,22,23]. For GSE109178 and GSE38417, the probe IDs were annotated using the platform GPL570, while GSE6011 was on the platform of GPL96. The expression dataset profile of GSE1004 was based on the GPL91 and GPL8300 platform. R package “limma” was used for background adjustment and quantile normalization. If a gene symbol corresponds to multiple probes, the average level of the expression value will be determined. Due to the small sample size could affect the power of statistical analysis and lead to inaccurate results, three datasets (GSE109178, GSE6011, and GSE38417) were integrated as a training set to expand the sample size. We selected an independent dataset GSE1004 to externally validate the gene expressions of key regulators. The Combat function of the “sav” R package was used to correct batch effects and principal component analysis (PCA) was introduced to evaluate the performance (Additional file 1: Figure S1).
Immune characteristics analysis for the microarrays datasets
Gene set enrichment analysis (GSEA) was used to explore the potential immunological pathways by GSEA software (version 4.1.0). Single-sample gene set enrichment analysis (ssGSEA) was conducted to assess the immunocyte fractions in DMD patients. The list of genes involved in gene-sets of infiltrating immune cells was obtained from the prior study . To identify the variation of biological processes between DMD and normal tissues, R package “GSVA” was introduced to run Gene set variation analysis (GSVA) enrichment analysis, and the latest version of immune response gene-sets was acquired from the platform MSigDB (http://software.broadinstitute/org/gsea/msigdb/). The Wilcox test was introduced to analyze the enrichment scores of immune response activity and immune cell abundance between muscle biopsy specimens from patients with DMD and non-DMD controls.
Isolation of muscle leukocytes and flow cytometry analysis
Four-week-old male C57BL/6 and mdx mice were euthanized via cervical dislocation, and the muscles from mouse limbs were harvested and rinsed in cold saline. Muscle tissues were then prepared for single-cell suspension by mesh rubbing method. Briefly, muscles were placed on a 150-mesh sieve, washed with saline three times, and a 25 mL small beaker was placed under the sieve. Then, the tissues were cut into pieces, rinsed with saline and collected in the beaker. The mixture was then filtered and a 300-mesh nylon sieve was used to remove cell debris followed with centrifuged at 157g for 5 min. We collected and re-suspended the pellet, layered it on an equal volume of Lymphocyte separation media (MultiSciences, Hangzhou, China), and centrifuged at 400g for 20 min. The interface of cells was collected, re-suspended in 4 ml of saline, and centrifuged at 157g for 5 min. The supernatants were discarded and the rest were re-suspended with saline. The following antibodies used for staining were purchased from MultiSciences (Hangzhou, China): APC-Cy7-anti-CD3, APC-anti-F4/80, FITC-anti-CD4, PE–anti-CD45, PE-Cy7-anti-CD11b, and PerCP-Cy5.5–anti-CD8. Optimal working dilutions were determined according to the relevant protocol. All antibodies were incubated for 30 min at 4 ℃. All flows were done using FACS ARIA II (BD Biosciences) and the data were analyzed using FlowJo 8.2.6 (Tree Star, Ashland, OR).
Histopathological and immunohistochemical (IHC) assay
The muscle biopsy specimens from DMD patients and non-DMD controls were freshly frozen in liquid nitrogen–cooled isopentane. The frozen muscle Sects. (8 μm) were stained with HE and pathological changes were observed under a light microscope. IHC assay was according to the previous manufacturer's suggestion . Briefly, the dry slides were preblocked in PBS containing 10% normal goat serum and incubated overnight with the primary antibodies for macrophages (rat monoclonal anti-mouse F4/80 antibody, Abcam, Cambridge, UK), CD4 positive T cells (mouse monoclonal antibody against CD4, Maxim, Fuzhou, China) or CD8 positive T cells (mouse monoclonal antibody against CD8, Maxim, Fuzhou, China), respectively. Then, the cells were rinsed and incubated with the appropriate secondary antibody (Proteintech, Wuhan, China) at 20 ℃ for 20 min. 3,3-diaminobenzedine tetrahydrocloride (Solarbio, Beijing, China) was used as chromogenic substrate. Lastly, the cells were counterstained with haematoxylin and mounted. The Ab binding was observed under a microscope.
Identification of m6A regulators
27 widely recognized m6A RNA methylation regulators were collected from published literatures. These regulators including 9 writers (CBLL1,METTL3, METTL5, PCIF1, RBM15, RBM15B, WTAP, ZC3H13, and ZCCHC4), 2 erasers (ALKBH5 and FTO), and 16 readers (ELAVL1, EIF3A, FMR1, G3BP1, G3BP2, HNRNPA2B1, HNRNPC, IGF2BP2, IGF2BP3, LRPPRC, PRRC2A, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3.). The online platform of String (cBioportal) (http://www.string-db.org/) and Cytoscape were utilized to evaluate protein–protein interaction (PPI). The correlation between m6A RNA methylation regulators was performed by R package “corrplot” (P < 0.05 as cut-off criteria).
Differential analysis of m6A regulators
Differentially expressed m6A regulators were performed by R package “limma”. Univariate logistic regression was introduced to determine m6A regulators in DMD patients (P < 0.05 as cut-off criteria). Least absolute shrinkage and selection operator (LASSO) regression was applied for minimizing the overfitting. Then the refined regulators were used to establish a predicting model. According to the coefficients obtained from the LASSO, the risk score equals the sum of coefficients and m6A regulator expression values. Receiver operating characteristic (ROC) curve analysis and the area under the ROC curve (AUC value) were finally used to evaluate the distinguishing performance.
Correlation analysis between m6A regulators and immune characteristics
Spearman correlation analyses were conducted to evaluate the relevance between m6A regulators and infiltrating immunocytes populations, immune response activity, and HLA gene expression. Heatmap was used for visualizing the results.
Unsupervised clustering for m6A regulators
By unsupervised clustering analysis, diverse m6A modification patterns were identified according to the expression profiles of m6A regulators. The cluster numbers and robustness were assessed by consensus clustering. We ran 1000 iterations of the above steps to guarantee the classification robustness with the R package “ConsensuClusterPlus”. The m6A modification patterns were further validated through PCA.
Biological pathway analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and HALLMARKS pathway were introduced to identify relevant enriched biological pathways in distinct m6A modification patterns. The expression matrix was transformed into the pathway activation score matrix through GSVA. Raw P values were corrected for multiple testing using the false discovery rate (FDR) and the thresholds were set at FDR < 0.05. Gene Ontology (GO) enrichment analysis was applied to access the major biological functions of m6A phenotype-associated genes by the R package “clusterProfiler” and adjusted P < 0.01 was considered as the cut-off criterion.
Identification of m6A regulator-mediated genes
To determine m6A regulator-mediated genes, R package “limma” was performed to identify differential expression genes (DEGs) between distinct m6A modification patterns. We overlapped the DEGs to determine the m6A phenotype-associated genes and visualized the result with Venn plot.
R (version 3.6.1) and SPSS (version 25.0) were introduced to perform data analysis and statistics. Student's t-test or Mann–Whitney U-test was carried out to compare differences between two independent groups. One-way ANOVA or Kruskal–Wallis tests were performed for the comparisons among three or more groups. Spearman correlation analysis was used to identify the relevance of gene expression. |R|> 0.25 and P < 0.05 were considered relevant and identified as statistically significant unless otherwise mentioned.
Characteristics analysis of immune microenvironment in DMD
The main immune-related biological processes and molecular functions associated with the pathogenesis of DMD were investigated by GSEA. The results suggested that many significant immune response-associated processes might be involved in the pathology of DMD, such as antigen progression and presentation, complement and coagulation cascades, and leukocyte transendothelial migration (Additional file 2: Figure S2A–F). Furthermore, we found 20/23 immune reaction related pathways significantly upregulated in DMD compared with non-DMD samples, indicating enhanced immune responses of muscle tissues in DMD (Additional file 3: Figure S3A and S3B, Tables S2). By utilizing ssGSEA method, we explored infiltrating immune cells difference between DMD and non-DMD groups and found that the extent of immune infiltration was significantly higher in DMD group (P < 0.05, Fig. 1A, Tables S3). Furthermore, we conducted FCM of the skeleton muscle in mdx and C57BL/6 mice to preliminarily validate the immune-cell infiltration status in DMD (Fig. 1B, C). As shown in Fig. 1D–H, the proportions of CD4+, CD8+ T cells and macrophages are significantly increased in mdx compared with the control group. Similarly, IHC staining for CD4+ T cell (CD4), CD8+ T cell (CD8) and macrophage (F4/80) using muscle samples of DMD patients and non-DMD controls validated the results (Fig. 1I). The staining signals for T cells and macrophages in DMD groups were notably higher than non-DMD control (Fig. 1J–L). In addition, we explored the HLA gene expression status and found that most of them were altered in DMD compared with non-DMD controls (Additional file 4: Figure S4, Tables S4).
The landscape of m6A regulators in DMD
To explore the potential role of m6A regulators in the immune environment of DMD, we evaluated the expression pattern of m6A regulators. 26 m6A regulators were selected for analysis in our work, including 9 writers, 16 readers, and 1 eraser (Additional file 7: Table S5). The results of PPI network indicated that m6A regulators had a tight association and functioned as a complex (Fig. 2A). Meanwhile, the correlation among 26 m6A regulators at transcription levels was analyzed and the results indicated a strong relationship exists between m6A regulators. Among them, G3BP2-YTHDF3 showed the most significant positive correlations with R = 0.69 (Fig. 2B). We further explored differential expression levels of 26 m6A regulators between DMD and non-DMD controls and 21 regulators were found to be significantly altered in the muscle tissues of DMD patients, including 19 down-regulated genes and 2 up-regulated genes (Fig. 2C and D). The expression values of 6 writers (CBLL1, ZC3H13, METTL5, RBM15, WTAP, and PCIF1), 12 readers (ELAVL1, HNRNPA2B1, LRPPRC, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, PRRC2A, G3BP1, EIF3A, and G3BP2), and 1 eraser (FTO) were reduced, whereas the expression values of reader FMR1 and writer RBM15B were significantly increased (P < 0.05). Among these differentially expressed regulators, PCIF1 showed the most statistically significant alteration and FMR1 showed the maximum fold-change in DMD (Additional file 7: Table S6). To elucidate the possible regulatory mechanism of m6A methylation modification, we analyzed the relationship between transcription factors and m6A regulators in DMD. Based on correlation coefficients greater than 0.6, 31 transcription factors associated with the m6A regulators were identified (Additional file 7: Table S7). As shown in Fig. 2E, there is a complex relation between m6A regulators and transcription factors. Among them, PRRC2A (reader), RBM15 (writer), and FMR1 (reader) were associated with diverse transcriptional factors; PRRC2A and RBM15 showed a significant positive correlation with most of factors, while FMR1 showed a negative correlation. In summary, our data identified aberrant expression levels of m6A methylation regulators and showed the complexity of gene regulation through m6A modification mechanisms in DMD patients.
Identification of key m6A regulators in DMD
The univariate logistic regression analysis was performed to determine critical m6A factors in DMD. Our result revealed that 22 m6A regulators were related to the development of DMD (Fig. 3A, Additional file 7: Table S8). Then, LASSO regression was introduced to avoid overfitting in the subsequent model construction (Fig. 3B, C). According to the optimum λ value, 7 genes (FMR1, FTO, G3BP1, IGF2BP3, LRPPRC, YTHDC1, and ZCCHC4) were selected as hub m6A regulators for DMD, which were then applied to construct a gene signature (Additional file 7: Table S9). Then we performed a logistic multifactor regression analysis (Fig. 3D) and calculated the diagnostic risk score of the gene signature to reveal its ability in distinguishing between normal and DMD samples. As is shown in Fig. 3E, the DMD group experienced a higher m6A risk score than the control group. The ROC curve analysis also suggested that the gene signature has a good performance in classifying the two groups (AUC = 1, Fig. 3F). Furthermore, the relationship between risk score and 26 m6A regulators in DMD samples was investigated. The risk score was negatively associated with most regulators, whereas positively linked to FRM1 and RBM15B (Fig. 3G). Moreover, a ROC curve of 7 m6A regulators was performed to estimate the accuracy of the candidate genes and the AUCs for these regulators ranged from 0.63 to 0.931 (Fig. 3H), indicating our results’ high accuracy. We selected an independent dataset GMS1004 from the GEO database to externally validate the gene expressions of key regulators. The result shows a similar tendency to the training set, which can prove the reliability of our analysis (Additional file 5: Figure S5).
The relevance between m6A regulators and immune characteristics in DMD
The relevance between immune cell infiltration and the expression values of m6A regulators was estimated through correlation analysis in the muscle samples of DMD and a significant association was found (Fig. 4A, Additional file 7: Table S10). For instance, activated CD8+ T cell abundance was positively correlated with IGF2BP3 (Fig. 4B), while activated CD4+ T cell was negatively correlated with PCIF1 (Fig. 4C). Similarly, we found the main immune-related pathways have also been linked to the expression values of m6A regulators in DMD suggesting that these immune related pathways and m6A regulators interact with each other or have a regulatory relationship (Fig. 4D, Additional file 7: Table S11). For instance, the TGF-β signaling pathway was positively associated with several m6A regulators, while the Toll-like receptor signaling pathway was negatively correlated with multiple m6A regulators. Moreover, we found the m6A reader, FMR1 and ELAVL1, were highly associated with many immune response gene sets. As seen in Figs. 4E and F, FMR1 was positively related to the TGF-β signaling pathway; in turn ELAVL1 was negatively related to cytokine receptor interaction. Besides, the relevance between m6A regulators and HLA expression was analyzed (Additional files 6 and 7 Figure S6A, Table S12). The result indicated that ZCCHC4 and HLA-DOB were the most positively correlated pair (Additional file 6: Figure S6B), but the most negatively were HNRNPA2B1 and HLA-DOA (Additional file 6: Figure S6C).
Consensus clustering of m6A regulators identified three types of patients with DMD
Consensus clustering was introduced to categorize patients with DMD into subgroups based on the expression levels of m6A regulators. With clustering stability increasing from k = 2 to k = 10, k = 3 was determined with appropriate clustering stability (Fig. 5A and B). Hence, DMD patients were clustered into three groups, including 9 samples in cluster A, 19 samples in cluster B and 28 samples in cluster C (Fig. 5C, Additional file 7: Table S13). PCA analysis further validated that the samples of DMD were separated into three non-overlapping clusters clearly (Fig. 5D). In addition, the expression differences of m6A regulators among the three cluster groups were evaluated and the distributions of m6A regulators' expression levels exhibit notable differences except for METTL3 and RBM15B (Fig. 5E, F).
Characteristics analysis of immune microenvironment in distinct m6A clusters
Infiltrated immunocytes were evaluated to characterize the immune infiltration among different m6A clusters. The proportion of 18/28 infiltrating immune cells was significantly heterogeneous in the different patterns (Fig. 6A). The abundance of most infiltrated immunocytes, including macrophages, activated CD4+ T cells, activated CD8+ T cells, and the natural killer cells were significantly higher in cluster C compared with cluster A or B. Immune response signaling pathways mediated by three clusters were also characterized (Fig. 6B). The result demonstrated that most of the immune pathways were activated in clusters B and C, while in a state of suppression in cluster A. In addition, the types of immune responses induced by cluster B and cluster C might be different. The immune reactions of ECM receptor interaction and cytokine receptor interaction were relatively more active in cluster C in contrast to TNF-α signaling via NF-κB which was stronger in cluster B. Moreover, the HLA gene expression showed a similar trend in these three patterns (Fig. 6C). Together, these data suggested the important role of m6A modification in shaping different immune microenvironments of DMD patients.
Biological characteristics of m6A modification clusters
To evaluate the biological functions in distinct m6A modification patterns, we utilized GSEA to perform pairwise comparisons of the HALLMARKS and KEGG pathways among the three clusters. According to the FDR < 0.05, representative of hallmark gene sets were enriched, and principally are oxidative phosphorylation and myc target. Additionally, significant pathways on KEGG gene sets were explored, including linoleic acid metabolism, RNA degradation, and pyruvate metabolism. According to the enrichment features of the two gene sets, cluster A was negatively correlated with inflammatory response and allograft rejection compared with cluster B or C, suggesting the low immunity (Fig. 7A–D). Additionally, the enriched pathways numbers were almost identical in cluster B and cluster C (Fig. 7E, F).
Identification of m6A regulators related genes in DMD
To elucidate the mechanisms of genes participated in m6A regulator mediated regulation, we investigated differentially expressed genes (DEGs) associated with the m6A phenotype among three clusters. In total, 225 DEGs were determined as m6A phenotype-associated genes (Fig. 8A and Additional file 7: Table S14). Then, we performed GO enrichment analysis based on these DEGs and the results are illustrated in Fig. 8B. The biological process (BP) analysis showed the process of protein catabolic, regulation of translation, and cytoplasmic translation. Cellular component (CC) analysis principally included the peptidase complex, endopeptidase complex, and proteasome complex. Molecular function (MF) analysis revealed ubiquitin-like protein ligase binding, ubiquitin protein ligase binding, and structural constituent of ribosome. Furthermore, the consensus clustering was performed based on the expression of genes associated with the m6A phenotype. By choosing a k value of 3, three different clusters of DMD patients were determined, among which, A, B, and C contained 11, 18, and 27 samples, respectively (Fig. 8C–E). Furthermore, the result of PCA showed that DMD patients in the three clusters were identifiable (Fig. 8F). We also explored the distribution of samples in different datasets, m6A cluster and m6A related gene cluster (Fig. 8G). The result revealed that the patients in m6A cluster A belong to the m6A related gene cluster A group. For patients in m6A cluster B, the majority of patients belong to the m6A related gene cluster C group, and the remaining patients belong to the A or B group. Additionally, patients in m6A cluster C are categorized into two distinct m6A related gene subcategories B and C (Additional file 7: Table S15). This further suggested the three distinct m6A modification patterns existed in DMD samples.
It is widely known that deficiency of dystrophin in DMD results in a series of symptoms, including progressive inflammatory response and muscle damage. Inflammation is the major factor that contributes to skeletal muscle fibrosis and ultimately results in progressive muscle wasting, functional disability, and reduced lifespan . Furthermore, the disordered immune microenvironment hampered the therapeutic response and clinical outcomes of many promising approaches aimed at restoring dystrophin (including stem cell transplantation, gene therapy, and exon skipping) [27, 28]. Thereby, it is biologically essential to have an in-depth understanding of the immunomodulatory mechanisms for developing new therapeutic strategies for DMD patients.
The present research first applied a comprehensive approach to characterizing the immune microenvironment in DMD skeletal muscle through bioinformatic analysis and experimental confirmation. By analyzing the expression data from DMD tissues and non-DMD muscle tissues by ssGSEA, we described the landscape of infiltrating immune cells in DMD skeleton muscle tissues and discovered the fraction of 26 immune cells was significantly altered. As the main infiltrating cells participated in dystrophin-deficient muscles , the increased infiltration of macrophages, CD4+ and CD8+ T cells was further validated by the FCM analysis in mdx and IHC in DMD patients. These results were in agreement with previous studies and suggested that both adaptive and innate immune cells may contribute to the pathogenesis of DMD [30, 31]. We also found the up-regulated immune-related pathways and altered HLA gene expression status in DMD skeleton muscle tissues, which further confirmed the complexity of immune microenvironmental changes in DMD.m6A RNA methylation is the most widespread pattern of post-transcriptional modification in eukaryotic. Previous studies have clearly shown that m6A can exert essential effects on regulating the immune system in a wide range of pathologies including tumorigenesis and viral infection [15, 17, 32]. However, to date, few studies have investigated the potential role of m6A modification in the immune microenvironment of DMD. We performed bioinformatics analyses to provide a picture of m6A modification in DMD immune microenvironment. The PPI network and expression correlation analyses revealed the close interactions among the m6A regulators which may help us to gain further insight into the regulatory mechanism of m6A modification. Of note, G3BP2 and YTHDF3 existed the most significant positive correlations. As one of the key components of stress granules (SGs), G3BP stress granule assembly factor 2 (G3BP2) is mainly mediated by the positive regulation of SG assembly and protein homooligomerization. Fu et al. confirmed the m6A-binding YTHDF proteins played an important role in SG formation. YTHDF1/3 depletion restricted SG formation and prevented the enrichment of mRNA signals in SGs . Therefore, we speculate the interaction between G3BP2 and YTHDF3 may mediate the progression of DMD by involving in related oxidative stress response pathways. Besides, we investigated the expression levels of 26 m6A regulators and most of them were altered in the DMD group compared with the non-DMD group, suggesting that m6A regulators may be relevant to the pathology of DMD.
Furthermore, we established a m6A related diagnostic signature including 7 genes (FMR1, FTO, G3BP1, IGF2BP3, LRPPRC, YTHDC1, and ZCCHC4). The ROC curve analyses revealed that the m6A signature has a good performance to discriminate between non-DMD and DMD. For m6A writer, zinc finger CCHC domain-containing protein 4 (also known as ZCCHC4) is mainly involved in 28S rRNA methylation  and relevant to the fate of core cytokines in inflammatory bowel diseases . For m6A reader, a recent study revealed that fragile X-linked mental retardation syndrome protein 1 (FMR1) knockout mice present with deficiencies in proinflammatory cytokine expression, specifically tumor necrosis factor-α expression and interleukin-6 in hippocampal . Additionally, the mutation of FMR1 in the drosophila model led to a decrease in bacterial phagocytosis . These shreds of evidence suggested that fragile X-linked mental retardation syndrome protein 1 (FMR1) could modulate the activity of immune system. It has been established that GTPase-activating protein SH3 domain-binding protein 1 (G3BP1, reader) can regulate the activation of the NF-κB pathway and type 1 interferon signaling, thereby affecting the immune response . As NF-κB signaling is regarded as a crucial signaling pathway involved in the chronic inflammation status of dystrophic muscle, we speculate that G3BP1 may contribute to cellular damage and progression of DMD and represent a potential therapeutic target for DMD. m6A reader insulin-like growth factor 2 mRNA binding protein 3 (IGF2BP3) may affect prognosis in hepatocellular carcinoma by modulating the TGF-β signaling pathway . LRPPRC, also named the leucine-rich PPR-motivated protein, is a member of the pentapeptide repeat (PPR) family. In antiviral immunity, deletion of LRPPRC expression results in increased activation of the IFN response . YTHDF1 (reader) belongs to the YTH domain family. Silencing of YTHDC1 led to increased expression of M1 phenotypic markers, enhanced production of proinflammatory cytokines, and promoted migration of microglial . For m6A eraser, fat mass and obesity associated protein (FTO) is an RNA demethylase and has been validated to participate in the regulation of muscle differentiation. New evidence is emerging that reduced FTO activity contributed to increased m6A methylation levels of IL-1β, IL-6 and TNF-α transcripts and aggravates inflammation in cardiomyocytes . Therefore, these 7 m6A regulators may contribute to the progression of immune response and further research is warranted to investigate the roles of these signatures in DMD.
Next, we characterized the potential mechanisms underlying the regulation of m6A modification in the immune microenvironment of muscle tissues in DMD to search the possible immunotherapeutic target. As expected, we found that most m6A regulators are closely associated with the infiltration of immune cells in DMD. For example, the abundance of IGF2BP3 was positively relevant to infiltration of activated CD8+ T cells and PCIF1 expression had a negative relevance to activated CD4+ T cell infiltration which is consistent with what has been previously observed in thyroid carcinoma but contrary to the findings in kidney renal clear cell carcinoma . Further research is warranted to delineate the involved biological processes. Similarly, there was a strong association between the m6A regulators and the activity of the immune pathways, suggesting a crucial role of m6A modification in the regulation of immune responses. Notably, TGF-β signaling pathway which is involved in chronic inflammatory response and fibrosis in DMD [44, 45], is associated with a variety of m6A regulators. In fact, it has been established that METTL3-METTL14-WTAP complex interacts with TGF-β pathway through the SMAD2/3 interactome . In addition, a recent study revealed that the m6A reader YTHDF3 can influence TGF-β signaling pathway by mediating peroxiredoxin 3 translation in liver fibrosis . Our data discovered the close relationship between the TGF-β signaling pathway and FMR1 in DMD which has not been reported previously and may provide a new insight into the regulation of this signaling pathway.
The role of epigenetic modifications has been increasingly appreciated based on its potentially relevant implications in identifying homogeneous groups of patients with different characters, which can advance our understanding of the pathophysiology and formulate individualized therapeutic strategies . More recently, molecular techniques such as genotyping chips and next generation sequencing (NGS) have enabled the rapid and cost-efficient studying of epitype . By analyzing the expression data from 56 DMD samples, we conducted consensus clustering and determined three DMD subtypes (clusters A/B/C) with diverse immune characteristics based on the expression of m6A regulators. We found that m6A cluster C presented the highest infiltration of immune cells and strongest immune responses than cluster A and cluster B. In addition, the immune related pathways affected by different m6A modification clusters varied greatly. The substantial differences exist in the immune microenvironment among the three clusters may lead to m6A different responses to therapy and have different outcomes. By identifying different m6A regulator-based expression patterns, it will be possible to develop more effective and targeted interventions to improve the prognosis of patients with DMD.
Although our work included a relatively large sample size by integrating GEO datasets to discover the role of m6A in the immune microenvironment of DMD, some limitations need to be considered. First, we investigated the immune cell infiltration through ssGSEA and chose the major infiltrating cell populations (macrophages, CD4+ and CD8+ T cells) for FCM and IHC validation. Further studies are still needed to characterize the infiltrating immunocytes and their exact mechanisms more thoroughly in DMD patients. Second, as our results are mainly based on the bioinformatic analysis of datasets, additional validation will likely need to be derived from experimental studies. In addition, the level of gene expression slightly differed between the training and validation dataset. This may be due to the inter-individual differences and smaller sample size in the GEO dataset, which inevitably affects the accuracy of results.
In this work, we characterized the overall landscape of the immune microenvironment in the skeletal muscle tissues of DMD and preliminary investigated the relevance between m6A regulators and DMD immune microenvironment. A diagnostic model involving seven m6A regulators was established with a well-performed risk score. Furthermore, three DMD subtypes (cluster A/B/C) were obtained with different immune microenvironmental characteristics through consensus clustering. The comprehensive analyses of the DMD m6A modification pattern may enhance our understanding of the immunomodulatory mechanisms in DMD and provide novel potential strategies for DMD therapy.
Availability of data and materials
Raw data in the study can be available from the corresponding author on reasonable request.
Duan D, Goemans N, Takeda S, Mercuri E, Aartsma-Rus A. Duchenne muscular dystrophy. Nat Rev Dis Primers. 2021;7(1):13.
Shieh PB. Emerging strategies in the treatment of Duchenne muscular dystrophy. Neurotherapeutics. 2018;15(4):840–8.
Hakim CH, Kumar SRP, Perez-Lopez DO, Wasala NB, Zhang D, Yue Y, et al. Cas9-specific immune responses compromise local and systemic AAV CRISPR therapy in multiple dystrophic canine models. Nat Commun. 2021;12(1):6769.
Nelson CE, Hakim CH, Ousterout DG, Thakore PI, Moreb EA, Castellanos Rivera RM, et al. In vivo genome editing improves muscle function in a mouse model of Duchenne muscular dystrophy. Science. 2016;351(6271):403–7.
Maffioletti SM, Noviello M, English K, Tedesco FS. Stem cell transplantation for muscular dystrophy: the challenge of immune response. Biomed Res Int. 2014;2014:964010.
Cordova G, Negroni E, Cabello-Verrugio C, Mouly V, Trollet C. Combined therapies for Duchenne muscular dystrophy to optimize treatment efficacy. Front Genet. 2018;9:114.
Rosenberg AS, Puig M, Nagaraju K, Hoffman EP, Villalta SA, Rao VA, et al. Immune-mediated pathology in Duchenne muscular dystrophy. Sci Transl Med. 2015;7(299):299rv4.
Tripodi L, Villa C, Molinaro D, Torrente Y, Farini A. The immune system in Duchenne muscular dystrophy pathogenesis. Biomedicines. 2021;9(10):1447.
Kumar A, Boriek AM. Mechanical stress activates the nuclear factor-kappaB pathway in skeletal muscle fibers: a possible role in Duchenne muscular dystrophy. FASEB J. 2003;17(3):386–96.
Spencer MJ, Montecino-Rodriguez E, Dorshkind K, Tidball JG. Helper (CD4(+)) and cytotoxic (CD8(+)) T cells promote the pathology of dystrophin-deficient muscle. Clin Immunol. 2001;98(2):235–43.
Porter JD, Merriam AP, Leahy P, Gong B, Khanna S. Dissection of temporal gene expression signatures of affected and spared muscle groups in dystrophin-deficient (mdx) mice. Hum Mol Genet. 2003;12(15):1813–21.
Zhang F, Liu H, Duan M, Wang G, Zhang Z, Wang Y, et al. Crosstalk among m(6)A RNA methylation, hypoxia and metabolic reprogramming in TME: from immunosuppressive microenvironment to clinical application. J Hematol Oncol. 2022;15(1):84.
Geula S, Moshitch-Moshkovitz S, Dominissini D, Mansour AA, Kol N, Salmon-Divon M, et al. Stem cells. m6A mRNA methylation facilitates resolution of naive pluripotency toward differentiation. Science. 2015;347(6225):1002–6.
Xiang Y, Laurent B, Hsu CH, Nachtergaele S, Lu Z, Sheng W, et al. RNA m(6)A methylation regulates the ultraviolet-induced DNA damage response. Nature. 2017;543(7646):573–6.
Deng X, Su R, Feng X, Wei M, Chen J. Role of N(6)-methyladenosine modification in cancer. Curr Opin Genet Dev. 2018;48:1–7.
Li H, Zhong Y, Cao G, Shi H, Liu Y, Li L, et al. METTL3 promotes cell cycle progression via m(6)A/YTHDF1-dependent regulation of CDC25B translation. Int J Biol Sci. 2022;18(8):3223–36.
Shulman Z, Stern-Ginossar N. The RNA modification N(6)-methyladenosine as a novel regulator of the immune system. Nat Immunol. 2020;21(5):501–12.
Fang C, He M, Li D, Xu Q. YTHDF2 mediates LPS-induced osteoclastogenesis and inflammatory response via the NF-kappaB and MAPK signaling pathways. Cell Signal. 2021;85:110060.
Liu Y, Liu Z, Tang H, Shen Y, Gong Z, Xie N, et al. The N(6)-methyladenosine (m(6)A)-forming enzyme METTL3 facilitates M1 macrophage polarization through the methylation of STAT1 mRNA. Am J Physiol Cell Physiol. 2019;317(4):C762–75.
Careccia G, Saclier M, Tirone M, Ruggieri E, Principi E, Raffaghello L, et al. Rebalancing expression of HMGB1 redox isoforms to counteract muscular dystrophy. Sci Transl Med. 2021;13(596):eaay8416.
Dang UJ, Ziemba M, Clemens PR, Hathout Y, Conklin LS, Investigators CV, et al. Serum biomarkers associated with baseline clinical severity in young steroid-naive Duchenne muscular dystrophy boys. Hum Mol Genet. 2020;29(15):2481–95.
Pescatori M, Broccolini A, Minetti C, Bertini E, Bruno C, D’Amico A, et al. Gene expression profiling in the early phases of DMD: a constant molecular signature characterizes DMD muscle from early postnatal life throughout disease progression. FASEB J. 2007;21(4):1210–26.
Khairallah RJ, Shi G, Sbrana F, Prosser BL, Borroto C, Mazaitis MJ, et al. Microtubules underlie dysfunction in duchenne muscular dystrophy. Sci Signal. 2012;5(236):ra56.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017;18(1):248–62.
Weller C, Zschuntzsch J, Makosch G, Metselaar JM, Klinker F, Klinge L, et al. Motor performance of young dystrophic mdx mice treated with long-circulating prednisolone liposomes. J Neurosci Res. 2012;90(5):1067–77.
Tidball JG, Welc SS, Wehling-Henricks M. Immunobiology of inherited muscular dystrophies. Compr Physiol. 2018;8(4):1313–56.
Charlesworth CT, Deshpande PS, Dever DP, Camarena J, Lemgart VT, Cromer MK, et al. Identification of preexisting adaptive immunity to Cas9 proteins in humans. Nat Med. 2019;25(2):249–54.
Mendell JR, Campbell K, Rodino-Klapac L, Sahenk Z, Shilling C, Lewis S, et al. Dystrophin immunity in Duchenne’s muscular dystrophy. N Engl J Med. 2010;363(15):1429–37.
Engel AG, Arahata K. Mononuclear cells in myopathies: quantitation of functionally distinct subsets, recognition of antigen-specific cell-mediated cytotoxicity in some diseases, and implications for the pathogenesis of the different inflammatory myopathies. Hum Pathol. 1986;17(7):704–21.
Mojumdar K, Liang F, Giordano C, Lemaire C, Danialou G, Okazaki T, et al. Inflammatory monocytes promote progression of Duchenne muscular dystrophy and can be therapeutically targeted via CCR2. EMBO Mol Med. 2014;6(11):1476–92.
Bhatnagar S, Kumar A. Therapeutic targeting of signaling pathways in muscular dystrophy. J Mol Med. 2010;88(2):155–66.
Liu C, Yang Z, Li R, Wu Y, Chi M, Gao S, et al. Potential roles of N6-methyladenosine (m6A) in immune cells. J Transl Med. 2021;19(1):251.
Fu Y, Zhuang X. m(6)A-binding YTHDF proteins promote stress granule formation. Nat Chem Biol. 2020;16(9):955–63.
Ren W, Lu J, Huang M, Gao L, Li D, Wang GG, et al. Structure and regulation of ZCCHC4 in m(6)A-methylation of 28S rRNA. Nat Commun. 2019;10(1):5042.
Nie K, Yi J, Yang Y, Deng M, Yang Y, Wang T, et al. A broad m6A modification landscape in inflammatory bowel disease. Front Cell Dev Biol. 2021;9:782636.
Hodges SL, Nolan SO, Taube JH, Lugo JN. Adult Fmr1 knockout mice present with deficiencies in hippocampal interleukin-6 and tumor necrosis factor-alpha expression. NeuroReport. 2017;28(18):1246–9.
O’Connor RM, Stone EF, Wayne CR, Marcinkevicius EV, Ulgherait M, Delventhal R, et al. A Drosophila model of Fragile X syndrome exhibits defects in phagocytosis by innate immune cells. J Cell Biol. 2017;216(3):595–605.
Deater M, Tamhankar M, Lloyd RE. TDRD3 is an antiviral restriction factor that promotes IFN signaling with G3BP1. PLoS Pathog. 2022;18(1):e1010249.
Chen CL, Tsukamoto H, Liu JC, Kashiwabara C, Feldman D, Sher L, et al. Reciprocal regulation by TLR4 and TGF-beta in tumor-initiating stem-like cells. J Clin Invest. 2013;123(7):2832–49.
Refolo G, Ciccosanti F, Di Rienzo M, Basulto Perdomo A, Romani M, Alonzi T, et al. Negative regulation of mitochondrial antiviral signaling protein-mediated antiviral signaling by the mitochondrial protein LRPPRC during hepatitis c virus infection. Hepatology. 2019;69(1):34–50.
Zhou H, Xu Z, Liao X, Tang S, Li N, Hou S. Low expression of YTH domain-containing 1 promotes microglial M1 polarization by reducing the stability of sirtuin 1 mRNA. Front Cell Neurosci. 2021;15:774305.
Dubey PK, Patil M, Singh S, Dubey S, Ahuja P, Verma SK, et al. Increased m6A-RNA methylation and FTO suppression is associated with myocardial inflammation and dysfunction during endotoxemia in mice. Mol Cell Biochem. 2022;477(1):129–41.
Jin MZ, Zhang YG, Jin WL, Wang XP. A pan-cancer analysis of the oncogenic and immunogenic role of m6A methyltransferase PCIF1. Front Oncol. 2021;11:753393.
Evans NP, Misyak SA, Robertson JL, Bassaganya-Riera J, Grange RW. Immune-mediated mechanisms potentially regulate the disease time-course of Duchenne muscular dystrophy and provide targets for therapeutic intervention. PM R. 2009;1(8):755–68.
Ceco E, Bogdanovich S, Gardner B, Miller T, DeJesus A, Earley JU, et al. Targeting latent TGFbeta release in muscular dystrophy. Sci Transl Med. 2014;6(259):259ra144.
Bertero A, Brown S, Madrigal P, Osnato A, Ortmann D, Yiangou L, et al. The SMAD2/3 interactome reveals that TGF beta controls m(6)A mRNA methylation in pluripotency. Nature. 2018;555(7695):256–9.
Sun R, Tian X, Li Y, Zhao Y, Wang Z, Hu Y, et al. The m6A reader YTHDF3-mediated PRDX3 translation alleviates liver fibrosis. Redox Biol. 2022;54:102378.
Giacopelli B, Zhao Q, Ruppert AS, Agyeman A, Weigel C, Wu YZ, et al. Developmental subtypes assessed by DNA methylation-iPLEX forecast the natural history of chronic lymphocytic leukemia. Blood. 2019;134(8):688–98.
Prager BC, Vasudevan HN, Dixit D, Bernatchez JA, Wu Q, Wallace LC, et al. The meningioma enhancer landscape delineates novel subgroups and drives druggable dependencies. Cancer Discov. 2020;10(11):1722–41.
We thank all study participants and laboratory staff in the Neurological Laboratory of Hebei Province. We appreciate Prof. Liang Liu of the Fourth Hospital of Hebei Medical University for his kindly technical guidance.
This research was supported by the National Natural Science Foundation of Hebei Province (H2021206223), and the medical research funding from Hebei Health and Family Planning Commission (20210038).
Ethics approval and consent to participate
The study procedures were in accordance with the ethical standards and gained the ethics approval from Research Ethics Committee of the Second Hospital of Hebei Medical University (approval number: 2022-P005; 2022-AE283).
Consent for publication
All the participants (parent, legal guardian, and/or child) signed written informed consents for the publication.
All authors report that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
: Figure S1. PCA plots for three datasets. A PCA plots for three datasets without processing. B PCA plots for three datasets with removing batch effects.
: Figure S2. Gene set enrichment analysis pathways involved in DMD. A antigen processing and presentation, B intestinal immune network for IGA production, C allograft rejection, D graft versus host disease, E complement and coagulation cascade, F leukocyte transendothelial migration.
: Figure S3. The activity differences of immune reaction gene-set between DMD patients and non-DMD controls. The overall landscape of immune reaction levels between DMD patients and non-DMD controls. B Violin diagrams showing the scores for the immune reaction levels of DMD patients and non-DMD controls. * P < 0.05; ** P < 0.01; *** P < 0.001; ns, no significance
: Figure S4. The box-plot shows the expression difference of HLA genes in DMD patients and non-DMD control. * P < 0.05; ** P < 0.01; *** P < 0.001; ns, no significance.
: Figure S5. The gene expressions of key m6A regulators were validated using an independent dataset. * P < 0.05; ** P < 0.01; *** P < 0.001; ns, no significance.
: Figure S6. The correlation between m6A regulators and HLA genes in DMD. A Heatmap showing the correlations between HLA genes and m6A regulators. B The relationship between HLA-DOB and ZCCHC4. C The relationship between HLA-DOA and HNRNPA2B1. The expression levels of genes are presented by a box plot on the right panel of B and C. * P < 0.05; ** P < 0.01; *** P < 0.001; ns, no significance.
: Table S1. Clinical background information of the human participants. Table S2. Difference on immune related gene sets between healthy and DMD samples. Table S3. Difference on immune cell infiltration fraction between healthy and DMD samples. Table S4. Difference on HLA gene expression between healthy and DMD samples. Table S5. m6A RNA methylation regulators. Table S6. Expression diversity of m6A regulators. Table S7. Correlations between m6A regulators and transcription factors. Table S8. Univariate logistic regression of m6A regulators. Table S9. LASSO regression coefficient. Table S10. Correlation between m6A regulators and immune related gene sets. Table S11. Correlations between m6A regulators and immune cell infiltration fraction. Table S12. Correlation between m6A regulators and HLA gene. Table S13. m6A modification patterns. Table S14. m6A related genes. Table S15. Change of sample attributes.
About this article
Cite this article
Han, X., Ji, G., Wang, N. et al. Comprehensive analysis of m6A regulators characterized by the immune microenvironment in Duchenne muscular dystrophy. J Transl Med 21, 459 (2023). https://doi.org/10.1186/s12967-023-04301-5
- Duchenne muscular dystrophy
- Immune microenvironment
- m6A regulators
- Diagnostic signature