An Integrative in silico Study to Discover Key Drivers in Pathogenicity of Focal and Segmental Glomerulosclerosis

Focal and segmental glomerulosclerosis (FSGS) is a clinical-pathologic syndrome characterized by segmental and focal lesions in glomeruli [1, 2]. Proteinuria and two microscopic features in renal biopsies, including segmental sclerotic lesions and podocyte foot process effacement are the main features of the FSGS [3]. A variety of factors can result in FSGS and usually, these heterogeneous factors are the basics of FSGS classification. The secondary type of FSGS can be caused by renal or nonrenal disorders, such as conditions that produce overactivity (hyperfiltration), stress, or high blood pressure, which impact the glomeruli. The genetic forms of FSGS are caused by an abnormal gene version that damages the glomeruli. Primary FSGS (or idiopathic) FSGS, on the other hand, has unknown origins. It is the most prevalent form of FSGS and is diagnosed when no secondary or genetic cause can be found [4, 5].

Glucocorticoids, immunologic drugs, and renin-angiotensin system blockers are the main currently used therapies for patients with FSGS; however, these treatments are not able to completely remove the disease manifestations and hinder its progression [6]. In addition, despite a huge number of studies, to date, the main underlying molecular mechanisms and the key drivers in the pathogenesis of FSGS are yet to be discovered [7]. Aiming to an efficient targeted therapy, it seems that more investigations are needed to shed a light on the disease’s underlying molecular mechanisms and to identify the key players in its pathogenicity.

Lately, a growing body of in silico and bioinformatics studies have analyzed high-throughput datasets aiming to understand the molecular basis of diseases. In this context, analyzing transcriptomics datasets could give a holistic view of expressional alterations in disease versus normal states [8-10]. Also, this practical strategy has become a popular one in the discovery of key genes as potential therapeutic targets or biomarkers in different diseases [11, 12].

Weighted gene co-expression network analysis (WGCNA) is one of the famous algorithms in the co-expression analysis of transcriptomics datasets. With applying a soft-threshold, WGCNA will introduce weight values to determine the interactions among genes. Along with building the weighted co-expression networks, the WGCNA algorithm is able to cluster co-expressed genes into individual modules and identify the most correlated module/s to a specific trait. Up to now, a huge number of investigations have applied the WGCNA in different biological processes like genetics, cancer, and brain imaging data analysis aiming for the discovery of novel biomarkers or therapeutic targets [13-17].

The objectives of this study were to analyze one FSGS-related transcriptomics dataset using the WGCNA algorithm to identify the most disease-correlated genes as potential drug targets and understanding the underlying pathogenic pathways in this pathological condition. In brief, after the analysis of the dataset, the intensity values of the identified differentially expressed genes (DEGs) in the array set were introduced to WGCNA in the R software environment. Afterward, the genes in the disease-correlated module were subjected to different functional analyses and key genes were identified by applying an integrative approach, considering the top genes in protein-protein interaction (PPI) network and a list ranked by module membership (kME) values of genes in the most disease-correlated module. Figure 1 shows the workflow of the present study.

