Analysis of dapagliflozin-induced expression profile of long noncoding RNAs in proximal tubular epithelial cells of diabetic kidney disease

2.1 Cell culture

HK-2 cells were obtained from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). HK-2 cells were cultured in DMEM (Hyclone, Los Angeles, CA, USA) supplemented with 10% fetal bovine serum (Gibco, Thermo Fisher, Waltham, MA, USA), 100 U/mL penicillin, and 0.1 mg/mL streptomycin at 37°C in a humidified environment supplemented with 5% CO₂. HK-2 cells were cultured with Ng (5.6 mmol/L) and Hg (30 mmol/L) for 24 h, and then were cultured with Ng, Hg, and Hg plus 2.5 μmol/L DAPA (DAPA, MCE, Shanghai, China) or 10 μmol/L canagliflozin (MCE, Shanghai, China) for 48 h.

2.2 Quantitative real-time polymerase chain reaction

Total RNA was extracted from HK-2 cells using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and reverse-transcribed in 20 μL reaction volume using PrimerScript RT Master Mix and random primers according to the manufacturer's protocol (Sangon Biotech, Shanghai, China). DNA amplification was performed on a 7500 Fast Real-Time PCR System (Quant Studio 5, Thermo Fisher, MA, USA). The sequences of the primers used for Real Time Quantitative Polymerase Chain Reaction (RT-qPCR) are listed in Table 1. The expression levels of lncRNAs were normalized to those of GAPDH and calculated using the 2^−ΔΔCt method.

2.3 Microarray analysis

Total RNA was extracted from HK-2 cells and RNA quantity and quality were determined using NanoDrop 2000 (Thermo Scientific, Shanghai, China). Cell transcriptome was analyzed using the Human LncRNA Microarray V4.0, which was designed for global profiling of human lncRNAs and mRNAs. All microarray analyses were performed by KangChen Bio-tech (Shanghai, China). Sample labeling and array hybridization were performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol (Agilent, Technology, Palo Alto, CA, USA).

2.4 Differential lncRNA and mRNA expression analyses

Agilent Feature Extraction software (version 11.0.1.1) was used to analyze the array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v11.5.1 software package (Agilent Technologies). Differentially expressed lncRNAs and mRNAs among the Ng, Hg, and DAPA groups were identified using scatter plots and volcano plots using the following criteria: fold change ≤2 and P < 0.05.

2.5 GO and KEGG pathway analysis

