Activation of complement C1q and C3 in glomeruli might accelerate the progression of diabetic nephropathy: Evidence from transcriptomic data and renal histopathology

ABSTRACT Aims/Introduction It is not unclear whether the complement system is involved in the pathogenesis of diabetic nephropathy (DN). We explored the role of the complement system in glomeruli from patients with DN using integrated transcriptomic bioinformatics analysis and renal histopathology. Materials and Methods Four datasets (GSE30528, GSE104948, GSE96804 and GSE99339) from the Gene Expression Omnibus database were integrated. We used a protein–protein interaction network and the Molecular Complex Detection App to obtain hub genes. Gene ontology and the Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were carried out to identify significant pathways. We also investigated the associations of C1q and C3 deposition on renal histopathology with clinical data, pathological parameters and renal survival in DN patients. Results We identified 47 up‐ and 48 downregulated genes associated with DN. C3, C1QB and C1QA were found to be complement‐related hub genes. The gene ontology and Kyoto Encyclopedia of Genes and Genomes analyses identified complement activation and humoral immune response as the significant oncology terms, with C1QB and C3 positioned at the center of the pathway. Regarding renal histopathology, patients with both C1q and C3 deposition had more severe glomerular classes. Multivariate Cox proportional hazards regression showed that the deposition of glomerular C1q and C3 was an independent risk factor for kidney failure. Patients with high C1q, C3 or C4d expression in glomeruli were more likely to progress to kidney failure, whereas glomerular mannose‐binding lectin was rare. Conclusions Complement activation is involved in the development of DN, and activation of the classical complement pathway in glomeruli might accelerate disease progression.


INTRODUCTION
Diabetic nephropathy (DN) is the leading cause of end-stage renal disease (ESRD) worldwide 1 . The estimated global prevalence of diabetes in 2019 was 9.3% (463 million), and is projected to increase to 10.9% (700 million) by 2045 2 . Nearly 30-40% of diabetes patients develop DN, and most patients with DN die of cardiovascular complications or infections before they progress to ESRD. Although various drugs are available, the existing drugs have limited therapeutic effect on DN. Therefore, it is necessary to explore its etiopathogenesis 3 . Pathophysiological studies have shown that the progression of DN is closely related to metabolic abnormalities, oxidative stress, hemodynamic changes and innate immunity 4,5 . As an important part of the innate immune system, the complement system has attracted much attention for its role in the development of DN [6][7][8][9] .
The complement system is the first barrier to accelerated inflammation, and it ameliorates the clearance of impaired cells and pathogenic microorganisms from the body through phagocytes and antibodies 10 . However, it remains unclear whether the complement system is involved in the pathogenesis of DN.
Omics analyses, especially transcriptomics and metabolomics, have been widely used to detect genes associated with the pathogenesis of DN 11 . Transcriptomic profiling using microarray chips provides important information about the molecular and cellular mechanisms of disease pathogenesis and progression 12 . However, independent microarray analyses might lead to a high false positive rate. The R software package, Robus-tRankAggreg (The R Foundation For Statistical Computing, Vienna, Austria), is helpful for ranking differentially expressed genes (DEGs) obtained from different platforms using a robust rank aggregation (RRA) approach 13 . Therefore, we integrated four transcriptomic microarray datasets from kidney tissue using RRA to determine whether complement-related genes are hub genes involved in DN. Then, we investigated the clinical and pathological significance of complement deposition in the renal histopathology of patients with DN to explore the potential role of the complement system in DN, which should provide new insights into potential targets for DN therapeutics.

