Rapid Identification of Potential Drugs for Diabetic Nephropathy Using Whole-Genome Expression Profiles of Glomeruli

Objective. To investigate potential drugs for diabetic nephropathy (DN) using whole-genome expression profiles and the Connectivity Map (CMAP). Methodology. Eighteen Chinese Han DN patients and six normal controls were included in this study. Whole-genome expression profiles of microdissected glomeruli were measured using the Affymetrix human U133 plus 2.0 chip. Differentially expressed genes (DEGs) between late stage and early stage DN samples and the CMAP database were used to identify potential drugs for DN using bioinformatics methods. Results. (1) A total of 1065 DEGs (FDR < 0.05 and fold change > 1.5) were found in late stage DN patients compared with early stage DN patients. (2) Piperlongumine, 15d-PGJ2 (15-delta prostaglandin J2), vorinostat, and trichostatin A were predicted to be the most promising potential drugs for DN, acting as NF-κB inhibitors, histone deacetylase inhibitors (HDACIs), PI3K pathway inhibitors, or PPARγ agonists, respectively. Conclusion. Using whole-genome expression profiles and the CMAP database, we rapidly predicted potential DN drugs, and therapeutic potential was confirmed by previously published studies. Animal experiments and clinical trials are needed to confirm both the safety and efficacy of these drugs in the treatment of DN.


Introduction
Diabetic nephropathy (DN), which is clinically characterized by proteinuria and morphological and ultrastructural changes in the kidney, is a serious complication of diabetes mellitus and is a major cause of end-stage renal disease worldwide. DN is a multifactorial progressive disease with complex pathogenesis, involving hyperglycemia, advanced glycation end products (AGEs), hemodynamic disorder, metabolic abnormalities, inflammatory factors, and oxidative stress [1]. Although our knowledge of DN is continuously increasing, no treatment strategies specifically target the pathogenesis of DN beyond controlling glucose levels, blood lipid levels, and high blood pressure [2]. As a result, the prognosis for most DN patients is poor, especially for those in late stages of the disease. The identification of potential drugs targeting the molecular pathogenesis of DN is critical for improving the prognosis and survival of patients with DN.
Whole-genome expression profiling is the simultaneous measurement of the expression of thousands of genes by microarray technology (or RNA-Seq) to create a global picture of tissue or cellular function. Comparing the wholegenome expression profiles of tissues (or cells) under physiologic and pathologic conditions may reveal the pathogenesis of DN. In addition to identifying differentially expressed and coexpressed genes, from which one can generate new hypotheses about the molecular mechanism of complex diseases, whole-genome expression data are also used to identify therapeutic drugs. The Connectivity Map (CMAP) database is a collection of genome-wide transcriptional expression data from cultured human cells treated with bioactive small molecules [3]. Lamb et al. findings showed that genomic signatures in the CMAP database can be used to identify potential new therapeutics, and signatures are often conserved across diverse cell types [3]. Therefore, the CMAP database can be used with whole-genome expression profiles to identify potential drugs for DN in the glomeruli of patients. Using the CMAP database, Zhong et al. predicted that the combination of an angiotensin-converting enzyme (ACE) inhibitor and a histone deacetylase inhibitor would maximally reverse the disease-associated expression of genes in a mouse model of HIV-associated nephropathy (Tg26 mice), and the renoprotective effect of the combined use of these inhibitors was proven in Tg26 mice [4]. It is feasible to utilize gene expression profiles of tissues under normal and physiopathological conditions to investigate potential therapeutic drugs based on bioinformatics methods. However, this kind of therapeutic drug identification in DN research is lacking. In our study, we utilized the gene expression profiles of microdissected glomeruli from DN patients and explored potential therapeutic drugs using in silico screening approaches.

