SLE-related public data were first downloaded from the Gene Expression Omnibus (GEO) database (19, 20), and seven datasets were collected: three datasets from PBMC samples [GSE121239 (GPL13158, Normal: 20, SLE: 292) (21, 22), GSE11907 (GPL96, Normal: 12, SLE: 110) (23), GSE81622 (GPL10558, Normal: 25, SLE: 30) (24)]; four datasets from whole blood samples [GSE65391 (GPL10558, Normal: 72, SLE: 924) (25), GSE100163 (GPL6884, Normal: 14, SLE: 55) (26–28), GSE45291 (GPL13158, Normal: 20, SLE: 292) (29), GSE49454 (GPL1261, Normal: 20, SLE: 157) (30)]. GEO2R was used for differential expression analysis of the GSE121239, GSE11907 and GSE81622 datasets. |log2 FC| > 0.5 and adj.P.Val <0.05 were screened as differentially expressed genes (DEGs). Information on the datasets was displayed in Supplementary File S1 , and the UMAP plots of the datasets were shown in Supplementary Figure S1 .

Five methods were used for the screening of biomarkers of SLE: WGCNA (31), SVM-RFE and RF (32, 33) based on DEGs results from the GSE121239 dataset, SLEDAI correlation analysis, and the difference analysis of two SLE datasets (GSE11907, GSE81622). Gene module classification and identification of DEGs was performed using the R package “WGCNA” (34, 35). First, the samples were clustered to assess whether there were any significant outliers. Second, the automatic network construction function was used to construct the co-expression network. A soft threshold β is calculated. Third, hierarchical clustering and dynamic tree cutting functions were used to detect modules. Finally, gene importance (GS) and module membership (MM) were calculated. Linking modules to clinical traits. The corresponding module extracted the corresponding module gene information for further analysis. The SVM classifier was constructed using the R package “e1071” and cross-validated to eliminate regression features. The “randomForest” package was used for the RF classification process. For correlation analysis, the Spearman correlation coefficient was applied, and GEO2R was used for DEGs analysis of the SLE dataset (GSE11907, GSE81622). Then, overlapping genes were selected from the above five methods for further analysis in this study. The specific screening conditions for the analytic methods were shown in Supplementary File S1 .

A total of 47 samples, including 23 SLE samples and 24 healthy controls, were collected from the First Affiliated Hospital of USTC (Anhui Provincial Hospital). This study was approved by the ethics committee of the First Affiliated Hospital of USTC, and informed consent forms were signed by the patients. Peripheral blood mononuclear cells (PBMCs) were isolated from EDTA anticoagulated by FicollPaque density gradient centrifugation. Then, the cells were cultured in RPMI-1640 medium. Cells were maintained in a humidified incubator at 37°C with 5% CO₂.