MATERIALS AND METHODS
Collection of gene expression datasets and platform records All datasets, derived from glomerular transcriptional profiles from DN patients, were downloaded from the Gene Expression Omnibus Database (https://www.ncbi.nlm.nih.gov/geo/). The gene expression profiling microarray data were eligible if they provided transcriptome comparison of glomeruli from kidneys with DN and control glomeruli tissues obtained from healthy living donors, and had more than three samples in each group. Finally, four datasets (GSE30528 14 , GSE104948 15 , GSE96804 16 and GSE99339 17 ) were included in the present study (Table 1).

Datasets pre-processing and identification of DEGs
We first used the "limma" package, a Bioconductor package of R (version 4.0.5), to normalize and select DEGs from the messenger ribonucleic acid expression microarray data. Then, we used the "RobustRankAggreg" package in R to integrate these four datasets with the RRA algorithm. Based on the criteria (|logFC (fold change)| >1 and adjusted P-value <0.05), the DEG frame results were finally established.

Functional enrichment analyses of DEGs
To explore the potential functions of the DEGs, we carried out gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the "clusterProfiler" package in R. The adjusted P-value <0.05 was considered significant. Using the "ggplot2" package in R, we visualized all data frame results.

Identification of hub genes
We used protein-protein interaction (PPI) 18 analysis to further explore the hub genes that were embedded in gene regulatory networks. First, the Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) database was used to search for the PPI network of DEGs with the highest confidence-score (>0.9). Then, the PPI network data, which were exported from STRING, were imported into Cytoscape 3.7.2 (Institute for Systems Biology, Seattle, WA, USA) for their visualization. Subsequently, we used the Molecular Complex Detection (MCODE), a clustering algorithm of Cytoscape, to carry out network analysis to obtain hub genes (the parameter criteria were as follows: degree cut-off = 2, k-score = 2, node score cut-off = 0.2 and max depth = 100). Consequently, the GO and KEGG pathway analysis of hub genes were carried out with the R package.

Clinical validation of DN patients
We retrospectively collected patients with type 1 or type 2 diabetes and nephropathy who underwent renal biopsy at the China-Japan Friendship Hospital, Beijing, China, between January 2013 and December 2020. All participants met the diagnostic criteria for diabetes proposed by the American Diabetes Association in 2021 19 and the criteria of DN followed by the Renal Pathology Society in 2010 20 . We recorded the classic histopathological characteristics, such as expansion of the mesangial matrix, glomerular hypertrophy, thickening of the glomerular basement membrane and nodular glomerulosclerosis. The exclusion criteria were as follows: (i) coincident nondiabetic renal disease, estimated glomerular filtration rate (eGFR) <15 mL/min/1.73 m 2 at baseline, inadequate renal tissue (<5 glomeruli) and follow-up interval <6 months or no follow-up information.

Clinical characteristics of DN patients
The baseline data of included patients were collected from the Electronic Medical Record System, including age, sex, duration of diabetes, 24-h urinary protein excretion, serum levels of albumin, creatinine, complement C1q (C1q), complement 3 (C3) and complement 4 (C4). The eGFR was calculated using the creatinine-based Chronic Kidney Disease Epidemiology Collaboration equation 21 . The follow-up interval was assessed from the time of renal biopsy to one of four end-points: ESRD, death, loss to follow up or study termination (October 2021). ESRD was defined as a need for maintenance dialysis or kidney transplantation. The primary end-point in the present study was the composite of a doubling of the baseline serum creatinine concentration or ESRD.

Renal histopathology
Renal biopsy specimens were evaluated using light microscopy, immunofluorescence staining and electron microscopy. Classification of DN and histological scoring were carried out according to the criteria of Tervaert et al 20 . The glomerular classification was as follows. Class I, isolated glomerular basement membrane thickening, and only mild, non-specific light microscopy changes that do not meet the criteria of class II through IV. Class II, mesangial expansion, mild (IIa) or severe (IIb): glomeruli classified as mild or severe mesangial expansion, but not meeting the criteria for class III or IV. Class III, nodular sclerosis (Kimmelstiel-Wilson lesions): at least one glomerulus with Kimmelstiel-Wilson nodules and <50% global glomerulosclerosis. Class IV, advanced diabetic glomerulosclerosis: >50% global glomerulosclerosis. The interstitial fibrosis and tubular atrophy (IFTA) scores were classified as follows: 0, absent; 1, <25%; 2, 25-50%; and 3, >50% of the total area. The interstitial inflammation was scored as follows: 0, absent; 1, infiltration only in areas related to IFTA; and 2, infiltration in areas without IFTA. Scores of vascular lesions were evaluated based on the presence of large-vessel arteriosclerosis and arteriolar hyalinosis. The intensity of immunofluorescence staining of the complement (C1q and C3) deposits in glomeruli was graded using a semiquantitative method on a scale of 0+ +++. When the presence of glomerular complement (C1q or C3) staining was detected at 1+ or a higher grade, it was taken as positive.
Immunohistochemical images were obtained from sections observed microscopically at 9200 magnification using an integrated digital camera system (Nikon, Tokyo, Japan). An average optical density was obtained by analysis of 10 different fields, and quantification was carried out semiquantitatively in a blinded manner. The Image-Pro Plus software 6.0 (Media Cybernetics, Bethesda, MD, USA) was used to analyze the average optical density of the stained areas.

