Identifying C1QB, ITGAM, and ITGB2 as potential diagnostic candidate genes for diabetic nephropathy using bioinformatics analysis

Introduction

Diabetic nephropathy (DN) is a serious complication caused by diabetic microangiopathy, which often progresses to end-stage renal disease (Singh, Winocour & Farrington, 2011). The pathological changes in DN included glomerular vascular injury, glomerulosclerosis, nodular lesion formation, renal function deterioration, and, ultimately, end-stage renal disease (Qi et al., 2017). Globally, 30 to 40% of diabetes cases have been reported to progress to DN. In spite of substantial advances in DN research in recent years, a large proportion of patients still irreversibly suffer from end-stage renal disease. Compared to other complications, metabolic disorders caused by DN make it more difficult to treat end-stage renal disease once they develop. Therefore, early diagnosis and timely intervention in DN are becoming more and more important (Nakagawa et al., 2011; Ruiz-Ortega et al., 2020).

The vast amounts of data generated by new technologies such as genome sequencing and microarray chips have expanded the ways in which data can be analyzed and interpreted for biological understanding and therapeutic advances (Bednar, 2000). However, most resent studies on DN have focused on exploring the expression of differential genes that affect glomerular pathological changes. These cannot fully reflect the pathogenesis of DN, because in the context of diabetes, renal tubular function plays an important role in regulating glomerular filtration (Vallon & Thomson, 2012). However, it is not appropriate to focus on the single transcription level of mRNA (Gao et al., 2021; Chen et al., 2022), and attention should be paid to other components in the transcriptional regulatory network.

MicroRNA (miRNA) are a class of endogenous small RNA with a length of about 20–24 nucleotides that play a variety of important regulatory roles in cells (Wang, Chen & Sen, 2016). MiRNA regulate gene expression mostly by degrading mRNA (Etheridge et al., 2011). In recent years, new progress has been made in the study of miRNA and human diabetes (Vasu et al., 2019).

Transcription factors (TFs) can activate or suppress gene transcription by binding to specific parts of chromatin. Previous studies mostly focused on the relationship between lncRNA-miRNA-mRNA in ceRNA networks (Guo et al., 2022), but TFs are also involved in complex transcriptional regulation. TF expression is closely related to the physiological and pathological state of cells, and their expression levels in tissues are specific in time and space. Therefore, changes in TF expression reflect changes in cell state and may lead to pathological processes (Lambert et al., 2018) such as diabetes (Polvani et al., 2016). However, in-depth research has revealed that TFs function through complex regulatory patterns, even though TFs in the same family may have a bidirectional effect on the progression of the same disease (Song et al., 2016; Rad et al., 2016). Therefore, it is crucial to accurately explore the role of both miRNA, and TFs in DN. Further understanding of the direct regulatory relationship between mRNA, miRNA and TF can provide new insights for personalized disease management, diagnosis and prognosis.

In this study, we selected the GSE30529 dataset containing DN renal tubule samples from the GEO database, used the robust multiple array averaging (RMA) method in R language to preprocess and normalize the data, and identified differentially expressed genes (DEGs) according to screening criteria. After that, gene set enrichment analysis (GSEA), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were applied to identify the biological processes, cellular components and molecular functions in which the DEGs were involved. Subsequently, a relevant protein-protein interaction (PPI) network was constructed, and the MCC algorithm of Cytohubba plug-in in Cytoscape software was used to identify the top 20 genes with the highest interaction scores as hub genes. Next, based on the GSE30122 dataset containing 69 samples, the ROC method was used to verify those genes. Subsequently, five Hub-genes with the most potential diagnostic values were selected. MiRNA and TF targets were predicted through four different online prediction websites, and the ‘nephroseq’ database was used to predict the correlation between genes and kidney function. The miRNA-mRNA-TF interaction was completed to elucidate the mechanism by which these genes and expression elements interact and cooperate to drive the occurrence and development of DN. Finally, the expression of hub genes in DN rats was further verified using qPCR. Altogether, these results have the potential to elucidate novel candidate genes associated with DN and provide a new insight into its molecular basis.

Materials & Methods

Microarray data acquisition