Patients.
The clinical study used a cross-sectional design. The control and DN kidney samples were obtained from leftover portions of diagnostic kidney biopsies. For the kidney biopsies, informed consent was obtained from the donors and patients. All of the participants provided written informed consent. The institutional review board of Jinling Hospital specifically approved this study [5].
A total of 18 DN patients diagnosed by renal biopsy were enrolled in the study. The baseline clinical characteristics of the DN patients are listed in Table 1.
For each biopsy specimen, light microscopy, immunofluorescence, and electron microscopy were routinely performed. Sections for light microscopy were stained with hematoxylin eosin, periodic acid-Schiff, Masson's trichrome, and periodic acid methenamine silver. All of the patients were categorized based on the pathologic classification of the Renal Pathology Society [5]. The glomerular classifications were as follows: class I, glomerular basement membrane thickening; class IIa, mild mesangial expansion; class IIb, severe mesangial expansion; class III, nodular sclerosis; and class IV, global glomerulosclerosis in >50% of glomeruli. Interstitial fibrosis and tubular atrophy (IFTA) were scored as follows: 0, absent; 1, <25%; 2, 25-50%; and 3, >50% of the total area. Interstitial inflammation was scored as follows: 0, absent; 1, inflammation only in relation to IFTA; and 2, inflammation in areas without IFTA. Arteriolar hyalinosis was scored as follows: 0, absent; 1, at least one area of arteriolar hyalinosis; and 2, more than one area of arteriolar hyalinosis. Arteriosclerosis was scored as follows: NA, absence of large vessels; 0, no intimal thickening; 1, intimal thickening less than thickness of media; and 2, intimal thickening greater than thickness of media. All of the specimens were scored by the same pathologist (Dr. Feng Xu) who was blinded to the clinical findings. In order to assess the reliability and reproducibility of the classification, biopsy slides were scored independently by another pathologist (Dr. Dandan Liang). The pathologic characteristics of the DN patients are listed in Table 1.
The DN patients were divided into 2 groups according to the following criteria: early stage DN group ( = 6), glomeruli isolated from the renal tissue of early stage DN patients who were identified with eGFR > 90 mL/min; late stage DN group ( = 12), glomeruli isolated from the renal tissue of late stage DN patients with eGFR between 15 mL/min and 60 mL/min. Control samples ( = 6) were obtained from living donor kidney biopsies. Control subjects were defined as having an eGFR of more than 90 mL/min, the absence of proteinuria, normal serum creatinine, and BUN.

Tissue Handling and Microdissection.
Tissues placed in RNALater (SIGMA, St. Louis, MO, USA) were manually microdissected at 4 ∘ C for glomeruli. In general, 10 glomeruli were collected from each biopsy tissue and were placed into cold RNeasy lysis buffer solution (Qiagen, Valencia, CA, USA).

RNA Extraction and Amplification.
Dissected glomeruli were homogenized, and RNA was prepared using RNeasy mini columns (Qiagen, Valencia, CA, USA), according to the manufacturer's instructions. RNA quality and quantity were determined using the Laboratory-on-Chip Total RNA Pico Kit Agilent Bioanalyzer. Samples without evidence of degradation were further amplified using the Ovation RNA amplification system kit (NuGEN, San Carlos, CA, USA).

Affymetrix Microarray Data and Preprocessing. The
Affymetrix microarray platform (Human U133 Plus 2.0) was used to produce the whole-genome gene expression profile data. Quality control and data processing were performed using R [6] and Bioconductor [7]. The CEL source files were processed into expression estimates, and background correction and quartile data normalization were performed using the RMA (robust multiarray average) algorithm [8].

Screening of Differentially Expressed Genes (DEGs).
The limma package [9] in R language was used to screen DEGs by pairwise comparison between groups. The statistical method implemented in the limma package is based on an approach called linear models. We used the method proposed by Benjamini-Hochberg (BH) for multiple testing correction. The adjusted values were the false discovery rates (FDR). The threshold criterion is a combination of FDR < 0.05 and fold change > 1.5. The DEGs between late stage and early stage DN samples were chosen for DN drug identification.

