Identification of m6A-Related Biomarkers in Systemic Lupus Erythematosus: A Bioinformation-Based Analysis

Introduction

Systemic lupus erythematosus (SLE) is characterized as a multifaceted and chronic autoimmune inflammatory disorder, showcasing a wide range of manifestations and presenting challenges in terms of treatment options.^1,2 Hence, the discovery of potential biomarkers in SLE via bioinformatics analysis holds immense clinical significance.

m6A, denoting methylation at the adenine base’s sixth position within RNA molecules, epitomizes the predominant modification in both coding and non-coding RNAs. As a pivotal post-transcriptional modulator, it significantly influences numerous physiological and pathological processes.³ Current studies suggest that m6A could substantially impact the progression of autoimmune diseases,⁴ indicating that m6A influences gene expression by regulating RNA metabolic processes such as splicing, stability, and translation, thereby playing a role in the pathogenesis of SLE.^5–8 METTL3 plays a key role in regulating the activation and differentiation of CD4⁺ T cells, and its downregulation in SLE patients may impact T cell-mediated immune responses.⁹ Conversely, the upregulation of circGARS, acting as a sponge for miR-19a, affects the expression of the m6A reader YTHDF2 and A20, consequently modulating the NF-κB pathway, potentially exacerbating SLE pathology.¹⁰ Additionally, the differential expression of m6A-modified long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) in SLE patients suggests their potential roles in the disease, possibly serving as new therapeutic targets for SLE.^11,12

Given the pervasive occurrence and significant influence of SLE, our research seeks to decipher the fundamental molecular mechanisms and formulate an accurate diagnostic model through merging transcriptome sequencing and machine learning algorithms. Furthermore, we investigated the causative link between MMP9 and IS utilizing an MR study. The insights garnered are anticipated to provide solid academic backing for the crafting of personalized SLE treatment and management plans, aiming ultimately at enhancing patient outcomes.

Materials and Methods

Data Collection and Processing

The workflow is illustrated in Figure 1. All the utilized data was sourced from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). We utilized microarray gene expression data from GSE50772, GSE61635, GSE72509 and GSE81622 to pinpoint key m6A regulators.¹³ The basic information of the above GEO datasets is shown in Supplementary Table S1. We applied limma-voom normalization to GSE72509 and then processed the datasets.¹⁴ The RobustRankAggreg package¹⁵ was preliminarily employed to identify differentially expressed genes. Before and after batch correction with sva package,¹⁶ an analysis using MsigDB C8 was conducted to screen and process the presence of non-PBMC cells and the correlated marker genes.^17,18

Figure 1 Flow chart of research design and analysis.

Identification of DEGs and Enrichment Analysis

In this study, the R limma package¹⁹ was utilized to pinpoint differentially expressed genes (DEGs) within the processed microarray data, with the filtering criteria established as: |log2 Fold Change (FC)| > 1 and Adjusted p value < 0.05. Genes that satisfied these conditions were designated as DEGs. Following this, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were undertaken to delve deeper into the functional annotations and pathways correlated with the pinpointed DEGs.²⁰

Identification of Potential m6A Regulators

To explore the influence of m6A regulatory genes in SLE, we performed a comprehensive analysis encompassing all samples from the combined dataset. Employing the linear modeling technique available in the R limma package, we discerned significant m6A regulatory genes showcasing statistically noteworthy disparities in expression levels between SLE samples and healthy controls (HCs) (p < 0.05). For a visual representation of these discrepancies, boxplots were crafted using the ggpubr package, and heatmaps were generated with the pheatmap package.

Formulation of RF, SVM, and Nomogram Models Related to m6A Regulators

Classifiers utilizing RF and SVM were developed to pinpoint characteristic m6A regulatory genes among significant genes and to project the occurrence of SLE. m6A regulatory genes exhibiting differential expression were discerned from datasets GSE50772, GS61635, GSE72509, and GSE81622. The precision of the models was evaluated through “boxplots of residual”, “reverse cumulative distribution of residual”, and “receiver operating characteristic (ROC)” curves. The RF model, crafted with the randomForest package²¹ in R, gauged the importance of pivotal m6A regulatory genes using the chosen optimal ntree. To forecast SLE incidence utilizing feature m6A regulatory genes, a nomogram model was devised utilizing the rms package²² in R. A calibration chart facilitated a comparison of predicted and actual values. Additionally, we executed a decision curve analysis (DCA) and plotted a clinical impact curve to evaluate the benefits of the model for patients.