The mRNA expression profile datasets were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). We uesd the keywords “Diabetic Nephropathy,” “Homo sapiens,” and “High throughput gene expression profile” as search criteria, and the datasets GSE30529 and GSE30122 fit the above screening conditions. Both GSE30529 and GSE30122 were located on the GPL571 platform ((HG-U133A_2) Affymetrix Human Genome U133A 2.0 Array). GSE30529 was used to screen for differential genes associated with DN, while GSE30122 was used to verify and explore the predictive value of genes. GSE30529 included 10 diabetic tubule samples and 12 control samples to investigate the influence of diabetes mellitus on DEGs. GSE30122 contained nine glomeruli and 10 tubule samples from kidneys with DN and 50 control glomeruli samples from the unaffected portion of tumor nephrectomies. The Series Matrix files for two datasets including gene expression counts were downloaded for downstream analysis. Gene ID conversion based on the genome “hg38” was performed on the count matrix, the differential expression values of the same gene were averaged, and the original data were log2 transformed and normalized. A workflow of our study is displayed in Fig. 1.

Results

Preprocessing analysis and identification of DEGs

DN and control samples from the GSE30529- GPL571 datasets consisted of 10 DN kidney tubules and 12 control samples. After GSE30529 was normalized, the box plot showed that samples had good distribution (Fig. S1A). We performed PCA to evaluate the repeatability of the expression matrix within GSE30529, and the results showed satisfactory repeatability based on the high proportion of variance explained by the first two principal components (PC1 and PC2) and the tight clustering of replicates in the PCA plot (Fig.S1B). Initially, we screened 19,820 genes in the microarray dataset GSE30529 after removing the duplicate probe. After setting criteria with the threshold of an adjusted absolute value of log2 (FC)> 1 and adj. P value < 0.05, 463 DEGs (340 increased and 123 reduced) were identified as DEGs. The DEGs are shown in Fig. 2A in volcano plots, and the top40 DEGs with the most distinct expression differences between DN and control samples are shown on the heat map (Fig. 2B).

Enrichment analysis

Functional and pathway enrichment analyzes were performed by the R software ‘clusterProfiler’ package. First, the expression profiles of all genes in the DN and control samples were cleared and uploaded to GSEA, and the C2 sub-collection CP gene collection was applied to conduct the GO enrichment analysis. Adj. p < 0.05 and Q < 0.25 were set as the filtering criterion for significant gene sets. The most enriched gene sets were related to mitosis (Fig. 3A; normalized enrichment score (NES) = 1.857, false discovery rate (FDR) = 0.017), immune responses (Figs. 3B–3C; NES = 1.928, FDR = 0.017; NES = 2.317, FDR = 0.017), angiogenesis (Fig. 3D; NES = 1.582, FDR = 0.017), and cytoskeletal dynamics and morphology pathways (NES = 1.629, FDR = 0.017), as determined by the GSEA results. At the same time, the gene sets related to immune responses were found to be among the top-ranked pathways in terms of enrichment score, indicating their potential importance in the context of diabetic nephropathy (Fig. 3F). Next, the results of GO, KEGG pathway, GSEA, and Reactome enrichment analyses suggested that the immune response, including the regulation of neutrophil activation, response to interferon-gamma, neutrophil mediated immunity, and neutrophil degranulation, in DN samples was more prominent than in the control samples (Fig. 4A). KEGG enrichment analysis revealed a strong association between the identified DEGs and pathways related to neutrophil functions, such as phagosome (KEGG ID: hsa04145, adjusted p-value < 0.01), complement and coagulation cascades (KEGG ID: hsa04610, adjusted p-value < 0.01), and cell adhesion molecules (KEGG ID: hsa04514, adjusted p-value < 0.01) (Fig. 4B). More details of the enriched GO/KEGG pathway are included in Table 1. In addition, GO and KEGG enrichment analyses of revealed that immune-related biological processes and pathways were predominantly enriched. When looking at an immune-related nephropathic disease like DN, it is crucial to understand its immune mechanisms.

Table 1:

GO and KEGG pathways.

Ontology	ID	Description	p.adjust
BP	GO:0042119	neutrophil activation	2.48E−15
BP	GO:0034341	response to interferon-gamma	1.36E−14
BP	GO:0002446	neutrophil mediated immunity	2.42E−14
BP	GO:0043312	neutrophil degranulation	2.90E−14
BP	GO:0002283	neutrophil activation involved in immune response	2.97E−14
CC	GO:0062023	collagen-containing extracellular matrix	2.30E−24
CC	GO:0034774	secretory granule lumen	1.41E−14
CC	GO:0060205	cytoplasmic vesicle lumen	4.40E−14
CC	GO:0031983	vesicle lumen	4.40E−14
CC	GO:0031091	platelet alpha granule	7.82E−14
MF	GO:0005539	glycosaminoglycan binding	9.45E−13
MF	GO:0008201	heparin binding	2.86E−12
MF	GO:1901681	sulfur compound binding	8.54E−10
MF	GO:0004866	endopeptidase inhibitor activity	8.95E−09
MF	GO:0030414	peptidase inhibitor activity	1.33E−08
KEGG	hsa04145	Phagosome	1.34E−12
KEGG	hsa05150	Staphylococcus aureus infection	2.83E−10
KEGG	hsa04610	Complement and coagulation cascades	1.33E−08
KEGG	hsa05140	Leishmaniasis	1.63E−08
KEGG	hsa04514	Cell adhesion molecules	1.08E−07
KEGG	hsa05152	Tuberculosis	1.34E−07
KEGG	hsa05323	Rheumatoid arthritis	1.40E−06
KEGG	hsa05144	Malaria	1.40E−06

