Weighted Gene Correlation Network Analysis (WGCNA) Detected Loss of MAGI2 Promotes Chronic Kidney Disease (CKD) by Podocyte Damage

Chronic kidney disease (CKD) is a worldwide public health problem, with increasing prevalence, poor outcomes and high treatment costs [1]. Currently, It is widely accepted that podocyte damage is not only associated with the degree of proteinuria, but also correlates with glomerulosclerosis development and renal function decline [2, 3]. Podocyte damage is one of the early-stage events in many human kidney diseases, including minimal change disease (MCD), focal segmental glomerulosclerosis (FSGS), diabetic kidney disease (DKD), membranous nephropathy (MN), lupus nephritis (LN) and other CKDs [4, 5]. The realization that CKD could be reversed by targeting this early-stage event has led to the current theory that novel therapeutics should inhibit podocyte loss. However, the pathogenesis of podocyte loss is complex and not well clarified, and the exact molecular mechanism of podocyte damage in CKD remains largely unexplored.

Traditional molecular biology, which focuses on assessing the functions of single genes or proteins, plays an important role in disease research at the molecular level. However, it is undeniable that traditional molecular biology can only explain the biological system locally, and a more comprehensive exploration of the whole biological system in diseases is very challenging.

Weighted gene correlation network analysis (WGCNA), a powerful method for constructing co-expression networks based on expression data, can reconstruct gene co-expression modules and summarize modules using module eigengenes (ME) and intramodular hub genes [6]. In this work, we utilized a computational strategy to perform a systematic assessment of differentially expressed genes (DEGs) based on 4 datasets, including the expression data from glomeruli of patients with 3 types of CKD (LN, DKD and IgAN) and the podocyte injury model induced by puromycin amino-nucleoside (PAN). Then, DGEs were used for WGCNA based on the GSE99339 dataset including 133 glomerular expression profiles of various CKD patients. WGCNA results combined with DAVID findings and podocyte standard genes (PSGs) defined by previous reports [7, 8] were used to identify podocyte injury-associated modules [9], and hub genes were detected by intramodular connectivity. Finally, the Nephroseq online tool (www.nephroseq.org Apr 2017, University of Michigan, Ann Arbor, MI) was used to assess the clinical significance of hub genes, and small interfering RNA (siRNA) was used to explore the role of hub genes in the physiological functions of podocytes.

The present study firstly provides a comprehensive understanding of podocyte injury-associated genes via a co-expression network, which may guide subsequent experimental studies on podocyte protection in CKD patients.

The 5 raw gene microarray expression data were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), and included GSE66107, GSE93798, GSE30528, GSE32591 and GSE99339. There were 3 normal podocyte specimens and 3 injured podocyte samples induced by puromycin amino-nucleoside (PAN) in the GSE66107 dataset. The remaining samples selected in this study were from the glomerular compartment of the human kidney, as follows. The GSE93798 dataset included samples from 20 IgA nephritis (IgAN) patients and 22 controls; GSE30528 encompassed specimens from 9 diabetic kidney disease (DKD) patients and 13 controls; since we only evaluated genome expression data of the glomeruli section, glomeruli samples from the 32 lupus nephritis (LN) and 14 glomeruli control patients were selected for further analysis, which only constituted a part of GSE32591. The GSE99339 dataset included 133 samples based on the GPL19184 platform, which consisted of various types of established CKDs such as LN, IgAN, membranous glomerulonephritis (MGN), hypertensive nephropathy (HT), focal segmental glomerular sclerosis (FSGS), tumorous nephrosis (TN), minimal change nephropathy (MCD), DKD, etc. GSE66107, GSE93798, GSE30528 and GSE32591 data were used for analyzing differential expressed genes (DEGs). Then, DEGs altered in more than one dataset were used in WGCNA to construct a co-expression network based on their normalized expression data in the GSE99339 dataset.