Construction of Weighted Gene Co-Expression Network Analysis

In our research, the R “WGCNA” package^23,24 was utilized to craft a co-expression network focusing on DEGs. To structure this network, a weighted adjacency matrix was formulated, ensuring the soft thresholding parameter accentuated strong gene associations and minimized weak ones. This matrix was subsequently converted into a topological overlap matrix (TOM) to gauge gene network connectivity. The TOM was instrumental in determining dissimilarity among neighboring genes, recognizing those with analogous expression patterns through average linkage hierarchical clustering. The ensuing clusters, distinguished by varied colors in the cluster tree, signified gene modules. By examining module interconnections and evaluating the interplay between gene modules and phenotypes, modules linked to clinical attributes were pinpointed.

Screening of Candidate Hub Genes

To predict SLE status, three machine learning methodologies were utilized: RF, LASSO, and SVM-RFE, denoting Random Forest, amalgamates multiple decision trees, providing predictions as an ensemble learning strategy, and notably manages classification and regression issues while showcasing resistance to noise and overfitting. Moreover, it provides a ranking of feature importance and adeptly manages absent data. LASSO, a regularization approach, introduces a penalty component to the linear regression’s loss function. It executes variable selection and coefficient shrinkage, effectively nullifying the coefficients of irrelevant variables to circumvent overfitting. SVM-RFE amalgamates SVM with backward elimination, training an SVM model with a comprehensive feature set, then iteratively discarding the feature with minimal weight until the desired feature quantity is attained. This strategy enables SVM-RFE to meticulously select the most pertinent features for classification, thereby enhancing the precision and efficacy of SVM.

Human Blood Specimen Collection

This study recruited SLE patients aged 18–65 years, including both newly diagnosed individuals in outpatient or inpatient settings and existing patients with exacerbations requiring hospitalization. Peripheral blood samples were collected from 32 diagnosed SLE patients and compared with samples from 20 healthy volunteers, who were enlisted from the Dermatology Departments of the First and Second hospital of Army Medical University. Disease activity was assessed using the SLEDAI score. Exclusion criteria included patients with controlled SLE after prolonged treatment (over 3 months), those currently on targeted medications such as rituximab, individuals with a lack of legal capacity, an inability to understand study risks and follow-up requirements, a history of substance abuse, and complex medical conditions like malignancies or blood disorders. All participants provided informed consent, in accordance with ethical standards. The separation of peripheral blood mononuclear cells (PBMCs) from peripheral blood was accomplished through Ficoll-Paque density gradient separation, with both PBMCs and serum subsequently being cryopreserved for future use. Compliance with the Declaration of Helsinki was ensured for all procedures involving human participants. Moreover, all experimental protocols received approval from the Ethics Committee at the First Affiliated Hospital of Army Medical University, Chongqing, China (approval no. KY2020317).

Real-Time Quantitative PCR Assay

RNA isolation was performed utilizing RNAiso Plus reagent (Takara, Shiga, Japan), and the generation of mRNA first-strand cDNA was executed with a PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara). Subsequent real-time PCR analysis was conducted employing TB Green Premix Ex Taq II (Takara) on a Bio-Rad real-time PCR apparatus (California, USA). GAPDH housekeeping gene primers (Sangon Biotech, Shanghai, China) were employed to assure normalization of expression data. Relative gene expression was evaluated through relative quantification, with the primer sequences being detailed in Table 1.

Table 1 Sequence of PCR Primers

MeRIP-qPCR