DOI: 10.7717/peerj.15437/table-1

Notes:

Annotation

BP

Biological Process

CC

Cellular Component

MF

Molecular Function

PPI network analysis, MCODE cluster modules and hub gene identification

We used STRING to construct the interaction network between proteins coded by DEGs, which consisted of 65 nodes and 183 edges, and the network was visualized by Cytoscape (Fig. 5A). We used the cytoHubba plugin to identify hub genes. The MCC algorithms identified 20 hub genes. Three gene cluster modules were identified with the MCODE plugin (Figs. 5B–5D), according to the filter criteria. Details of these twenty hub genes identified by cytoHubba are given in Table 2. These genes, including PSMB8, IRF8, IRF9, MX1, BST2, GBP2, IFITM2, ISG20, IFITM3, HLA-B, HLA-G, HLA-DRA, HLA-DPA1, ITGB2, ITGAM, HLA-DPB1, PTPRC, CCL2, HLA-DMA, and C1QB, were the most crucial in the PPI network and played an important role in the pathogenesis of DN.

Table 2:

Twenty hub genes identifed by cytoHubba.

Gene symbol	Description	log2FC	Q value	Regulation
PSMB8	Proteasome subunit beta type-8	1.342044794	0.000164732	UP
IRF8	Interferon regulatory factor 8	1.543495504	0.00054193	UP
IRF9	Interferon regulatory factor 9	1.015989833	0.000995947	UP
MX1	Interferon-induced GTP-binding protein Mx1	1.194190613	0.004106856	UP
BST2	Bone marrow stromal antigen 2	1.16619317	0.037461717	UP
GBP2	Guanylate-binding protein 2	1.447589391	0.003980602	UP
IFITM2	Interferon-induced transmembrane protein 2	1.216171959	0.020762203	UP
ISG20	Interferon-stimulated gene 20 kDa protein	1.018021414	0.00758014	UP
IFITM3	Interferon-induced transmembrane protein 3	1.241954263	0.020591992	UP
HLA-B	HLA class I histocompatibility antigen, B-7 alpha chain	1.499097739	0.000194636	UP
HLA-G	HLA class I histocompatibility antigen, alpha chain G	1.118671026	0.000373508	UP
HLA-DRA	Major histocompatibility complex, class ii, dr alpha	1.777805059	0.000258662	UP
HLA-DPA1	Major histocompatibility complex, class ii, dp alpha 1	2.468768075	0.000131723	UP
ITGB2	Integrin beta-2	1.538936861	0.000113105	UP
ITGAM	Integrin alpha-M	1.308463099	0.0000864	UP
HLA-DPB1	HLA class II histocompatibility antigen, DP beta 1 chain	1.727551628	0.00104215	UP
PTPRC	Receptor-type tyrosine-protein phosphatase C	2.737544685	0.0000112	UP
CCL2	C-C motif chemokine 2	1.772954596	0.000270229	UP
HLA-DMA	Major histocompatibility complex, class ii, dm alpha	1.526072847	0.000553678	UP
C1QB	Complement C1q subcomponent subunit B	2.021596465	0.001339804	UP

DOI: 10.7717/peerj.15437/table-2

Notes:

Annotation

FC

fold change

Q value

adjust P-value

Verification of the 20 hub genes by the GSE30122 database

The 20 hub genes were verified using the GPL17586- GSE30122 dataset, which consisted of 41 DN samples and 20 control samples. We used the R software ‘ggplot2’ package and ‘ggpubr’ package to construct box plots and perform Mann–Whitney U test statistical analysis. In line with our predictions, we found that the expression levels of the 20 hub genes in the DN group were different (Table 3), except for BST2 (Mann–Whitney U test, p = 0.753), IFITM2 (Mann–Whitney U test, p = 0.145), IFITM3 (Mann–Whitney U test, p = 0.161), and MX1 (Mann–Whitney U test, p = 0.054), which were higher than those in the control group (Fig. 6A). The expression levels of the hub genes in the control samples were decreased compared with those in the DN samples (Figs. 6B–6N).