Gene ontology term enrichments were performed to explore the functions of the differentially expressed mRNAs using the online database DAVID (http://david.abcc.ncifcrf.gov/) [16]. KEGG pathway analyses were performed to further explore the functions of the differentially expressed mRNAs according to DAVID (http://david.abcc.ncifcrf.gov/). All terms and pathways were considered significant with an enrichment score of >2.0 and P < 0.05.

2.6 Protein–protein interaction (PPI) network analysis

To construct the PPI network, we selected ten mRNAs that were significantly differentially expressed between the Ng and Hg groups, Hg and DAPA groups, using the Search Tool for the Retrieval of Interacting Genes (STRING, http://string.embl.de/) [17]. The cut-off criteria based on interaction score was a combined score of ≥0.9 (0.9 indicates the highest confidence). Subsequently, KEGG pathway analysis was performed using the DAVID online database (http://david.abcc.ncifcrf.gov/).

2.7 lncRNA–mRNA and lncRNA–miRNA–mRNA co-expression (CE) networks

The lncRNA–mRNA CE network was built using LncTar database (http://www.cuilab.cn/lnctar) [18] and expression correlations between them were measured using Spearman correlation coefficient (threshold criteria was set at Spearman correlation coefficient R > 0.8, either positive or negative). Statistical significance was set at P < 0.05. Eighteen lncRNAs were selected to generate different networks using the Cytoscape software [19].

The lncRNA–miRNA–mRNA CE network was constructed based on the correlation between the lncRNAs and miRNAs. The miRNA-binding sites were predicted using miRcode (http://www.mircode.org/) and miRanda (http://www.microrna.org/), and the miRNA–mRNA interactions were predicted using TargetScan (http://www.targetscan.org/). Based on the established CE data, lncRNA–miRNA interaction networks were mapped using Cytoscape software.

2.8 Statistical analysis

All values are expressed as mean ± standard deviation (SD) derived from at least three experiments, with similar results. Student t-test and one-way analysis of variance were performed to determine the statistical significance of differences between groups. All statistical analyses were performed using GraphPad Prism 6.0. Statistical significance was set at P < 0.05.

3.1 Differentially expressed lncRNAs and mRNAs in HK-2 cells cultured in Ng, Hg, and Hg plus DAPA

The microarray probes detected thousands of transcripts in HK-2 cells cultured with Ng, Hg, and DAPA. Scatter and volcano plots were generated to assess the variation in the expression level of lncRNAs and mRNAs between the Ng, Hg, and DAPA groups (Figure 1). A total of 7734 lncRNAs and 3427 mRNAs were identified. Among them, expression levels of 3656 and 3105 lncRNAs were up- and down-regulated, respectively, in the Hg group compared to those of Ng group, and expression levels of 1106 and 2056 mRNAs were up- and down-regulated (fold change ≥2.0, P < 0.05). Meanwhile, expression levels of 714 and 259 lncRNAs were up- and down-regulated, and those of 138 and 127 mRNAs were up- and down-regulated in the DAPA group compared to the Hg group, suggesting that these lncRNAs and mRNAs may participate in the biological processes of HK-2 cells. The values representing expression of the ten most strongly up-regulated and down-regulated lncRNAs induced by Hg with or without DAPA treatment are presented in Tables 2 and 3, respectively.

Figure 1

Differentially expressed lncRNAs and mRNAs in HK-2 cells cultured in in Ng, Hg, and Hg plus DAPA. (A, B) Scatter plots of the differentially expressed lncRNAs between Ng and Hg groups (A); DAPA and Hg groups (B). (C, D) Scatter plots of the differentially expressed mRNAs between Ng and Hg groups (C); DAPA and Hg groups (D). (E, F) Volcano plots of the differentially expressed lncRNAs between Ng and Hg groups (E); DAPA and Hg groups (F). (G, H) Volcano plots of the differentially expressed mRNAs between Ng and Hg groups (G); DAPA and Hg groups (H). (Fold-changes ≥2.0 and P < 0.05). DAPA, dapagliflozin; Hg, high glucose; Ng, normal glucose.

3.2 Signature lncRNA expression in the Hg and Hg + DAPA groups

The previous results suggested that DAPA plays a key role in the physiological function of HK-2 cells, highlighting the need to investigate the characteristics of DAPA-induced lncRNAs including their source, classification, length distribution, and chromosome distribution. The results suggested that natural antisense lncRNAs represented the largest category (26%) of all the differentially expressed lncRNAs, followed by intronic antisense lncRNAs (23%) and bidirectional lncRNAs (16%) (Figure 2A). We also found that most lncRNAs were derived from RNA-sequence sources (Figure 2B). These lncRNAs were 401–800 nt and 1201–1600 nt long (Figure 2C) and chromosome distribution analysis revealed that DAPA-associated lncRNAs were located on different chromosomes (Figure 2D).

Figure 2

Signature lncRNA expression in the Hg and Hg + DAPA groups. (A) The numbers of lncRNAs from various databases. (B) The classification of lncRNAs. (C) The distribution of differentially expressed lncRNAs based on the length of nuclear acids. (D) The number of reads mapped to each chromosome. DAPA, dapagliflozin; Hg, high glucose.

3.3 Gene Ontology (GO) and KEGG pathway enrichment analysis

Although lncRNAs do not encode proteins, they regulate the expression of neighboring and overlapping protein-coding genes, thereby performing their functions. To investigate the functions of the differentially expressed mRNAs in HK-2 cells, GO terms enrichment analysis was performed. The results showed that the mRNAs whose expression levels were up-regulated between the Ng and Hg group were associated with biological processes such as anatomical structure development, development process, single-organism developmental process, cell surface receptor signaling pathway, and cell differentiation (Figure 3A). The transcripts with down-regulated expression levels were mainly associated with the cell cycle, mitotic cell cycle, and chromosome organization (Figure 3B). For the Hg and DAPA groups, the mRNAs with up-regulated expression were mainly associated with the regulation of collagen biosynthetic process, regulation of collagen metabolic processes, and response to hepatocyte growth factor (Figure 3C); whereas mRNAs with down-regulated expression levels were mainly involved in the regulation of apoptosis, apoptotic signaling pathway, positive regulation of response to stimulus, extrinsic apoptotic signaling pathway, and vasculature development (Figure 3D). We also performed KEGG pathway enrichment analysis and found ten pathways associated with up-regulated mRNAs and ten pathways related to down-regulated mRNAs in the Hg group compared with the Ng group. The up-regulated mRNAs were mainly related to the type I diabetes mellitus signaling (Figure 3E), while the down-regulated mRNAs mainly enriched the cell cycle signaling (Figure 3F). Significantly, KEGG pathway analyses suggested that the differentially expressed genes between the Hg and DAPA groups were related to transcriptional dysregulation observed in cancer, DNA replication, and the Notch signaling pathway (Figure 3G). Genes with down-regulated expression were enriched in the pathways related to central carbon metabolism in cancer, alanine, aspartate, and glutamate metabolism, cytokine receptor interaction, ErbB signaling pathway, HIF-1: hypoxia inducible factor-1 (HIF-1) signaling pathway, MAPK: mitogen-activated protein kinase (MAPK) signaling pathway, p53 signaling pathway, P13K-Akt signaling pathway, and focal adhesion (Figure 3H). Taken together, these results suggest that the impaired expression of genes and lncRNAs leads to dysregulation of pathways that might be strongly related to diabetes mellitus or its complications, and provide key clues to explore the relationship between DKD and DAPA treatment in the future.

Figure 3

GO and KEGG pathway enrichment analysis of differentially expressed lncRNA corresponding genes. (A, B) GO annotations for biological process category of the up-regulated (A) and down-regulated (B) genes between the Ng and Hg groups. (C, D) GO annotations for biological process category of the up-regulated (C) and down-regulated (D) genes between the DAPA and Hg groups. (E, F) KEGG analysis of the up-regulated (E) and down-regulated genes (F) between the Ng and Hg groups. (G, H) KEGG analysis of the up-regulated (G) and down-regulated genes (H) between the DAPA and Hg groups. Enrichment score values were calculated as −log10 (P-Values). DAPA, dapagliflozin; GO, gene ontology; Hg, high glucose; KEGG, Kyoto Encyclopedia of Genes and Genomes; Ng, normal glucose.

3.4 PPI network construction and KEGG pathway analysis

To further investigate the cellular processes related to the differentially expressed genes between the Ng, Hg, and DAPA groups, the PPI networks of the encoded proteins were generated using STRING database. The interaction score (≥0.9) of the differentially expressed genes was used to evaluate the protein interactions, and the interacting proteins were subjected to KEGG pathway enrichment analysis. We established the PPI network (Figure 4A) between the Hg and Ng groups and found that proteins encoded by differentially expressed genes were mainly involved in metabolic pathways and peroxisome proliferator-activated receptor (PPAR) signaling pathway (Figure 4B). We also constructed the PPI network (Figure 4C) for the proteins encoded by differentially expressed genes between the DAPA and Hg groups and found that cell cycle and cancer pathways play an important role in HK-2 cells subjected to DAPA treatment (Figure 4D). Overall, these data demonstrated that most of the differentially expressed genes were mainly associated with metabolism and cell cycle signaling.

Figure 4

PPI network construction and KEGG pathway analysis. (A, B) PPI network (A) and KEGG pathway analysis (B) for differentially expressed mRNA between the Ng and Hg groups. (C, D) PPI network (C) and KEGG pathway analysis (D) for differentially expressed mRNA between the Hg and DAPA groups. Each protein is presented as a network node, and the edges represent meaningful protein–protein associations. CDB, curated data bases; CE, co-expression; DAPA: dapagliflozin; ED, experimentally determined; GC-O, gene co-occurrence; GF, gene fusions; GN, gene neighborhood; Hg: high glucose; KEGG, Kyoto Encyclopedia of Genes and Genomes; Ng: normal glucose; PH, protein homology; PPI, protein–protein interaction; TM, text mining.

3.5 CE of lncRNA-mRNA pairs and their functional prediction

To date, the functions of most lncRNAs have not been elucidated, and the prediction of most lncRNAs functions is partly based on the annotations of co-expressed mRNAs. Based on the results of the previous analyses, we selected six lncRNAs that were significantly differentially expressed in the Ng, Hg, and DAPA groups to build CNC networks according to the degree of correlation. The networks indicated that one lncRNA could regulate multiple genes in different ways, and that one gene could also be regulated by multiple lncRNAs (Figure 5A, C). The differentially expressed mRNAs were implicated in several biological processes, such as metabolic pathways, protein processing in the endoplasmic reticulum, MAPK signaling, and the Advanced glycation end products-receptor For Advanced Glycation End(AGE-RAGE) signaling pathway in diabetic complications (Figure 5B, D). The network constructed with the lncRNAs and mRNAs differentially expressed between the Hg and Ng groups showed that lncRNA ENST00000514869 (up-regulated expression) and lncRNA ENST00000607068 (down-regulated expression) were positively correlated with the transcripts encoding zinc finger protein (ZNF), C-X-C motif chemokine ligand 2 (CXCL2), VEGFA, and TGF-β (Figure 5A). Analysis of Hg and DAPA groups revealed that lncRNA ENST00000608684 (up-regulated expression) and the lncRNA T127088 (down-regulated expression) correlated with the transcripts encoding vascular endothelial growth factor A (VEGFA), C-C motif chemokine ligand 2 (CCL2), and bolA family member 2B (BOLA2B) (Figure 5C). Thus, this approach allowed us to not only deduce the lncRNAs related to the mRNAs of interest but also assess the effect of DAPA on HK-2 cells.

Figure 5

Co-expression of lncRNA–mRNA pairs and their functional prediction. (A) LncRNAs and the potential cis- and trans-regulated target genes located nearby are presented in the network between the Ng and Hg groups. (B) KEGG analysis for targeted genes of differentially expressed lncRNAs between the Ng and Hg groups. (C) LncRNAs and the potential cis- and trans-regulated target genes located nearby are presented in the network between the Hg and DAPA groups. (D) KEGG analysis for targeted genes of differentially expressed lncRNAs between the Hg and DAPA groups. Red color and blue color of diamond nodes represent up and down-regulated lncRNAs, respectively. The square nodes of yellow and green represent up- and down-regulated mRNAs, respectively. CE, co-expression; DAPA: dapagliflozin. Hg: high glucose; KEGG, Kyoto Encyclopedia of Genes and Genomes; Ng: normal glucose.

We hypothesized that the DAPA-induced lncRNAs play a critical role in ameliorating HK-2 cell injury. To evaluate this, we constructed a Venn diagram of the lncRNAs and mRNAs differentially expressed between the Hg and Ng groups and the Hg and DAPA groups. These analyses revealed 89 lncRNAs whose expression levels underwent significant changes simultaneously in the three groups (fold change ≥2.5, P < 0.05) (Figure 6A, B). Among them, we selected six lncRNAs according to their fold-changes, and constructed a CNC network using their related mRNAs. To detect significant modules in this CNC network, we used cytoHubba [20] in the Cytoscape to determine key genes based on their properties (degree, edge-percolated component, maximum neighborhood component, density of maximum neighborhood component, maximal clique centrality, and six centralities) (Figure 6C, E). We further performed KEGG enrichment analysis and found that these lncRNA-related mRNAs mainly participated in the following pathways: microRNAs in cancer pathway, MAPK signaling pathway, oxidative phosphorylation, neuroactive ligand–receptor interaction, and transcriptional mis-regulation, among other biological signaling pathways (Figure 6D, F). Next, we verified the changes in the expression profile of six lncRNAs identified in HK-2 cells treated with or without DAPA using RT-qPCR. Expression levels of NR_047558, ENST00000605056, and uc031qjg.1 were significantly increased in the Hg group and decreased after DAPA treatment. Meanwhile, expression levels of NR_027051 and ENST00000606424 were decreased in the Hg group and increased after DAPA treatment, which was consistent with the results of the microarray analysis (Figure 6G). Additionally, we explore the role of canagliflozin (CANA), another class of SGLT2 inhibitor, in HK-2 cells induced by Hg. Interestingly, the expression levels of NR_027051 and ENST00000606424 were decreased in the Hg group but there were no changes in the CANA group. Meanwhile, NR_047558, uc031qjg.1, and ENST00000505796 were significantly increased in the Hg group but decreased after CANA treatment (Figure 6H). Collectively, these results suggest that SGLT2 inhibitors may reverse HK-2 cell damage induced by Hg, and thus, provide new therapeutic strategies for DKD.

Figure 6

The CE network of lncRNA–mRNA between the Ng, Hg, and DAPA groups. (A, B) Venn diagram of overlapping differentially expressed lncRNAs reversed by DAPA. (C) LncRNA–mRNA network and the hub genes which were down-regulated in the Hg group but up-regulated in the DAPA group. (D) KEGG analysis for targeted genes which were down-regulated in the Hg group but up-regulated in the DAPA group. (E) LncRNA–mRNA network and the hub genes which were up-regulated in the Hg group but down-regulated in the DAPA group. (F) KEGG analysis for targeted genes which were up-regulated in the Hg group but down-regulated in the DAPA group. (G, H) Validation of the expression of significant lncRNAs by qRT-PCR. CANA, canagliflozin; CE, co-expression; DAPA: dapagliflozin; Hg: high glucose; KEGG, Kyoto Encyclopedia of Genes and Genomes; Ng: normal glucose; SD, standard deviation. The data displayed in the histograms are presented as the means ± SD. ^*P < 0.05 vs. Ng; ^#P < 0.05 vs. Hg.

3.6 LncRNA–miRNA–mRNA CE network construction

CeRNAs regulate other transcripts by competing for shared miRNA response elements (MREs). Thus, we constructed a ceRNA network (lncRNA–miRNA and miRNA–mRNA CE networks) for HK-2 cells to demonstrate the regulatory relationships among lncRNAs, mRNAs, and miRNAs. Four differentially expressed mRNAs among the Ng, Hg, and DAPA groups, sharing common MRE sites in several lncRNAs, were selected (Figure 7). Some miRNAs are reportedly associated with diabetes mellitus. For example, miR-27a protects HK-2 cells against Hg-induced apoptosis and autophagy [21]. In contrast, repression of miR-29a increases collagen IV production in HK-2 cells [22]. In our study, these miRNAs and mRNAs indirectly predicted the partial functions of lncRNAs. We believe that further evaluation of different RNA interactions would provide additional insights into the mechanisms underlying therapeutic effects of DAPA in DKD.

Figure 7

lncRNA–miRNA–mRNA CE network construction. Green nodes represent lncRNAs, red nodes represent miRNAs, and blue nodes represent mRNAs. CE, co-expression; lncRNA, long noncoding RNAs.