All procedures were conducted following the guidelines outlined in the EpiGentek MeRIP kit manual (EpiGentek, P-9018). Total RNA was extracted from cells or tissues and DNase I treatment was performed to remove any DNA contamination. The RNA concentration and quality were assessed using NanoDrop. Around 400 ng of total RNA was subjected to cleavage and removal of non-target m6A sequences using the Cleavage Enzyme Mix and Cleavage Buffer, incubated at 37 °C for 30 min. The cleaved RNA fragments were subsequently immunoprecipitated with 5 µg of anti-m6A antibody and Immunoprecipitation Buffer, with incubation at 25 °C for 1 h. Immunoprecipitation with Affinity Beads followed at 4 °C for 5 h. Affinity Beads were washed thoroughly with Wash Buffer to remove non-specifically bound RNA fragments and impurities. The bound m6A RNA fragments were eluted from the Affinity Beads using Elution Buffer, incubated at 55 °C for 15 min. The eluted m6A RNA fragments were then purified with RNA Binding Beads and RNA Purification Solution, incubated at room temperature for 10 min, followed by elution from the RNA Binding Beads with Elution Solution at room temperature for 2 min. The concentration and quality of the purified m6A RNA fragments were measured using NanoDrop to ensure there was no degradation or contamination. The expression levels of genes of interest within the m6A RNA fragments were analyzed using RT-qPCR and compared with input RNA or non-immunoprecipitated (IgG) controls.

Functional Enrichment Analysis of MMP9

GSVA analysis, a variant of the GSEA algorithms, explored variations in biological pathways between different pattern clusters based on enrichment scores. The GSVA R Package²⁵ was utilized to conduct functional enrichment analysis on SLE samples within the merged dataset, yielding enrichment pathways. Analysis was performed using “c2.cp.kegg.Hs.symbols.gmt”, with an adjusted p value < 0.05 indicating a statistical difference between distinct clusters. For immune function exploration, samples from the control group were excluded, and repetitive samples underwent average processing, discarding rows with an average expression value of zero. Immune function gene sets were sourced from the “immune.gmt file”, and samples were analyzed using the ssGSEA method to determine the enrichment level of immune function gene sets in each one. Raw scores for immune function were normalized between 0 and 1 for subsequent analysis. Samples were bifurcated into two groups based on MMP9 expression level, categorizing samples above the median as “High” and those below as “Low”. Box plots illustrated the score distribution of various immune functions in the two sample groups, with each plot representing a specific immune function, the box delineating the interquartile range, and the lines depicting the full score range.

Mendelian Randomization

All the Data were sourced from an accessible database, and a two-sample Mendelian Randomization (MR) was deployed to explore any causal links between the hub gene and the susceptibility to SLE, defining SNPs as instrumental variables (IVs). The data for the hub gene was extracted from publicly accessible Genome-Wide Association Study (GWAS) resources. We focused on the m6A-related gene, MMP9, and IS as an illustrative complication of SLE for MR analysis. And specific datasets were employed for both the exposure and outcome variables. The detailed information on the MMP9 exposure dataset is accessible at the GWAS catalog with the dataset ID ebi-a-GCST90010240, available at https://gwas.mrcieu.ac.uk/datasets/ebi-a-GCST90010240/. Correspondingly, the dataset for the outcome variable, Ischemic Stroke, is listed under the ID ebi-a-GCST006908, which can be found at https://gwas.mrcieu.ac.uk/datasets/ebi-a-GCST006908/. These datasets were meticulously chosen to enhance the precision and relevance of our MR analysis. The analysis of MR was facilitated using the TwoSampleMR package²⁶ in R. Data IDs for the designated exposure factor and outcome were pinpointed. Considering the limited number of SNPs associated with the specific gene at the standard GWAS threshold (5e-08), we opted for a p value threshold of 5e-05 to ensure an adequate number of SNPs as IVs. All selected SNPs had an R² greater than 0.01 and an F-statistic greater than 10, indicating their sufficient strength as effective IVs. Furthermore, the robustness of the results was validated through MR Egger regression analysis. The corresponding SNP information from the GWAS data for the outcome was retrieved using the list of instrumental SNPs for the exposure factor. The next step involved merging the GWAS data for the exposure factor and outcome, aligning effect alleles to ensure a uniform direction of effect for the same SNP in both datasets. MR analysis was then performed using various methods including “IVW” (Inverse Variance Weighting), “MR Egger”, “Weighted”, “Weighted Median”, and “Simple Mode” to gauge the causal effect between the exposure factor and the outcome, producing confidence intervals and p values in the process. We proceeded with heterogeneity analysis and pleiotropy tests and created scatter plots, forest plots, and funnel plots for visual representation of the MR results. To conclude, a leave-one-out sensitivity analysis was executed to evaluate the stability of the findings.