Statistical analysis
We used SPSS 21.0 (SPSS Inc.; www.spss.com) to analyze clinical and pathological data; graphs were created using GraphPad Prism 8.0 (GraphPad Software Inc.; www.graphpad.com). Continuous variables were represented as meanstandard deviation. Categorical variables were summarized as numbers (percentages). Student's t-test or one-way ANOVA was used for continuous variables. The v 2 -test or Fisher's exact test was used for categorical variables. We used the linear regression analysis to calculate the slopes of the rates of change in eGFR. The results were also analyzed by calculating the difference between the baseline and the final eGFR and dividing by time. Cox proportional hazard models were carried out to assess the hazard ratio (HR) and 95% confidence interval (CI) for parameters associated with the definitive renal outcome. We first used univariate log-rank tests to estimate the association of baseline parameters with the definitive renal outcome. The variables that had P-values <0.1 by univariate analysis or were clinically significant were subsequently entered into a final multivariate Cox regression model. Differences with a two-sided P-value <0.05 was considered statistically significant.

Construction of DEGs by integrated analysis
The flowchart of the present study is shown in Figure 1. By integrating four microarray datasets (GSE30528, GSE104948, GSE96804 and GSE99339), we identified 95 DEGs, including 47 upregulated genes and 48 downregulated genes, which were the most significant genes that were associated with DN (Table S1). The top 20 significantly upregulated and downregulated genes ranked according to the fold change value were shown in Figure 2a. GO analysis showed that upregulated genes were mainly enriched in humoral immune response (9.69E-10), extracellular matrix organization (1.58E-9), extracellular structure organization (1.62E-9) and regulation of complement activation (9.74E-9), whereas downregulated genes were mainly enriched in striated muscle tissue development (3.25E-8), muscle organ development (4.97E-8), muscle tissue development (5.20E-8) and carbohydrate biosynthetic process (9.08E-7; Figures S1 and S2). We then carried out KEGG pathway enrichment analysis of upregulated genes, and found that the top three enriched pathways were complement and coagulation cascades (3.47E-07), pertussis (5.25E-06) and extracellular matrix-receptor interaction (1.08E-05; Figure 2b). In contrast, the top three pathways that were related to downregulated genes, were mainly focused in HIF-1 signaling pathway (0.000567), growth hormone synthesis, secretion and action (0.000788), and FoxO signaling pathway (0.00113; Figure 2c).

PPI analysis and identification of hub genes
The PPI network of the aforementioned 95 DEGs was built according to the STRING database and was then ranked with topological analysis on nodes. The exported data was visualized by Cytoscape software ( Figure S3). Four discrete subclusters were obtained through MCODE plug-in by Cytoscape software, and the top 24 hub genes ranked based on the degree score was shown in Table S2.
GO and KEGG enrichment analysis of hub genes GO analysis of hub genes were found to be associated with 96 biological process terms. Figure 3a shows the top 10 terms of GO analysis. We found that the hub genes were mainly enriched in extracellular organization (2.76E-10), complement activation (1.87E-07) and humoral immune response (4.01E-07). To further explore the pathogenesis of DN, we again carried out the KEGG pathway analysis. Figure 3b presents the results of the top five terms of the KEGG enrichment analyses, and most of the hub genes were linked to complement and coagulation cascades (6.29E-07), extracellular matrix-receptor interaction (7.49E-07), focal adhesion (2.12E-06) and PI3K-protein kinase B signaling pathway (3.79E-06). Based the above GO and KEGG analyses results, we could find that some hub genes (C1QA, C1QB and C3), which were involved in the complement activation process, did play a vital role in the pathogenesis of DN.