The FSGS dataset GSE129973 related to the glomeruli tissue of FSGS patients with 20 control and 20 case samples was downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). In order to detect and remove possible outliers, principal component analysis (PCA) was performed prior to the main analysis. Using the Network analyst online tool (http://www.networkanalyst.ca), samples were subjected to normalization steps and after removing the dataset probes related to multiple genes, Linear Model for Microarray Analysis (Limma) was selected for the analysis procedure and significant DEGs were identified based on false discovery rate (FDR) cutoff <0.049.

The WGCNA algorithm was used in the R programming environment to find co-expressed modules. After introducing the intensity values related to the identified DEGs to the WGCNA algorithm and sample clustering step, the “pickSoftThreshold” feature was used to find the best soft-threshold power from 1 to 20. The degree of independence and mean connectivity values were used to determine an effective soft threshold. To obtain the co-expression modules, several steps were completed, including the construction of an adjacency matrix, the estimation of a topological overlap matrix, module construction, and dynamic branch cutting with a merging threshold of 0.25. Following parameters were designated for identification of the co-expression modules: “soft-threshold power = 6, minModuleSize = 30, mergeCutHeight = 0.25.” To describe the adjacencies among genes in the identified modules and to verify the module’s division accuracy, a heatmap based on topological overlap matrix dissimilarity was created. After the assignment of samples as control or FSGS, the most correlated module was found by calculating the association between the co-expression modules and these two states (phenotypes).

Among all the DEGs, only the ones with |log2 fold change (FC)| ≥ 1 were regarded as seeds for the construction of a PPI network. The network was built using the STRING database (confidence >0.4) [18], and Cytoscape was used for network visualizations. After topological analysis of the constructed network, DEGs with a degree of centrality ≥5 were considered as hub genes.

Functional enrichment analyses, including gene ontology (GO) and pathway enrichment analyses, were performed on the isolated genes from the module of interest, and the selected DEGs based on log2 FC cutoff: |log2 FC| ≥ 1. The GO and Reactome pathway enrichments were performed using Cytoscape software (version 3.8.2) [19] and the CluGO module (version 2.5.7) [20] (significant enrichment threshold was set at p < 0.05).

The selection of candidate hub genes was performed using both the kME value of genes in the selected module, as well as the degree connectivity of the genes in the constructed PPI network among the selected DEGs based on log2 FC. The kME is simply defined as the correlation of a gene with the module eigengene [21]. As mentioned above, after the construction of a PPI network considering DEGs with |log2 FC| ≥ 1, the genes with a degree of connectivity ≥5 were chosen as hub genes in the network analysis strategy. In addition, considering the kME and gene significance (GS) values (kME > 0.8, GS > 0.6), the top genes of the selected module were extracted and considered as hub genes. Using a Venn diagram, the intersection genes of hub genes in the module of interest and hub genes in the PPI network were identified as the key genes in FSGS disease.

GSE121233 and GSE104948, two other FSGS array sets, were used to verify the expression pattern of the identified hub genes. Both of the chosen datasets included microarray data from FSGS and healthy individuals’ glomeruli tissues. GSE121233 included data from 5 healthy and 5 FSGS individuals and GSE104948 included data from 18 healthy and 10 FSGS patients. The same protocol was used to analyze these two datasets as it was for the main dataset GSE129973.

To ensure the precision of the core analysis, multiple preprocessing steps, including PCA and normalization procedures were completed before the dataset analysis. The resulted PCA plot showed no outlier among all the control or FSGS samples (Fig. 2a). Additionally, quantile normalization was performed to ensure that the expression distributions of each sample are identical over the entire dataset. According to the FDR cutoff, 4,843 significant DEGs (2,125 downregulated and 2,718 upregulated) were detected and extracted for further analysis. A volcano plot and a heatmap representing top DEGs based on log2 FC and top 50 DEGs based on FDR value are shown in Figure 2b and c. The results of the analyzed dataset, including DEGs’ names and values, are provided in an online supplementary File (for all online suppl. material, see www.karger.com/doi/10.1159/000524133) (online suppl. file 1 – DEGs and network analysis results).

Figure 3a and b shows the results of sample clustering and the scale independence and mean connectivity plots. The obtained results of the WGCNA algorithm like cluster dendrogram, network heatmap plot, trait heatmap, and mean connectivity platform of the selected module are depicted in Figure 4a–f. According to the results of scale independence and mean connectivity, a soft threshold of 6 was selected as the best power to get an approximate scale-free topology. Furthermore, genes were clustered into 9 separate co-expression modules after hierarchical clustering and module merging steps (black, blue, brown, cyan, green-yellow, magenta, purple, magenta, and tan) (Fig. 4a, b). Figure 4c depicts a network heatmap plot displaying the module division accuracy. The produced heatmap depicted the topological overlap adjacency among genes across modules, and as can be noticed, genes in the same module have a higher correlation. The eigengene adjacency heatmap and clustering dendrogram also revealed the division of the established modules into two clusters (Fig. 4d). Names and kME values of the clustered genes in co-expressed modules are provided (online suppl. file 2 – GS and kME values).

The most disease-relevant modules were defined by calculating the associations of each module and two phenotypes (healthy and disease states). The black module showed the highest positive correlation with the disease state (r = 0.87; p = 2E−18) (Fig. 4e). Magenta (r = 0.76; p = 1E−8) and purple (r = 0.67; p = 3E−6) were other modules showing a positive correlation with the disease state. Other modules, including blue, brown, cyan, green-yellow, and salmon showed negative correlations with the disease state. Scatter plot representing GS versus kME values of all genes in the black module is shown in Figure 4f. A high correlation between GS and kME values in the black module (r = 0.83, p < 1e−200) implying that the genes are associated with both the phenotype (disease state) and the module eigengene.

The results of functional enrichment analyses on the black module’s genes are shown in Figure 5. According to the Reactome pathway analysis results, the genes were mostly enriched in vesicle-mediated transport, immune system, and cell cycle control. Biological process terms also confirmed the Reactome analysis results, where the primarily enriched terms were related to vesicle-mediated transport, intracellular protein transport, and cell cycle regulation. Other top enriched terms of both categories were protein modifications and ephrin receptor signaling pathway. The GO terms of molecular function showed that the black module genes were mainly enriched in the RNA binding, cadherin binding, and kinase binding. About the cellular component terms, the top enriched terms were nucleus, cytoskeleton, focal adhesion, Golgi membrane, lysosome, and endocytic vesicles.

A PPI network with 104 nodes and 302 edges was constructed with the participation of DEGs with |log2 FC| ≥ 1 (Fig. 2e) (online suppl. file 2). In the constructed PPI network, 45 DEGs with a degree cutoff ≥5 were considered as hub genes. GO enrichment analysis was performed for all the DEGs with |log2 FC| ≥ 1 and in terms of biological process, genes were mostly enriched in inflammatory responses, extracellular matrix organization, platelet degranulation, and response to lipopolysaccharides. Integrin binding and heme binding were the top enriched molecular function terms and in the case of cellular component terms, genes were mostly enriched in extracellular vesicles (EVs) and extracellular places.

In terms of hub gene selection, 45 DEGs were identified as hub genes in the constructed PPI network and were listed for further key gene selection steps (Table 1). In the WGCNA part, all the DEGs were subjected to the analysis and after clustering, and analyzing the module-trait relationships, the black module’s genes were recognized as the most disease-correlated genes. In this step, 193 DEGs were selected as hub genes according to both GS and kME values in the black module. Considering both the lists of hub genes coming from the WGCNA study, as well as topological analysis of the PPI network, two genes, including bone morphogenic protein 2 (BMP-2) and collagen type IV alpha 1 (COL4A1), were recognized as key genes with a potential role in the pathogenicity of FSGS (Fig. 6a). The expression profiles of the identified key genes in other FSGS-related datasets (GSE121233 and GSE104948) confirmed the upregulation patterns of these two genes with close log2 FC values (Fig. 6b).

The progression to ESRD is most common in patients with FSGS among glomerular diseases and the current treatments are not efficient enough to impede such progress [22]. According to research, up to 56% of patients who receive a kidney transplant will experience FSGS recurrence in the transplanted kidney [23, 24]. Corticosteroid and immunosuppressive therapy are widely used to treat FSGS; however, both drug types have proven to have a limited response in resistance and FSGS recurrence patients [25]. Thus, there is a crucial need for more efficient therapies to manage FSGS. The application of high-throughput transcriptomics techniques will pave the way for better understanding the disease underlying molecular mechanisms and the discovery of novel molecules for diagnosis and targeted-based treatment. The present experiment applied an integrative bioinformatics analysis on a transcriptomics dataset coming from FSGS patients in order to identify the most disease-correlated genes as potential drug targets and clarify the involved molecular pathways of the disease. After the analysis of the transcriptomics dataset and identification of DEGs, two approaches, including WGCNA and network analysis were followed to reach the most disease-correlated co-expression modules and top genes in a constructed PPI network.

According to the obtained results, BMP-2 and COL4A1 were two identified key genes spotted by the two applied strategies in the present study on the GSE129973 dataset, which was confirmed by two other FSGS transcription datasets, including GSE121233 and GSE104948. However, in articles related to these datasets, the importance of overexpression of BMP-2 and COL4A1 has not been considered and discussed [26, 27]. However, the results of our analysis emphasized the role of these two important genes in FSGS. COL4A1 encodes the alpha 1 chain of type IV collagen, which is a key structural component of the extracellular matrix and is mainly expressed in the brain, muscles, kidneys, and eyes [28, 29].

It has been shown that COL4A mutations are associated with idiopathic FSGS and predispose some patients to FSGS [30-32]; probably because abnormal expression of type IV collagen can damage the renal parenchyma cell [31]. In addition, increased levels of the transforming growth factor-beta (TGF-ß) signaling have also been shown to induce aberrant collagen IV isotype (4A1 and 4A2) synthesis in the glomerular basement membrane of the kidneys as a condition that is closely linked with fibrosis and sclerosis [33]. In order to identify new biomarkers and therapeutic targets for podocytopathies such as FSGS, one experiment analyzed a podocyte-specific mRNA expression profile in two different models of FSGS (Actn4−/−podocytes) at 2, 6, and 44 weeks of age. As a result, other types of collagen, including collagen VIII alpha 1 and collagen III alpha 1 were recognized as two of the 25 most highly upregulated genes in the FSGS model at 6 weeks compared with control group [34].

Regarding the abnormal expression of COL4A in other chronic kidney diseases where FSGS is seen as a histologic pattern, the same overexpression of collagens is seen as a marker of interstitial fibrosis and glomerular sclerosis [35, 36]. During fibrosis, scar tissue formation in the interstitial space is the consequence of increased extracellular matrix and elevated accumulation of ECM components [37]. A study by Holderied et al. [38] has also shown a link between the thickness of Bowman’s capsule and the enhanced levels of collagen type IV (all six α-chains) as a result of activated parietal epithelial cells in human diabetic nephropathy.

BMP-2 is a low molecular weight glycoprotein that belongs to TGF-β superfamily [39]. It has been shown that BMP-2 acts as a novel fibrosis antagonizing cytokine partly by downregulating type I TGF-ß receptors and Smads [40]. Generally, it is clear that the lack of ECM homeostasis regulation can cause FSGS. Based on a study conducted on the ECM proteome changes in a murine FSGS model, some of the BMP signaling inhibitors like sclerostin were among the downregulated genes in the models [41]. Elevated levels of the Bmp-2 gene were also shown in previous experiments analyzing the gene expression profiles of podocytes, mesangial cells, and endothelial cells of the Cd2ap and bigenic Cd2ap+/−, Fyn−/− mutant mouse models of FSGS [42, 43].

The dysregulation of BMP-2 has also been seen in other chronic kidney diseases because of the key role of TGF-β superfamily in these diseases. In addition, the released EVs from the damaged endothelial cells of the CKD patients were shown to have a high content of BMP-2 and calcium, which are the mediators of calcification and osteogenic differentiation of vascular smooth muscle cells. Based on research, stimulation of human endothelial cells by TNF-α could result in nuclear translocation and activation of NF-κB and consequently induction of BMP-2 expression in cells [44].

According to the GO and pathway enrichment analysis on the most correlated co-expression module of black, most of the genes were acting in extracellular places and vesicles, in the immune systems pathways, and in vesicle-mediated transport pathways. Involvement of the immune systems pathways has previously been shown in different CKDs and FSGS [8, 45]. Despite the well-known role of EVs in cell waste excretion, they also were shown to be engaged in different physiological and pathological conditions by participating in intercellular communications [46, 47]. EVs were also considered as main participants in various functions affecting the ECM components such as their synthesis, organization, degradation, and calcification [48, 49]. Some studies support the notion that fibrosis is an EV-regulated pathology in multiple organs such as kidneys. In one experiment, in vivo inhibition of exosome secretion in CKD animal models ameliorated renal fibrosis and decreased the expression of fibronectin, αSMA, and collagen I [50]. Munkonda et al. [51] revealed that in vitro addition of EVs from human differentiated podocytes to cultured human proximal tubule epithelial cells increase p38 and Smad3 phosphorylation and activate TGF-βR1, resulting in enhanced expression of fibronectin and collagen IV. On the other hand, EVs can also serve as biomarkers since they comprise modified cargo components in diseased conditions when compared with a healthy state [52, 53]. As well as, within the kidney, the origination of EVs could be endothelial cells, podocytes, or blood cells and they can be found within the urine, tissue, or circulation. Therefore, they can be a source of biomarker detection in FSGS [52].

One limitation of the present study was about the type of FSGS disease which was unknown to us (not mentioned in the analyzed dataset) and might affect our final results. In addition, although COL4A and BMP-2 can be introduced as key involved genes in the FSGS pathogenesis because they are mainly involved in fibrosis and sclerosis, we cannot say for sure whether they are disease-specific or related to CKD in general. It seems that more investigations having other control types from other CKDs are needed to verify our result.

The results of the present study added and highlighted some information about the underlying processes and involved molecules in the pathogenesis of FSGS. Two hub molecules (BMP-2 and COL4A) and some pathways (vesicle-mediated transport) were recognized as potential players in the pathogenesis of FSGS. The potential of these targets to stop the pathologic processes especially in CKDs like FSGS is something that needs vigorous investigations.

We thank members of the Regenerative Medicine Research Center for their help in some parts of the bioinformatic analyses steps.

An Ethics statement was not required for this study type, no human or animal subjects or materials were used. However, the protocol of this experiment was reviewed by Ethics Committee of Isfahan University of Medical Sciences.

The authors declare that they have no conflict of interest.

The authors received no funds for this work.

Amir Roointan and Alieh Gholaminejad participated in design, analysis of dataset using the WGCNA algorithm, interpretation of data, and drafting the manuscript. Maryam Ghaeidamini was participated in dataset selection, data analysis, as well as preparing the manuscript. Jesus Simal-Gandara helped in interpretation of the analyzed data and preparing the manuscript draft. All authors reviewed the manuscript.

The analyzed dataset by the present study is available in the GEO repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129973).