Statistical Analysis

Disparities in gene expression between individuals with SLE and HCs were analyzed using Student’s t-test. The discriminative power of key genes was assessed through ROC curve analysis and AUC computation using the pROC package²⁷ in R. Statistical evaluations were performed employing R software (version 4.2.1) and GraphPad Prism (version 9.0.0).

Results

Identification of DEGs Between SLE Patients and HCs

In a retrospective analysis of four microarrays - GSE50772, GSE61635, GSE72509, and GSE81622 - a collective total of 289 SLE patient samples and 93 HC samples were evaluated. Utilizing the limma R software package post batch-effect mitigation, 171 DEGs were pinpointed from the metadata: with 86 genes being notably up-regulated and 85 genes distinctly down-regulated in SLE patient samples relative to HCs (Figure 2A and B).

Figure 2 Genes differentially expressed between the two cohorts. (A) Heat map illustrating the differential gene expression analysis from the amalgamated dataset following the mitigation of batch effects. (B) Volcano plot depicting the differential analysis results of the data set. (C) GO enrichment analysis of potential hub genes. (D) KEGG pathways related to candidate hub genes.

Functional Enrichment Analysis

The clusterProfiler package²⁸ in R was utilized to execute GO and KEGG analyses on the 171 identified DEGs. The findings revealed that the DEGs primarily participated in the response to virus, vesicle lumen, and olfactory receptor activity (Figure 2C). Moreover, KEGG analysis unveiled several notable signaling pathways, including Olfactory transduction, Coronavirus disease (COVID-19), and Influenza A (Figure 2D).

Identification of m6A Regulators in SLE

Seventeen m6A regulatory genes were identified, encompassing six writers, nine readers, and two erasers. Disparities in the expression levels of these seventeen genes between SLE patients and HCs are showcased (Figure 3A). Twelve genes exhibited statistically significant expression deviations (p < 0.05). Specifically, YTHDF3, IGFBP2, and IGFBP3 were upregulated in SLE, while METTL14, ZC3H13, RBM15B, YTHDF1, YTHDF2, HNRNPC, IGF2BP1, FTO, and ALKBH5 were downregulated. The heatmap for these twelve pivotal genes is depicted (Figure 3B).

Figure 3 The m6A regulatory genes in samples from SLE and HC. (A) The boxplot illustrating differential expression of 17 m6A regulatory genes between SLE patients and HC. (B) Heat map displaying differential expression analysis of the combined dataset following the removal of batch effects. (*p < 0.05, ***p < 0.001).

Construction of RF, SVM and Nomogram Models

The RF model showcased commendable performance, highlighted by a diminished residual value in comparison to the Support Vector Machine (SVM) model (Figure 4A and B). This superior performance of the RF model can be ascribed to its proficiency in managing intricate interactions and non-linear associations among variables, as well as its adeptness at mitigating overfitting and bias by averaging outcomes across decision trees. Thus, we opted for the RF model for ensuing analyses (Figure 4C). Upon ranking the notable m6A regulatory genes based on their importance scores, it was observed that all genes held scores exceeding 5 (Figure 4G), making them eligible for inclusion within the nomogram model (Figure 4H). In the nomogram, each gene was scored individually, and the comprehensive score to predict SLE incidence was computed by tallying these scores. The calibration curves (Figure 4D) illustrated a tight congruence between the solid and dotted lines, signaling the model was well-calibrated. Moreover, the decision curve analysis (DCA) (Figure 4E), demonstrated a notable divergence of the red line, representative of the m6A genes, from the baseline gray and black lines. This divergence underscores the potential utility of m6A genes in contributing to the prognostic assessment of SLE. Complementing this, the clinical impact curve (Figure 4F) provided additional support for the relevance of these genes, indicating their significance in the broader context of SLE pathophysiology.