Raw expression data were preprocessed with R v3.4.1 (https://www.r-project.org/). The preprocessing of gene expression profile data was performed using the Robust Multi-array Average (RMA) algorithm in affy package of Bioconductor (http://www.bioconductor.org/), including background correction, quantile normalization and pro summarization. To minimize false positives and maintain high expression levels of DEGs for downstream analysis, we only kept the expressed probes with the criterion of “present (P)” in more than 50% of all samples for each dataset, which was identified by the mas5calls function in the affy package [10]. Finally, a total of 26357, 11118, 11296, 7194 and 7114 gene probes were identified in the GSE66107, GSE93798, GSE30528, GSE32591 and GSE99339 datasets for subsequent preprocessing, respectively. As the GSE99339 dataset included several sample batches, the empirical Bayes algorithm (sva package) [11] was used to correct the batch effect of gene expression value in GSE99339. The detailed information for the samples included in this study is summarized in Table 1.

The Linear Models for Microarray Data package in Bioconductor was applied to identify DEGs by comparing expression values between groups. The corresponding p values of gene symbols after T-test were defined as adjusted p-values; adjusted p-value < 0.05 and |fold change(FC)| ≥ 1.5 were selected as cut-off criteria for DEGs.

To explore the interactions between genes, a system biology approach, WGCNA, which converts co-expression measures into connection weight or topology overlap measures, was applied for gene co-expression network construction [12]. The co-expression methodology was typically used for exploring the correlation between gene expression levels. Genes involved in the same pathway or functional compound tend to show similar expression patterns [13]. Therefore, the construction of a gene co-expression network facilitates the identification of genes with similar biological functions [14]. In this work, all DEGs altered in more than one datasets were selected from GSE99339 and entered to construct weighted co-expression modules using the WGCNA package in the R language.

Ju et al. [7, 8] have developed an iterative algorithm (in silico Nano dissection), which is a machine learning based approach to statistically define cell type-specific genes using high-throughput data. By applying the algorithm to kidney microarray datasets, it was cross-validated with previously reported data and predicted cell type-positive standard genes for mesangial cells and podocytes in kidney diseases. According to the results, podocyte-positive standard genes (PSGs) included ACTN4, NPHS1, NPHS2, SYNPO, PTPRO, NES, DDN, LRRC7, WT1, TCF21, MAGI1, MAGI2, RAB3A, KIRREL, KIRREL2, KIRREL3, CD2AP, CASK, IQGAP1, TJP1, TRPC6, NCK2, DES, PDPN, PODXL, LMX1B, FAT1, CDH3, PLCE1, CLIC5, CDH13, MME, UCHL1, CD80, Sv2b, PALLD, FYN, PLAUR, DAG1, AGRN, EFNB1, EZR, FOXC2, MAFB, UTRN, MYOC, BASP1, SULF1, SCEL, and CX3CL1.

DAVID (Database for Annotation, Visualization and Integrated Discovery) (http://david.abcc.ncifcrf. gov/) [15] is a gene functional classification tool that integrates a set of functional annotation tools for investigators to analyze biological functions in massive amounts of genes. GO (Gene Ontology) function and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses of all modules were performed based on DAVID.

The immortalized mouse podocytes MPC5 (provided by Dr. Peter Mundel, Harvard University) have been described in detail previously [16]. Podocytes were maintained in RPMI-1640 (Gibco, USA) supplemented with 5% FBS. To propagate podocytes, recombinant mouse γ-IFN (50 U/ml, R&D Systems, Minneapolis, MN) was added to the culture medium for incubation at 33°C. To induce differentiation, podocytes were maintained at 37°C for 14 days, and used in experiments. Adriamycin (ADR) is widely used to assess renal diseases by inducing nephrotic syndrome in vivo and podocyte injury in vitro. Therefore, ADR (0.25µg/ml TargetMol, USA) was used to cause podocyte injury, by incubation at 0, 24 and 36h, respectively. Then, cells of each group were harvested at 48h after differentiation. MAGI2 siRNA and scramble siRNA control were obtained from Applied Biological Materials (ABM) Inc. Transfection was performed with Lipofectamine RNAiMAX reagent (Invitrogen), according to the manufacturer’s protocol. Forty-eight hours after transfection, total RNA and protein were extracted. All experiments were performed in triplicate.

Total RNA from podocytes was extracted using TRIzol reagent (TaKaRa Bio, Japan). The gene expression levels of MAGI2, podocin (NPHS2), CD2-associated protein (CD2AP), synaptopodin (SYNPO) were assessed by quantitative real-time reverse transcription polymerase chain reaction (qRT-PCR) using the RT-PCR and SYBR Green kit (TaKaRa Bio, Japan) on a LightCycler 480 Real Time PCR System. The primers used were: MAGI2, 5′- TAAAGACACGACGAACCCCA-3′ (forward) and 5′- AAGTTGGGTCTGGGCCTATC′ (reverse); NPHS2, 5′- CCCAGGTACACAGAGCACAA -3′ (forward) and 5′- CTCACGTCACAACCTTCAG-3′ (reverse); CD2AP, 5′- AGCTGCTTTCCTAACTCCGT -3′ (forward) and 5′- CGCATGGCCTTCTCTTCTTC -3′ (reverse); SYNPO, 5′-AGTGAGGAAGAGGAAGTGCC -3′ (forward) and 5′- GAGGCTAGTGTGGTGAGTGT -3′ (reverse). Conditions for quantitative PCR were 37°C for 15 min and 85 °C for 5 s for RT, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. GADPH was used as an endogenous control for data normalization.

Immunoblot was used to assess the protein amounts of MAGI2, NPHS2, CD2AP and SYNPO. Briefly, cells were lysed with RIPA buffer (Thermo Fisher Scientific, USA), and Western blot was performed by standard procedures. The primary antibodies used were goat anti-MAGI2 polyclonal antibody (R&D Systems, USA), rabbit anti-NPHS2 polyclonal antibody (Proteintech Group, China), rabbit anti-CD2AP polyclonal antibody (Santa Cruz Biotechnology, USA), and rabbit anti-SYNPO polyclonal antibody (Proteintech, CHINA). Anti-goat and anti-rabbit horseradish peroxidase-conjugated (HRP) antibodies (Proteintech, China) were used as secondary antibodies. Immunoreactive bands were detected by enhanced chemiluminescence (ECL; Beyotime Biotechnology, China).

On coverslips, podocytes were fixed with 4% paraformaldehyde, followed by permeabilization with 0.3% Triton X-100. For F-actin staining, fixed and permeabilized cells were incubated with fluorescein isothiocyanate (FITC) -phalloidin (Yeasen Bio, Shanghai, China). For detecting the NPHS2 protein, after incubation with primary antibody for two hours at room temperature, the slides were incubated with Alexa Fluor® 647 Goat Anti-Rabbit IgG (Invitrogen). Analysis was performed by confocal laser-scanning microscopy (Zeiss Lsm510 Meta, Jena, Germany). Photomicrographs of podocytes stained with various antibodies were selected randomly and analyzed in a blinded manner.

Based on the above cut-off criteria, a total of 2969 up-regulated and 3185 down-regulated DEGs were obtained in injured podocyte samples induced by PAN compared with controls (in GSE66107). Meanwhile, there are 1678 up-regulated and 2204 down-regulated DEGs in glomeruli from IgAN patients compared with controls (in GSE93798). A total of 344 up-regulated and 1078 down-regulated DEGs were identified in DKD patients compared with controls (in GSE30528); 473 up-regulated and 367 down-regulated DEGs were identified in LN patients compared with controls (in GSE30528). As shown in Supplementary Fig. 1A-D (for all supplementary material see www.karger.com/10.1159/000495205), hierarchical cluster analysis of the data revealed that DEGs could be used to distinguish injured podocytes or most CKD samples from respective controls except for three samples included in the GSE93798 dataset (Supplementary Fig. S1A-D). The details of the top 50 up-regulated and 50 down-regulated DEGs in each dataset are provided in Supplementary Table S1-4. The majority of DEGs (87.5% up-regulated and 82.0% down-regulated) of PCGs were detected in only one dataset (Supplementary Fig. S1E, F). Meanwhile, only 15 DEGs (4 up-regulated and 11 down-regulated) showed differential expression in all four gene expression profiles, including CDKN1C, MYO5C, C1orf21, KBTBD11, NR3C2, FZD2, MAGI2, LPL, ZFPM2, USP46, PRKAR2B, TRIM38, DESI1, CTSS and ITGB2. Since these 15 DEGs were differently expressed not only in injured podocytes but also in glomeruli of CKD patients, they were considered biochemical makers of CKD, especially in case of podocyte injury.

Weighted gene correlation network analysis (WGCNA) has been successfully applied in cancer-related studies [6, 17]. In this study, we applied WGCNA to explore the potential mechanisms of podocyte damage using DEGs obtained in the above analysis. A total of 9217 DEGs explored were altered in more than one datasets, and 1260 DEGs which were overexpressed in one dataset but downregulated in another compared to corresponding controls were removed. Therefore, 7957 DEGs remained for analysis. Consequently, 4031 DEGs found in the processed data of GSE99339 were selected for further investigation. Fig. 1A shows the cluster of 133 glomeruli samples from various CKD patients, with no abnormal sample. Then, a co-expression network was constructed with the WGCNA package in the R software. To achieve a scale-free topology (R²> 0.9), we set β = 6 in Fig. 1B, then converted the pairwise correlation into an adjacency matrix of connection strengths through the soft-thresholding approach (connection strength = |correlation|^β). A dissimilarity matrix based on topological overlap measure (TOM) was used to identify gene modules with a dynamic tree-cutting algorithm [6]. All modules were assigned corresponding colors (Fig. 1C). In this study, 4031 DGEs were ultimately categorized into 12 modules with sizes ranging from 52 to 1224 genes. Each module was assigned a unique color: turquoise, blue, brown, yellow, green, red, black, pink, magenta, purple green-yellow and tan represented 1214, 708, 526, 235, 211, 195, 100, 87, 70, 58, 57 and 52 genes, respectively. Of all genes, 518 were not assigned to any of the 12 modules, and designated as grey (Table 2, Fig. 1C). A complete list of genes in each module is shown in Supplementary Table S5.

Module eigengenes (ME) were used to represent each module, which was determined by the first principal component. First, we analyzed the relationship between the modules and pathological parameters by calculating the correlation value of ME for each module (Supplementary Figure S2). The results showed no specific link between single modules and serum creatinine (Scr) or estimated glomerular filtration rate (eGFR). This could be because of the small sample size, and altered gene expression in glomerular damage is an early stage event in CKD, while pathological changes of cardinal features such as tubular cell atrophy and interstitial injury, as common cardinal features in CKD, mainly lead to decreased eGFR. Meanwhile, several modules were significantly associated with certain clinic parameters, e.g. the turquoise (P=0.045) module, which was negatively associated with body mass index (BMI), while significant positive correlation was found between the yellow module and mean blood pressure (P=0.008).

To further determine the biological functions of these genes in the 12 modules, DAVID [9] was used to mine the modules’ biological significance, including GO biological process (BP) terms and KEGG pathways. Supplementary Table S6 lists the top 30 significantly enriched GO BP terms (p< 0.05) for each module. Fig. 2 displays three representative terms for each module. Supplementary Table S7 lists complete KEGG pathway terms (p< 0.05) for each module. Table 3 lists the top three KEGG pathways with p< 0.05 for each module. According to major biological processes and the number of PSGs included in various modules, the latter were roughly parceled into two sections, including podocyte-specific and nonspecific injury associated modules, as follows.

For instance, the yellow and brown modules. The yellow module, which contained 251 genes and only two PCGs, and was positively associated with blood pressure (r=0.23, P=0.008, Supplementary Figure S2), showed functional enrichment in inflammatory response, immune response, cellular response to lipopolysaccharide and apoptotic process, and so on (Fig. 2, Supplementary Table S6). Pathway analysis further revealed that genes in this module were enriched in signal transduction pathways such as NF-kappa B [18], TNF [19],Toll-like receptor [20] and NOD-like receptor [21] signaling pathways (Table 3, Supplementary Table S7). These GO terms or pathways were shown to be associated with podocyte damage and blood pressure regulation [22, 23] in multiple reports. However, these pathways are not podocyte-specific, and also involved in the functions of other kidney cells, including mesangial [24] and endothelial [25] cells. Similar observations were made for the brown module, which was composed of 533 genes and no PSG, and negatively correlated with blood pressure (r =-0.2, P=0.021, Supplementary Figure S2). The BP terms of the brown module were associated with the functions of oxidation-reduction, metabolic process and fatty acid beta-oxidation, and their disorders were revealed to be connected with a predisposition to cell injury (Fig. 2, Supplementary Table S6). In addition, pathway analysis showed genes in the brown module were mainly enriched in Metabolic pathways, Carbon metabolism, and Glyoxylate and dicarboxylate metabolism, which are also associated with the regulation of cell metabolism (Table 3, Supplementary Table S7). Therefore, it is reasonable to infer that genes in the yellow and brown modules may play roles in maintaining survival of not only podocytes, but also glomerular endothelial cells and mesangial cells.

According to major biological processes and KEGG pathway analysis, the green and red modules stood out. There were 211 DEGs in the green module, including 10 PSGs and 6 co.DEGs (Table 2), most than in the other modules. as shown in Fig. 2 and Supplementary Table S6, genes in the green module were significantly enriched in processes such as extracellular matrix disassembly, glomerulus development, and glomerular visceral epithelial cell development. Significant KEGG pathways included several signaling pathways relevant to podocyte function such as Adherens junction and the Rap1 signaling pathway (Table 3, Supplementary Table S7). The results were confirmed and visualized in a functionally grouped network created by the Cytoscape software v3.5.1 combined with ClueGO plug-in (v2.3.5) (Supplementary Fig. 3). Similarly, 211 genes in the red module, which contained 4 PSGs and 1 co.DEGs, also showed significant enrichment in functions involved in response to endoplasmic reticulum stress, cellular response to mechanical stimulus and regulation of focal adhesion assembly. Most of these GO terms in the red and green modules have been reported to associate with podocyte damage in CKD patients.

To explore the relationship between modules, clustering analysis for ME was applied. The results showed the green and red modules were very close, suggesting that genes in the red and green modules played similar roles in maintaining podocyte function. Especially, the green module included more PSGs and co.DEGs than other modules, so we suspect that genes included in the green module were more likely to specifically participate in triggering the injury mechanism of podocytes compared with those of other modules. Therefore, the green module was defined as podocyte-injury associated module in the present study.

Highly connected hub nodes are central to the network’s architecture [6] and studies suggested that genes more centralized in the network are more likely to be key drivers to respective cellular functions compared with peripheral genes [26]. Therefore, more important nodes can be obtained by identifying hub nodes. In this study, the top 5% of nodes with highest connectivity were defined as hub genes in each module [27]. Fig. 3A-D are the networks of hub genes in the brown, yellow, green and red modules, only displaying connections with weight (w) above threshold values of 0.3, 0.15, 0.05, and 0.05, respectively. Clearly, hub genes (red nodes) exhibited high connectivity with neighboring genes whose functions are consistent with GO and KEGG findings. Taking the brown module as an example, which showed functional enrichment in metabolic and oxidation-reduction processes, some hub genes have been also reported to participate in similar processes, e.g. the SLC13 Family (SLC13A3) [28, 29], SLC22 Family (SLC22A6) [30] or MSRA that plays a protective role in the progression of UUO-induced kidney fibrosis via suppression of fibrotic responses caused by oxidative stress [31]. Similarly, some hub genes in the yellow module have been identified as pro-inflammatory factors in a variety of human diseases, including CTSS and ITGB2 [32], consistent with the BP term or KEGG pathway enriched in the yellow module. Meanwhile, some well-known inflammatory factors also had high connectivity to these hub genes, e.g. IL1B, CXCR4 and CASP1. These observations further supported the notion that proposed roles for hub genes of unknown functions may be inferred from clusters of genes similarly expressed across many biological conditions.

For the red module, one of the hub genes was ARHGEF12, which was highly connected to the well-known podocyte-associated gene CD2AP, which was reported to be enriched in glomerular expression and necessary for normal glomerular filtration barrier function by spatiotemporal regulation of RhoA activity [33, 34]. The most important module was the green, which was defined as podocyte injury-associated module according to the above analysis. As shown in Fig. 3C, it included more PSGs and co.DEGs compared with all other modules. In this module, 11 hub genes exhibited high connectivity with neighboring genes (FGF1, AIF1, CR1, NPSH2, PLCE1, etc.), whose functions are involved in stabilizing the intracellular environment and maintaining podocyte survival by balancing anabolic metabolism and catabolism. Specifically, NPSH2, a podocyte-specific gene that encodes a protein that plays a role in regulating glomerular permeability was found; mutations in this gene cause steroid-resistant nephrotic syndrome [35]. MAGI2, one of the hub genes in the green module, particularly attracted our attention; it was also defined as one of the PSGs and significantly downregulated in injured podocytes and the glomeruli of IgA, DKD or LN patients by gene expression data analysis. This strongly suggested that MAGI2 not only could be well suited as a CKD biomarker, especially in case of podocyte damage, but also contributed to the process of occurrence and development of podocyte damage.

To further explore the relationship between CKD and MAGI2 in patients, Nephroseq (www.nephroseq.org Apr 2018, University of Michigan, Ann Arbor, MI) data were analyzed in various CKD patients. As shown in Fig. 4A-F, the expression levels of MAGI2 were reduced in the glomeruli of LN (Fig. 4A), IgAN (Fig. 4B) and DN (also referred to as DKD, as above) (Fig. 4D) patients, in agreement with the above expression profile data, as well as in vasculitis (Fig. 4C) and focal segmental glomerulosclerosis (FSGS) samples (Fig. 4E, F) compared with normal controls. Meanwhile, the expression levels of MAGI2 were decreased in LN patients with WHO class III (Fig. 4G) and IV (Fig. 4H), compared with WHO class II cases. In addition, MAGI2 was downregulated in kidney transplant patients with acute rejection compared with cases without rejection (Fig. 4i). What’s more, the expression levels of glomerular MAGI2 were significantly positively correlated with the glomerular filtration rate (GFR) in patients with nephritic syndrome (NS) (Fig. 4J) or diabetic nephropathy (DN) (Fig. 4K), and negatively associated with urinary protein content in LN patients (Fig. 4L). These results suggested that reduced expression of glomerular MAGI2 may participate in the pathogenesis of kidney damage and CKD progression.

To confirm the role of the MAGI2 gene in the pathogenesis of podocyte damage, conditionally immortalized transgenic mouse podocytes (MPC5) were used. These cells grow into a cobblestone morphology at the permissive temperature of 33°C and are maintained at 37°C for 14 days before experimentation to induce differentiation. First, we assessed MAGI2 expression in undifferentiated podocytes or injured podocytes induced by Adriamycin (ADR), which was often used in nephropathic models. As shown in Fig. 5A and C, MAGI2 protein and mRNA levels were increased in differentiated podocytes compared with the undifferentiated counterparts. Meanwhile, in differentiated podocytes, exposure to ADR caused a substantial loss of MAGI2 at the mRNA and protein levels (Fig. 5A and C). To determine if MAGI2 influences cytoskeletal remodeling in podocytes, as an early event in CKD development that drives a shortening and simplification of podocyte foot processes, we knocked down MAGI2 in differentiated podocytes by siRNA interference. Western blot and real-time PCR confirmed the knockdown efficiency of MAGI2 siRNA (50 nM) at 48h. In addition to MAGI2 reduction, the mRNA and protein levels of NSPH2, SYNPO and CD2AP, which have fundamental functions in basic biochemical processes in podocytes [36] (Fig. 5B and D), were also decreased after MAGI2 silencing. Similar results of down-regulated NSPH2 are shown in Fig. 5E, which depicts confirmation by immunofluorescence. In addition, F-actin had overtly recovered arrangement by confocal fluorescence microscopy and FITC-phalloidin staining. In the scramble siRNA group, F-actin was characterized by highly ordered, parallel, contractile actin filament bundles, and disorganized, with reorganized, short, branched actin filaments filled in the cytoplasm of podocytes treated with MAGI2 siRNA 50 nM for 48h (Fig. 5E). What’s more, flow-cytometric analysis showed that apoptosis was increased after MAGI2 silencing (Fig. 5F and G). These results suggested that MAGI2 inhibition promoted cytoskeletal rearrangement and injury in podocyte cultures.

As a global health problem, CKD affects 10-16% of the adult population in Asia, Europe, and the USA [37-39], and may further develop into kidney failure. Over the past ten years, accumulating evidence indicates that proteinuria is associated with renal function decline, and severe proteinuria can individually predict all-cause mortality in CKD patients, independent of numerous confounding factors [40]. It is widely accepted that effective interventions to reduce proteinuria in CKD patients should constitute the top priority for future studies [41]. Although emerging evidence suggests that podocytes serve as the final barrier to urinary protein loss through the maintenance of foot-processes and an interposed slit-diaphragm (SD) [42], current understanding of the mechanism of podocyte injury is very limited and remains to be elucidated.

Gene expression analysis by microarrays has become one of the most widely used high-throughput methods for gathering genome-wide functional data. In this study, we first analyzed differential gene expression in four datasets, including typical podocyte injury models induced by PAN and microarray data of glomeruli from three common CKD types, to provide insights into the possible mechanisms of podocyte damage and facilitate the early prediction of podocyte injury in CKD patients. As a result, 15 DEGs (co.DEGs) were altered significantly in the four datasets, including CDKN1C, MYO5C, C1orf21, KBTBD11, NR3C2, FZD2, MAGI2, LPL, ZFPM2, USP46, PRKAR2B, TRIM38, DESI1, CTSS and ITGB2. Identification of these dysregulated genes helps predict patient outcomes and guides treatments. Meanwhile, some of these genes have been previously reported. For example, CDKN1C, also known as P57, is very important in cell cycle regulation, and cells deficient in this gene exhibit high rates of cell growth and decreased differentiation, and CDKN1C/p57 is down-regulated in many proliferation diseases, including Beckwith-Wiedemann Syndrome and autosomal dominant polycystic kidney disease (ADPKD) [43, 44]. In addition, previous findings reveal that damaged podocytes may inhibit p57 protein expression, but activate a cyclin D1-independent cell-cycle mechanism to promote glomerular epithelial cell (GEC) proliferation in the cellular lesions of FSGS [45]. The NR3C2 gene encodes the mineralocorticoid receptor that mediates aldosterone effects on salt and water balance within restricted target cells. Mutations in NR3C2 cause autosomal dominant pseudo hypoaldosteronism type I [46], a disorder characterized by urinary salt wasting. A recent study showed that NR3C2 is downregulated in stable and chronic rejection (CR) patients compared with kidney-transplant recipients with “operational tolerance” [47]. CTSS, also known as cathepsin S, is a lysosomal cysteine proteinase that may participate in the degradation of antigenic proteins to peptides for presentation on MHC class II molecules; in addition, circulating cathepsin S levels are associated with GFR decline in mice and humans [48], and specific inhibition of cathepsin S could diminish the effects of cardiovascular disease (CVD), especially in patients with CKD [49]. Still, only few co.DEGs have been reported in association with CKD to date, especially for their roles in podocyte damage. Thus, the unreported co.DEGs could provide helpful information for possible biomarkers in further studies.

It is well accepted that co-expressed genes are more likely to be co-regulated and functionally related [50]. Therefore, identifying co-expressed protein-coding genes can help assign the functions of disease-causing gene [51]. This has been successfully applied to assess protein-coding genes, e.g. distinguishing dysfunctional regulatory subnetworks and identifying candidate biomarkers in multiple diseases [52]. However, few studies have used this tool to assess genes related to podocyte damage. The objectives of the current study were to gain molecular insights into CKD progression and identify candidate gene groups involved in podocyte damage. Finally, we obtained 12 modules with sizes ranging from 52 to 1,214 genes. Next, we revealed that the green module containing 211 DEGs likely played key roles in the multistep development of podocyte damage, covering a wide range of functions in glomerulus development, extracellular matrix disassembly, and glomerular visceral epithelial cell development, and also comprised more co.DEGs and PSGs compared with the other modules. Meanwhile, multiple genes in the green module have been reported to be involved directly in podocyte lesion, e.g. MME [53], CR1 [54], NPHS2 and IQGAP2 [55]. These results imply that genes in the green module may play important roles in multiple CKDs, particularly participating in podocyte damage. Therefore, it is reasonable to define the green module as a podocyte injury-associated module; this provides a comprehensive analysis of podocyte injury-associated genes based on the co-expression network, and may guide subsequent experimental studies on podocyte protection.

It should be noted that membrane associated guanylate kinase, WW and PDZ domain containing 2 (MAGI2) was differentially expressed in all 4 datasets, and also defined as the one of hub genes in the green module according to average connectivity. MAGI2 was also one of the PSGs, also known as AIP1 or SSCAM. Its encoded protein is characterized by two WW domains, a guanylate kinase-like domain, and multiple PDZ domains, with structural similarity to the membrane-associated guanylate kinase homologue (MAGUK) family. As a tight junction protein in epithelial tissues, MAGI2 has a critical role in maintaining the functional structure of the slit diaphragm, and is essential in the functioning of the kidney filtration barrier [56]; as one of the WT1 target genes, it is required for podocyte development [57]. More importantly, MAGI2 null mice show progressive proteinuria as early as 2 weeks postnatally [58]. In the present study, glomerular MAGI2 expression and its clinical implication were examined in many types of CKD using Nephroseq data, and the results showed significantly downregulated MAGI2 in CKDs, corroborating previous findings. In addition, there were negative and positive associations of MAGI2 expression with proteinuria and eGFR, respectively. Podocytes express several specific cytoskeletal proteins, including the slit diaphragm (SD)-localized protein NPHS2 (podocin) and actin-associated proteins SYNPO (synaptopodin) and CD2AP, which participate in the maintenance of slit diaphragm and foot process integrity [59-61]. Finally, we used MAGI2 knockdown in podocytes to examine how this gene affects glomerular function and CKD progression; the results indicated that MAGI2 silencing in podocytes led to significantly reduced expression levels of NPHS2, CD2AP and SYNPO, and caused distinct cytoskeletal rearrangement, which is widely considered an early event in the development of CKD. This observation could at least in part explain the clinical significance of MAGI2 in CKD, as mentioned above. These results imply that the MAGI2 gene may constitute a candidate biomarker in CKD diagnosis and serve as a therapeutic gene target for podocyte protection in CKD patients.

This study presented a novel approach using WGCNA to explore gene changes during the CKD process, especially when podocyte damage occurs. According to the network constructed by WGCNA, 12 modules were identified, and four were selected for detailed analysis. The green and red modules were identified as podocyte injury-associated modules, which may provide new insights into pathways and genes underlying CKD and podocyte damage. Meanwhile, as a hub gene of the green module, MAGI2 was associated with CKD progression, with its loss promoting podocyte cytoskeletal rearrangement. Overall, WGCNA is an efficient approach to systems biology. By using this procedure, the key node gene MAGI2 was identified, and found to be involved in the regulation of cytoskeletal rearrangement in podocytes, with its loss predisposing to proteinuria and CKD.

This work was supported by the National Natural Science Foundation of China (Grant No. 81600549, to Wan-Peng Wang), the natural science research project of Huai’an city (Grant No. HAB201737, to Wan-peng Wang); National Natural Science Foundation of China (Grant No. 81700586, to Jian-Xiao Shen); the Foundation Research Project of Jiangsu Province (The Natural Science Fund No. BK20160707 to Zhi Zuo); Shanghai Science Health Bureau(No. 20154Y0003, to Xing-hua Shao); Chinese medical association special foundation for medical science(No.15020110599 to Xing-hua Shao). All experiments of this study were authorized by the Ethics Committee of LianShui County People’s Hospital.

The authors declare that they have no competing interests.

Z. Zuo, J.-X. Shen and Y. Pan contributed equally to this work.