Public microarray data containing clinical information on CKD and normal kidney tissues was obtained from the NCBI GEO GSE104948 and GSE116626 datasets. The GSE104948 data set (RNA was extracted from the glomerular compartment), including 50 CKD kidney samples and 18 normal samples, was based on the Affymetrix Human Genome U133 Plus 2.0 Array of the GPL22945 platform. The GSE116626 data set (RNA was extracted from archival formalin-fixed paraffin-embedded kidney biopsy samples), including 74 CKD kidney samples and 7 normal samples, was based on the Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip of the GPL14951 platform. To merge the multiple datasets, we utilized the inSilicoMerging (Taminau et al., 2012) R package to process the datasets. In addition, we used the Johnson et al. method (Johnson et al., 2007) to remove group effects. In total, 124 CKD samples and 25 normal samples of tissues were included in the follow-up analysis of this study.

Here, differential analysis was conducted utilizing the limma R package (Ritchie et al., 2015) to get genes that differ between the CKD group and the control group. The statistical criterion for screening RNA expression was | fold-change (FC) | > 1.5 and p-value <.05. Based on the org. Hs.eg.db R package, the KEGG rest API, and the Molecular Signatures Database, we obtained gene annotations for GO, KEGG, and C5, respectively. Then, we performed function analysis utilizing the ClusterProfiler R package to get the DEGs enrichment results. p-value <.05 was statistically significant. The maximum gene set is 5000 and the minimum gene set is 5.

Using gene expression profiles, we calculated the mean absolute deviation (MAD) for each gene and excluded the 50% of DEGs with the lowest mean absolute deviation. In addition, we used the R package WGCNA to remove outlier DEGs and probes to construct a scale-free co-expression network. Specifically, first the Pearson’s correlation matrices and average linkage method were both performed for all pair-wise Genes. Then, a weighted adjacency matrix was constructed using a power function A_mn = |C_mn|^β (C_mn = Pearson’s correlation between Gene_m and Gene_n; A_mn = adjacency between Gene m and Gene n). β was a soft-thresholding parameter that could emphasize strong correlations between Genes and penalize weak correlations. After choosing the power of 6, the adjacency was transformed into a topological overlap matrix (TOM), which could measure the network connectivity of a Gene defined as the sum of its adjacency with all other Genes for network Gene ration, and the corresponding dissimilarity (1-TOM) was calculated. To classify Genes with similar expression profiles into Gene modules, average linkage hierarchical clustering was conducted according to the TOM-based dissimilarity measure with a minimum size (Gene group) of 10 for the Genes dendrogram. To further analyze the module, we calculated the dissimilarity of module eigen Genes, chose a cut line for module dendrogram, and merged some modules. A total of eight co-expression modules were obtained by merging the modules with a distance less than 0.25. Lastly, GS and MM were calculated according to correlations between gene expressions with clinical subtype and correlations between gene expression and module feature vector, respectively. In the clinically significant module, 16 highly connective genes were screened as key genes according to the cut-off criteria [(MM) > 0.8 and (GS) > 0.1].

The PPI networks for modules with very robust filtering conditions (score >0.7) were analyzed using the STRING database. Cytoscape Software (version 3.8.2) was utilized to visualized these PPI networks. The main functional modules were analyzed using Cytoscape’s Molecular Complex Detection Technology (MCODE) plug-in. Selection criteria are defined as follows: K Core = 2, Cut Grade = 2, Maximum Depth = 100, Cut Node Score = 0.3. Cytoscape’s plugin cytoHubba uses the MCC (Maximum Clique Centrality) algorithm to score each node gene. The pivot genes were screened using the top 5 nodal genes of each algorithm’s MCC score. Predictions of gene function and mapping genes with comparable effects were generated by GeneMANIA, a website for constructing PPI networks. Some of the bioinformatics methods employed by network integration algorithms are physical interactions, co-expression, co-localization, gene enrichment analysis, genetic interactions, and locus prediction. In this study, we used GeneMANIA to identify PPI networks of eigengenes.

To verify the predicted value of the screende hub genes, we constructed a logistics model using PROC in the R package (version 3.6.3) and used the GGPLOT2 package to visualize the results. The diagnostic value of the identified biomarkers was assessed by the area under the ROC curve (AUC, AUC was between 0.5 and 1). The closer the AUC is to 1, the more effective the diagnosis is. In addition, we performed a controlled reliability analysis using the RNA expression datasets GSE93798 and GSE104066 as validation sets. GSE93798 includes 20 CKD samples and 22 control samples (RNA was extracted from the glomerular compartment). GSE104066 includes 70 CKD samples and 6 control samples (RNA was extracted from the glomerular compartment).