Figure 4 The RF model and SVM model. (A) Boxplots of residuals for the RF model. (B) Residual distribution for RF and SVM models. (C) Error levels for SLE groups, HC groups, and all samples. (D) Calibration curves for the nomogram model. (E) Decision curve showing feasibility of the model. (F) Clinical impact curve demonstrating model’s potential clinical benefit. (G) Bubble plot illustrating gene importance. (H) Nomogram model predicting incidence based on m6A regulatory genes.

Identification and Validation of Potential Biomarkers

The LASSO was employed to scrutinize the 171 DEGs for potential biomarkers of SLE, revealing 34 genes of potential diagnostic relevance (Figure 5A). The SVM-RFE algorithm also analyzed the 171 DEGs and yielded 39 genes (Figure 5B–D). Moreover, Random Forrest identified 31 feature genes (Figure 5E and F). The WGCNA was constructed based on the DEGs. A total of 19 modules containing 348 hub genes were identified (Figure 6A–C). Only four genes, EPSTI1, USP18, HP, MMP9 were common to all approaches (Figure 6D), suggesting their importance in SLE diagnosis. We executed ROC analysis of the four crucial genes and the results showed that EPSTI1 (AUC=0.888), USP18 (AUC=0.870), HP (AUC=0.817), MMP9 (AUC=0.753) (Figure 6E).

Figure 5 Selection of potential biomarkers for SLE. (A) Screening the tuning feature in LASSO model. (B) RMSE (Root Mean Squared Error) is a metric used for evaluating the performance of the SVM-RFE model. A smaller RMSE value indicates a better predictive capability of the model. (C and D) The feature selection in SVM-RFE model aims to identify the optimal subset of features to enhance the model’s generalization capability and prediction accuracy. (E) The performance of the random forest model varies with the number of trees, enabling the optimization of model settings to improve prediction accuracy while avoiding overfitting or underfitting. (F) Bubble plot of gene importance obtained from the random forest tree model.

Figure 6 Selection of diagnosis marker candidates for SLE. (A) Gene modules associated with SLE. (B) Clustering of all genes in the merged dataset based on a topological overlap matrix (1-TOM), yielding co-expression modules depicted in various colors. (C) A scatter plot highlighted the yellow module’s robust positive correlation with SLE in the merged dataset. (D) A Venn diagram showcased four diagnostic biomarkers shared among LASSO, SVM-RFE, RF, and WGCNA. (E) ROC curves were employed to assess the diagnostic relevance of hub genes in SLE.

MMP9 Associated with m6A Modification in SLE

The expression levels of EPSTI1, USP18, HP and MMP9 were all significantly up-regulated in SLE compared to HCs (Figure 7A and B). The increased methylation level of MMP9 RNA was confirmed by MeRIP qPCR while the others showed no significance (Figure 7C).

Figure 7 The m6A modification level and gene expression level of hub genes. (A) Violin plots illustrating the expression profiles of four hub genes. (B) Validation of the expression levels of four hub genes in SLE using RT-qPCR. (C) Verification of the m6A modification levels of four crucial genes in SLE through MeRIP-qPCR. (*p < 0.05, **p < 0.01, ***p < 0.001, **** p < 0.0001). Abbreviation: ns, no significance.

To characterize the immune cell composition of the two cohorts, CIBERSORT was employed, revealing statistically significant differences in macrophages and CD8⁺ T cells (Figure 8A and B). This finding suggests a correlation between MMP9 expression and immune cell infiltration, implying that MMP9 may play a role in regulating immune cells in SLE patients (Figure 8A–C).

Figure 8 MMP9 is related to immunocyte infiltration levels. (A) Correlation map illustrating the relationship between MMP9 expression and macrophages. (B) Correlation map illustrating the relationship between MMP9 expression and CD8+ T cells. (C) Lollipop plot visualizing the correlation between MMP9 and immune cells. Text highlighted in red indicates p values that are less than 0.05, denoting statistical significance.