Validation of Microarray Expression Data.
The relative mRNA levels of 10 genes were validated in new selected glomerular samples. The clinical and pathologic characteristics of these DN patients are listed in Table 2. The processes used for patient screening, tissue handling, microdissection, and total RNA were performed as previously described. The mRNA levels of the target genes were analyzed by quantitative real-time RT-PCR (qRT-PCR) using the Applied Biosystems6 7900HT Fast Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). The qRT-PCR results were normalized to 18S ribosomal RNA using the 2 −ΔΔCT method [10], and significance was set to < 0.05.  between the late and early stage DN samples, and the upregulated probes and the downregulated probes were saved in GRP files. The CMAP website provides a KS statistic algorithm; we uploaded the GRP files as the query gene signature, which was then compared to each rank-ordered list to determine whether upregulated query genes appeared near the top of the list and downregulated query genes near the bottom (positive connectivity) or vice versa (negative connectivity), yielding a connectivity score ranging from 1 to −1. A high negative connectivity score indicated that the corresponding perturbagen reversed the expression of the query signature and might have the potential to treat DN. [4]. The IDs of the probe table downloaded from the CMAP database were converted to Entrez gene symbols using the Affymetrix lookup table associated with the platform. If more than one probe ID corresponded to the same gene, the gene rank was condensed to the median of the probe ranks. Then, we extracted the 500 top and bottom DEGs. Potential DN drugs should reverse the DEGs of DN; genes upregulated in DN should be downregulated by a potential DN drug and vice versa. The DEGs (FDR < 0.05 and fold change > 1.5) between the late and early groups were used to calculate the reversing scores using formula (1). The perturbagen with the highest reversing scores might have the potential to treat DN. Consider

Drug Identification Using a Matching Algorithm
where Score indicates the reversing score for a drug in th experiment in the CMAP database; ∩ indicates the intersection between two sets; up indicates a list of upregulated genes during disease; down indicates a list of downregulated genes during disease; up indicates genes upregulated by a drug in the th experiment in the CMAP database; and down indicates genes downregulated by a drug in the th experiment in the CMAP database.

Functional Enrichment Analysis.
The functional enrichment analysis of the screened DEGs and the genes reversed by potential drugs was performed via the GeneAnswers package in R language [11]. The GeneAnswers package functionally categorizes genes based on Fisher's exact test with annotation libraries of the gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG).

DEGs in the Glomeruli of DN Patients.
We separately compared the samples in the 2 stages of DN with the controls. A total of 105 DEGs were identified between the early stage DN samples and the controls, including 54 upregulated genes and 51 downregulated genes; 2572 DEGs were identified between the late stage DN samples and the controls, including 1626 upregulated genes and 946 downregulated genes.
Only 105 DEGs were identified between the early stage DN samples and the controls, and the enriched GO categories were primarily involved in "response to stimulus." These results were in accordance with the mild pathological changes in the glomeruli from early stage DN patients.
In contrast, 1065 DEGs were identified between the late and early stages of DN samples, including 815 upregulated genes and 250 downregulated genes. The heat map in Figure 1 shows the expression levels of the top 100 regulated genes across the 24 samples. As shown in Figure 2, the most enriched GO categories of the DEGs were "extracellular matrix," "protein binding," "cell adhesion," and "immune system process." The most enriched KEGG pathways were "ECM-receptor interaction," "complement and coagulation cascades," "focal adhesion," "cytokine-cytokine receptor interaction," and "PI3K-Akt signaling pathway." These GO categories and KEGG pathways are closely related to DN progression [12,13].
As shown in Figure 3, qRT-PCR analysis was performed to confirm the degree and direction of the expression changes in 10 genes. All 10 genes assayed by qRT-PCR were found to be significantly differentially expressed in the microarray analysis between the late and early stage DN samples. As determined by qRT-PCR analysis, 10 out of the 10 selected genes demonstrated a change in expression in the same direction (i.e., up-or downregulated) (Figure 3(a)). Similarly, the direction of the change in gene expression determined by qRT-PCR analysis agreed with the directions obtained for the genes that were found to be significantly differentially expressed by microarray analysis between the early stage DN samples and the controls (5 out of 5 genes, Figure 3(b)) and between the late stage DN samples and the controls (10 out of 10 genes, Figure 3(c)).