Table 3:

Statistical analysis of Hub-gene between control and DN groups.

The t test was used when the observed variables were close to normal distribution in each group. When samples did not meet the requirements of normality test, Mann-Whitney U test (Wilcoxon Rank Sum test) was used.

GENE	t/T	Δ	95% CI	p
BST2	0.316	0.118	−0.627–0.863	0.753
C1QB	80	1.98	1.263–2.555	0.000
CCL2	5.128	1.075	0.657–1.494	0.000
GBP2	3.746	0.966	0.43–1.503	0.001
HLA-B	3.507	0.807	0.348–1.266	0.000
HLA-DMA	6.062	1.135	0.761–1.509	0.000
HLA-DPA1	6.375	1.408	0.967–1.849	0.000
HLA-DPB1	3.605	1.052	0.449–1.654	0.001
HLA-DRA	160	1.645	1.012–1.97	0.000
HLA-G	3.23	0.406	0.155–0.656	0.002
IFITM2	366	0.485	−0.203–1.174	0.145
IFITM3	1.418	0.388	−0.158–0.935	0.161
ISG20	211	0.686	0.324–1.175	0.000
IRF8	6.141	1.301	0.878–1.724	0.000
IRF9	4.487	0.462	0.257–0.668	0.000
ITGAM	4.841	1.038	0.592–1.485	0.000
ITGB2	5.155	1.296	0.773–1.819	0.000
MX1	331	0.459	−0.008–0.866	0.054
PSMB8	3.665	0.727	0.331–1.124	0.000
PTPRC	159	1.272	0.671–1.969	0.000

DOI: 10.7717/peerj.15437/table-3

Notes:

The t test was used when the observed variables were close to normal distribution in each group. When samples did not meet the requirements of normality test, Mann–Whitney U test (Wilcoxon Rank Sum test) was used.

The ROC curves of the 20 hub genes in DN and control groups

The R package ‘pROC’ was used to analyze the 20 hub gene expression profiles between the DN and control groups and paint the ROC curves. The AUC helped to evaluate sensitivity and specificity, and describe the original effectiveness of the diagnostic tests. Based on the statistical results, most of these hub genes had noteworthy diagnostic values in the DN group. Of these, CIQB had the largest AUC (AUC: 0.911). The ROC curves of genes with an AUC greater than 0.85 are shown in Fig. 7. We screened hub genes with AUC>0.85 to identify better candidate genes based on their good diagnostic performance in the DN group. C1QB, ITGB2, HLA-DPA1, ITGAM, and IRF8 levels were upregulated in DN with statistical significance after filtering, as determined by a p-value < 0.05 for upregulated genes. The specific p-values for each gene were: C1QB (Mann–Whitney U test, p < 0.01), ITGB2 (Mann–Whitney U test, p < 0.01), HLA-DPA1 (unpaired t-test, t = 6.375, p < 0.01), ITGAM (unpaired t-test, t = 4.841, p < 0.01), and IRF8 (unpaired t-test, t = 6.141, p < 0.01) (Table 4). Accordingly, we speculate that C1QB, ITGB2, HLA-DPA1, ITGAM, and IRF8 could be potential candidate genes for DN diagnosis based on our current results.

Table 4:

The evaluation of the hub-gene as prognostic predictor.

gene	AUC	95% CI	Optimal cut-off value	sensitivity	specificity
BST2	0.508	0.360–0.655	−0.801	0.38	0.778
C1QB	0.911	0.837–0.985	0.244	0.82	0.833
CCL2	0.838	0.731–0.945	0.163	0.72	0.833
GBP2	0.829	0.686–0.971	0.292	0.92	0.722
HLA-B	0.761	0.614–0.909	0.379	0.78	0.722
HLA-DMA	0.846	0.719–0.972	0.449	0.84	0.833
HLA-DPA1	0.884	0.795–0.974	0.388	0.84	0.833
HLA-DPB1	0.781	0.639–0.923	0.323	0.84	0.667
HLA-DRA	0.823	0.674–0.973	1.175	0.98	0.722
HLA-G	0.749	0.604–0.894	0.424	0.88	0.611
IFITM2	0.608	0.449–0.767	0.962	0.92	0.333
IFITM3	0.587	0.432–0.741	−1.27	0.22	1.000
ISG20	0.771	0.631–0.912	0.645	0.94	0.556
IRF8	0.852	0.749–0.955	0.149	0.72	0.833
IRF9	0.818	0.697–0.938	0.012	0.66	0.889
ITGAM	0.864	0.751–0.978	0.496	0.92	0.722
ITGB2	0.891	0.804–0.978	0.719	0.94	0.722
MX1	0.636	0.477–0.795	0.54	0.84	0.444
PSMB8	0.749	0.604–0.894	0.206	0.74	0.722
PTPRC	0.827	0.696–0.957	0.54	0.86	0.722

