Identification of Fibrotic Biomarkers Associated with Macrophages in Diabetic Nephropathy

Background

Diabetic nephropathy (DN) is the main cause of end-stage renal disease. The annual mortality rate of diabetic patients with end-stage renal disease is approximately 20%, which is higher than that of many solid cancers, such as prostate, breast, or renal cell carcinoma [1]. Renal fibrosis is an important pathological feature of kidney injury. The entire kidney’s histological structure can be affected by fibrosis, specifically termed glomerulosclerosis, interstitial fibrosis, and perivascular fibrosis of vasculature [2]. As fibrosis progresses, kidney function is gradually lost, leading to renal failure. Currently, there are no good therapeutic means to cure or reverse renal fibrosis. The kidney has limited regenerative potential, especially the glomeruli, which renders treatment of renal fibrosis more difficult. The lack of diagnostic approaches to assess renal fibrosis is one of the factors affecting treatment [3]. Therefore, validating new biomarkers is of great significance for improving the prognosis of patients with DN.

Both the innate and adaptive immune systems play key roles in the development of DN [4]. In recent years, the role of macrophages in DN fibrosis has attracted attention [5]. Increased reactive oxygen species have been identified as the single unifying upstream event of DN [6], which mediates the infiltration of macrophages and the production of inflammatory factors, while the release of inflammatory mediators can further promote the production of reactive oxygen species [7]. It is suggested that there is crosstalk between oxidative stress and macrophages recruitment, and the two influence each other to promote the renal failure in patients with DN. However, few studies have been conducted from the perspective of combining macrophage recruitment with oxidative stress.

In the present study, we performed a systematic bioinformatics profiling via integrating oxidative stress as well as immune infiltration to identify the biomarkers for DN. Our results provide key oxidative stress-related genes associated with macrophages in DN that may reflect the progression of renal fibrosis.

Material and Methods

DATA SOURCE:

In this study, we obtained a gene expression microarray from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo). The microarray data were obtained from the platform GPL17586 with the dataset GSE96804 (DN=41, control=20). We obtained 464 oxidative stress-related genes from the Molecular Signatures Database (MSigDB).

IMMUNE AND STROMAL CELL ANALYSES:

We applied the “xCell” package [8] to determine the enrichment score for each cell type and forecast the levels of infiltrating immune cells in the glomeruli of DN patients. xCell is a cell enrichment tool taking advantage of single-sample gene set enrichment analysis that calculates the scores for 64 cell types, including 34 immune cells, 30 stromal cells, and other cells. It reduces dependencies between closely related cell types through a spillover compensation technique.

WEIGHTED GENE CO-EXPRESSION NETWORK ANALYSIS AND MODULE GENE SELECTION:

Weighted gene co-expression network analysis (WGCNA) was performed via the R package “WGCNA” to identify the modules with the highest relevance to immune cells in patients with DN [9]. For the first step, hierarchical clustering was implemented on the specimens to uncover the outlier specimens and eliminate the abnormal samples. After that, to establish a scale-free network, the optimal “soft” power was determined relying on the function pickSoft Threshold. For the next step, the correlation matrix was constructed by the WGCNA software package. The optimal soft threshold was chosen to convert the correlation matrix into an adjacency matrix, and then transformed to a topological overlap matrix. Genes with similar expression patterns were categorized into gene modules, which are represented by branches of the clustering tree and different colors. Then, the correlation coefficients between differentially infiltrating immune cells and each module were counted. The modules with high linkage were deemed as candidates linked to differentially infiltrating immune cells and chosen for subsequent analyses. For the selected candidate module, we set |module membership| >0.80 and |gene significance| >0.20 as the screening standard for deciding key genes in the candidate module.

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED IMMUNE-RELATED OXIDATIVE STRESS GENES:

The differentially expressed genes (DEGs) between DN and normal specimens in the GSE96804 dataset were identified deploying the “limma” R package. The parameters |log2FC| >0.585 and adjusted P<0.05 were set as the criterion for DEGs. On the side, the heatmap and volcano plot were produced with R package “pheatmap” and “ggplot2”. The DEGs were intersected with 464 oxidative stress-related genes to obtain differentially expressed genes related to oxidative stress (DEOSGs). Subsequently, the crossover of DEOSGs and genes in key modules was done through the “Venn Diagram” R package and named as DEIOSGs. Also, the Spearman correlation coefficient between DEIOSGs and immune cells was counted via the “corrplot” in R package, and the functional comparability of DEIOSGs was analyzed via the “GOSemSim” package.