Functional Analysis of Different Molecular Subtypes

GSVA enrichment analysis found that gene sets pertaining to MMP9 grouping predominantly encompass pathways of immunological significance and metabolic relevance (Figure 9A). In addition, the immunological functional analysis boxplot reveals that in various immune cell types or functions, such as aDCs (activated dendritic cells). A significant difference in scores is evident between the high and low gene expression groups. (Figure 9B).

Figure 9 Functional enrichment analysis of MMP9. (A) Pathway enrichment analysis and a bar graph representation of DEGs in different subtypes were performed using GSVA. (B) Box plot analysis of immune functions illustrates the distribution of scores for various immune functions, based on the grouping of samples by the expression levels (high and low) of the target gene MMP9. Each box plot represents a specific immune function, with the box itself depicting the interquartile range of the scores and the whiskers extending to show the full range of the scores. (*p < 0.05, **p < 0.01, ***p < 0.001).

MMP9 is Causally Associated with Ischemic Stroke

Ischemic stroke is one of the significant complications associated with SLE.^7,8 Recent studies on MMP9 and ischemic stroke has provided valuable insights into the complex interplay between MMP9 and various aspects of ischemic stroke, including blood brain barrier (BBB) dysfunction, cognitive impairment, therapeutic targets and pathophysiological changes.^9–11 The SNPs characteristics of MMP9 and Ischemic stroke were considered robust instrumental variables (Supplementary Table S2). Figures 10A and C illustrate the individual genetic variations’ causal effects on SLE. Utilizing the IVW method, we identified a significant association between MMP9 and IS, with an odds ratio (OR) of 1.0134 (95% CI = 1.0004–1.0266, p = 0.0440). However, the MR-Egger method did not reveal a statistically significant effect (OR = 1.0173, 95% CI = 0.9951–1.0400, p =0.1314) (Supplementary Table S3). The funnel plot representing the causal effect displayed approximate symmetry (Figure 10B), and the MR Egger regression’s intercept indicated no horizontal pleiotropy (p = 0.6725), further affirming the absence of biased causal effects (Supplementary Table S4). The heterogeneity tests yielded a Q statistic of 72.8263 with 70 degrees of freedom (p = 0.3852) for MR Egger and a Q statistic of 73.0138 with 71 degrees of freedom (p = 0.4117) for the IVW method, indicating no significant heterogeneity in the causal effects (Supplementary Table S5). Figure 10D presents the systematic MR analysis conducted on the remaining SNPs after each SNP was excluded. The consistent results suggest that the causal significance remains intact for all SNPs, indicating the absence of a dominant SNP influencing MMP9 levels and IS. These findings validate the previous MR results.

Figure 10 Two sample Mendelian randomization study results. (A) Scatter plot depicting the causal impact of MMP9 on the susceptibility to ischemic stroke. (B) Funnel plots utilized to assess the overall heterogeneity of MR estimates for the influence of MMP9 on ischemic stroke. (C) Forest plot showcasing the causal effect of individual SNPs on the risk of ischemic stroke. (D) Leave-one-out plot employed to visualize the causal effect of MMP9 on ischemic stroke risk when excluding one SNP at a time.

Discussion

SLE is a multifaceted autoimmune disease characterized by a wide array of clinical manifestations and constrained therapeutic strategies.^1,2 The exploration of alternative mechanisms contributing to SLE is of paramount importance. Recent investigations have underscored the profound influence of epigenetics, specifically m6A RNA modification, in the etiology of SLE.^7,8,12

Within the immune system’s framework, m6A plays an indispensable role in orchestrating the activation, differentiation and function of both innate and adaptive immune cells, including T cells,^29–31 B cells³² and dendritic cells,^33,34 etc. Moreover, m6A is implicated in DNA repair and immune tolerance, processes intrinsically linked to the genesis of autoantigens and autoimmune responses.³⁵ With the advent of cutting-edge molecular biology techniques, the scrutiny of m6A-associated genes in SLE’s pathogenesis has emerged as a focal research domain, targeting the discovery of innovative and efficacious therapeutic interventions.