DOI: 10.7717/peerj.15437/table-4

Notes:

AUC

area under the ROC curve

CI

confidence interval

The co-regulatory network of miRNA-mRNA-TF

After searching the databases (miRWalk3.0, TargetScan, MiRDB, and starBase) for scores > 0.95, we identified 46 miRNAs acting on a 3′ UTR region of the genes as potential biological prognostic predictor. Subsequently, 11 TFs (HIF1A, KLF5, MBD1, RFX5, RFXANK, RFXAP, RUNX1, SP1, SPI1, STAT1, and WT1) were screened using the online database TRRUST (Table 5). Finally, using the results of the screening, we predicted 46 miRNA–mRNA pairs and 13 TF-mRNA pairs and then integrated the predicted pairs to describe a TF–miRNA–mRNA regulatory network constructed using Cytoscape (Fig. 8). Among these regulatory networks, we noted that IRF8 had the highest number of miRNA connections, which formed 25 miRNA-mRNA pairs in total. TFs with the highest concentration were found on ITGB2 including HIF1A, KLF5, RUNX1, SP1, and SPI1. It is noteworthy that of the five genes, only C1QB had no associated TFs based on the predicted outcomes.

Table 5:

Target gene co-regulatory network.

Target gene	TF	miRNA
ITGB2	none	hsa-mir-26b-5p hsa-mir-146a-5p hsa-mir-335-5p hsa-mir-148b-5p hsa-mir-17-3p hsa-mir-22-5p hsa-mir-3677-5p hsa-mir-182-5p hsa-mir-27a-3p
ITGAM	ARID1B POLR2A IRF1	hsa-mir-191-5p hsa-mir-26a-5p hsa-mir-372-3p
IRF8	ZNF143 NFIC TCF7 ZNF76 WRNIP1 MTA2 EBF1 CBFB TBX21 EED EZH2 POU2F2	hsa-mir-130a-3p hsa-mir-130b-3p hsa-mir-146a-3p hsa-mir-181a-5p hsa-mir-181c-5p hsa-mir-186-5p hsa-mir-218-5p hsa-mir-3605-3p hsa-mir-365a-3p hsa-mir-365b-3p hsa-mir-664a-3p hsa-mir-101-3p hsa-mir-146a-5p
HLA-DPA1	ATF1 WRNIP1 RFXANK HDGF CBFB ZNF382 STAT1 MAZ RFX5 PML CEBPB	hsa-mir-1343-3p hsa-let-7a-5p hsa-let-7b-5p hsa-let-7f-5p hsa-let-7g-5p hsa-let-7i-5p hsa-mir-1179 hsa-mir-146a-5p hsa-mir-146b-5p hsa-mir-148a-5p hsa-mir-148b-5p hsa-mir-21-3p hsa-mir-320a hsa-mir-320b hsa-mir-320c hsa-mir-320d hsa-mir-342-5p hsa-mir-3679-5p hsa-mir-374a-5p hsa-mir-374b-5p hsa-mir-589-5p hsa-mir-7-5p hsa-mir-98-5p hsa-mir-126-3p hsa-mir-129-2-3p hsa-mir-200c-3p hsa-mir-214-3p hsa-mir-26a-5p
C1QB	NR2F1 YBX1 ADNP	hsa-mir-26b-5p hsa-mir-124-3p hsa-mir-129-2-3p hsa-mir-146a-5p

DOI: 10.7717/peerj.15437/table-5

Notes:

TF

Transcription factors

miRNA

MicroRNA

Hub gene expression was positively correlated with the degree of kidney injury

In order to fully explore the clinical value of the selected five hub genes, we screened two gene expression datasets containing 186 and 60 CKD patients, respectively (Saint-Mezard et al., 2009; Ju et al., 2015). Correlation analysis between biological predictors and GFR was carried out in the Nephroseq V5 database (Figs. 9A–9E). All five biological predictors were negatively correlated with GFR, with the following Pearson correlation coefficients and p-values: C1QB (r =−0.601, p-value < 0.01), ITGB2 (r =−0.506, p-value < 0.01), HLA-DPA1 (r = −0.589, p-value < 0.01), ITGAM (r = −0.560, p-value < 0.01), and IRF8 (r =−0.501, p-value < 0.01). Thus, a higher expression of those hub genes indicated worse renal function in patients with CKD, which may play a role in kidney deterioration and damage in patients with CKD.

