Complement C7 is Specifically Expressed in Mesangial Cells and is a Potential Diagnostic Biomarker for Diabetic Nephropathy and is Regulated by miR-494-3p and miR-574-5p

Introduction

Diabetic nephropathy (DN) is a common microvascular complication that is both the most serious complication of diabetes (DM) and the most common cause of end-stage renal disease.¹ The major pathological features of DN include glomerular vascular injury, glomerular sclerosis, the formation of nodular lesions and the deterioration of renal function, which ultimately leads to end-stage renal disease.^1,2 It has been reported that 30–40% of diabetes cases develop into diabetic nephropathy worldwide.^3,4 Due the complexity of metabolic disorders, once diabetes has developed into end-stage renal disease, it is often more difficult to treat than other kidney diseases; thus, timely prevention and treatment is of great significance to delay diabetic nephropathy.

Weighted gene coexpression network analysis (WGCNA) was used to construct a free-scale gene coexpression network and is widely used to analyze large-scale data sets and identify highly associated gene modules.^5,6 Related modules have been used to identify candidate biomarkers or therapeutic targets in many diseases, such as depressive order⁷ and chronic kidney disease.⁸ Single-cell RNA sequencing (scRNA-seq) has emerged as an indispensable tool to dissect cellular heterogeneity and offers enormous potential for de novo discovery. Single-cell transcriptomic atlases provide unprecedented resolution to reveal complex cellular events and deepen our understanding of biological systems.⁹ Although RNA sequencing (RNA-seq) has been used in several studies of DN, scRNA-seq can pinpoint the cell types in which specific genes are differentially expressed. A recent scRNA-seq study performed on the patient’s kidney cortex revealed strong angiogenic signature changes in glomerular cell types, proximal convoluted tubules, distal convoluted tubules, and principal cells.¹⁰

The complement system is a mainstay of systemic innate immunity comprising several interacting components.¹¹ C3a-mediated proinflammatory and prefibrotic responses in rats with type 2 diabetes mellitus (T2DM) aggravate renal damage,¹² while inhibition of complement C5 could attenuate glomerular mesangial proliferation and urinary protein excretion in rats,¹³ suggesting that complement system activation is associated with the pathogenesis of DN. MicroRNAs (miRNAs) are an abundant class of noncoding RNAs that are reported to be important for many biological processes through the regulation of gene expression, including in DN.¹⁴ Studies have revealed that some miRNAs could be regulators of immune and inflammatory responses and are potential therapeutic targets in DN.^15,16 However, there have been few studies on genetic variation in miRNA and complement levels in diabetic nephropathy.

In the present study, we selected the GSE30528, GSE96804, GSE9006, GSE131882 and GSE161884 datasets from the GEO database to identify a potential biomarker gene for the early diagnosis of DN and to pinpoint the cell types that specifically express the biomarker. First, we performed WGCNA, GO analysis, KEGG analysis and GSEA to identify several candidate genes. Then, qRT-PCR and single-cell RNA sequencing reanalysis were used to confirm and explore the results. The miRNA targets of C7 were predicted by online databases and further validated by qRT-PCR and luciferase reporter gene assays.

Materials and Methods

Collection of Clinical Samples

Serum from total of 31 DN patients and 30 healthy donors was collected at Chu Hsien-I Memorial Hospital (Tianjin, China). All the samples were stored at −80°C until RNA extraction was performed. The research protocol for this study was approved by the Ethics Committee of Chu Hsien-I Memorial Hospital. Informed consent was obtained from all the participants. We adhered to the principles of the 1983 Declaration of Helsinki for all experiments.

We selected 31 patients (16 males and 15 females) with diabetic nephropathy who were hospitalized at the Chu Hsien-I Memorial Hospital (Tianjin, China) from January 2020 to December 2020. The patients had an average disease course of 10.4±5.5 years and an average body mass index (BMI) of 25.6±3.8. Patients with acute cardiovascular and cerebrovascular events, infections, stress states, other endocrine and metabolic diseases, recent ketoacidosis, or hyperosmolar nonketotic coma were excluded. The 2012 K-DOQI standard of the American Kidney Disease Foundation was adopted to diabetic nephropathy disease. Patients with extensive albuminuria or microalbuminuria with diabetic retinopathy or type 1 diabetes with a disease course of more than 10 years with microalbuminuria can be diagnosed with DK. According to UACR staging, there were 14 patients with microalbuminuria (UACR 30~300 mg/g) and 17 patients with massive albuminuria (UACR>300 mg/g). In addition, 30 healthy people who received a physical examination at the center during the same period served as a healthy control (HC) group.