KYOTO ENCYCLOPEDIA OF GENES AND GENOMES AND GENE ONTOLOGY ENRICHMENT ANALYSES:

To reveal the potential biological functions and underlying mechanisms of DEIOSGs, we used the R package “cluster Profiler” to implement Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses [10]. The GO terms biological processes, cellular components, molecular functions, and KEGG pathways with adjusted P<0.05 were considered statistically significant.

CONSTRUCTION OF PROTEIN–PROTEIN INTERACTION NETWORK AND SCREENING OF HUB GENES:

To assess protein functional relationships, all DEIOSGs were linked to the Search Tools for the Retrieval of Interacting Genes database [11] (confidence=0.15), and the results were then imported into Cytoscape, which was used to visualize the protein–protein interaction (PPI) network. The plug-in in the Cytoscape software, molecular complex detection, was used to identify significant gene clusters. The nodes with the most interactions with neighboring nodes were considered as the key nodes. According to the importance of key nodes, the top 3 are selected as hub genes. Since the third of the top 3 had 2 genes in our research, 4 genes were actually obtained.

RECEIVER OPERATING CHARACTERISTIC CURVE AND EXPRESSION ANALYSIS:

Receiver operating characteristic (ROC) curves of hub genes were drawn to prove their accuracy. Specifically, the “pROC” package was used for the ROC curve profile. Relying on the area under the ROC curve (AUC) result, the hub genes with AUC >0.7 were considered as promising biomarkers for DN diagnosis [12]. The expression hub genes between DN and control specimens were visualized in the boxplots created by the “ggplot2” R package. In addition, the 4 hub genes were validated in an external dataset, GSE30528.

IN VITRO ANALYSES:

The mouse podocyte cell line, MPC5, was cultured in RPMI-1640 culture media containing 10% fetal bovine serum, 100 U/mL penicillin, and 100 μg/mL streptomycin. To induce a hyperglycemic condition, MPC5 cells were incubated in high glucose (HG; 30 mM) for 72 h, in comparison with normal glucose (5.5 mM). The expression of hub genes was measured by quantitative real-time polymerase chain reaction, according to an earlier approach [13]. All primers were used in the synthesis (GenePharma, Shanghai, China), and sequences are listed in Table 1.

Reducing sugars can react non-enzymatically with the amino groups of proteins to initiate a complex series of rearrangements and dehydrations, and then to produce a class of irreversibly cross-linked moieties called advanced glycation end products (AGEs) [14]. The AGEs can bind to AGEs receptors on the surface of macrophages to induce excessive inflammation and oxidative stress [15]. The human monocytic leukemia cell line, THP-1, was cultured in RPMI-1640 medium with 10% fetal bovine serum. Differentiation of THP-1 macrophages was induced with 100 ng/mL phorbol 12-myristate 13-acetate, followed by treatment with 200 μg/mL AGEs (Abcam, USA) for 24 h. The expression of hub genes was determined by quantitative real-time polymerase chain reaction. Tumor necrosis factor alpha (TNF-α; PT512, Beyotime Biotechnology, Shanghai, China) and superoxide dismutase (A001-3-1, Jiancheng Bioengineering Institute, Nanjing, China) were detected using enzyme-linked immunosorbent assay on the cell culture supernatant, following the manufacturer’s instructions.

STATISTICAL ANALYSIS:

GraphPad Prism 8.0 (GraphPad Software, CA, USA) and R software (version 4.2.1) were used for statistical analyses. Correlation tests were performed using Spearman analysis. The diagnostic value of hub genes was evaluated with ROC curve analysis. Hub genes were analyzed for correlation with clinical features via Pearson analysis. An unpaired t test was performed for comparisons between 2 groups. The results were considered statistically significant when the P values were less than 0.05.

Results

IMMUNE INFILTRATING CELL ANALYSIS:

Studies have indicated that immune infiltration plays an important role in the development of DN [4]. We further explored xCell to uncover the cell types latently implicated in glomerular injury in the diabetic kidney (Figure 1). The enrichment scores of 64 cell types were obtained for each sample. Box-plot showed that 18 types of immune cells were increased in DN patients, including dendritic cells, CD4+ T cells, macrophages, M1 macrophages, and T helper 2 cells; in contrast, 3 types of immune cells, including neutrophils, natural killer cells, and T helper 1 cells, decreased. The results suggested that different immune cells might play different roles in the development of DN.

CONSTRUCTION OF CO-EXPRESSION NETWORKS:

WGCNA was used to identify the function of immune-related genes in DN. We chose β=9 (scale-free R2=0.86) as the “soft” threshold based on the scale independence and average connectivity (Figure 2A). Accordingly, 10 gene co-expression modules were generated, which are presented in different colors (Figure 2B). Proceeding further, via relevance analysis of the traits (immune cells) and modules, we detected that the yellow module was highly and positively linked to dendritic cells (Cor=0.89, P=1.2e–21), macrophages (Cor=0.81, P=2.6e–15) and M1 macrophages (Cor=0.75, P=3.2e–12) (Figure 2C). Therefore, the yellow module, containing 1000 key genes, was selected as the critical module linked to immune infiltrating cells for further analysis. The significant correlations between gene significance for macrophages/M1 macrophages and module membership are presented in the yellow module (Figure 2D).

IDENTIFICATION OF DEGS:

Differential expression analysis was implemented in the GSE96804 dataset. In total, 1507 DEGs (|log2FC| >0.585 and adjusted P<0.05) were derived, of which 568 DEGs were upregulated and 939 were downregulated. The volcano plot of the DEGs is shown in Figure 3A. The top 5 genes with the most distinct differences were MIR4521, MIR4256, MIR1256, VTRNA1-3, and SPDYE8P among the upregulated genes, and FOS, DUSP1, CXCR2, PDK4, and HSD17B14 among the downregulated genes. Meanwhile, we used the heatmap to show the top 20 genes with the most obvious differential expression between the DN and control samples (Figure 3B).

ENRICHMENT FUNCTION ANALYSIS OF DEGS:

To explore the role of DEGs in DN occurrence and development, GO and KEGG enrichment analyses were performed. In the biological processes analysis, DEGs mainly participated in small molecule catabolic process, urogenital system development, and organic acid catabolic process. In the cellular components analysis, DEGs significantly participated in collagen-containing extracellular matrix (ECM), apical part of cell, apical plasma membrane, mitochondrial matrix and peroxisomal matrix. The molecular functions result suggested that DEGs were notably linked to glycosaminoglycan binding, ECM structural constituent, and oxidoreductase activity (Figure 3C). KEGG analysis showed that the DEOSGs were primarily enriched in focal adhesion, carbon metabolism, protein digestion and absorption, AGE-RAGE signaling pathway in diabetic complications, and ECM-receptor interaction (Figure 3D).

IDENTIFICATION OF DEIOSGS AND FUNCTIONAL ENRICHMENT ANALYSIS:

A total of 464 oxidative stress-related genes were obtained from MSigDB. From the intersection of the 1507 DEGs, oxidative stress-related genes, and genes in the yellow module, 8 DEIOSGs were obtained, as shown in the Venn diagram in Figure 4A. Spearman analysis revealed significant positive correlations among the 8 DEIOSGs. Meanwhile, correlation analysis showed a significant positive correlation between the 8 DEIOSGs and the infiltration of fibroblasts, dendritic cells, macrophages, and M1 macrophages (Figure 4B, 4C). Functional similarity analysis suggested that the 8 DEIOSGs had certain similarities in function (Figure 4D).

To further realize the functions of the 8 DEIOSGs in DN, the GO terms as shown in Figure 5A–5C were used. In the biological processes analysis, DEIOSGs were mainly linked to response to oxidative stress, anatomical structure formation, circulatory system development, and ECM organization. In the cellular components analysis, DEIOSGs were notably linked in the ECM, collagen-containing ECM, myofibril, and death-inducing signaling complex. Molecular functions results suggested that DEIOSGs were notably linked to endopeptidase activity, heme binding, protease binding, interleukin-1 receptor antagonist activity, and cysteine-type endopeptidase activity involved in the apoptotic signaling pathway. KEGG analysis showed that DEIOSGs were primarily enriched in the AGE-RAGE signaling pathway in diabetic complications and the TNF signaling pathway (Figure 5D).

EXPRESSION ANALYSIS, ROC CURVE ANALYSIS, AND HUB GENE VALIDATION FOR DN DIAGNOSIS:

The 4 most important genes were selected by the PPI network (Figure 6A), namely MMP2, COL1A1, HIF-1α, and CASP3, which were significantly upregulated in the DN samples compared with in the control samples in the GSE96804 dataset (Figure 6B). To investigate the validity of these 4 genes as the diagnostic biomarkers for DN, the creation of hub single-gene ROC curves was run (Figure 6C). Notably, the AUC value for CASP3 was 0.948, which was the highest. The AUC values for COL1A1, MMP2, and HIF-1α were 0.924, 0.896, and 0.759, respectively. We next validated the 4 hub genes using an external dataset, GSE30528, as shown in Figure 6D. Except for COL1A1, the AUC areas of all other genes were >0.7 in the GSE30528 dataset. The AUC values for MMP2, CASP3, and HIF-1α were 0.966, 0.803, and 0.795, respectively. The above-mentioned results illustrated that the 3 genes had good diagnostic value.