Validation of five hub genes’ expression with qRT-PCR

To further assess the expression of these hub genes, a total of three DN model rats and three healthy controls were enrolled as a validation cohort. We first measured serum creatinine, BUN, and albumin levels, and calculated the urinary protein/creatinine ratio in the DN rat group and control group. We found that the level of kidney injury markers, such as serum creatinine and BUN, in the DN group was higher than in the control group (unpaired t-test, t = 3.391, df = 4, p = 0.0275, r = 0.861). Meanwhile, the DN group had a higher urinary protein/creatinine ratio (unpaired t-test, t = 17.23, df = 16, p < 0.001, r = 0.974), suggesting more significant proteinuria (Fig. 10A). RT-PCR technology was then used to confirm the differential expression levels from participant serum and kidney samples. Consistent with the microarray data, ITGAM (unpaired t-test, t = 10.35, df = 4, p < 0.001, r = 0.981) and C1QB (unpaired t-test, t = 13.60, df = 4, p < 0.001, r = 0.989) expression were significantly upregulated in both serum and kidney samples (Figs. 10B and 10C) between DN and control groups. At the same time, ITGB2 was highly expressed in the kidney samples in the disease group (unpaired t-test, t = 5.512, df = 4, p = 0.0053, r = 0.940). However, there was no difference in the levels of HLA-DPA1. Notably, the expression level of IRF8 in the DN group was lower than in the control group, which was contrary to our prediction.

Discussion

DN is usually diagnosed on the basis of increased urine albumin-creatinine ratio (UACR) or decreased estimated glomerular filtration rate (eGFR) when a patient has no other CKD. However, DN diagnosis relies heavily on assumptions and updated specific DN markers are urgently needed, especially for early-stage disease (Forbes & Cooper, 2013; Ioannou, 2017). Therefore, in this study, based on the GSE30529 microarray data taken from the GEO database, we screened DEGs between the DN and control group and identified 463 DEGs (340 increased genes and 123 decreased genes).

According to our GSEA, the expression levels of the immune inflammatory response related signaling pathways, such as the interleukin and neutrophil degranulation pathways were upregulated in the DN group. This may reveal the cause of glomerular vascular injury, glomerulosclerosis, nodular lesions, and deterioration of renal function. The GO analysis concluded that DEGs were involved in neutrophil activation, neutrophil-mediated immunity, neutrophil degranulation, neutrophil activation, and response to interferon-gamma, which indicates that the level of immune response and regulation in DN group was more obvious than that in the control. In addition, the expression of DEGs may act on the molecular functions of glycosaminoglycan binding, heparin binding, sulfur compound binding, endopeptidase inhibitor activity, and peptidase inhibitor activity, indirectly demonstrating the complexity of the DN’s pathogenesis. The enriched pathways identified by KEGG analysis included phagosome, staphylococcus aureus infection, complement and coagulation cascades, leishmaniasis, cell adhesion molecules, tuberculosis, rheumatoid arthritis, and malaria. Thus, based on GSEA, GO, and KEGG enrichment analysis results, DN showed strong immune activation and signal transduction, which is in agreement with previous studies (Wada & Makino, 2013).

We used the GSE30122 microarray dataset, including 19 DN and 50 control samples for verification after constructing a PPI network to identify 20 hub genes. Statistical analysis of hub gene expression levels explained that the expression level of most hub genes in DN was higher than in the control group, and there was a statistically significant difference (P < 0.05). In the ROC curve analysis, the AUC of C1QB, HLA-DPA1, IRF8, ITGAM, and ITGB genes were higher than 0.85, indicating a higher DN predictive value. Therefore, based on the five genes screened above, we constructed a miRNA-mRNA-TF regulatory network using the online database to try to interpret the pathogenesis of DN in the process of gene expression. We then evaluated the association of selected genes with kidney injury using the Nephroseq V5 database. At the same time, we further verified the expression levels in both serum and kidney samples of these five genes using the qPCR test and a DN rat model, and the results confirmed our previous prediction.