Potential DN Drugs Predicted by the KS Statistic Algorithm.
To explore the potential drugs targeting the molecular mechanisms of DN, we used the DEGs between the late and early stages of DN. These genes were enriched for their specific contribution to nephropathy because genes that are differentially regulated in human diabetes per se, in the absence of nephropathy, were excluded by this strategy. After uploading 812 upregulated and 188 downregulated probe IDs to the CMAP database, the top 20 drug perturbations that most strongly reversed the DRGs are listed in Table 3. Among these drugs, MG-132 [14] and MG-262 [15] are proteasome inhibitors, piperlongumine inhibits PI3K/Akt/mTOR signaling [16] and NF-B activity [17], 15d-PGJ2 (15-delta prostaglandin J2) activates PPAR [18], and vorinostat and trichostatin A are histone deacetylase inhibitors (HDACIs).

Potential DN Drugs Predicted by the Matching Algorithm.
We further utilized a matching algorithm and the DEGs between early and late stage DN samples to calculate the reversing score for each drug in the CMAP database. The top 20 drug perturbations that had the highest scores are listed in Table 4. We clustered the drugs by the similarity of reversed DEGs, and the drugs with similar pharmacological characteristics were clustered together ( Figure 4). For example, HDACIs, including vorinostat, trichostatin A, and valproic acid, were clustered together. Among these drugs, piperlongumine, 15d-PGJ2, vorinostat, and trichostatin A were also identified by the KS statistic algorithm. Valproic acid and parthenolide are also histone deacetylase inhibitors [19], resveratrol inhibits cAMP phosphodiesterase [20], and LY-294002 inhibits PI3K [21].

