Figure 1 shows the overall study design. We downloaded series matrix files of GSE99967 from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo) database to identify the hub genes related to LN. The platform of the GSE99967 dataset was GPL21970 (Affymetrix Human Gene 2.0 ST Array). After obtaining the data, one sample GSM2666765 was removed for its inconsistent expression data (Supplementary Figure 1). The remaining whole peripheral blood samples from 29 LN patients and 16 normal controls were kept for further analysis (11). The expression data had been log₂ transformed.

The R package “WGCNA” (7) was utilized to construct the co-expression network of all genes in the GSE99967 dataset. Briefly, we calculated Pearson's correlation to construct a pairwise coefficient matrix and then transformed into a weighted adjacency matrix with an appropriate soft-thresholding power. Next, we converted the adjacency matrix into a topological overlap matrix (TOM), and used hierarchical clustering and dynamic tree cut algorithm to classify genes on TOM-based dissimilarity. Then we calculated the correlations of module eigengene (ME), which represents the first principal component of a module, and merged similar modules with highly correlated eigengenes. In addition, a module-trait relationship was calculated based on the correlation between the gene module and clinical trait. Finally, gene significance (GS) which refers to the correlation between gene expression and the clinical trait was used to identify clinical-related modules.

To investigate the potential biological function of genes in clinical-related modules, we performed GO enrichment analysis using Database for Annotation, Visualization and Integrated Discovery DAVID (https://david.ncifcrf.gov/). GO terms with p < 0.05 were considered statistically significant. The R package “ggplot2” was utilized to show results in the bubble chart.

The GEO2R online tools (12) were used to identify differentially expressed genes (DEGs) between LN patients and normal controls. Genes with adjusted p < 0.05 and | log₂ fold change (FC)| ≥ 1 were considered as DEGs. The heatmap and volcano plot of all DEGs were made by the R package “pheatmap” and “ggplot2,” respectively.

We used CytoHubba (13), a Cytoscape software (14) plugin to select the hub genes in the key module. The top ten genes with higher maximal clique centrality (MCC) values were screened. Among these ten genes, five DEGs were identified as hub genes. We investigated the GS and module membership (MM) of the selected hub genes to verify their reasonability.

To validate five hub genes, we downloaded another two datasets GSE72798 (15) and GSE32591 (16) from the GEO database. The dataset GSE72798 consists of 10 LN patients and 10 normal controls and the RNA was extracted from blood samples. In the GSE32591 dataset, the total RNA was extracted from kidney biopsy of 32 LN patients and 14 normal controls.

To further verify selected genes, we recruited 78 LN patients, 67 healthy controls (HCs), and 25 patients with IgA nephropathy from the First Affiliated Hospital, College of Medicine, Zhejiang University. Patients or the public provided written informed consent to participate in this study. All LN patients were confirmed by kidney biopsy, and the Systemic Lupus Erythematosus Disease Activity Index (SLEDAI) was calculated, which contains 24 items reflecting disease activity. Active LN was defined as biopsy-proven active nephritis or SLEDAI ≥10 and at least two renal elements of the SLEDAI; inactive LN was defined as biopsy-proven pure class V or VI nephritis or SLEDAI <10.

Serum was separated from blood samples and stored at −80°C. Peripheral blood mononuclear cells (PBMCs) were isolated by human PBMCs separation medium (tbdscience, Tianjin, China).

Total RNA was extracted from PBMCs samples using TRIzol® (Invitrogen, CA, USA) and then reversely transcribed into cDNA by the PrimeScript™ II Reverse Transcriptase (Takara, Shiga, Japan). According to the manufacturer's instructions, reverse transcription was conducted at 37°C for 15 min and 85°C for 5 s. qRT-PCR was performed by SYBR Green on CFX96™ Real-Time PCR Detection Systems (Bio-rad, CA, USA). All primer sequences used were shown in Supplementary Table 1. The conditions of qRT-PCR were as follows: 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 31 s, with dissociation at 95°C for 15 s, 60°C for 1 min and 95°C for 15 s. The average threshold cycle (Ct) was calculated for each transcript, which was performed in triplicate. The relative expression of each mRNA was normalized by GAPDH and analyzed using the 2^−ΔΔCT method [ΔΔCT= (CT of gene) – (CT of GAPDH) – (CT of HC)].

The serum levels of IFI44 and IFIT3 were assessed with ELISA kits (G-Biosciences, CA, USA for IFI44, and MyBioSource, CA, USA for IFIT3). Briefly, diluted serum samples were incubated in the microplate. The membrane was then incubated with biotin-labeled antibody after washing twice. Similarly, the membrane was incubated with avidin-labeled horseradish peroxidase (HRP) and then with substrated solution (TMB). Finally, the absorbance values at 450 nm of samples were read, and the concentration of IFI44 and IFIT3 were calculated using standard curves run on each ELISA plate.

Statistical analyses were performed using SPSS 20 (IBM) and GraphPad Prism 8.0 (GraphPad Software, CA, USA). The unpaired t-test was used to analyze differences between two groups. Results were expressed as the mean ± standard deviation (SD). The receiver operator characteristic (ROC) curves were drawn and area under curve (AUC) was calculated to assess the diagnostic capability of biomarkers. Logistic regression analysis was performed to create a combined score. P < 0.05 was considered statistically significant.

Samples in the GSE99967 dataset were divided into two groups (29 LN samples and 16 normal controls). Expression data of all genes from GSE99967 were used to conduct WGCNA (Figure 2). Based on scale-free topology model fit index and mean connectivity, the soft-threshold power was set as four (Figure 2A). Then, we merged similar modules to acquire 12 co-expression modules by setting the cut height of module eigengenes as 0.5 (Figures 2B,C). Figure 2D showed the cluster dendrogram and adjacency heatmap of eigengenes, which means that 12 modules were primarily separated into two clusters.

From the heatmap of module-trait relationships, we found that multiple modules were associated with LN (Figure 3A). Then, we calculated the gene significance of all genes in 12 co-expression modules (Figure 3B). The results showed that the darkgreen module was most significantly related to LN, followed by the purple module. The correlation between MM and GS in these two modules are shown in Figure 3C, which mean that MM was highly correlated with GS for LN in both darkgreen module and purple module (cor = 0.77, p = 5.8e−76, and cor = 0.41, p = 5.2e−15, respectively).

To investigate the potential biological function of genes in the above two modules, we performed GO enrichment analysis. The top 20 significant GO biological process terms of each module were shown in Figures 3D,E. We found that genes in the purple module were mainly enriched in defense response to virus, type I interferon signaling pathway, immune response, and inflammatory response, which are highly associated with LN. However, genes in the darkgreen module were mostly distributed in terms like cell diversion and proliferation, mitotic nuclear division, and DNA replication, which mainly play roles in metabolism. Thus, we defined the purple module as the key module related to LN and performed further analysis.

Hub genes should be DEGs between LN patients and normal controls. In dataset GSE99967, a total of 137 DEGs were identified, including 84 up-regulated genes and 53 down-regulated genes (Supplementary Table 2). The heatmap and volcano plot of all DEGs are shown in Supplementary Figures 2, 3.

To identify hub genes in the key module, we calculated MCC values of genes in the purple module and extracted the top 10 genes with higher scores (Supplementary Figure 4). Among these ten genes, 5 DEGs were identified as hub genes, including interferon induced protein 44 (IFI44), interferon induced protein with tetratricopeptide repeats three (IFIT3), HECT and RLD domain containing E3 ubiquitin protein ligase 5 (HERC5), radical S-adenosyl methionine domain containing two (RSAD2), and DExD/H-box helicase 60 (DDX60). The GS and MM value of these genes were all > 0.4 and 0.9, respectively, indicating that hub genes were significantly associated with the key module and LN trait.

As expected, the expression levels of hub genes including IFI44, IFIT3, HERC5, RSAD2, and DDX60 were significantly upregulated in LN samples from the GSE99967 dataset (Figure 4A). For verifying hub genes, we obtained another two datasets GSE72798 and GSE32591, and analyzed the expression levels of the above five genes between LN patients and normal controls (Figures 4B,C). Except for DDX60 in GSE72798, the other four hub genes were markedly increased in LN samples.

In addition, to further validate these four hub genes, we collected nine blood samples from four healthy controls and five LN patients (patient characteristics are shown in Supplementary Table 3) to perform qRT-PCR. The results showed that compared with HCs, the mRNA levels of IFI44, IFIT3, HERC5, and RSAD2 were all significantly elevated in LN patients (Supplementary Figure 5), meaning that they are expected to be potential biomarkers in identifying LN.

Among five hub genes, we selected type I interferon-response genes IFI44 and IFIT3 for further analysis. To detect the serum levels of IFI44 and IFIT3, 124 subjects, including 51 healthy controls and 73 LN patients (patient characteristics are shown in Supplementary Table 4) were used to perform ELISA assays (Figure 5). As expected, IFI44 and IFIT3 were markedly upregulated in LN patients compared with the healthy controls (Figure 5A). Furthermore, ROC curves were used to calculate the sensitivity and specificity of these two biomarkers for identifying LN patients. The AUC values of IFI44 and IFIT3 were 0.811 and 0.758, respectively (Figure 5B, Supplementary Table 5). Then, we combined these two biomarkers using logistic regression analysis (Supplementary Figure 6A, Supplementary Table 6). However, the AUC value of combined score (AUC = 0.811) was the same as IFI44 (Supplementary Figure 6B), indicating that IFI44 had the highest power to distinguish between the two groups. Additionally, we found that IFI44 was significantly elevated in active LN compared with inactive LN patients, while IFIT3 failed to discriminate active LN from inactive ones (Figure 5C). And the AUC of IFI44 in differentiating active LN from inactive LN patients was 0.697, which was similar to clinical indicators such as anti-dsDNA (AUC = 0.600), C3 (AUC = 0.665), and C4 (AUC = 0.673) (Figure 5D; Supplementary Table 7).

To examine whether serum IFI44 and IFIT3 can reflect the pathological class of kidney biopsy, we classified 73 LN patients into five groups according to LN class. The results showed that they were both elevated in class III(±V), class IV(±V), and class V compared with HCs (Supplementary Figure 7), meaning that both IFI44 and IFIT3 can indicate the prognosis of LN.

Moreover, we identified the specificity of IFI44 for LN. An independent cohort of 37 subjects, including 12 healthy controls and 25 patients with IgA nephropathy (patient characteristics are shown in Supplementary Table 8), were used for investigation. Interestingly, there was no statistical difference between HCs and patients with IgA nephropathy (Figure 5E), indicating that IFI44 is specific for LN.