Diagnostic biomarkers serve as vital instruments for discerning pathological alterations in patients. Contemporary research is concentrated on authenticating biomarkers explicitly correlated to SLE.³⁶ Zhao et al indicated that the methylation status of the IFI44L promoter serves as a potent discriminator between patients with SLE, healthy individuals and those with other autoimmune disorders.³⁷ Recent studies employed bioinformatics to pinpoint DEGs and appraise their diagnostic accuracy through ROC curve analysis, accentuating IFI44 as a propitious diagnostic biomarker.^38–41 Concurrently, other inquiries utilizing DEG and WGCNA analyses have identified immune-related pathways associated with SLE, with IFI27 manifesting potential diagnostic capabilities.³⁸ Subsequent research has unveiled a positive association between TNFSF13B and OAS1 in SLE and metabolic syndrome, elucidating insights into comorbidity and prospective combined treatments for SLE with metabolic syndrome.⁴² Additional studies have laid the foundation for probing key genes (SPI1, MMP9, C1QA, CX3CR1, MNDA) linked to atherosclerosis in SLE patients.⁴³ Furthermore, bioinformatics explorations have discerned DEGs predominantly connected to diverse viral infections, with HMGB1, ITGB2 and CREB5 as central genes.⁴⁴

Our findings denote significant enrichment of DEGs in the olfactory transduction pathway, while the majority of other genes are affiliated with viral infections. The olfactory transduction pathway encompasses odor molecule binding, ion channel opening, neuron depolarization and signal transmission, coupled with sensitivity to environmental fluctuations. Recent research has correlated olfactory transduction with SLE, including alterations in olfactory function in SLE patients and its association with olfactory dysfunction.⁴⁵ Integrative models of N6-methyladenosine immunity in SLE have demonstrated a positive correlation between olfactory transduction and the pathway involving neuroactive ligand-receptor interactions.⁵ Additionally, longitudinal studies on SLE and systemic sclerosis have associated olfactory dysfunction in SLE with variables such as age, disease activity, heightened anxiety and depression symptoms, suggesting potential diagnostic and therapeutic implications.⁴⁶ Nevertheless, comprehensive research is requisite to fully comprehend the exact relationship and underlying mechanisms between the olfactory transduction pathway and SLE.

In the pursuit of potential biomarkers, we intersected key genes derived from Random Forest, SVM-RFE, LASSO and WGCNA analyses. MeRIP-qPCR validation corroborated the association of m6A modification with MMP9 among these genes. Additionally, the RM2target database predicts interactions between MMP9 and ALKBH5, FTO, METTL3, METTL14 and YTHDF2, substantiated by in vitro experimental evidence.⁴⁷

MMP9, an enzyme proficient in degrading extracellular matrix and immune complexes, is augmented in SLE patients and intimately associated with disease activity and organ damage.⁴⁸ The pathogenesis of MMP9 encompasses multifarious mechanisms, including self-antigen degradation, immune complex clearance, cell apoptosis regulation, cell adhesion and migration modulation, angiogenesis participation and inflammatory responses.^49–51 MMP9 also facilitates the release of NETs, which can instigate autoimmunity and tissue impairment.⁵¹ Additionally, the escalation of MMP9 expression in the kidneys of SLE patients is correlated with lupus nephritis development and progression.^52,53 Moreover, MMP9 is connected to atherosclerosis and myocardial injury in the cardiovascular system of SLE patients.⁴³ MMP9’s ramifications extend to the nervous system, with associations to neuro-psychiatric symptoms through blood-brain barrier degradation and immune cell infiltration.^54–56 These multifaceted impacts accentuate MMP9’s significance as a prospective therapeutic target for managing immune dysregulation, tissue damage and organ complications in SLE.