Paraffin-embedded kidney tissue sections from Chaohu Hospital of Anhui Medical University, including CKD patient group (n = 3) and normal control group (n = 3), were obtained according to Institutional Review Board-approved protocols, and informed consent forms were signed by the patients. The expression of IL10RA (1:50, ZENBIO, China), CD45 (1:50, ZENBIO, China), CTSS (1:50, ZENBIO, China), and C1QA (1:50, Affinity, USA) were detected according to the instructions of the immunohistochemistry kit (ZSBIO, China).

CIBERSORT (Chen et al., 2018) transforms normalized gene expression matrices into immune cell invasiveness components to estimate the relative frequency of immune invasion and is a 1,000 permutation deconvolution algorithm. Then build a histogram to display the 22 types of content. A correlation heatmap of immune cell infiltration in each sample was developed to visualize correlations between immune cell subtypes. Additionally, differential analysis between CKD and normal tissue immune cells was also visualized by a violin plot. Importantly, associations between the identified biomarkers and the level of infiltrating immune cells were explored and visualized by dot-bar graphs using Spearman’s rank correlation analysis.

We utilized the GSEA analysis (Subramanian et al., 2005) to explore regulatory target genes, biological process (BP) GO terms, KEGG pathways, and Human phenotype Ontology in which the selected IL10RA might be involved in CKD. The samples were divided into low expression group (<50%) and high expression group (≥50%) according to the expression level of IL10RA. Datasets in the Molecular Signatures Database, including c3.all.v7.4. symbols.gmt, c5. go.bp.v7.4.symbols.gmt, c2.cp.kegg.v7.4.symbols.gmt, and c5.hpo.v7.4.symbols.gmt, were served as reference gene sets. Statistical significance was assessed by comparing the enrichment score with the enrichment results generated from 1000 random permutations of the gene set to obtain p-value, and p < .05 was considered significant for GSEA analysis using default parameters.

To explore drug-gene interactions, existing or/and potentially relevant drug substances were identified using the STITCH database (Szklarczyk et al., 2016). The PubChem database (Kim et al., 2021) and the PDB database (Karuppasamy et al., 2020) were used to obtain the molecular structures of ligands and target proteins. Docking simulations and visualization were performed through PyMOL software (Nguyen et al., 2020) and AutoDock Vina (Lam and Siu, 2017).

All data were processed and analyzed using R software. The Mann-Whitney U test (Wilcoxon rank sum test) was used to analyze differences between two groups of continuous non-normal variables. A possible correlation between two variables was detected by the Pearson correlation coefficient. p < .05 considered the difference to be statistically significant.

After merging the GSE104948 and GSE116626 datasets, we removed batch-to-batch variance from the matrix of gene expression (Supplementary File S1, S2). In Figure 2A, the box diagram shows that the sample distribution of each data set is quite different before the batch effect is removed, revealing that the batch effect exists. The sample distributions of the two datasets tend to be consistent after excluding the batch-to-batch variance, and the medians are on the same straight line. Figure 2B depicts UMAP results for multiple datasets with different colors representing different datasets before showing batch deletion. As shown, the two datasets do not intersect with each other and are independent of each other. After removing the batch-to-batch variance, the sample distributions between datasets tend to be consistent. From the density map in Figure 2C, we can observe that there is a great difference in the sample distribution of each data set before excluding the batch effect. The sample distributions between the datasets tend to be consistent after eliminating the batch effect.

After preprocessing the data with R software, we extracted the DEGs in the gene expression matrix. Under the criteria of p-value <.05 and | fold-change (FC) | >1.5, 657 genes were identified as DEGs, with 521 genes up-regulated and 136 genes down-regulated (Supplementary File S3). Figures 3A, B show a volcano plot of DEGs and a heatmap of the top 50 DEGs. Next, Human phenotype ontology, GO and KEGG signaling pathway enrichment analyses were performed to dissect the biological functions and signaling pathways involved in 657 selected DEGs (Supplementary Files S4–S6).

The top 10 results of Human phenotype Ontology show that nephritis, membranoproliferative glomerulonephritis, and impaired oxidative burst were significantly enriched (Figure 3C), which indicates the reliability of our data. More importantly, the top 10 GO analysis shows that a large number of biological processes related to immune and inflammatory responses are significantly enriched, including cell activation, immune response, immune system process, leukocyte activation, and myeloid leukocyte activation (Figure 3D). In terms of KEGG Pathway, complement and coagulation cascades, ECM-receptor interaction, and Fc gamma R-mediated phagocytosis are significantly enriched (Figure 3E). The results above strongly suggest that autoimmunity and inflammation play essential roles in the development process of CKD.