Complement C1q B chain (C1QB) is a protein coding gene. C1q associates with the proenzymes C1r and C1s to yield C1, the first component of the serum complement system (Reid, 2018). In the blood plasma, the complement system is a proteolytic cascade that mediates innate immunity, a non-specific defense mechanism to pathogens (Daugan et al., 2017). The main consequences of complement activation are the opsonization of pathogens, recruitment of inflammatory and immunocompetent cells, and direct killing of pathogens. Although complement cascades are an important part of the innate immune system, uncontrolled activation can lead to serious disease. One recent review found that excessive complement activation in atypical hemolytic uremic syndrome resulted in renal failure if not treated (Willows, Brown & Sheerin, 2020). Experimental and clinical evidence suggests that several components of the complement system are involved in the pathogenesis of DN (Flyvbjerg, 2017). The deposition of membrane attack complex formation (MAC) in the complement system is believed to be closely related to the pathogenesis of advanced DN (Fearn & Sheerin, 2015). Lu et al.’s (2022) study suggested that C1 and C2 could be two distinct immune-associated genes in the pathogenesis of DN. Sun et al. (2019) found an association between C1q and C3c complement deposition on renal histopathology and more severe kidney damage in DN patients. Our results provide evidence that C1QB is highly expressed in DN models and also emphasize the roles of C1 and local complement activation in DN, while the expression of C2 and C3 related genes is worthy of further investigation.

The ITGAM gene encodes the integrin alpha M chain, and ITGB2 encodes an integrin beta chain that combines with multiple different alpha chains to form different integrin heterodimers. Integrins are integral cell–surface proteins that participate in cell adhesion as well as cell–surface mediated signaling (Siegers, 2018). Integrin ITGAM/ITGB2 is a receptor for the iC3b fragment of the third complement component, which is involved in the composition of the complement system that enhances macrophage apoptotic neutrophil phagocytosis (Kristóf et al., 2013). ITGB2 was identified as a central gene in the complement cascade pathway that is negatively correlated with glomerular filtration rate (GFR; Xu et al., 2021). Complement activation contributes to MAC formation, which drives neutrophil activation and endothelial injury (Riedl et al., 2014). Meanwhile, integrin ITGAM/ITGB2 can recognize P1 and P2 peptides of fibrinogen gamma chain and is also a receptor for coagulation factor Xa (FXa). Fibrinogen is a soluble glycoprotein that plays an important role in the coagulation system (Mosesson, 2005) and inflammation (Davalos & Akassoglou, 2012). There is increasing evidence showing that serum fibrinogen, an acute phase marker of inflammation, plays an important role in regulating the inflammatory response and aggravating the progression of renal disease by binding to receptors expressed on the surface of different cells (Wang et al., 2017). Elevated plasma fibrinogen levels have been shown to be an indicator of a proinflammatory state and are closely associated with diabetes (Le et al., 2008). In addition, FX is mainly synthesized in the liver (Borensztajn, Peppelenbosch & Spek, 2008). However, there is evidence that FX extrahepatic expression by macrophages contributes to the development of asthma (Shinagawa et al., 2007). FXa activates PAR-1 and PAR-2 to enhance inflammation (Jose & Manuel, 2020). DM directly increased FX expression in macrophages (Oe et al., 2016). Local FX synthesis in macrophages infiltrating the kidney may contribute to DN progression. In our study, although these two genes were consistent as previously predicted and showed differences in expression between different groups of kidney tissue samples, we found that the level of differences in expression between groups of ITGAM was more significant and statistically maintained in blood samples. Previous studies did not examine the gene expression levels in specific samples in order to focus on differences in expression between the two (Hu et al., 2021). Our validation test results suggests that ITGAM can be a better predictor of DN than ITGB2. Blood samples are more accessible to obtain than kidney samples in clinical practice, which further increases the value of testing ITGAM to predict DN.

HLA-DPA1 is a member of the class II alpha chain HLA mimicry. Class II is a heterodimer consisting of an α chain (DPA) and β chain (DPB), both anchored to the cell membrane. It plays an important role in the immune system and manifests itself as peptides produced by extracellular proteins. Class II molecules are expressed in antigen presenting cells (APC) such as B lymphocytes, dendritic cells, and macrophages. Previous studies have shown that HLA-DPA1, the closest centromere gene expressed in HLA-DO α, may increase the risk of diabetes and diabetic kidney disease (DKD), which is a complication of diabetes (Varney et al., 2010). This gene also plays a role in the rheumatoid arthritis (RA) and systemic lupus erythematosus (SLE) related pathways (Ma et al., 2017). Consistent with Ma et al.’s (2017) research, the conclusion obtained from our preliminary analysis based on the microarray suggested that HLA-DPA1 was highly expressed in the disease group, while our q PCR verification results suggested that there was no statistical difference between the groups. These results show that further research is needed to determine whether HLA-DPA1 in DN can affect pathologic processes, similar to those in RA and SLE.