Discussion
Current clinical strategies to treat DN focus on the intensification of glycemic control and the control of blood pressure and blood fat. Renoprotective drugs based on the molecular pathogenesis of DN are unavailable because the molecular pathogenesis of DN is complicated. Recently, high-throughput transcriptome technology has been used to explore the molecular pathogenesis of complex diseases [22,23], and drug screening methods based on transcriptome data have attracted increasing attention [4,24].
The most common kidney lesions in people with diabetes are those that affect the glomeruli, and DN is characterized by morphological and ultrastructural changes in the kidney, including expansion of the molecular matrix and loss of the charge barrier on the glomerular basement membrane [25]. The gene expression profiles of glomeruli microdissected from DN biopsy samples will provide an excellent opportunity to explore the molecular mechanisms of this complex disease and to identify potential drugs.
In this study, we obtained the whole-genome transcriptome profiles of glomeruli from DN patients and normal controls. The glomeruli contain podocytes, endothelial cells, and mesangial cells, whereas the CMAP database was built upon four types of cultured human cell lines [3]. Because the signatures of drugs are often conserved across diverse cell types [3], we can utilize the CMAP database to identify potential drugs for DN.
There were only 105 DEGs between early stage DN samples and controls, mainly involving "response to stimulus." These results indicated that there was minimal gene expression change involving molecular pathogenesis in the early stage of DN. The DEGs between the late and early stages  of DN samples were related to "extracellular matrix," "cell adhesion," "immune system process," "ECM-receptor interaction," "complement and coagulation cascades," "cytokinecytokine receptor interaction," and "PI3K-Akt signaling pathway" (Figure 2). These functional categories and pathways have been widely related to the pathogenesis of DN [13]. Coexpression network analysis and association analysis indicated that some DEGs between DN samples and controls exhibit no correlation with the progression or prognosis of DN (data not shown). These genes, which may be differentially regulated in human diabetes per se, were excluded from the comparison of the late and early stage DN samples. Therefore, it is more suitable to use the DEGs between the late and early stages to identify potential therapeutic drugs for DN treatment.
A KS statistic algorithm and a matching algorithm were applied in this study. However, using either algorithm, some drugs among the top 20 are unlikely to have renal protection on the basis of prior clinical and pharmacological knowledge. Therefore, the results based on the two different algorithms were combined to enhance the reliability of the potential therapeutic drugs, which provided a good foundation for the in vitro and in vivo studies.
Piperlongumine, 15d-PGJ2, vorinostat, and trichostatin A were identified by both algorithms. The molecular mechanisms of the 4 drugs include inhibition of NF-B activity, histone deacetylase (HDAC) activity, PI3K-Akt signaling pathway, and the activation of PPAR . The transcription factor NF-B is induced by various cell stress-associated stimuli, including growth factors, vasoactive agents, cytokines, and oxidative stress. NF-B in turn controls the regulation of genes encoding proteins involved in immune and inflammatory responses. The activation and nuclear translocation of NF-B in human DN have been demonstrated in the intrinsic cells of the kidney [26]. Upregulation of HDACs has been reported in the kidneys of patients with DN as well as in type 1 and type 2 in vivo animal models of diabetes [27]. HDACIs have anti-inflammatory and antifibrotic effects in the kidney and may prove to be a novel class of agents in the treatment of diabetic nephropathy [27]. The PI3K-Akt signaling pathway is directly related to cellular proliferation, migration, differentiation, and survival. There are many  Perturbagen name: the name given to a perturbagen; instance ID: the ID uniquely identifying each instance; score: the reversing score calculated using formula (1); reversed: number of genes in the query signature, the expression of these genes was reversed by the corresponding perturbagen; aggravated: number of genes in the query signature, the expression of these genes was aggravated by the corresponding perturbagen.
known factors that enhance the PI3K-Akt pathway, including insulin, EGF, and IGF-1. In rat mesangial cells and db/db mice, high-glucose decreased the expression of MNSOD via the PI3K-Akt signaling pathway and further aggravated oxidative stress [28]. The nuclear receptor PPAR is located in all three types of glomerular cells, with prominent expression in podocytes. PPAR agonists, which have emerged as promising candidates for treating DN, are effective in delaying and even preventing disease progression. Piperlongumine can inhibit both NF-B activation and the PI3K-Akt signaling pathway [16,17] and increase mRNA levels of PPAR 2 [29]. Piperlongumine is a main component of the root of Piper longum, a plant used by some Indian tribes to treat diabetes, digestive disorders, and obesity [30]. In streptozotocin-(STZ-) induced diabetic rats, Piper longum root aqueous extract can significantly decrease fasting blood glucose levels and protect liver and kidney function [30].
Vorinostat and trichostatin A belong to the same group of HDACIs based on their chemical structures (hydroxamic acid). They are broad inhibitors of HDAC activity and inhibit class I and class II enzymes [33]. Vorinostat is FDA approved for use against refractory cutaneous T cell lymphoma [34,35]. HDACIs have beneficial effects in diabetic nephropathy. In cultured proximal tubule cells, vorinostat treatment reduced EGFR protein and mRNA and attenuated cellular proliferation [36]. Daily treatment of diabetic rats with vorinostat for 4 weeks blunted renal growth and glomerular hypertrophy [36]. In STZ-induced diabetic mice, long-term administration of vorinostat decreased albuminuria, mesangial collagen IV deposition, and oxidative-nitrosative stress through an eNOS-dependent mechanism [37]. In STZ-induced diabetic rats, trichostatin A prevented extracellular matrix accumulation and epithelial-to-mesenchymal transition in diabetic kidney [38].
In addition to these four drugs, we have identified potential therapeutic drugs to treat DN. Proteasome inhibitors, including MG-132 and MG-262, have anti-inflammatory and antifibrotic effects [15]. MG132 alleviates kidney damage by inhibiting Smad7 ubiquitin degradation, SnoN degradation, and TGF-activation in STZ-induced DN rats [39,40]. Valproic acid, another histone deacetylase inhibitor, has beneficial effects on proteinuria, glomerulosclerosis, and renal inflammation in Adriamycin-induced nephropathic mice [41]. LY-294002, a PI3K inhibitor, prevented the quantitative and distributional changes of CD2AP induced by high-glucose and advanced glycosylation end products in mouse podocytes [42]. In db/db mice, LY-294002 decreased the levels of phosphorylated Akt and phosphorylated FoxO3a, increased the level of MnSOD expression, and further decreased oxidative stress [28]. Resveratrol, a phosphodiesterase inhibitor, improved diabetic nephropathy in several animal models of types 1 and 2 diabetes through its antioxidative effects resulting from direct radical scavenging or modulation of antioxidant enzymes [43].