Based on the screened 657 DEGs expression profile, WGCNA was performed to identify the major modules most associated with CKD (Supplementary Files S7–S10). Eight modules were identified after merging strong association modules with a cluster height limit of 0.25 (Figure 4A). The module feature vector clustering was investigated next, and the results revealed the distance between them (Figure 4B). Furthermore, the correlations between modules and clinical symptoms were explored. The red module (r = 0.45, p = 1.0e-8), the turquoise module (r = 0.52, p = 1.6e-11), the black module (r = 0.50, p = 6.9e-11), and the blue module (r = 0.53, p = 5.8e-12) are positively correlated with CKD, while the pink module (r = −0.61, p = 1.2e-16), the brown module (r = −0.41, p = 1.6e-7), the magenta module (r = −0.36, p = 5.5e-6), and the grey module (r = −0.58, p = 5.7e-15) are negatively correlated with CKD (Figure 4C).

We performed functional enrichment to explore more about the biological functions of the DEGs in eight modules (Supplementary Files S11, S12). The results of GO and KEGG analysis revealed that DEGs in the turquoise module were linked to a large number of biological processes and pathways related to autoimmunity, inflammation, and infection. GO enrichment analysis showed that turquoise module DEG genes have leukocyte activation involved in immune response, cell activation involved in immune response, leukocyte degranulation, and neutrophil activation (Figure 4D). KEGG analysis was associated with Chemokine signaling pathway, Natural killer cell mediated cytotoxicity, Complement and coagulation cascades, and Viral protein interaction with cytokine and cytokine receptor (Figure 4E). According to GS > 0.8 and MM > 0.1, 16 genes in the turquoise module are identified as key genes (MS4A6A, RAC2, GPR65, LYZ, MYO1F, PYCARD, LCP1, CTSS, AOAH, IL10RA, CD53, EVI2A, C1QA, NCF2, PTPRC, MS4A4A).

To further discover CKD-related hub genes and their mechanisms, we mapped the above 16 key genes with high expression in the turquoise module of the CKD group, uploaded them to the online STRING database, and constructed a PPI network (Supplementary File S13). A PPI network with 15 nodes and 43 edges was realized (Figure 5A). Among the 15 nodes, the top 4 genes with a high binding degree were found by Cytoscape (version 3.8.2) MCODE and MCC calculation methods. These genes, which were identified to play hub roles in CKD, are listed as follows: IL10RA, CD45, CTSS, and C1QA (Figure 5B). The specific information of the hub genes is shown in Table 1.

Next, we explored the co-expression networks and potential functions of hub genes according to the GeneMANIA database (Figure 5C) (Supplementary File S14). They revealed the sophisticated PPI networks with the protein domains of 0.60%, pathway of 1.88%, genetic interactions of 2.87%, co-localization of 3.64%, predicted of 5.37%, co-expression of 8.01%, physical interactions of 77.64%. Function analysis indicated that they are mainly related to a variety of immune and inflammatory pathways, including humoral immune response mediated by circulation, B cell mediated immunity, complement activation, adaptive immune response, humoral immune response, interleukin-8 production, and receptor signaling pathway via JAK-STAT, revealing their essential role in contributing to CKD pathogenesis.

We conducted ROC analysis to study the relationship between hub gene expression and the prognosis of CKD patients (Supplementary File S15). An AUC greater than 0.800 is considered to have excellent specificity and sensitivity for the diagnosis of CKD. As shown in Figure 6A, the AUC value of IL10RA was 0.821 (95% CI: 0.730-0.911), CD45 was 0.836 (95% CI: 0.739-0.933), CTSS was 0.861 (95% CI: 0.773-0.949), and C1QA was 0.836 (95% CI: 0.749-0.923). More importantly, the combination of all 4 hub genes is 0.881 (95% CI: 0.795-0.966). The results showed that IL10RA, CD45, CTSS, and C1QA have high diagnostic value.

Furthermore, two new CKD-related datasets, including GSE93798 (Figure 6B) and GSE104066 (Figure 6C), validated the above four hub DEGs (Supplementary File S16). Through verification, the mRNA expression of each hub gene was significantly overexpressed in CKD, compared with the control. The validation results fully support the assumption that IL10RA, CD45, CTSS, and C1QA may be diagnostic markers of CKD.