IN VITRO ANALYSES:

In MPC5 cells, compared with the control group, the HG group had increased relative mRNA expression of MMP2 (P=0.0009), CASP3 (P=0.0049), and HIF-1α (P=0.001) (Figure 7A). Further, MMP2 (P=0.0064), CASP3 (P=0.0045), and HIF-1α (P=0.0015) were also highly expressed in THP-1 macrophages treated with AGEs (Figure 7B). Enzyme-linked immunosorbent assay showed decreased superoxide dismutase activity and increased TNF-α levels after AGE stimulation in THP-1 cells (Figure 7C, 7D).

Discussion

As one of the common microvascular complications of diabetes, DN is a primary contributor of end-stage renal disease, which imposes a huge medical and financial burden. Renal fibrosis is the best predictor of disease progression, and currently, the treatment methods are very limited. Macrophage infiltration and oxidative stress perform vital roles in renal fibrosis and are inseparable [16]. Some studies have identified biomarkers related to immune infiltration in DN [17], but few studies have focused on both mechanisms.

In this study, 8 DEIOSGs were identified. The enrichment results showed that the DEIOSGs were mainly involved in oxidative stress response, anatomical structure formation, circulatory system development, and ECM organization. These factors play an important role in diabetic renal fibrosis. Excessive ECM deposition is an established hallmark of renal fibrosis [18]. KEGG analysis revealed that the DEIOSGs were primarily enriched in the AGE-RAGE signaling pathway in diabetic complications and TNF and MAPK signaling pathways, consistent with the current understanding of fibrogenic pathways stimulated by HG [19]. Earlier studies showed that the nonsteroidal anti-inflammatory drug-activated gene-1, or growth differentiation factor-15, significantly reduced glycogen deposition and renal fibrosis via targeting AGE-RAGE-mediated inflammation pathways [20]. Neuregulin 4 ameliorated tubulointerstitial fibrosis and attenuated the expression of AGEs in DN through TNF-R1 signaling [21]. Multiple drugs suppressed HG-induced renal fibrosis and inflammation by inhibiting the MAPK pathway [22,23].

Further, we noted a huge difference in immune cell levels between DN and control specimens. Among the 18 types of immune cells increased in patients with DN, dendritic cells, CD4+ T cells, macrophages, M1 macrophages, and T helper 2 cells were included. The results suggest that both the innate and adaptive immune systems play a significant role in the pathogenesis of DN. Previous studies have shown that there is a high proportion of macrophages, which cause podocyte damage through multiple pathways, in patients with DN [24]. M1 macrophages and M2 macrophages are closely related to inflammatory responses, which play important roles in the process of fibrosis [25]. M1 is involved in the initiation phase of tissue damage and is related to proinflammatory functions. In contrast, the M2 phenotype is involved in the repair phase by producing anti-inflammatory cytokines, growth factor, and proangiogenic cytokines [26]. In general, in the process of kidney fibrosis, first, resident immune cells produce chemoattractant substances that enhance inflammatory responses by recruiting more M1; after that, the change in tissue environment shapes the macrophage phenotype toward M2, with functional properties that meet the tissues needing to address the danger [5]. The time course of the different phenotypes is not clear. Studies have shown that reduction of macrophage (M1 or M2) recruitment by the use of agents able to reduce its accumulation in diabetic kidneys improves renal function and fibrosis [27,28]. In our study, all DEIOSGs were positively correlated with macrophages and M1 macrophages. This may be associated with different stages of the patients. Reactive oxygen species are essential for the induction and maintenance of the M1 type of macrophage polarization [29]. Therefore, the results suggested that all 8 genes might play a part in kidney fibrosis and likely influence the development of DN through the recruitment or polarization of macrophages.