Recent investigations have elucidated the multifaceted connection between SLE and Ischemic Stroke. These studies encompass the correlation between the genetic predisposition to SLE and the concomitant risk of IS, the manifestation of neurological complications within the SLE patient population and the methodologies employed in the therapeutic intervention of IS in the context of SLE.^57–60 MMP9’s activation during the acute phase can precipitate blood-brain barrier disruption, brain edema formation and hemorrhagic transformation.^56,61,62 Studies intimate a correlation between MMP9 levels and the risk of severe brain edema, hemorrhagic transformation and adverse functional outcomes.⁶³ SNPs in the MMP9 gene may also influence the occurrence and prognosis of IS. Some studies have found that the genetic polymorphism of MMP9 is related to the risk and prognosis of IS.^64–66 For example, the MMP9 rs17576 polymorphism is associated with an increased risk of IS in the Han Hakka population and interacts with the MMP12 rs660599 polymorphism, further increasing the risk of IS.⁶⁷ In addition to the genetic level, MMP9 is also closely related to IS at the protein level. Some studies have found that the level of MMP9 in the serum or plasma of IS patients is elevated and is related to the severity and adverse outcomes of IS.^68–70 Other studies have found that there is an elevation of MMP9 in choroid plexus epithelial cells (CPECs) of IS patients, leading to the degradation of tight junction proteins (TJs) and disruption of the blood-cerebrospinal fluid barrier, thus promotes the infiltration of immune cells and neuronal death.⁷¹ In summary, MMP9 has a complex and important relationship with IS and can affect the occurrence and development of IS from multiple aspects. Some therapeutic strategies targeting MMP9 have been validated in animal models or clinical trials, such as using MMP inhibitors or neutralizing antibodies to reduce the activity or expression of MMP9.⁶² However, these strategies still have some limitations and uncertainties and need further optimization and validation. In our study, a two-sample MR analysis employing GWAS data of MMP9 and IS identified a causal association between MMP9 expression levels and an augmented risk of IS, furnishing more reliable outcomes compared to conventional observational studies.^72,73

In our exploration of the connections between MMP9, m6A modifications, IS, and SLE, we have unveiled a complex web of interrelated factors and influences among these elements. The role of MMP9 in SLE is particularly crucial, as it is not only involved in immune cell migration and inflammatory responses but is also closely related to the prognosis of ischemic stroke. Furthermore, specific genetic variations in MMP9 have been linked to an increased risk of IS, a common complication in SLE, underscoring a potential connection between MMP9 and the co-occurrence of these two conditions.

m6A modifications are pivotal in SLE, particularly through their impact on the expression of hub genes, including MMP9. The m6A modification of MMP9 plays a critical role, indirectly influencing the occurrence of IS in patients with SLE. Changes in this epigenetic mechanism might modulate inflammatory responses and vascular stability, thereby impacting the risk of IS.

Our research suggests that alterations in m6A modifications and the upregulation of MMP9 within the context of SLE could jointly affect inflammatory processes, and the upregulation of MMP9 might be associated with the risk of IS in SLE patients. These insights provide a new perspective on understanding SLE and its complications, offering valuable guidance for future therapeutic strategies. By considering these factors collectively, our analysis not only deepens the understanding of the links between SLE and IS but also directs future research and clinical practice.

Despite the significant discoveries made in our study, there are several limitations. Firstly, there is a scarcity of GSE datasets directly relevant to the pathogenic mechanisms of m6A modification in SLE. This limits our ability to deeply explore the role of m6A modification in SLE. Secondly, the relatively small sample size used in our study may impact the generalizability and reliability of our findings. Additionally, potential biases in sample selection suggest that our results might not be fully applicable to a wider population. Moreover, our machine learning models were mainly used for feature selection and preliminary biomarker screening, rather than for complex predictive or classification tasks. Finally, while advanced bioinformatics techniques and statistical methods were employed, our analysis was confined to gene expression data. Future studies incorporating other types of data, such as protein expression and clinical data, could provide more comprehensive insights.

Conclusion

In conclusion, the m6A regulatory genes may exert a significant influence on the pathogenesis of SLE. In the context of SLE, the expression of MMP9 is subject to regulation by m6A modification and exhibits a significant correlation with ischemic stroke associated with SLE. This observation may contribute to the advancement of pre-symptomatic diagnostic methodologies and furnish additional insights into the underlying molecular mechanisms, thereby enhancing the comprehension of risk genes pertinent to SLE.