Cell Line

Human mesangial cells were purchased from American Type Culture Collection (ATCC) and cultured in Dulbecco’s modified Eagle’s medium (DMEM) media supplemented with 10% fetal bovine serum (Invitrogen), penicillin (100 U/mL), and streptomycin (100 μg/mL). For the experiments, cells were serum starved for 24 h before supplementing the media with 0,5,10 mM glucose.

Data Collection

The expression profiles of RNAs were screened from the National Center of Biotechnology Information Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/gds/). The GSE30528¹⁷ dataset is composed of 22 mRNA expression data of human kidney samples, including 9 DN patients and 13 healthy controls. For dataset GSE96804,¹⁸ there are 61mRNA expression data from human kidney samples including 41 DN patients and 20 healthy controls. The R software package was used to process the downloaded files and to convert and reject the unqualified data. The data were calibrated, standardized, and log2 transformed. For dataset GSE131882,¹⁰ there are three controls and three DNs’ single-cell sequencing data from human kidney samples.

Differential Expression Analysis

The differently expressed mRNAs between Diabetic Human Kidney and Control Human Kidney groups were screened using Limma package, with the criterion of |log ₂(fold change [FC]| > 1.5 and p-value < 0.05).

Weighted Gene Co-Expression Network Analysis

The expression profile of 22 samples was obtained from the GSE30528 dataset. After validation, these data were used to construct a co-expression network using the WGCNA package in R (version 3.6.0). The WGCNA methodology analysis was performed as a previous report.¹⁹

Gene Ontology and Pathway Enrichment Analysis

To assess the function of DEGs in DN, Gene Ontology annotation and Kyoto Encyclopedia of Genes and Genomes pathway analyses were performed by using the cluster profiler package of R software.²⁰ P-value < 0.05 was set as the cut-off criterion.

Gene Set Enrichment Analysis

GSEA was performed to examine critical pathways represented under different conditions.²¹ The estimated expression levels of all identified genes were applied to GSEA and then ES values were calculated according to the ranked-ordered gene list. As the predefined gene sets, the Kyoto Encyclopedia of Genes and Genomes (KEGG) information deposited in NCBI database was analyzed using an in-house script.

Validation with qRT- PCR

QRT-PCR was used to confirm mRNA expression. Total RNA obtained from participants’ serum with Trizol (Invitrogen, Carlsbad, CA, USA). RNA purity was measured using the NanoDrop Spectrophotometer (Nanodrop™ Thermo Fisher Scientific, Vantaa, Finland). Then Power SYBR Green (Takara) was used to perform qRT-PCR. β-actin was applied as an internal reference. Primer sequences are summarized as follow, β-actin (F: 5′-CTCCATCCTGGCCTCGCTGT-3′, R: 5′-GCTGTCACCTTCACCGTTCC-3′), C3 (F: 5′-GGGGAGTCCCATGTACTCTATC-3′, R: 5′-GGAAGTCGTGGACAGTAACAG-3′), C7 (F: 5′-TGTAAAACGACGGCCAGT-3′, R: 5′-CAGGAAACAGCTATGACC-3′). All PCR reactions were conducted in triplicate. Relative expression was calculated using the 2^−ΔΔCt method. miRNAs were synthesized by Genepharm Technologies (Shanghai, China). Sequences of the miRNAs used in the study are as follows: miR-494-3p (5′-UGAAACAUACACGGGAAACCUC-3′), miR-574-5p (5′-UGAGUGUGUGUGUGUGAGUGUGU-3′), miR-6769b-5p (5′- UGGUGGGUGGGGAGGAGAAGUGC-3′) and miR-18a-5p (5′- UAAGGUGCAUCUAGUGCAGAUAG −3′). miRNAs were transfected using Lipofectamine 2000 reagent.

