The landscape of immune cell infiltration in the glomerulus of diabetic nephropathy: evidence based on bioinformatics
BMC Nephrology volume 23, Article number: 303 (2022)
Increasing evidence suggests that immune cell infiltration contributes to the pathogenesis and progression of diabetic nephropathy (DN). We aim to unveil the immune infiltration pattern in the glomerulus of DN and provide potential targets for immunotherapy.
Infiltrating percentage of 22 types of immune cell in the glomerulus tissues were estimated by the CIBERSORT algorithm based on three transcriptome datasets mined from the GEO database. Differentially expressed genes (DEGs) were identified by the “limma” package. Then immune-related DEGs were identified by intersecting DEGs with immune-related genes (downloaded from Immport database). The protein–protein interactions of Immune-related DEGs were explored using the STRING database and visualized by Cytoscape. The enrichment analyses for KEGG pathways and GO terms were carried out by the gene set enrichment analysis (GSEA) method.
11 types of immune cell were revealed to be significantly altered in the glomerulus tissues of DN (Up: B cells memory, T cells gamma delta, NK cells activated, Macrophages.M1, Macrophages M2, Dendritic cells resting, Mast cells resting; Down: B cells naive, NK cells resting, Mast cells activated, Neutrophils). Several pathways related to immune, autophagy and metabolic process were significantly activated. Moreover, 6 hub genes with a medium to strong correlation with renal function (eGFR) were identified (SERPINA3, LTF, C3, PTGDS, EGF and ALB).
In the glomerulus of DN, the immune infiltration pattern changed significantly. A complicated and tightly regulated network of immune cells exists in the pathological of DN. The hub genes identified here will facilitate the development of immunotherapy.
Diabetic nephropathy (DN) is an extremely common complication of diabetes mellitus, affecting up to 30% of individuals with type 1 diabetes mellitus and 40% of individuals with type 2 diabetes mellitus . Over the last two decades, no new drugs have been approved to specifically prevent DN or to improve kidney function . Current strategies to treat this disease mainly focus on intensification control of glycemic, blood pressure and other risk factors . But a significant proportion of patients still progress to end stage renal failure overtime with a well control of risk factors, which highlight the urgent need to identify effective therapies.
Recently, DN is increasingly considered as an inflammatory disease, with immune modulation being involved in every stage. Immune cells infiltrated in kidney tissues, such as macrophages, T cells, and mast cells produce various inflammatory cytokines, metalloproteinases, and growth factors, which modulate the local response and increase inflammation [4,5,6,7,8].
Although there is substantial evidence for the vital roles of macrophages, T cells, B cells, and mast cells, the effect of other immune cells including neutrophils, NK cells, and dendritic cells are less clear. Furthermore, the immune system contains various types of immune cell, which constitute a complex regulatory network, potential intercorrelations of them in the pathological process of DN are less clear. Therefore, an integrated analysis of immune infiltration, differentially expressed genes, and aberrant pathways activation in kidney tissues is essential.
Here, we used the CIBERSORT algorithm to analyze the abundance of 22 kinds of immune cells in the glomerulus samples of DN based on three microarray datasets mined from the Gene Expression Omnibus (GEO) . Then an integrated bioinformatics analysis was conducted to screen immune-related differentially expressed genes(IR-DEGs)and explore aberrant pathways activation. Besides, we evaluated the correlation between IR-DEGs with immune infiltration cells and clinical indicator. Our research provides a potential value for immunotherapy of DN.
Materials and methods
We systematically searched publicly available DN gene expression datasets on Gene Expression Omnibus (GEO) database using the following terms: (“DN” or “DKD” or “diabetic nephropathy” or “diabetic kidney disease”) and “Homo sapiens” [porgn:_txid9606] and (“Expression profiling by array”[Filter]).
Datasets included in this study have to meet the following criteria: (1) The organism of samples is Homo sapiens. (2) The sample type is glomerular tissue. (3) gene expression was detected by microarray. (4) dataset including DN and NC (negative control) samples (5) The dataset contains more than 3 samples per group. Finally, three datasets were included in this research (GSE96804, GSE30528, and GSE47183). Database searches were conducted in February 2022. The detailed workflow of this study is shown in Fig. 1.
The probe ID was converted to a gene symbol using the platform annotation files. If a gene matched multiple probes, the maximum value of multiple probes expression was applied. To avoid sample heterogeneity arising from different technical platforms and dissimilar methods of processing data, batch effect removal and normalization were conducted by the “limma” package [10, 11]. Then Principal component analysis (PCA) was performed to Confirm homogeneity after data pre-processing.
Immune cell infiltration profiling
A deconvolution algorithm—CIBERSORT—was used to estimate the relative proportion of 22 types of infiltrating immune cell in glomerular tissue by identifying 547 immune cell-related genes .
The normalized gene expression data were uploaded to the CIBERSORT web portal (http://CIBERSORT.stanford.edu/) with a parameter of 1000 permutation and a relative mode. P < 0.05 was considered to be statistically significant. The Pearson correlation coefficient between each immune cell was calculated and visualized with a correlation heatmap by the “corrplot” package.
Differentially expressed gene analysis
The “limma” package was used to reveal the differentially expressed genes (DEGs) between DN and NC samples with a threshold of |log2FC (fold change) |> 1 and adjusted P-value < 0.05 . The visualization of DEGs was conducted by the “ComplexHeatmap” and “ggplot2” packages.
Gene set enrichment analysis (GSEA)
The enrichment analyses for Gene Ontology (GO) terms and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways were carried out by the gene set enrichment analysis (GSEA) method with the “clusterProfiler” package.
Protein–protein interaction (PPI) analysis of immune-related DEGs
The immune-related gene list was downloaded from the ImmPort Database (https://www.immport.org/shared/home). Then immune-related differentially expressed genes (Immune-related DEGs) were identified by intersecting DEGs with immune-related genes. The protein–protein interactions of Immune-related DEGs were estimated by the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org). Then genes with a confidence score of interaction more than 0.95 were visualized by Cytoscape (version 3.8.2). The “Cytohubba” plug-in was used to identify the highest linkage hub genes in the network and the parameters used were as follows: the top 10 nodes, ranked by degree. .
Correlation analysis of Hub genes with infiltrated immune cells and clinical indicator
Correlation analysis between hub genes and immune cells were calculated and visualized with the “corrplot” package. Then the hub genes which has a medium to strong correlation with immune cells were conducted correlation analysis with renal function indicator—estimated Glomerular Filtration Rate (eGFR), and visualized with the “ggpubr” package.The eGFR and corresponding gene expression matrix used here were mined from the Nephroseq database (includes transcriptome data of 9 DN patients and 13 NC volunteers glomerular tissue).
All statistical analysis was performed on R software (version 4.0.2). An independent t-test was performed to compare continuous variables between two groups, p-value < 0.05 was considered statistically significant. Correlations were assessed by the Pearson correlation coefficient (0.1 <| r |< 0.4 represents a weak correlation; 0.4 <| r |< 0.7 represents a medium correlation; and | r |> 0.7 represents a strong correlation).
Data integration and removal of batch effects
First, the GEO series of GSE96804 (41 DN and 20 Negative Controls—NC samples), GSE30528 (9 DN and 13 NC samples), and GSE47183 (14 DN and 14 NC samples) were acquired from the GEO database and merged as one big dataset with a total of 111 samples (64 DN and 47 NC samples). The second series (GSE30528) employs excellent Negative Controls (donor biopsies), while the NC samples of the other two come from the unaffected portion of tumor nephrectomies.
Then, the batch effect was evaluated and visualized by a box plot and a PCA score plot (Fig. 2A-B). Then batch effect removement and data normalization were conducted with the “limma” package. At last, the homogeneity between GEO series was confirmed (Fig. 2C-D).
Infiltration proportion of immune cells changed in DN group
The proportion of 22 types of immune cells in the glomerulus tissues were estimated by the CIBERSORT algorithm. As shown in Fig. 3A, the fraction of immune cells varied significantly among the samples and groups. Due to the limitations of the CIBERSORT algorithm, the infiltration of some low abundance immune cell subsets wasn’t fully revealed. The violin plot shows that 11 types of cells changed significantly in DN. Among which, 7 types were up-regulated and 4 types were down-regulated (Fig. 3B). These results demonstrated that there were different infiltration patterns of immune cells in DN glomerular tissues.
To evaluate whether different types of immune cells interact with each other, correlation analysis was conducted. And we found multiple pairs of positively or negatively related immune cells. Among them, T cells follicular helper and B cells memory show the most synergistic effect (r = 0.51, p < 0.001); meanwhile, Mast cells activated and Mast cells resting show the most competitive effect (r =—0.56, p < 0.001) (Fig. 4).
Identification of differentially expressed genes(DEGs)
In differential expression analysis (64 DN vs. 47 NC), a total of 143 differentially expressed genes were identified, of which 58 were upregulated and 85 were downregulated in DN. To visualize the DEGs more intuitively, we plotted a volcano plot and a heatmap plot (Fig. 5A-B). The expression details of the top 10 upregulated and top 10 downregulated genes are listed in Table 1.
Gene set enrichment analysis (GSEA)
To explore the possible pathways and gene sets associated with immune functions, we performed GSEA. Results show that several pathways related to immune, autophagy, diabetes, and Kidney disease were significantly activated in the DN group, including Leukocyte transendothelial migration, Chemokine signaling pathway, FoxO and PI3K − Akt signaling pathway etc. (Fig. 6A-B).
To explore the immune-related gene function in the pathological of DN, GO biological processes (GO-BP) analysis was conducted using GSEA. Eight representative immune-related and eight metabolic-related biological processes were shown in Fig. 6C-D. These results show that immune and metabolic-related processes were activated aberrantly in the glomerulus of DN.
PPI network construction of immune-related DEGs
Thirty-three immune-related differentially expressed genes were obtained by intersecting DEGs with the lists of immune-related genes (IRGs) obtained from the ImmPort database (Fig. 7A). To get a further understanding of the regulation effects among these genes, a PPI (Protein–protein interaction) was constructed by the “Cytoscape” (Fig. 7B), and ten hub genes were identified by the “Cytohubba” plug-in (Fig. 7C).
Correlation analysis between hub genes and infiltrated immune cells
To understand the molecular mechanisms of alteration in infiltrating immune cells in detail, the correlations between immune related hub genes and immune cells need to be unveiled. So, we calculated Pearson's correlation coefficient for hub genes with immune cells. Figure 8 demonstrates that 9 of the 10 hub genes have moderate to high correlations with differential infiltrated cells(r > 0.4), Especially FOS with neutrophils (r = 0.7). However, none of the hub genes shows a moderate or higher association with non-differentially infiltrating immune cells. This suggests that these genes may be associated with the retention of immune cells in the glomerulus.
Exploration of the association between hub genes and clinical indicator
Furthermore, we explored the correlations between the 9 hub genes and the vital clinical renal function indicator—eGFR—based on data mining from the Nephroseq database (version: v5). The data show that 3 genes are negatively correlated with eGFR levels (SERPINA3, LTF, C3), 3 genes are positively correlated with it (PTGDS, EGF, ALB), but the other 3 did not meet the statistical threshold (Fig. 9).
With the advancement of bioinformatics, efficient algorithms have accelerated the transition of various omics big data to new therapeutic targets. As one of the most robust algorithms for analyzing immune cell infiltration—CIBERSORT—has been used in multiple tissues [13,14,15,16]. In this study, we used it to explore the landscape of immune infiltration in the glomerulus of DN.
In summary, we found that a variety types of infiltrated immune cell have altered during the pathogenesis of DN. Correlation analysis shows that the immune response may act as a complicated and tightly regulated network. In addition, we also screened out 6 immune-related hub genes which not only have correlations with differently infiltrated immune cells but also have a medium to strong correlation with eGFR (the most important renal function indicators).
Serpin Family A Member 3 (SERPINA3) was revealed to associated with acute immune responses in human renal allograft, and could be a novel biomarker for timely diagnosis [17, 18]. And Andrea Snchez-Navarro found urinary SerpinA3 could detect renal fibrosis and inflammation, with a particular potential for the early detection of AKI to CKD transition . Lactotransferrin (LTF) was reported could be delivered to the Intracerebral hemorrhage affected brain by infiltrating polymorphonuclear neutrophils, and may assist in hematoma detoxification . Qi Zhao found LTF could regulate the immune microenvironment of prostate cancer through JAK/STAT3 Pathway . Nick Wlazlo revealed that Complement C3 (C3) were longitudinally associated with insulin resistance in type 2 diabetes over a 7-year follow-up period, and may reflect progression of metabolic dysregulation . A diabetes control and complications Trial found PTGDS (prostaglandin-H2 D-isomerase), strongly and positively associated with the albumin excretion rate in Type 1 Diabetes Mellitus . Boris B. Betz performed urinary peptidomics analysis in a rodent DN model (Cyp1a1mRen2, renin-dependent hypertension and diabetes mellitus synergized to produce massive albuminuria), and found that urinary epidermal growth factor (EGF) was more than twofold reduced in rats with DN in comparison with controls . Hypochlorite modified albumin (ALB) is an active compound that is formed during the reaction between proteins, and was found to be much higher in concentration and to act as a mediator of oxidative stress and inflammation in patients with uremia or diabetes mellitus . These results show that the 6 hub genes might act as vital roles in the pathological of DN.
It is of interest that neutrophils decreased drastically in DN patients, which is in contrast to the study of an STZ (streptozotocin) induced diabetic mouse model . There are several possible reasons for this opposite conclusion: (1) the impact of intervention drugs may impair neutrophils activity. SCIENCE TRANSLATIONAL MEDICINE reported that angiotensin-converting enzyme inhibitors (ACEI) intervention significantly reduced neutrophils activity compared to treated with an (angiotensin receptor blocker) ARB or with no drug in vitro and in vivo. What's more, ACE knockout neutrophils showed decreased survival signaling and increased apoptosis . (2) the negative feedback regulation of the immune microenvironment; (3) the sample size of this study is still not enough, resulting in abnormal results.
There are some related bioinformatics researches based on DN transcriptome data, but these literatures mainly focus on screening DEGs and exploring related pathways [28,29,30,31]. To the best of our knowledge, there is no research report about using transcriptome data to reveal the landscape of immune cell infiltration in DN renal tissue. In addition, the hub genes we reported here were not only highly correlated with renal function, but more importantly, were highly correlated with differentially infiltrated immune cells. In short, we focused on immune infiltration and screened out 6 hub genes may associate with DN immune infiltration.
Nevertheless, some limitations still existed in our study: first, the infiltration of immune cells was estimated by algorithm, rather than based on more direct evidence, such as cell markers. Second, although batch effect removal and normalization were done, the homogeneity of samples between different datasets is still not fully guaranteed. For example, the DN samples may come from patients with very different pathological stages, or with other co-morbidities. And some NC samples come from tumor nephrectomies, which could have changes related to inflammation or alterations in vasculature. So, the results of immune cells differential infiltration stated here is far from definitive. Third, the correlation analysis of hub genes with eGFR is done based on a relatively small sample data (only 22 samples), so this may lead to over- or underestimation of its importance. And the exact mechanism of hub genes involved in immune infiltration had not been investigated through biological experiment in the present study. So, there is an urgent need to conduct in vivo experiments to confirm the results with higher-level evidence.
In conclusion, our work provides a global picture of the immune environment in the glomerulus of DN. The potential therapeutic targets identified in the present study may provide new insightful for clinical treatment.
Availability of data and materials
All raw data are available in GEO datasets (GSE96804, GSE30528, and GSE1009).
Differentially expressed genes
Gene Expression Omnibus
Kyoto Encyclopedia of Genes and Genomes
Principal component analysis
Gene set enrichment analysis
Search Tool for the Retrieval of Interacting Genes
Estimated glomerular filtration rate
Reutens AT. Epidemiology of diabetic kidney disease. Med Clin North Am. 2013;97(1):1–18.
Breyer MD, Susztak K. The next generation of therapeutics for chronic kidney disease. Nat Rev Drug Discov. 2016;15(8):568–88.
Fineberg D, Jandeleit-Dahm KA, Cooper ME. Diabetic nephropathy: diagnosis and treatment. Nat Rev Endocrinol. 2013;9(12):713–23.
Guiteras R, Sola A, Flaquer M, Manonelles A, Hotter G, Cruzado JM. Exploring macrophage cell therapy on Diabetic Kidney Disease. J Cell Mol Med. 2019;23(2):841–51.
Bending JJ, Lobo-Yeo A, Vergani D, Viberti GC. Proteinuria and activated T-lymphocytes in diabetic nephropathy. Diabetes. 1988;37(5):507–11.
Xiao X, Ma B, Dong B, Zhao P, Tai N, Chen L, Wong FS, Wen L. Cellular and humoral immune responses in the early stages of diabetic nephropathy in NOD mice. J Autoimmun. 2009;32(2):85–93.
Okon K, Stachura J. Increased mast cell density in renal interstitium is correlated with relative interstitial volume, serum creatinine and urea especially in diabetic nephropathy but also in primary glomerulonephritis. Pol J Pathol. 2007;58(3):193–7.
Galkina E, Ley K. Leukocyte recruitment and vascular injury in diabetic nephropathy. J Am Soc Nephrol. 2006;17(2):368–77.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.
Benito M, Parker J, Du Q, Wu J, Xiang D, Perou CM, Marron JS. Adjustment of systematic microarray data biases. Bioinformatics. 2004;20(1):105–14.
Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. https://doi.org/10.1093/nar/gkv007.
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4):S11.
Yuan WH, Xie QQ, Wang KP, Shen W, Feng XF, Liu Z, Shi JT, Zhang XB, Zhang K, Deng YJ, et al. Screening of osteoarthritis diagnostic markers based on immune-related genes and immune infiltration. Sci Rep. 2021;11(1):7032.
Tan L, Xu Q, Shi R, Zhang G. Bioinformatics analysis reveals the landscape of immune cell infiltration and immune-related pathways participating in the progression of carotid atherosclerotic plaques. Artif Cells Nanomed Biotechnol. 2021;49(1):96–107.
Wang L, Wei Q, Zhang M, Chen L, Li Z, Zhou C, He M, Wei M, Zhao L. Identification of the prognostic value of immune gene signature and infiltrating immune cells for esophageal cancer patients. Int Immunopharmacol. 2020;87: 106795.
Wang D, Xiong T, Yu W, Liu B, Wang J, Xiao K, She Q. Predicting the Key genes involved in aortic valve calcification through integrated bioinformatics analysis. Front Genet. 2021;12:650213.
Wang X, Zu Q, Lu J, Zhang L, Zhu Q, Sun X, Dong J. Effects of donor-recipient age difference in renal transplantation, an Investigation on renal function and fluid proteome. Clin Interv Aging. 2021;16:1457–70.
Yang B, Sylvius N, Luo J, Yang C, Da Z, Crotty C, Nicholson ML. Identifying biomarkers from transcriptomic signatures in renal allograft biopsies using deceased and living donors. Front Immunol. 2021;12: 657860.
Sanchez-Navarro A, Mejia-Vilet JM, Perez-Villalva R, Carrillo-Perez DL, Marquina-Castillo B, Gamba G, Bobadilla NA. SerpinA3 in the Early Recognition of Acute Kidney Injury to Chronic Kidney Disease (CKD) transition in the rat and its Potentiality in the Recognition of Patients with CKD. Sci Rep. 2019;9(1):10350.
Zhao X, Ting SM, Sun G, Roy-O’Reilly M, Mobley AS, Bautista Garrido J, Zheng X, Obertas L, Jung JE, Kruzel M, et al. Beneficial role of neutrophils through function of lactoferrin after intracerebral hemorrhage. Stroke. 2018;49(5):1241–7.
Zhao Q, Cheng Y, Xiong Y. LTF Regulates the immune microenvironment of prostate cancer through JAK/STAT3 pathway. Front Oncol. 2021;11:692117.
Wlazlo N, van Greevenbroek MM, Ferreira I, Feskens EJ, van der Kallen CJ, Schalkwijk CG, Bravenboer B, Stehouwer CD. Complement factor 3 is associated with insulin resistance and with incident type 2 diabetes over a 7-year follow-up period: the CODAM Study. Diabetes Care. 2014;37(7):1900–9.
Shao B, Zelnick LR, Wimberger J, Himmelfarb J, Brunzell J, Davidson WS, Snell-Bergeon JK, Bornfeldt KE, de Boer IH, Heinecke JW. Albuminuria, the high-density lipoprotein proteome, and coronary artery calcification in type 1 diabetes mellitus. Arterioscler Thromb Vasc Biol. 2019;39(7):1483–91.
Betz BB, Jenks SJ, Cronshaw AD, Lamont DJ, Cairns C, Manning JR, Goddard J, Webb DJ, Mullins JJ, Hughes J, et al. Urinary peptidomics in a rodent model of diabetic nephropathy highlights epidermal growth factor as a biomarker for renal deterioration in patients with type 2 diabetes. Kidney Int. 2016;89(5):1125–35.
Tang DD, Niu HX, Peng FF, Long HB, Liu ZR, Zhao H, Chen YH. Hypochlorite-modified albumin upregulates ICAM-1 expression via a MAPK-NF-kappaB signaling cascade: protective effects of apocynin. Oxid Med Cell Longev. 2016;2016:1852340.
Mohamed R, Jayakumar C, Ranganathan PV, Ganapathy V, Ramesh G. Kidney proximal tubular epithelial-specific overexpression of netrin-1 suppresses inflammation and albuminuria through suppression of COX-2-mediated PGE2 production in streptozotocin-induced diabetic mice. Am J Pathol. 2012;181(6):1991–2002.
Cao DY, Giani JF, Veiras LC, Bernstein EA, Okwan-Duodu D, Ahmed F, Bresee C, Tourtellotte WG, Karumanchi SA, Bernstein KE, et al. An ACE inhibitor reduces bactericidal activity of human neutrophils in vitro and impairs mouse neutrophil activity in vivo. Sci Transl Med. 2021;13(604):eabj2138.
Tang YL, Dong XY, Zeng ZG, Feng Z. Gene expression-based analysis identified NTNG1 and HGF as biomarkers for diabetic kidney disease. Medicine (Baltimore). 2020;99(1):e18596.
Zeng M, Liu J, Yang W, et al. Identification of key biomarkers in diabetic nephropathy via bioinformatic analysis [published online ahead of print, 2018 Nov 28]. J Cell Biochem. 2018. https://doi.org/10.1002/jcb.28155.
Geng XD, Wang WW, Feng Z, Liu R, Cheng XL, Shen WJ, Dong ZY, Cai GY, Chen XM, Hong Q, et al. Identification of key genes and pathways in diabetic nephropathy by bioinformatics analysis. J Diabetes Investig. 2019;10(4):972–84.
Xu B, Wang L, Zhan H, Zhao L, Wang Y, Shen M, Xu K, Li L, Luo X, Zhou S, et al. Investigation of the mechanism of complement system in diabetic nephropathy via bioinformatics analysis. J Diabetes Res. 2021;2021:5546199.
This work was financially supported by Shiyuan Jin Academic Experience Inheritance Studio Project, Guangdong Provincial Hospital of Traditional Chinese Medicine [2018 No. 7] and Xiaotao Wang Academic Experience Inheritance Studio Project, Guangdong Provincial Hospital of Traditional Chinese Medicine [2018 No. 7].
Ethics approval and consent to participate
All methods were performed in accordance with the relevant guidelines and regulations.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
ZHOU, W., LIU, Y., Hu, Q. et al. The landscape of immune cell infiltration in the glomerulus of diabetic nephropathy: evidence based on bioinformatics. BMC Nephrol 23, 303 (2022). https://doi.org/10.1186/s12882-022-02906-4