The interferon (IFN) regulatory factor (IRF) family specifically binds to the upstream regulatory region of type I IFN and IFN-inducible MHC class I genes (Salem et al., 2014). IRF family proteins bind to IFN-stimulated response elements (ISRE) and adjust the expression of genes, namely IFN-α and IFN-β. In diabetic mice, the expression of IRF8 is closely related to the activation of microglia and the regulation of the inflammatory response, thereby influencing the occurrence of diabetic retinopathy (Liu, Xu & Gao, 2021). However, the role of IRF8 in DN, which is also a microvascular complication of diabetes, still needs to be further investigated. As a TF, IRF8 has a complex regulatory relationship with other TFs. STAT1 binds to the IRF promoter region in colon cancer cells, while MBD1 inhibits IRF8 expression (McGough et al., 2008), which is another indication evidence that WT1 in leukemia and IRF8 are anticorrelated (Vidovic et al., 2010). Therefore, it is necessary to further explore the interaction between IRF8 and other transcription related factors in DN. In addition, the results from the online database indicated an inverse relationship between IRF8 and renal function, while qPCR indicated a low expression level of IRF8 in the DN group. These differences may be caused by the different regulatory effects of IRF8 as a TF in different species.

Notably, in the miRNA-mRNA-TF networks, all of the 11 predicted TFs have been proven to be key modulators and potential therapeutic targets for various inflammatory diseases. The TF SPI1 could regulate the expressions of both ITGAM and ITGB2.SPI1 encodes an ETS-domain TF that activates gene expression during myeloid and B-lymphoid cell development (Wittwer, Marti-Jaun & Hersberger, 2006). This was consistent with our analysis which showed that SPI1 expression was significantly increased in DN, suggesting that DN progression may be related to bone marrow cell activation and B-lymphoid cell development.

MicroRNA is one of the most important endogenous epigenetic factors that can inhibit the post-transcriptional gene expression of target genes. Recent studies have revealed that miRNAs may be regulators of immune and inflammatory responses and are potential therapeutic targets in DN (Kato & Natarajan, 2015; Zhou et al., 2021). Among the miRNAs predicted by the online databases, some have been shown to be closely related to the occurrence and development of diabetes. For example, hsa-miR-424-5p binds to PD-1 signaling molecules, stimulates the immune response through the mTORC signaling pathway, and is involved in the pathogenesis of type 1 diabetes (Wang et al., 2020). There is evidence that these miRNAs also act as protective factors in other diseases. Mir-5581-3p plays an inhibitory role in the progression of hepatocellular carcinoma (HCC) by regulating the expression of cardiolipin synthase 1 (CRLS1; Yin et al., 2019). The lower expression level of mir-8485 is associated with significantly lower overall survival rate of oral squamous cell carcinoma (OSCC; Gholizadeh et al., 2020). Mir-8485 can also alleviate cardiomyocyte injury in chronic heart failure (CHF) by targeting TP53INP1 (Luo et al., 2022). However, whether the predicted miRNAs can participate in regulating cell viability, apoptosis, inflammatory responses in DN, and disease progression still deserves further investigation.

Although this is the first study to construct a miRNA-mRNA-TF regulatory network to explore potential predictors of DN through two microarray datasets, this study still had several limitations. First, due to the small number of data sets integrated, only data sets related to DN genes were obtained from NCBI database, so other undiscovered DEGs still need to be explored based on multiple databases. Second, although identified hub genes have been reported to be closely related to the regulation of immune responses, there is still little evidence that they can clearly regulate immune responses in the DN disease process (Li et al., 2021). The observed high expression of immune-related genes in DN may be due to the enrichment of immune cells in DN tissue. Thus, more precise and detailed studies at the cell or animal level, such as single cell sequence studies, are needed to further validate the potential factors identified in our bioinformatics findings. Third, our qPCR results were based on rat models, but there are differences in gene expression between species. Lastly, we recognize that the number of laboratory animals used in our study was limited. The small sample size may have led to inconclusive results, and validation with larger cohorts of animals is needed to strengthen the evidence supporting our bioinformatics findings. Therefore, our prediction results need to be carried out across large clinical cohort studies or randomized controlled studies to evaluate its actual clinical value.

Conclusion

In summary, we conducted comprehensive analysis and validation based on multiple microarray datasets, and the results showed that C1QB, ITGAM, and ITGB may be potential candidate genes that are upregulated in DN. These genes have remarkable diagnostic properties. The predicted miRNA-mRNA-TF pathway network was built based on online databases. Our results at the transcriptional level may operate as reference when exploring the pathogenesis of DN in the future. It is important that these interactions be treated with caution given the complexity of miRNA and TF crosstalk.

Supplemental Information