The 3 hub genes (MMP2, CASP3, and HIF-1α) were selected by the PPI network and validated in an external dataset. MMP2 is involved in ECM remodeling and glomerular basement membrane structure destruction [30]. Our results indicated that MMP2 was a key gene participating renal fibrosis, which is consistent with previous studies [31]. The oxidation state of the cysteinyl sulfydryl in the propeptide residue of MMP2 is important for its activity [32]. Oxidative stress can contribute to the increase in MMP2 activity, possibly by S-glutathiolation, and this can result in calponin-1 loss and maladaptive vascular remodeling in early hypertension [33]. It was proposed that macrophage recruitment to and activation at overgrowing tissue is a multi-step process requiring reactive oxygen species- and JNK-mediated MMP2 upregulation [34]. The MMP2 knockout mice showed improvement in renal fibrosis and inhibition of macrophage infiltration [35]. However, other studies have shown that MMP2 degrades a variety of ECM proteins, and its elevation is a compensatory mechanism in DN [36]. Combination therapy of antifibrotic drugs stimulated MMP2 expression or activity and suppressed macrophage infiltration in unilateral ureteral obstruction models [37,38]. The conflictions may be related to the different models and stages of kidney disease, which need to be confirmed by further studies.

Oxidative stress and inflammation are important factors of programmed cell death. Apoptosis and pyroptosis, leading to the reduction of parenchymal cells, are different forms of programmed cell death in renal fibrosis. Caspases are crucial mediators of programmed cell death. Among them, CASP3 is the most abundant of the executioner caspases [39]. CASP3 is involved in several signaling pathways leading to renal cell death and is indispensable for chromatin condensation and DNA fragmentation [40,41]. Apoptotic cells or cell residues have a strong chemotactic effect on macrophages. Studies have shown that TNFα/Casp3/Gasdermin E-mediated pyroptosis is responsible for the initiation of ureteral obstruction-induced renal tubule injury, and CASP3 deletion alleviated renal tubule damage and inflammation and finally prevented the development of kidney fibrosis [42]. CASP3 may influence the progression of DN by participating in macrophage recruitment.

HIF-1α (the star molecule in diabetes research) is an oxygen-sensitive monitoring and regulatory protein, although its specific role in DN remains unclear. DN is characterized by renal hypoxia, oxidative stress, and defective nutrient deprivation signaling, which cause both HIF-1α activation and HIF-2α suppression [43]. The shift in the balance of the 2 inducible isoforms promotes proinflammatory and profibrotic pathways in glomerular and renal tubular cells [44]. A key role of HIF-1α in response to hypoxia is to increase oxygen delivery by promoting angiogenesis, which is accompanied by the production of proinflammatory cytokines and activation of profibrotic mechanisms [45]. In recent years, sodium-dependent glucose cotransporters-2 inhibitors and glucagon-like peptide 1 receptor agonists have been given serious attention for their renal protective effects. Interestingly, the renal protective effects of both drugs may be related to the SIRT1-HIF-1α pathway [46,47]. HIF-1α is the downstream molecule of Sirt1 under anoxic conditions; inhibition of HIF-1α can reduce inflammation and fibrosis in rat mesangial cells cultured with HG [48]. Recent studies showed that the decrease in HIF-1α hydroxylation in diabetic renal macrophages increased its stability and expression. The upregulation of HIF-1α promoted the glycolysis of macrophages, thus affecting the polarization of macrophages [49]. Our results indicated that HIF-1α was positively related to M1 macrophages. In vitro studies confirmed that HIF-1α was highly expressed in HG-treated podocytes and AGE-treated THP-1 macrophages. In addition, the secretion of TNF-α increased and the activity of superoxide dismutase decreased in AGE-stimulated THP-1 cells. These findings indicate that HIF-1α might play important roles in the polarization of macrophages in DN fibrosis through regulating oxidative stress and inflammatory cytokines secretion. Notably, HIF proline hydroxylase inhibitors have been approved for use in renal anemia, the application of which can contribute to the recovery of HIF-1 function, thereby reducing the overproduction of reactive oxygen species and protecting against kidney injury [50]. Therefore, the long-term effects of HIF-1α on chronic kidney disease may be influenced by many factors, including the primary disease and disease stage, and need to be further investigated.

The strength of this study is the use of bioinformatics analysis and experiments to explore key oxidative stress-related genes associated with macrophages in DN, which is beneficial to better understand the pathogenesis of renal fibrosis. The hub genes have predictive value for the progression of renal fibrosis. Therapeutic strategies targeting macrophages are promising. However, this study still had a few limitations. First, validation of the hub genes in in vivo models or clinical samples is needed to confirm the current conclusions. Second, the hub genes may be involved in the progression of fibrosis by regulating oxidative stress and macrophage recruitment or polarization. Further experiments are needed to explore the exact underlying mechanisms and their interactions with other factors.

Conclusions

Our results provide new biomarkers for DN reflecting the progression of renal fibrosis. These genes were closely associated with oxidative stress and macrophage infiltration. As the interaction between oxidative stress and macrophage recruitment or polarization is complex, further research is needed to identify effective strategies.