RNA was extracted using total RNA isolation reagent (Epizyme, Shanghai, China). A cDNA synthesis kit (Bioer, Hangzhou, China) was used for reverse transcription of total RNA. GAPDH (forward: 5’-GTCTCCTCTGACTTCAACAGCG-3’, reverse: 5’-ACCACCCTGTTGCTGTAGCCAA-3’), MX2 (forward: 5’-CACCGAGCTAGAGCTTCAGGA-3’, reverse: 5’-CCGGGAAGGTCAATGATGGT-3’), and GAPDH was subsequently processed as internal references. Relative expression was calculated by the 2^-ΔΔCt method. P <0.05 was considered statistically significant.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using the R package “clusterProfiler” (36) to explore the function and pathways of the screened genes. Gene set enrichment analysis (GSEA) (https://www.gsea-msigdb.org/gsea/index.jsp) was used to observe changes in pathways and biological processes in gene expression datasets. Additionally, the single-sample GSEA (ssGSEA) algorithm was utilized for the calculation of core pathway activity scores for each sample based on gene expression levels.

Four algorithms, ssGSEA, CIBERSORT, xCell and ImmuCellAI, were applied to perform immune cell infiltration analysis. The ssGSEA was performed using the R package “GSVA” to calculate the relative infiltration levels of 28 immune cell types (37). The xCell algorithm was implemented using the R package “xCell” to calculate the fraction of 64 immune cells in the sample. The CIBERSORT algorithm was applied to calculate the proportion of 22 immune cells in all samples. The standardized gene matrix was also used to estimate the abundance of immune cells in all samples through the ImmuCellAI website (http://bioinfo.life.hust.edu.cn/ImmuCellAI/#!/) (38). The results of immune infiltration were visualized by the R package “ggplot2” and the bioinformatics online website (http://www.bioinformatics.com.cn/). Violin plots were drawn to elucidate the expression levels of individual immune infiltrating cells among different subgroups. Correlation heatmaps were drawn to illustrate Spearman correlations of biomarkers or pathways with immune infiltrating cells. The CellMark database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/index.jsp) was applied to query neutrophil-related marker genes (39). Based on the HPA database (40, 41) (https://www.proteinatlas.org/), transcriptional data and single-cell data were used to analyze the expression of biomarkers in immune cells.

The STRING database (https://cn.string-db.org/) was utilized to analyze the protein interaction relationships of biomarkers (42). The Cistrome DB database (http://cistrome.org/db/#/) was used to find transcription factors of the genes corresponding to these associated proteins (43). The DREIMT database (http://www.dreimt.org/) was applied to predict potential therapeutic drugs for SLE (44).

All data processing and analysis were performed by R software (version 4.1.0) and GraphPad Prism (version 8.0.2). Differences between two groups of continuous non-normally distributed variables were analyzed using the Mann‒Whitney U test (Wilcoxon rank sum test). The Kruskal‒Wallis test was used for comparison of variables among multiple groups. P <0.05 was considered to be statistically significant.

To identify potential biomarkers of SLE, five methods were applied. In the WGCNA algorithm, the sample clustering map was shown in Supplementary Figure S2 . The soft threshold power was chosen as 30 when 0.85 was used as the correlation coefficient threshold ( Figure 2A ). The WGCNA network clustering analysis divided 2512 DEGs (GSE121239) into 12 modules, and then we performed correlation analysis on all modules and found eight modules closely associated with disease occurrence, with a total of 2181 genes ( Figure 2B ). The clustering dendrogram of DEGs associated with SLE was shown in Figure 2C . In total, 1646 genes were identified as key candidate biomarkers for SLE by the SVM-RFE algorithm ( Figure 2D ). In addition, 394 genes were screened as potential biomarkers using the RF algorithm ( Figure 2E ). Based on the correlation with SLEDAI, 176 potential key biomarkers were identified from the DEGs. In addition, 3337 and 435 DEGs were obtained from the GSE11907 and GSE81622 datasets, respectively ( Figures 3A, B ). We finally obtained the unique key biomarker MX2 by the above algorithms ( Figure 3C ).

As shown in Figures 3D, E , we found high MX2 expression in SLE based on multiple PBMC and whole blood datasets. Subsequently, the same result was validated by RT‒qPCR ( Figure 3F ). Further exploration of the clinical value revealed that the expression of MX2 was positively correlated with SLE severity and SLEDAI (r = 0.38, P < 0.05, Figures 4A, B ). ROC curves showed the probability of MX2 as a valuable biomarker with AUCs of 0.958 and 0.976, respectively ( Figures 4C, D ), with high accuracy.

To determine the detailed mechanisms of MX2 in SLE, we performed a multifaceted analysis and validation. enrichKEGG analysis showed that MX2 was mainly enriched in the NOD-like receptor signaling pathway, and the same results were deduced from the GSE11907 and GSE81622 datasets ( Figures 5A, B and Supplementary Figure S3C ). The GSEA of the above three datasets also identified the NOD-like receptor signaling pathway as an important association pathway of MX2 ( Figures 5C, D and Supplementary Figure S3D ). To further verify the relationship between MX2 and the NOD-like receptor signaling pathway, we performed correlation analysis between MX2 and genes of this pathway. The results showed that MX2 was correlated with most genes of the pathway ( Figure 6A , Supplementary Figures S4A, B ). Among them, OAS1, OAS2, OAS3, STAT1, and IRF9 showed a high correlation with MX2 in all three datasets ( Figure S4C ). Next, the ssGSEA algorithm was used to score the NOD-like receptor signaling pathway for the SLE samples. Based on the score, the correlation between MX2 and the NOD-like receptor signaling pathway was calculated ( Figure 6B , Supplementary Figure S4D ), and the results showed that MX2 was clearly correlated with the NOD-like receptor signaling pathway. There were significant differences in pathway scores between the high and low MX2 expression subgroups ( Figure 6C , Supplementary Figure S4E ) and similarly significant differences in MX2 expression between the high and low pathway score subgroups ( Figure 6D , Supplementary Figure S4F ).

Four immune infiltration analysis methods were employed to explore the relationship between MX2 and immune cells. Based on the results of ssGSEA, 28 differences in immune cells between the normal and SLE groups were shown in Figures 7A, C , and differences in immune cells between the high and low MX2 expression subgroups were shown in Figures 7B, D . Correlation analysis between MX2 and immune cells showed that MX2 was associated with activated dendritic cells, central memory CD8 T cells, eosinophils, natural killer T cells, neutrophils, regulatory T cells and type 2 T helper cells with a high correlation ( Figure 7E ). Further validation analysis of these seven immune cells by CIBERSORT, xCell and ImmuCellAI algorithms revealed that the relationship between MX2 and activated dendritic cells and neutrophils also remained consistent in these three algorithms, and the trend of neutrophils was much more stable than that of activated dendritic cells ( Figure 7F ). The results of immune cells infiltration analysis of all CIBERSORT, xCell and ImmuCellAI were presented in Supplementary Figures S5 , S6A–C . Based on the CellMark database, the correlation of MX2 with neutrophil marker genes was analyzed ( Figure 7G ), and the results showed that MX2 correlated with all the neutrophil marker genes. In immune cell analysis using the HPA database, HPA and Monaco transcriptome data showed that MX2 was predominantly expressed in neutrophils ( Figures 8A, B ), and gene expression clustering data showed that MX2 was part of the cluster 50 (neutrophil-cellular senescence) gene cluster with a confidence level of 0.98129 ( Supplementary Figure S6D ). Then, in single-cell analysis, the flow cytometry data showed that MX2 in PBMCs was mainly expressed in neutrophils ( Figure 8C ), and the gene expression clustering data showed that MX2 was part of the cluster of cluster 19 (monocytes and neutrophils-innate immune response)-related genes with a confidence level of 1 ( Supplementary Figure S6E ).

To explore the relationship between MX2, the NOD-like receptor signaling pathway and neutrophils, based on four immune infiltration analysis algorithms and correlation analysis, we found that the NOD-like receptor signaling pathway was significantly associated with neutrophils ( Figure 9A ). The difference in neutrophils between the high and low pathway groups was shown in Figure 9B . The results showed that neutrophil scores were significantly higher in the high NOD-like receptor signaling pathway scoring group in the GSE121239 and GSE11907 datasets. The same trend was observed in the GSE81622 dataset. The results of the correlation of the NOD-like receptor signaling pathway with all immune cells of the four algorithms were shown in Supplementary Figures S7A–D .

To determine the regulatory mechanism of MX2, we used the STRING database for PPI analysis to predict the proteins closely related to MX2, including HERC5, IFIT2, IFIT3, IRF7, MX2, OAS1, OAS2, OAS3, OASL and RSAD2 ( Figure 10A ). Besides, Cistrome DB database was applied to explore these transcription factors (TFs) of the corresponding genes of related proteins ( Figure 10B , Supplementary Figure S8 ). The results showed that the TFs regulating MX2 mainly included CUX1, TCF7L1 and SPI1, and the “gene-transcription factor” regulatory network was established by Cytoscape, which contains 132 nodes and 200 edges ( Figure 10C ).

A total of 24 genes, including MX2 and its highly related NOD-like receptor signaling pathway genes (r >0.5), were used as pathogenic targets for drug prediction, and the distribution of the screened drugs was shown in Figure 11 , of which 652 drugs had |tau| scores >99, with detailed scores in Supplementary File S6 .