To verify the expression of IL10RA, CD45, CTSS, and C1QA in CKD, we treated kidney tissues with IHC and found that IL10RA, CD45, CTSS, and C1QA were all highly expressed in the CKD tissues compared with the control subjects, which is consistent with our bioinformatics prediction (Figures 7A–H).

To examine differences in immune patterns between CKD and normal tissues, the matrix of gene expression estimated the infiltration ratio of 22 immune cells using the CIBERSORT method (Supplementary Files S17–S19). In each sample, a histogram depicted the composition of 22 types of immune cells (Figure 8A). Colors on every histogram exhibit the percentages of immune cells, with a sum of 1 for each sample. The results indicated that naive B cells (137), neutrophils (136), M1 macrophages (133), M2 macrophages (131), and resting CD4 memory T cells (129) were the most abundant immuno-infiltrating cells in all 149 samples. In the following study, 22 kinds of immune cells in CKD samples were evaluated for their correlation (Figure 8B). The correlation heat map of 22 immune cells showed that naive B cells (GBD 2016 Causes of Death Collaborators, 2017), regulatory T cells (Dekker et al., 2017), monocytes (Zimmermann et al., 1999), M2 macrophages (Zoccali et al., 2006), and activated NK cells (Honda et al., 2006) are associated with most immune cells. However, activated CD4 memory T cells (Holle et al., 2022), M0 macrophages (Holle et al., 2022), CD8 T cells (Quon et al., 2011), resting CD4 memory T cells (Quon et al., 2011), and eosinophils (Quon et al., 2011) are only associated with a few immune cells. Violin plots of the difference in immune cell infiltration indicated that, compared with the normal control sample, gamma delta T cells, activated NK cells, monocytes, M0 macrophages, and M1 macrophages infiltrated more, while naive B cells, regulatory T cells and activated dendritic cells infiltrated less (Figure 8C).

We further explored whether there is a potential correlation between immune cell abundance and hub gene expression using Pearson’s correlation analysis (Supplementary File S20). As shown in Figures 9A–D, a total of seven immune cell populations that are related to all four core genes, of which naive B cells, resting memory CD4 T cells, regulatory T cells, and activated dendritic cells were statistically negatively with IL10RA, CD45, CTSS, and C1QA, while gamma delta T cells, monocytes, M0 macrophages, and M1 macrophages were positively correlated with them, suggesting they may play essential roles in CKD development.

Since IL10RA plays an important role in immune infiltration, and the log2FC of IL10RA is the largest of the central genes, we performed an IL10RA analysis using the GSEA method to gain insight into the biological processes and predict the potential signal pathways of IL10RA expression in CKD (Supplementary Files S21–S24). The top 10 results of MSigDB C5 Human phenotype Ontology showed IL10RA was involved in abnormalities in various immune cells, including abnormal leukocyte, abnormal granulocyte, abnormal neutrophil, abnormal myeloid leukocyte morphology, and abnormal lymphocyte physiology (Figure 10A). In addition, high levels of IL10RA may affect several manipulated downstream potential genes, including MAML1, PEA3, BACH2, ELK1, ZNF597, ETS2, MIR92A, and TEL2 (Figure 10B).

More importantly, the top 10 results of MSigDB C5 GO biological processes showed IL10RA regulated a large number of immune system responses, including antigenic processing and presentation of polypeptide antigen, immune response regulating signal pathway, regulation of response to biotic stimulus, T cell receptor signaling pathway, myeloid leukocyte mediated immunity, I kappaB kinase NF kappaB signaling, and response to interferon gamma (Figure 10C). Meanwhile, MSigDB C2 KEGG gene sets found that in addition to B cell receptor signaling pathway, T cell receptor signaling pathway, and Fc gamma R-mediated phagocytosis, IL10RA was also related to a large number of inflammation-related pathways, including chemokine signaling pathway, apoptosis, Nod like receptor signaling pathway, Toll like receptor signaling pathway, and cell adhesion molecules cams. Thus, the above results suggest that the immune and inflammatory responses play essential roles in IL10RA contributing to the CKD pathogenesis (Figure 10D).

Searching for targeted drugs for IL10RA provides a new strategy for potential drug therapy for CGN. Based on the STITCH database, we obtained 4 small molecular drugs, including chitin, selenomethioni, leupeptin, and isosorbide din (Supplementary File S25). Then, the above 4 bioactive compound ligands were docked with the protein IL10RA to evaluate the binding potential. As shown in Figures 11A–D, the docking 3D model of protein IL10RA and four small molecules drugs with the firmest binding, showing their potential to alleviate or even reverse CKD development (Supplementary File S26).