Single-Cell RNA Sequencing Reanalysis

The Single-nucleus sequencing data from GSE131882 dataset was downloaded from GEO database. Single nucleus RNA sequencing of three early human diabetic kidney samples and three controls were included.¹⁰ The single-cell data was visualized by Kidney Integrative Transcriptomics (K.I.T.) database (http://humphreyslab.com/SingleCell/), which was created by Ben Humphrey’s lab of Washington University. The database consists of RNA-seq and snATAC-seq to mouse and human kidney and lung.

Regulatory Predictive Analysis

The miRwalk tool²² and TargetScan 7.2 (http://www.targetscan.org)²³ were used to predict the candidate regulatory target miRNAs of C7.

Luciferase Reporter Gene Assay

The Mutant (MUT) 3ʹUTR and wild-type (WT) 3ʹUTR sequences of C7 were amplified and cloned into a pGL3-Basic reporter vector (Promega). MES cells from miR-NC, miR-494-3p and miR-574-5p groups were placed in 24-well plates and transfected with pGL3-MUT and pGL3-WT plasmids (Promega). After cultured for 48 h at 37°C, cells were collected. Firefly and Renilla luciferase activities were measured using Luciferase Assay System (Promega). Firefly luciferase activity was normalized to Renilla luciferase activity.

Statistical Analysis

Data analysis was performed using GraphPad prism 6.0 (GraphPad software, San Diego, CA, USA). The data are expressed as means ± SEMs. Unpaired t-test was used to distinguish the differences between two groups. P< 0.05 was considered as a statistically significant difference with 2-tailed tests. Receiver operating characteristic (ROC) curve was used to analyze the diagnosis value, including the area under the curve (AUC), sensitivity, and specificity.

Results

Weighted Coexpression Network Construction and Key Module Identification

To determine the co-expression network most highly associated with diabetic nephropathy progression, WGCNA was performed on the GSE30528 dataset. A total of 22 samples with clinical data were included in the coexpression analysis (Figure 1A). In this study, a power of β =17 (scale-free R² = 0.85) was selected as the soft-thresholding parameter to ensure a scale-free network (Figure 1B and C). Then, Pearson correlation coefficients were calculated for pairwise genes to yield a similarity matrix, which was transformed into an adjacency matrix using the threshold and power values listed above. A total of 10 modules were identified via average linkage hierarchical clustering (Figure 1D). The green module was found to have the highest association with DN development (green, r=0.71, p=2e⁻⁴) (Figure 1E), and this module was selected as the significant module for further analysis. The correlations of these 10 modules are shown in Figure 1F.

Figure 1 Co-expression network construction of 22 diabetic nephropathy samples of GSE30528 dataset via WGCNA. (A) The clustering was based on the expression data of differentially expressed genes between tumor samples and normal samples in DN. (B) Relationship between scale-free topology model and soft-thresholds (powers). (C) Relationship between the mean connectivity and various soft-thresholds (powers). (D) Dendrogram of all differentially expressed genes clustered based on a dissimilarity measure (1-TOM). (E) Heatmap of the correlation between module eigengenes and clinical traits. (F) The eigengene adjacency heatmap of the modules.

Functional Enrichment Analysis

There were 35 DEGs in the green module (Table S1). The relative expression of these genes in normal and DN samples is shown in Figure 2A. To identify the biological function and network of the 35 DEGs, the protein–protein interaction was assessed by the STRING database (Figure 2B). We found that the expression of 34 genes were upregulated and only 1 gene had downregulated expression. Then, KEGG pathway enrichment analysis was used to explore the functions of the DEGs. The top 10 most significantly enriched pathways are shown in Figure 2C, in which C3, C1qrs and C6, 7, 8, and 9 were all upregulated. This may suggest that the “complement and coagulation cascades” pathway is activated in DN (Figure 2D).

Figure 2 The functional analysis of DEGs in GSE30528 dataset. (A) Heatmap of the DEGs in green module by WGCNA of GSE30528 dataset. (B) The protein-protein interaction network of the DEGs in green module. (C) The top 10 significant KEGG pathways. (D) The pathway map of Complement and coagulation cascades. The red squares represent upregulated genes.

Validation by Another GEO Dataset

To confirm that the “complement and coagulation cascades” pathway is significantly changed and may play a key role in DN development, another GEO dataset (GSE96804) was used. First, 126 DEGs (|log2FC|>1.5, p <0.05) were identified. We annotated the DEG functions using gene ontology (GO) analysis (Figure S1A–D). Then, KEGG pathway analysis based on the DEGs was performed (Figure 3A). As a result, we found that the “complement and coagulation cascades” pathway was significantly enriched. In addition, GSEA of the whole genes also found that this pathway was promoted in DN (p = 0.0039) (Figure 3B and C). All DEGs identified in the pathway were upregulated, including KNG, C3, CR1, Clqrs and C6, 7, 8, and 9. These results were consistent with our previous analysis of the GSE30528 dataset. Moreover, we compared the common DEGs between the GSE30528 green module and the overall GSE96804 dataset. There were nine common DEGs, including C3 and C7 (Figure 3D). Furthermore, the expression of these two genes was upregulated (Figure S2A–F) in both the GSE30528 and GSE96804 datasets, which agreed with the KEGG analysis.

Figure 3 Complement and coagulation cascades was significantly activated in DN. (A) The top 10 significant KEGG pathways in GSE96804 dataset. (B) GSEA analysis of GSE96804 dataset. (C) The pathway map of Complement and coagulation cascades. The red squares represent upregulated genes. (D) The Venn map of common DEGs between GSE30528 green module and GSE96804 dataset.

Validation of C3 and C7 Expression with qRT-PCR and ROC Curve Analysis

To further assess the expression of C3 and C7, a total of 31 DN patients and 30 healthy controls were enrolled as a validation cohort. qRT-PCR technology was used to confirm the differential expression levels from participant serum samples. Consistent with the microarray data, C7 expression was significantly upregulated (Figure 4A) between DN and HC patients. However, there was no difference in the levels of C3 (Figure 4B). To assess the potential value of increased C7 expression for DN diagnosis, we further performed ROC curve analysis. We found that the ROC curve of C7 expression showed a significant distinguishing efficiency with an AUC value of 0.865 (95% CI: 0.774–0.957, ***p <0.001) (Figure 4C), with the best cutoff value of 3.4, a sensitivity of 90.3% and a specificity of 73.3%, which indicated that C7 could be a candidate biomarker for DN diagnosis.

Figure 4 RT-qPCR analysis of C7 and C3. (A) Validation the expression level of C7 (**p< 0.01). (B) Validation the expression level of C3 (p>0.05). (C) The receiver operating characteristic curve of C7 in distinguishing DN. Abbreviations: DN, diabetic nephropathy; HC, healthy control.

Validation of Increased C7 Levels in a Single-Cell RNA Sequencing Reanalysis

To explore complex cellular events and deepen our understanding of biological systems, we reanalyzed the expression of common DEGs between the GSE30528 and GSE96804 datasets in different cell groups in DN. The kidney cortex was sampled from three nondiabetic controls and 3 DN patients following nephrectomy in the GSE131882 dataset. As shown in Figure 5A, a total of 20 DEGs were differentially expressed in 9 cell groups. Notably, the upregulated DEGs C7, TNC, CDH6, and FN1 were specifically expressed in renal cells. We found that C7 expression was mainly significantly upregulated in mesangial (MES) cells. In addition, a total of 23,980 nuclei passed filters with an average of 2541 genes and 6894 unique molecular identifiers per nucleus. By using the Kidney Interactive Transcriptomics database (http://humphreyslab.com/SingleCell/), eleven kidney cell types were identified by unsupervised clustering and expression analysis of lineage-specific markers following batch correction (Figure 5B). We also found that C7 expression was mainly significantly upregulated in the MES cells (Figure 5C–E). This suggested that C7 was mainly secreted by MES cells. In addition, the expression proportion of C7 in the MES cells of DN samples was much higher than that of the MES cells of the control samples. Similarly, we performed the same analysis of C3; however, C3 was not significantly enriched in any cell type (Figure S3A–D). Thus, the single-cell RNA sequencing reanalysis indicates that C3 may not be a biomarker of DN. Furthermore, we reanalyzed the GSE9006 dataset and found that C7 was not significantly expressed in type 1 and 2 diabetic samples (Figure S4A–B). These results revealed that C7 is a potential MES cell-specific biomarker of DN.

Figure 5 The scRNA-seq analysis validation of C7. (A) The expression of common DEGs between GSE30528 and GSE96804 datasets in different types of cell groups in DN were predicted. PCT, proximal convoluted tubule; CFH, complement factor H; LOH, loop of Henle; DCT, distal convoluted tubule; CT, connecting tubule; CD, collecting duct; PC, principal cell; IC, intercalated cell; PODO, podocyte; ENDO, endothelium; MES, mesangial cell; LEUK, leukocyte. The yellow diamond represents the cell type, red represents up-regulation, blue represents down-regulation, and the circle size is based on Degree. (B) DN and control samples were integrated into a single dataset. (C) The distribution of C7 in different cells between healthy controls and DN. (D) The relative expression level of C7 in identified cell types between healthy controls and DN. (E) The percent expression of C7 in identified cell types between healthy controls and DN.

MiR-494-3p and miR-574-5p are Regulators of C7

To investigate the underlying regulatory mechanism of C7, we identified the differentially expressed miRNAs in DN based on the GSE161884 dataset (Figure 6A and B). Then, we predicted the target miRNA that regulates C7 using TargetScan (http://www.targetscan.org/) and miRwalk (http://mirwalk.umm.uni-heidelberg.de/) software. The common differentially expressed miRNAs and predicted target miRNAs were identified, including miR-494-3p, miR-574-5p, miR-6769b-5p and miR-18a-5p (Figure 6C). In the next step, we verified the expression of these four miRNAs in MES cells. As shown in Figure 6D, the relative expression of C7 mRNA in MES cells was dependent on glucose concentration. In addition, the expression of miR-494-3p and miR-574-5p significantly decreased as glucose concentration increased. Hence, we focused on these two miRNAs to determine their regulatory relationship with the relative expression of C7 in MES cells. The results suggest that miR-494-3p and miR-574-5p mimic can downregulate the expression of C7 (Figure 6E). A luciferase reporter assay confirmed that C7 could bind to miR-494-3p and miR-574-5p (Figure 6F).

Figure 6 MiR-494-3p and miR-574-5p are regulators of C7. (A) and (B) The heatmap and volcano plot of differentially expressed miRNAs. (C) The common miRNAs compared with miRNAs predicted by Targetscan and miRwalk. (D) The relative expression of C7 and different miRNAs in MES cell treated with different glucose. (E) The relative C7 expression in DN cell treated with miR-494-3p mimics and miR-574-5p mimics. (F) Detection of the binding of miR-494-3p and miR-574-5p with C7 by luciferase reporter gene assay. *p < 0.05, **p< 0.001. Abbreviation: N.S., no significant.

Discussion

To date, several factors, including urinary albumin levels, have been considered for assessing DN. However, the diagnosis of DN greatly depends on assumptions, and newer specific markers for DN, especially for early-stage disease, are urgently needed.²⁴ To date, several studies have revealed that C7 can be used as a potential biomarker of DN,^17,25 but there is a lack of comprehensive analysis based on multiple data sets and a lack of further exploration at the cellular level. WGCNA is a commonly used bioinformatics analysis tool used to identify the key modules and genes associated with specific clinical traits.²⁶ Single-cell RNA sequencing is a precise way to detect molecular changes in diseases, provides deep insights into cell-to-cell variation in the tumor and microenvironment, allows high-resolution dissection of the pathogenic mechanisms of diseases, and affords potential clinical utilities.²⁷ This approach has been increasingly used in disease research.⁹

In this study, we performed a WGCNA, a single-cell sequencing analysis and other analyses to detect a new potential biomarker for DN diagnosis. To further verify C7 as a novel biomarker of DN with the ability to distinguish it from other nephropathies, four datasets (GSE30528, GSE96804, GSE9006, and GSE131882 datasets) with different control samples were selected. Through a series of bioinformatics analyses, we found that complement and coagulation cascades were significantly enriched in DN. This may suggest that the complement cascade plays an important role in DN development. The complement system is a proteolytic cascade in the blood plasma and a mediator of innate immunity, a nonspecific defense mechanism against pathogens.²⁸ There are three pathways of complement activation: the classical pathway, the lectin pathway, and the alternative pathway.²⁹ All these pathways generate crucial enzymatic activity that, in turn, drives the activation of effector complement molecules. The main consequences of complement activation are the opsonization of pathogens, the recruitment of inflammatory and immunocompetent cells, and the direct killing of pathogens.³⁰ In cancer, the complement system plays a major and complex role in killing antibody-coated tumor cells, supporting local chronic inflammation, or hampering antitumor T cell responses favoring tumor progression.³¹ The complement system is also involved in various autoimmune diseases, such as primary biliary cholangitis (PBC), primary sclerosing cholangitis (PSC) and autoimmune hepatitis (AIH).³² While the complement cascade is an important component of the innate immune system, uncontrolled activation can cause severe disease. In a recent review, excessive complement activation in atypical hemolytic uremic syndrome caused renal failure if untreated.³³

C7 encodes a serum glycoprotein that forms a membrane attack complex (MAC) together with complement components C5b, C6, C8, and C9 as part of the terminal complement pathway of the innate immune system.³⁴ This protein initiates membrane attack complex formation by binding the C5b-C6 subcomplex and inserts into the phospholipid bilayer, serving as a membrane anchor.³⁵ A previous study showed that C7 might be a potential tumor suppressor.³⁶ C7 is also a critical component of the terminal pathway of complement activation.³⁵ C7 is expressed in endothelial cells as a trap for assembling MACs. MAC deposition has been reported in smooth muscle,³⁷ dermatomyositis microvasculature^38,39 and in the mechanism of advanced diabetic kidney disease.⁴⁰ Thus, MAC formation will increase the expression level of C7. In our study, qRT-PCR validation and ROC curve analysis were applied to evaluate the diagnostic value of C7. ScRNA-seq analysis revealed that C7 was specifically highly expressed in the MES cells in DN. These findings showed that in the early stages of diabetes, there is increased expression of C7 in DN. C7 is released into circulation and is reflected as high levels of C7 in the blood. These findings are supported by those of other studies, in which vascular deposition of complement proteins has been reported.^17,25 For example, Falk⁴¹ reported complement deposition in the kidneys of diabetic patients. Sircar et al²⁵ reported that C7 expression was upregulated in DN, but did not describe a specific mechanism. Although the upregulation of C7 expression in DN has been reported, the specific cell type that expresses C7 is unknown. Single-cell sequencing has paved the way for further research on C7. C7 is released into circulation and is reflected as high levels of C7 in the blood. Our results showed that as a biomarker, C7 was not significantly expressed in type 1 and 2 diabetic samples, which was in accordance with a previous study.⁴¹ This is the first study to reveal that elevated C7 gene and gene products in MES cells is a potential biomarker for early diabetic nephropathy.

MiRNAs are noncoding RNAs that can inhibit the expression of RNA through translation inhibition or messenger RNA degradation.^42,43 In recent years, studies have revealed that some miRNAs could be regulators of immune and inflammatory responses and are potential therapeutic targets in DN.^15,16 Single nucleotide polymorphisms of miRNA-related genes can affect the production, expression, and function of mature miRNAs, thus affecting disease susceptibility and the immune response. We identified the differentially expressed miRNAs in DN using the GSE161884 dataset. Bioinformatics retrieval and luciferase reporter assays confirmed that miR-494-3p and miR-574-5p could be regulators of C7 in DN. C7 is the direct target of miR-494-3p and miR-574-5p, and these miRNAs can negatively regulate its expression level. This is the first study to reveal that the elevated C7 gene expression level in MES is regulated by miR-494-3p and miR-574-5p in early diabetic nephropathy.

Taken together, our findings indicate that C7 might represent a novel biomarker that is specifically expressed in MES cells in diabetic nephropathy that can be used for diagnosis and that C7 expression levels are regulated by miR-494-3p and miR-574-5p. However, several questions need to be addressed in future studies, such as determining the functions of C7 at the cellular level and in animal models and investigating and why C7 expression is upregulated in MES cells in DN.