Comparison of the clinical and pathological characteristics of DN patients
This study enrolled 151 patients with pathologically diagnosed DN ( Figure S4). C3 or C1q deposits were found in glomerular capillary walls, the mesangium, Bowman's capsule and the tubular basement membrane. The grouping of C3 or C1q deposits in the present study was based on the presence of deposition in the glomeruli. On direct immunofluorescence   microscopy, C1q and C3 co-deposits (C1q+ C3+), only C1q or C3 deposits (C1q+ C3-/C1q-C3+), and no complement deposits (C1q-C3-) were found in 24 of 151 (15.89%), 55 of 151 (36.42%) and 72 of 151 (47.68%) patients, respectively. Patients with C1q+ C3+ deposition in glomeruli had significantly lower serum albumin (P < 0.01) and higher urinary protein excretion (P < 0.05) levels compared with those with only one or no complement deposition in glomeruli. Furthermore, patients with C1q+ C3+ deposition in glomeruli had significantly more severe glomerular classes (P < 0.05; Table 2).

Survival analysis according to C1q and C3 deposition in DN patients
Of the 151 DN patients, ESRD or doubling of the baseline serum creatinine level was noted in 90 cases. Patients who died before they reached the composite renal outcome were excluded from the present study. The mean follow up of this cohort was 23 -19 months. As shown in Figure 4a, patients with C1q+ C3+ deposition in glomeruli had worse renal survival compared with those with C1q-C3-or with C1q+ C3-/C1q-C3+ (P < 0.05). The event-free survival probability was significantly lower in DN patients with than in those without C3 deposition in glomeruli (P < 0.05 Figure 4b), as well as in those with than in those without C1q deposition in glomeruli (P < 0.001, Figure 4c). These results showed that glomerular C1q and C3 deposition were associated with poorer renal survival.
Immunohistochemical analyses of C3, C1q, C4d and MBL in kidney tissues C3 deposition in the kidney indicates activation of the common complement pathway, so further research is required on the complement pathways involved in DN. We carried out a preliminarily immunohistochemical analysis of complement components in kidney tissues from 54 DN patients who were followed for 1 and 2 years. Of the 54 patients, seven were followed up for only 1 year. Figure 6 shows immunohistochemical staining for C1q, C3, C4d and MBL in glomeruli (Figure 6a-d), and the relationship between the expression of Transcriptional levels of C1QA and C3 in different types of kidney diseases We then further explored the gene expression levels of complement components in DN when compared with other kidney diseases. The transcriptional data were downloaded from Nephroseq v5 online platform (https://www.nephroseq.org/). The levels of gene expression of C1QA and C3 in DN glomeruli were comparable to those in glomeruli from immunoglobulin A nephropathy or lupus nephritis, but even significantly higher than those from membranous nephropathy and minimal change disease (Figure 7).

DISCUSSION
In the present study, we identified DEGs associated with DN from four independent datasets using RRA, which overcomes noise in different individual studies. A total of 95 overlapping DEGs, including 47 up-and 48 downregulated genes, were identified. The top 24 hub genes ranked based on the degree score by MCODE were selected for GO and KEGG enrichment analyses to investigate their molecular mechanisms. In the GO analysis, most hub genes were enriched in extracellular matrix organization (2.76E-10), complement activation (1.87E-07) and humoral immune response (4.01E-07). The C1QA, C1QB, VSIG4, C3 and CLU hub genes were found to be involved in both complement activation and humoral immune response, which is consistent with previous studies [22][23][24] . In KEGG pathway enrichment analysis, we found that most of the top 10 hub genes were linked to the complement and coagulation cascades, supporting the potential role of the complement system in the development of DN. Five of the top 24 hub genes were involved in complement system function. The hub gene, CLU (the abbreviation of clusterin), a kind of complement regulatory protein, combines with C5b-7 and hinders the production of the membrane attack complex 25 . Clusterin was upregulated in both the urine and glomeruli of DN patients 23,26 , and it could attenuate the development of renal fibrosis 27 . The B7 family-related protein, VSIG4, acts as a complement receptor for C3 28 , but it remains unclear how VSIG4 interacts with the complement system to affect the course of DN. C3 not only plays a pivotal role in the classical/lectin complement pathway, but also in the alternative complement pathway as the central component in the complement cascade 29 .
Next, we investigated the renal histopathology of 151 patients with DN. C1q and C3 staining was routinely applied by direct immunofluorescence. We compared the clinical and histological characteristics according to whether there was C3 and C1q deposition in the glomeruli of patients with DN, and found that patients with both glomerular C3 and C1q deposition had lower serum albumin levels, higher 24-h urinary protein excretion, and more severe renal pathologic lesions compared with those without glomerular C3 and C1q deposition. Sun et al. 7 found an association between C1q and C3c complement deposition on renal histopathology and more severe kidney damage in DN patients. These results emphasize the potential role of local complement activation in DN. Several studies have also shown the vital role of complement in the pathogenesis of DN [30][31][32][33] . For instance, Pelletier et al. 31 found that renal C1q was strongly associated with inflammation, fibrosis, proteinuria and the rate of renal function decline in patients with DN, and renal C1q deposition mainly in glomeruli hili and arterioles was correlated with the occurrence of DN 30 . Angeletti et al. 32 suggested that loss of C3 convertase regulator prompted podocyte injury and glomerulosclerosis. BTBR ob/ob mice, which show renal C3 deposition and increased expression of the C3a and C3a receptor, developed glomerular injury, podocyte loss and albuminuria 33 . Consistently, a wide range of inflammatory   processes was confirmed in glomeruli from DN mice by single-cell sequencing 34 . Complement activation is a primary pathogenic mechanism in numerous inflammatory disorders 35 .
Here, we found that the rate of kidney function decline was significantly higher in patients with both C3 and C1q deposition than in patients without either C3 or C1q deposition in the glomeruli, and that glomerular C1q+ C3+ deposition was an independent risk factor for renal survival. This suggests that complement activation plays a role in DN progression. The deposition of C3 in glomeruli indicates activation of the common complement pathway, and it is necessary to determine which complement pathway is activated in DN. Therefore, we examined the immunohistochemical expression of C1q (classical pathway), MBL (lectin pathway), C4d and C3 in glomeruli. C4d is a C4 cleavage product that binds covalently to the target tissue, and can arise from both the classical and MBL pathways. Those with high baseline C1q, C3 or C4d expression in glomeruli were more likely to progress to ESRD, whereas glomerular MBL was rare in the present cohort. These findings suggest that local classical complement pathway in glomeruli is involved in the progression of DN. Therefore, suppression of the complement system might be a novel therapy for DN 36 . Further study of the role of the complement system in DN might help to initiate new therapeutic approaches 37 .
The present study had several limitations. First, how the complement system plays a role in DN and its mechanism requires further research. Second, complement deposition in renal tubules and vessels was not evaluated, but it might also be involved in the development of DN. Further research is required to verify the present results and elucidate the precise mechanism of complement in the development and progression of DN.
In conclusion, local complement system activation in glomeruli is involved in DN, and C1q and C3 deposition in glomeruli might be associated with worse renal survival in patients with DN. A new method of regulating the complement system might lead to the development of new drugs to prevent or slow the progression of DN.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of the article.
Table S1 | The 47 upregulated differential expressed genes and 48 downregulated differentially expressed genes identified through integrating four transcriptomics microarray datasets.