Transcriptome Analysis Reveal Candidate Genes and Pathways Responses to Lactate Dehydrogenase Inhibition (Oxamate) in Hyperglycemic Human Renal Proximal Epithelial Tubular Cells.

Diabetic kidney disease (DKD) is the leading cause of both chronic kidney disease (CKD) and end-stage renal disease (ESRD). Previous studies showed that oxamate could regulate glycemic homeostasis and impacted mitochondria respiration in a hyperglycemia-dependent manner in the rat proximal tubular cells. To explore the transcriptome gene expression profiling of kidney tissues in human renal proximal epithelial tubular cell line (HK-2), we treated HK-2 cells with high D-glucose (HG) for 7 days before the addition of 40 mM oxamate for a further 24 hours in the presence of HG in this study. Afterwards, we identified 3,884 differentially expressed (DE) genes based on adjusted P-value ≤ 0.05 and investigated gene relationships based on weighted gene co-expression network analysis (WGCNA). After qRT-PCR validations, MAP1LC3A, MAP1LC3B (P-value < 0.01) and BECN1 were found to show relatively higher expression levels in the treated groups than the control groups, while PGC1α (P-value < 0.05) showed the lower expressions. Accordingly, enrichment analyses of GO terms and KEGG pathways showed that several pathways [e.g., lysosome pathway (hsa04142) and p53 signaling pathway (hsa04115)] may be involved in the response of HK-2 cells to oxamate. Moreover, via WGCNA, we identified two modules: both the turquoise and blue modules were enriched in pathways associated with lysosome. However, the p53 signaling pathway was only found using all 3,884 DE genes. Furthermore, the key hub genes IGFBP3 (adjusted P-value = 1.34×10-75 and log2(FC) = 2.64) interacted with 6 up-regulated and 12 down-regulated DE genes in the network that were enriched in the p53 signaling pathway. This is the first study reporting co-expression patterns of a gene network after lactate dehydrogenase inhibition in HK-2 cells. Our results may contribute to our understanding of the underlying molecular mechanism of in vitro reprogramming under hyperglycemic stress that orchestrates the survival and functions of HK-2 cells.

Introduction

Diabetic kidney disease (DKD) is the leading cause of chronic kidney disease (CKD) and end-stage renal disease (ESRD) worldwide. Globally, due to the increasing incidence of diabetes and aging, DKD continues to increase stably and imposes the highest burden and the strongest correlation with mortality in patients with diabetes (1, 2). However, the pathological pathways that precipitate the development and/or progression of DKD remain to be fully elucidated. In the kidney, proximal tubular epithelial cells (PTECs), which represent ~90% of the outer kidney cortex and ~50% of the entire kidney mass, present more mitochondria than other renal cell types, because they reabsorb 80% of the filtrate that passes through the glomerulus, including large amounts of water, electrolytes, glucose, and other small molecules (3). Increased glucose reabsorption of PTECs under hyperglycemic stress triggers the initial period of hyperfiltration and tubule-glomerular feedback, exacerbates hypoxia and drives a poorly understood dominos of renal injuries including atrophy of proximal tubular, atubular glomeruli and tubulointerstitial fibrosis that associate with culminating in renal failure and further contributes to the emergence of proximal tubulocentric research (4, 5).

Recently, both in vivo and in vitro tracer studies have demonstrated an increased glycolysis-driven lactate production, or pyruvate-to-lactate production, in early renal changes associated with diabetes in rat kidney and PTECs (6). Lactic acid is the sink for three-carbon compounds, and the hub of glycolysis, mitochondrial energy metabolism and gluconeogenesis in PTECs (7). Pyruvate-to-lactate process reflects the redox buffer of NADH/NAD+ ratio across cells (7). Exploration of the pathophysiology involved mechanisms of lactate metabolic adaptation to reshape the field of metabolic reprogramming in hyperglycemic PTECs could offer a comprehensive update on strategies targeting renal tubules.

Oxamate is an isosteric and isoelectronic analogue of pyruvate, i.e., a lactate dehydrogenase (LDH) inhibitor wildly used in tumor cells or in tumor cell energy metabolism and apoptosis research. Recent studies found that oxamate could regulate glycemic homeostasis from both central and peripheral tissues (8). The previous results showed that oxamate, in rat proximal tubular cells line NRK-52E, impacted mitochondria respiration in a hyperglycemia-dependent manner (6). The flexibility of mitochondrial respiration adapting to nutrient stress is dependent on mitophagy and mitochondrial biogenesis; expressions of the mitophagy or mitochondrial biogenesis related genes and their encoded proteins might be involved, e.g., BCL2, BECN1, MAP1LC3A, MAP1LC3B and PGC1α. Therefore, we aim to explore the transcriptomics of human PTECs cultured with high D-glucose (HG) combined with oxamate in vitro to uncover the unexplored in vitro reprogramming under hyperglycemic stress and the underlying mechanism for this reprogramming that orchestrates the survival and functions of HK-2 cells (human cortex proximal tubular immortalized cell line), and the therapeutic potentials of targeting the pathways to reprogram the DKD.

The weighted gene co-expression network analysis (WGCNA) is widely used for gene co-expression networks that are constructed by genes with the significant co-expression relationships (9). It shows the co-expressed genes with the similar expression patterns across samples that are controlled by the same transcriptional regulatory programs (10, 11), which has been used in the integrated meta- and microRNA- analysis (12, 13); therefore it can be potentially used in the co-expression analysis of HK-2 cells after oxamate treatment. In this study, three independent replicates of human HK-2 cells with oxamate treatment and three independent replicates of control groups were used. The main objectives of our study are as follows:

1. Conduct a genome-wide transcriptome study on HK-2 cells in the absence and presence of oxamate treatment under HG condition to reveal differentially expressed (DE) genes.

2. Co-expressions analysis for HK-2 cells and construct the weighted gene co-expression networks.

3. Perform the enrichment analysis and identify the potential biological functions for the DE genes on HK-2 cells.

4. Reveal the candidate biomarkers to affect protein levels on HK-2 cells in the network.

Materials and Methods

HK-2 Cell Culture

HK-2 cells, the human renal proximal tubular epithelial cell line, were obtained from Beijing Beina Chuanglian Biotechnology Research Institute. HK-2 cells were cultured in Dulbeccos modified eagles medium (DMEM) with the low glucose (Sigma-Aldrich Inc, St. Louis, USA) supplemented with 10% fetal bovine serum (Sigma-Aldrich, St. Louis, USA), 10 units/ml penicillin and 10 mg/ml streptomycin. They were maintained in the continuous culture at 37°C in a humidified atmosphere (5% CO2) in an incubator. Growth medium was changed every second day, and cells were sub-cultured until further measurements at 80% colony confluency.

Oxamate could induce apoptosis in a dose-dependent manner (14); therefore, the xCELLigence real-time cell analyzer (RTCA) S16 system (Agilent Technologies, Santa Clara, USA) was used to monitor the dynamic real-time cell viability to explore a tolerable dose of oxamate. We exposed HK-2 cells to 5.5 mM low D-glucose (LG) medium or 25 mM HG medium with different concentrations of oxamate (i.e., 0 mM, 20 mM, 40 mM, 60 mM, 80 mM and 100 mM) for 24-hour intervention. One day prior cells were seeded at a density of 1×103 cells per well, respectively, in the 16-well assay plate in medium containing 5.5 mM (LG) or 25 mM (HG) D-glucose. The medium was replaced by 5.5 mM LG medium or 25 mM HG medium with different concentrations of oxamate (i.e., 0 mM, 20 mM, 40 mM, 60 mM, 80 mM and 100 mM). Accordingly, we normalized the starting point 2 hours later after introducing the experimental variables and allowed the cells to grow for the oxamate intervening 24 hours via monitoring every 30 minutes. Here, each experiment was carried out in triplicate.

Western Blot Analysis

The radioimmunoprecipitation assay lysis buffer (RIPA buffer) and phenylmethylsulfnoyl fluoride (PMSF) (Servicebio, Wuhan, China) were used to extract proteins at 4°C from the HK-2 cells treated with HG (25 mM) combined with or without different concentrations of oxamate (e.g., 0 mM, 20 mM, 40 mM, 80 mM) for 24 hours. We used BCA protein assay kit (Servicebio, Wuhan, China) to measure and adjust the protein concentrations. The proteins were denatured with 5 × SDS loading buffer (Servicebio, Wuhan, China) at 98°C for 15 min. Then, the prepared proteins were separated by SDS-polyacrylamide gel electrophoresis and concreted with 5% concentrated gel. Polyvinylidene fluoride membranes were used to transfer the proteins at 25 Voltage overnight. We used 5% skim milk solution to block the protein, incubated the membranes with the specific primary antibodies and the secondary antibodies, and visualized the membranes by ECL Plus reagents (Servicebio, Wuhan, China). The primary antibodies [Anti-PGC1α (ab54481), Anti-CASP3 (9665), Anti-BCL2 (ab59348), Anti-BAX (GB11690), Anti-BECN1 (bs-1353R), Anti-MAP1LC3 (14600-1-ap) and Anti-β-actin (GB12001)] and all secondary antibodies were purchased from abcom (Cambridge, MA, USA), Proteintech (Wuhan, China), BIOSS (Boston, Massachusetts, USA) and Servicebio (Wuhan, China), respectively. All primary antibodies were included at the ratio of 1: 3000. AlphaEaseFC software (Alpha Innotech, Miami, USA) was used to calculate the grayscale value of the proteins, and the protein images were processed by the Adobe software.

RNA Sequencing, Read Quality Control and Alignment to Reference Genome

After HK-2 cell culture exposed to different concentrations of oxamate and western blot analysis, 40 mM was chosen as the appropriate oxamate treatment concentration. Thus, six replicates of cells were firstly exposed to HG treatments for seven days, then three replicates of them were exposed to 40 mM oxamate (case group), while the other three replicates were still under the same HG condition but 40 mM oxamate expose (control group). Unfortunately, one replicate of case group failed in RNA sample preparation, so five replicates of cells were finally chosen for RNA sequencing in this study. A total amount of 2 μg RNA per sample was used for RNA sample preparation. The sequencing libraries were generated using NEBNext Library Prep Kit (NEB, USA) for Illumina following the manufacturer’s recommendations. Then, they were sequenced on the Illumina Hiseq X ten platform to generate the 150 bp paired-end reads.

In order to guarantee the data quality, raw data was filtered for the clean reads by removing the contaminated reads for adapters (> 5 bp adapter sequences), low quality reads (Phred quality value <= 19 more than 15%) and reads with Ns (Ns > 5%). Afterwards, clean reads were aligned to the human reference genome GRCh38.p13 (Genome Reference Consortium Human Build 38) using HISAT2 software (version 2.1.0) (15) that uses a modified BWT algorithm to convert reference genome to index for faster speed and fewer resources.

FPKM Calculation and Differentially Expressed Gene Analysis

We used HTSeq software (version 0.6.0) (16) to calculate the read counts for each gene. Then, the fragments per kilobase million mapped reads (FPKM) () was calculated to estimate the expression levels of genes in each sample, following the formula

where F is the number of fragments in a certain sample that is assigned to a certain gene, N is the total number of mapped reads in the certain sample and L is the length of the certain gene. FPKM could eliminate the effect of sequencing depth and gene length on gene expression levels, so it was most commonly used in the previous studies (17, 18).

DE gene analysis between two groups was performed using R package DESeq2 (version 1.30.0) that was designed for differential gene expression analysis between two groups with biological replicate samples under the theoretical basis following the hypothesis of negative binomial distribution (19–21). DESeq2 estimates the expression level of each gene for each sample by the linear regression model to consider the genes of the same expression levels that share the similarity deviations or own the expression characteristics. It calculates the P-value by Wald test and corrects the multiple testings by Benjamini–Hochberg (BH) procedure to achieve the adjusted P-value. Here, genes with thresholds of adjusted P-value ≤ 0.05 are identified as DE genes. In addition, fold change (FC) of each gene was calculated based on the averaged FPKM value for case group and control group; thus, the up‐regulated and down‐regulated genes were defined when the log2(FC) values were positive and negative, respectively.

All genes were visualized in a volcano plot by the R function plot based on FCs and adjusted P-values. The DE genes were clustered in a heat map by the R function heatmap based on the transforms of log10(FPKM + 1) values, where the transforms were scaled for each gene and clustered using the default complete hierarchical clustering method.

Differentially Expressed Gene Validation by qRT-PCR Experiment

The approach of qualitative reverse transcription polymerase chain reaction (qRT-PCR) was applied to validate five genes that are B-cell lymphoma 2 apoptosis regulator (BCL2), beclin1 (BECN1), microtubule-associated proteins 1A/1B light chain 3A (MAP1LC3A), microtubule-associated proteins 1A/1B light chain 3B (MAP1LC3B) and peroxisome proliferator-activated receptor-g coactivator 1α (PGC1α), using the re-cultured and re-treated cells. The selected genes might be involved in mitophagy and mitochondrial biogenesis for the flexibility of mitochondrial respiration adapting to nutrient stress. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was chosen as the internal gene for qRT-PCR experiments and 2×SYBR Green qPCR Master Mix (High ROX) (Servicebio, Wuhan, China) was used to perform qRT-PCR experiments. The method of 2−ΔΔCt was used to calculate the relative gene expression levels. All the primers for qRT-PCR experiments are shown in the .

Table 1

The primers used in the qRT-PCR experiments.

Gene name	Forward Primer	Reverse Primer	Length
BCL2	GAGGAAGTCCAATGTCCAGCC	GCATCCCAGCCTCCGTTATC	157
BECN1	ACATGAGCGAGTTGGTCAAGATC	CCCAGTGACCTTCAGTCTTCG	162
MAP1LC3A	ACATGAGCGAGTTGGTCAAGATC	ACATGAGCGAGTTGGTCAAGATC	141
MAP1LC3B	GTTGGCACAAACGCAGGGTA	ACACTGCTGCTTTCCGTAACAA	305
PGC1α	ACCCAGAACCATGCAAATCACA	ACCCAGAACCATGCAAATCACA	166
GAPDH	ACCCAGAACCATGCAAATCACA	ACCCAGAACCATGCAAATCACA	168

Gene Co-Expression Network of Differentially Expressed Genes and Their Associations With Oxamate Treatment

The gene co-expression network was constructed by R package WGCNA (9) for the similarity measurement between the gene expression profiles by Pearson correlation coefficients of matrix. It transformed the similarity matrix into an adjacency matrix (A) raised to a β exponent (soft threshold) based on the free-scale topology model. A total of 20,377 genes were filtered out of 29,483 genes based on the median absolute deviation (MAD) ≥ 0.01 in this study. The β power parameter (soft threshold) was equal to 8 when the R2 of the free-scale topology was equal to 0.8 ().

In the network construction, we set the minimum module size equal to 30 for detection with the unsigned TOM type. Additionally, we set the dendrogram cut height for modules merging at 0.25 by clustering module eigengenes using the dissimilarity, so the modules whose eigengenes are highly correlated above 0.75 would be merged on each branch.

Module association between the module eigengenes (MEs) (i.e., the first principle component to represent the overall expression level of the module) and the oxamate treatment status (i.e., 0, 0, 0, 1, 1 for the control and case groups, respectively) was calculated for the relevant module identifications. We calculated the module significance (MS) to evaluate the correlations. Normally, modules with the highest MS score are considered as the key modules (9), thus MS genes in the association analysis (P-value < 0.1) were assigned for functional enrichment analysis.

Gene Ontology and Pathway Enrichment Analysis

Both Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis were conducted to investigate the important enrichments for the associations between the identified DE genes and gene-related biological functions. We used R package clusterProfiler (version 3.6) (22) to test the statistical enrichments of GO terms and KEGG pathways using DE genes (adjusted P-value < 0.05) and DE genes in the key modules with high MS scores, where the GO and pathway enrichments were both based on over-representation analysis (ORA).

Protein-Protein Interactive Analysis of Differentially Expressed Genes

In this study, protein-protein interactions (PPIs) were predicted for the top DE genes based on the STRING database (https://string-db.org). Blastx software (version 2.2.28) (https://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to align the target gene sequences to the selected reference protein sequences using the default settings and the protein networks were built according to the known interactions of human species. Afterwards, we used Cytoscape software (version 3.5.1) to visualize the networks of PPIs (23).

Results

HK-2 Cell Culture and Tolerable Dose Exploration

Comparing with the cells cultured without oxamate (0 mM concentration), the cells treated with 20 mM or 40 mM oxamate retained the proliferation and viability either under LG or HG conditions. However, the cell viability progressively decreased after treating with even higher concentrations of oxamate (i.e., 60 mM, 80mM and 100mM) for 24 hours in both LG medium () and HG medium (). It was suggested that oxamate affected the cell viability in a dose-dependent manner in either LG or HG culture condition within 24 hours ().

Figure 1

(A) Cell viabilities in low D-glucose (LG) with different concentrations of oxamate (LGOXA-0mM, LGOXA-20mM, LGOXA-40mM, LGOXA-60mM, LGOXA-80mM and LGOXA-100mM) for 24 hours. (B) Cell viabilities in high D-glucose (HG) with different concentrations of oxamate (HGOXA-0mM, HGOXA-20mM, HGOXA-40mM, HGOXA-60mM, HGOXA-80mM and HGOXA-100mM) for 24 hours.

Further western blot assay was used to detect the protein expressions of PGC1α, CASP3, BCL2, BAX, BECN1 and MAP1LC3 (). The results showed that oxamate significantly decreased the protein expression levels of PGC1α and BAX under the culture of LGOXA-40mM comparing with the cells cultured without oxamate (0 mM concentration), whereas their expressions remained stable under the treatment of HGOXA-40mM (). Culturing HK-2 cells in 25 mM HG medium with 40 mM oxamate for 24 hours is suggested to prompt a diabetic state affiliated metabolic reprogramming without disrupting cell growth ability.

Figure 2

The protein expression levels of peroxisome proliferator-activated receptor-g coactivator 1α (PGC1α), caspase 3 (CASP3), B-cell lymphoma 2 apoptosis regulator (BCL2), BCL2 associated X apoptosis regulator (BAX), beclin1 (BECN1) and microtubule-associated proteins 1A/1B light chain 3 (MAP1LC3) in HK-2 cells treated with low D-glucose (LG) with different concentrations of oxamate (LGOXA-0mM, LGOXA-20mM, LGOXA-40mM and LGOXA-80mM) and high D-glucose (HG) with different concentrations of oxamate (HGOXA-0mM, HGOXA-20mM, HGOXA-40mM and HGOXA-80mM) for 24 hours. (A) Western blot gels. (B) Relative grey values. * indicates P-value < 0.05 after student’s t-test.

Statistics of Alignment to Human Reference Genome

On average, 46,017,884 (97.14%) clean reads were retained from 47,372,938 raw reads after quality control by removing 251,576 low-quality reads, 110,726 ploy-N reads and 992,752 adapter polluted reads. Afterwards, 44,710,680 of them were uniquely mapped to the human reference genome with mapping rate of 97.16% ().

Table 2

Statistics of quality control and alignment to human reference genome (GRCh38.p13).

Sample	Raw read number	Low-quality read number (%)	Ploy-N read number (%)	Adapter polluted read number (%)	Clean read number (%)	Mapped read (%)
HG1	48,909,654	219,180 (0.45)	70,820 (0.14)	732,722 (1.50)	47,886,932 (97.91)	46,597,334 (97.31)
HG2	47,359,384	337,124 (0.71)	83,178 (0.18)	675,864 (1.43)	46,263,218 (97.69)	44,866,400 (96.98)
HG3	47,542,838	192,338 (0.41)	261,400 (0.55)	2,126,872 (4.47)	44,962,228 (94.57)	43,678,133 (97.14)
HGOXA2	46,419,566	243,330 (0.52)	60,204 (0.13)	698,242 (1.50)	45,417,790 (97.84)	44,197,980 (97.31)
HGOXA3	46,633,246	265,908 (0.57)	78,028 (0.17)	730,058 (1.57)	45,559,252 (97.70)	44,213,554 (97.05)
Mean	47,372,938	251,576 (0.53)	110,726 (0.23)	992,752 (2.09)	46,017,884 (97.14)	44,710,680 (97.16)
SD	980,639	55,115 (0.12)	84,669 (0.18)	634,430 (1.33)	1,144,377 (1.44)	1,135,766 (0.15)

% indicates the percentage. HGOXA2 and HGOXA3 indicate the case group that were treated with 40 mM oxamate in 25 mM HG medium. HG1, HG2 and HG3 indicate the control groups that were only treated with 25 mM HG medium without oxamate.

Differentially Expressed Gene in Up/Down Regulation Status and Their Experimental Validations

In this study, we identified 3,884 DE genes (adjusted P-value ≤ 0.05) including 1,664 up-regulated DE genes and 2,220 down-regulated DE genes. All the details of 3,884 DE genes with gene names, read counts for each gene along different samples, log2(FC) values, P-values, adjusted P-values and regulation status are listed in . The top ten up-regulated DE genes and down-regulated DE genes were shown in . We found that the most up-regulated DE gene was SAT1 (adjusted P-value = 5.57×10-83 and log2(FC) = 1.75), while the most down-regulated DE gene was DYNC1H1 (adjusted P-value = 3.66×10-41 and log2(FC) = -1.14). In addition, the log2(FC) values varied from -1.68 to 3.16 among the top ten up-regulated DE genes and ten down-regulated DE genes ().

Table 3

Top ten up-regulated and ten down-regulated differentially expressed genes.

Ensembl gene ID	Gene name	Full description	Log2(FC)	Adjusted P-value	Regulation status
ENSG00000130066	SAT1	Spermidine/spermine N1-acetyltransferase 1	1.75	5.57×10-83	Up-regulated
ENSG00000146674	IGFBP3	Insulin like growth factor binding protein 3	2.64	1.34×10-75	Up-regulated
ENSG00000165272	AQP3	Aquaporin 3 (Gill blood group)	1.97	1.40×10-59	Up-regulated
ENSG00000126709	IFI6	Interferon alpha inducible protein 6	2.04	1.15×10-53	Up-regulated
ENSG00000163220	S100A9	S100 calcium binding protein A9	3.16	5.50×10-51	Up-regulated
ENSG00000187134	AKR1C1	Aldo-keto reductase family 1 member C1	2.12	8.45×10-50	Up-regulated
ENSG00000019186	CYP24A1	Cytochrome P450 family 24 subfamily A member 1	1.57	5.53×10-49	Up-regulated
ENSG00000169715	MT1E	Metallothionein 1E	1.42	1.07×10-48	Up-regulated
ENSG00000198886	MT-ND4	Mitochondrially encoded NADH dehydrogenase 4	1.04	2.15×10-44	Up-regulated
ENSG00000151632	AKR1C2	Aldo-keto reductase family 1 member C2	2.82	9.21×10-43	Up-regulated
ENSG00000197102	DYNC1H1	Dynein cytoplasmic 1 heavy chain 1	-1.14	3.66×10-41	Down-regulated
ENSG00000167548	KMT2D	Lysine methyltransferase 2D	-1.68	9.12×10-38	Down-regulated
ENSG00000109971	HSPA8	Heat shock protein family A (Hsp70) member 8	-1.03	1.58×10-37	Down-regulated
ENSG00000148773	MKI67	Marker of proliferation Ki-67	-1.38	6.53×10-37	Down-regulated
ENSG00000243156	MICAL3	Microtubule associated monooxygenase, calponin and LIM domain containing 3	-1.39	1.46×10-35	Down-regulated
ENSG00000149503	INCENP	Inner centromere protein	-1.36	5.14×10-33	Down-regulated
ENSG00000178209	PLEC	Plectin	-0.92	6.87×10-33	Down-regulated
ENSG00000127481	UBR4	Ubiquitin protein ligase E3 component n-recognin 4	-1.10	1.11×10-30	Down-regulated
ENSG00000113810	SMC4	Structural maintenance of chromosomes 4	-1.43	3.91×10-29	Down-regulated
ENSG00000130723	PRRC2B	Proline rich coiled-coil 2B	-1.10	1.44×10-27	Down-regulated

displayed the obvious division between up-regulated DE genes and down-regulated DE genes based on log2(FC) values. The range of log2(FC) values reached to -5 for down-regulated DE genes and 5 for up-regulated DE genes (). After the clustering of transformed FPKM values, we found that the DE genes with high expression levels (red color) in two samples (HGOXA2 and HGOXA3) from case group were clustered together, and vice versa (low expression levels in yellow color); thus, the gene expression levels showed an apparent partition between two treatment groups (). After the qRT-PCR experiment validations, we found that MAP1LC3A, MAP1LC3B (P-value < 0.01) and BECN1 had relatively higher expression levels in the 40 mM oxamate treated groups than the control groups, while PGC1α (P-value < 0.05) showed the lower expressions ().

Figure 3

Differentially expressed (DE) analysis results. (A) Volcano plots for all genes based on fold change (FC) values and adjusted P-values. Note: Genes with thresholds of adjusted P-value ≤ 0.05 are identified as DE genes. (B) Heatmap plots for all DE genes. Note: The color from yellow to red indicates the gene expression levels from low to high after the transforms of log10(FPKM + 1) values. HGOXA2 and HGOXA3 indicate the case group that were treated with 40 mM oxamate in 25 mM HG medium. HG1, HG2 and HG3 indicate the control groups that were only treated with 25 mM HG medium without oxamate. (C) Relative expression levels of the selected genes validated by qRT-PCR experiments. Note: HGOXA-40mM indicates the case group that was treated with 40 mM oxamate in 25 mM HG medium. ** and * indicate the P-value < 0.01 and P-value < 0.05, respectively, after the student’s t-test.

Co-Expression Network of Differentially Expressed Genes and the Key Modules Associated With Oxamate Treatment

Using 3,884 DE genes for the sample clustering, we found the samples with oxamate treatment were clustered together (). The eigengene adjacency heatmap indicated that these modules of DE genes could be clustered further together into groups, where treatment status was grouped with blue pattern (). The weighted DE gene network construction was visualized in the topological overlap matrix (TOM) clusters that showed a high level of overlap densities among the two clusters (). Obviously, DE genes were grouped into 2 modules that had similar co-expressions using the average linkage hierarchical clustering algorithm (), where 3,087 DE genes were grouped into turquoise module as the key module, followed by 797 DE genes into blue module. The module-trait relationship results revealed that oxamate treatment had high correlations with the turquoise module (correlation coefficient = -0.97, P-value < 0.01) and blue module (correlation coefficient = 0.96, P-value < 0.01) ().

Figure 4

Co-expression network results of differentially expressed (DE) genes. (A) Sample cluster heatmap with oxamate treatment for DE genes. (B) Eigengene dendrogram and adjacency heatmap of different co-expression modules for DE genes. (C) Weighted DE gene network heatmap of the co-expression interactions based on the topological overlap matrix (TOM) dissimilarity. Note: The gene dendrogram and module assignment are shown along the left side and the top, where the axe colors indicate the different modules. The color intensity inside the heatmap represents the overlap degree, where light color represents low overlap and darker red color represents higher overlap. (D) Module-trait relationship heatmap between oxamate treatment and control groups for DE genes. Each row indicates module eigengenes (the first principle component of the module) with the correlation coefficients (P-values in the brackets), where red color represents positive correlation and blue color represents negative correlation.

Significant Enrichments of GO Term and Pathway

Based on 3,884 DE genes, a total of 1,556 significant GO terms were derived for three domains including 998 biological process (BP), 336 cellular component (CC) and 222 molecular function (MF) ontologies (). Based on the regulation status of DE genes, 430 up-regulated and 1126 down-regulated significant GO terms (), and 22 up-regulated and 58 down-regulated significant pathways () were achieved respectively.

Figure 5

Number of the significant up-regulated and down-regulated (A) GO terms and (B) pathways.

We found that the three most significant GO terms were chromosome segregation (GO:0007059, adjusted P-value = 2.95×10-30) with 115 DE genes, mitotic nuclear division (GO:0140014, adjusted P-value = 2.33×10-28) with 104 DE genes and DNA replication (GO:0006260, adjusted P-value = 2.33×10-28) with 102 DE genes in the BP ontology of down-regulated category (). Similarly, the three most significant pathways were also in the down-regulated category including RNA transport (hsa03013, adjusted P-value = 1.43×10-9) with 55 DE genes, cell cycle (hsa04110, adjusted P-value = 2.04×10-9) with 42 DE genes and spliceosome (hsa03040, adjusted P-value = 1.27×10-6) with 42 DE genes (). Based on 3,087 DE genes in the turquoise module after WGCNA analysis, the three most significant GO terms and pathways () were similar to the enrichment results based on 3,884 DE genes, except the longevity regulating pathway - multiple species (hsa04213, adjusted P-value = 3.46×10-6) with 22 DE genes () only in the turquoise module. However, enrichment results of 797 DE genes in the blue module were quite different from those in the turquoise module (), where only RNA transport was still one of the three most significant pathways ().

Figure 6

Significant (A) GO terms and (B) pathways for differentially expressed (DE) genes (n = 3,884) in the down-regulated and up-regulated categories.

Figure 7

Significant (A) GO terms and (B) pathways for differentially expressed (DE) genes in the turquoise module (n = 3,087) in the down-regulated and up-regulated categories.

Figure 8

Significant (A) GO terms and (B) pathways for differentially expressed (DE) genes in the blue module (n = 797) in the down-regulated and up-regulated categories.

In the up-regulated category, protein localization to endoplasmic reticulum (GO:0070972, adjusted P-value = 6.26×10-18, DE gene = 51), neutrophil mediated immunity (GO:0002446, adjusted P-value = 1.14×10-17, DE gene = 102) and cotranslational protein targeting to membrane (GO:0006613, adjusted P-value = 1.36×10-17, DE gene = 43) in the BP ontology were the most significant GO terms based on 3,884 DE genes (). Meanwhile, lysosome (hsa04142, adjusted P-value = 2.30×10-13, DE gene = 45), ribosome (hsa03010, adjusted P-value = 3.54×10-11, DE gene = 47) and coronavirus disease - COVID-19 (hsa05171, adjusted P-value = 5.07×10-8, DE gene = 53) were the most significant pathways (). Moreover, we found that p53 signaling pathway with 16 DE genes (hsa04115, adjusted P-value = 1.73×10-2) was also involved in the top significant pathways (). Based on 3,087 DE genes in the turquoise module, the three most significant GO terms and pathways () were similar to the enrichment results based on 3,884 DE genes, but enrichment results of 797 DE genes in the blue module were still different ().

Protein-Protein Interaction Networks of Differentially Expressed Genes

The PPI network results were visualized using the top eight up-regulated DE genes (SAT1, IGFBP3, AQP3, IFI6, S100A9, AKR1C1, CYP24A1 and AKR1C2) and the top six down-regulated DE genes (DYNC1H1, KMT2D, INCENP, SMC4, MKI67 and PRRC2B). It showed that most DE genes had the PPIs with at least two other DE genes (). For example, the up-regulated IGFBP3 had the PPIs with the other 6 up-regulated DE genes (CSF2, MT1F, RAD54B, MKI67, CCNE2 and DDIT3) and 12 down-regulated DE genes (STEAP1, ANGPTL4, FOS, PRR15, HTRA1, LCN2, C1R, MUC1, SAA1, NDUFA4L2, MX1 and C1S) (). The down-regulated KMT2D had the PPIs with the other 13 up-regulated DE genes (RP11-313J2.1, ZNF699, CCNE2, ZNF488, ZNF770, RAD548, SZT2, EP400, MICAL3, TRRAP, DYNC1H1, ZIC4 and PHC1P1) and 6 down-regulated DE genes (HIST2H2AA3, ZSCAN16, HIST1H1C, ZNF880, HIST1H2AC and IFIT1) ().

Figure 9

Protein-protein interactions (PPIs) network of 14 top DE genes including (A) 8 top up-regulated DE genes and (B) 6 top down-regulated DE genes.

Discussion

Proximal tubules present a very high density of mitochondria required for energy consumption (24). Mitochondria occupy about 16.3% relative volume of the human cross-sectioned proximal convoluted tubules (S1 and S2 combined) (25). Mitochondrial dysfunction, with persistent energy depletion and deficiency of oxygen utilization efficiency, is the main driver in the progression of hyperglycemic proximal tubules (24). The process of selective removal of dysfunctional and depolarized mitochondria, known as mitophagy, interfaces and coordinates with mitochondrial biogenesis to maintain mitochondrial homeostasis. Previous studies showed that oxamate could impact mitochondria respiration in a hyperglycemia-dependent manner in the rat proximal tubular cells (6). Here, our results further suggested that under diabetic state, this metabolic reprogramming caused by oxamate was related to its impact on the resilience in both mitophagy and mitochondrial biogenesis, in which MAP1LC3 (MAP1LC3A and MAP1LC3B), BECN1 and PGC1α were involved.

Oxamate Features on Cell Viability/Survival and Mitochondrial Biogenesis

In mammals, lactate and pyruvate as two sinks for three-carbon compounds, enable the crosstalk between glycolytic flux and mitochondrial respiratory function, and serve as a circulating redox buffer that equilibrates the NADH/NAD+ ratio across cells and tissues together (7). For a long time, oxamate had been considered as a classic LDH enzyme inhibitor. Later, the kinetic modeling results revealed the multiple sites targeting of oxamate (14). Besides LDH inhibition function, oxamate strongly blocked pryruvate kinase (PYK) and enolase (ENO) activities, and slightly inhibited hexosephosphate isomerase (HPI), aldolase (ALD) and glucose-6-phosphate dehydrogenas (Glc6PDH) activities (14). Oxamate effect on cell fate in a dose-dependent manner depended on its cumulative influences on multiple combinations of metabolic enzymes and metabolites. Heretofore, we found that treating hyperglycemic NRK-52E with a relatively high concentration of oxamate (e.g., 100mM) significantly reduced the cell proliferations and survivals (6). To avoid the bias related to cell death caused by excessive drug concentration, we explored the dose-dependent manner of oxamate on cell viability/survival and mitochondrial biogenesis. Finally, cells treated with 40mM oxamate under HG condition that retained cell viability/survival is suggested to abrogate the PGC1α inhibiting properties caused by oxamate (, ); thus, 40mM oxamate is considered as a plausible concentration for our transcriptomic study.

Differentially Expressed and Co-Expressed Genes After Oxamate Treatment

Oxamate was found to downregulate the expression of IL6 and upregulate the expression of PDE3B in the skeletal muscle of db/db mice (26). However, IL6 (adjusted P-value = 2.43×10-2 and log2(FC) = 0.54) was identified as one up-regulated DE gene, whereas PDE3B (adjusted P-value = 1.80×10-2 and log2(FC) = -0.52) as one down-regulated DE gene in our study (). As a B-cell differentiation factor, interleukin 6 (IL6) is a multifunctional cytokine to regulate the hematopoiesis, immune and acute-phase responses, and inflammations of interleukin 1 (IL1), tumor necrosis factor alpha (TNF-α), and lipopolysaccharide (LPS) as the stimuli (27–30). PDE3B with PDE3A from the PDE3 family hydrolyze cAMP and cGMP, and its isoforms are in higher expressions than PDE3A in tissues that regulate energy homeostasis, including adipose tissue, liver and pancreatic β cells (31, 32). Ahmad et al. (33) demonstrated that the role of PDE3B regulated NLRP3 inflammasome in modulating inflammatory responses to contribute to a reduced inflammatory state in adipose tissue. In addition, our study found SAT1 (Spermidine/spermine N1-acetyltransferase 1) was the most up-regulated DE gene (adjusted P-value = 5.57×10-83 and log2(FC) = 1.75) (). It is reported that SAT1 is an important transporter involved in sulfate homeostasis as an anion exchanger and regulates the expression of the hepatocellular sulfate by glyoxylate that could be a metabolic link between liver and kidney (34).

The HIF inhibiting properties were removed using oxamate in the hyperglycemic rat proximal tubular cells (6). Studies have been reported that oxamate could induce autophagy via downregulation HIF-1α through inhibiting the Akt-mTOR signaling pathway in cancer cells (35, 36). BECN1 is a haploinsufficient tumor suppressor gene as a core component of the class III phosphatidylinositol 3-kinase that is involved in autophagosome formation and vesicular trafficking (37). It has been suggested that BECN1 was required for mitophagy completion, and suppressed mitophagy by BECN1 deficiency could cause aberrant mitochondria quality control in adipocyte to cause lipodystrophy and metabolic dysregulation (38). The qRT-PCR experiment results of MAP1LC3A (adjusted P-value = 1.30×10-3 and log2(FC) = 0.71) and MAP1LC3B (adjusted P-value = 8.47×10-4 and log2(FC) = 0.35) showed higher expression levels in the 40 mM oxamate treated groups than the non-treated groups after validations (P-value < 0.01 after student’s t-test) (), which is consistent with the transcriptome results in the up-regulated status (). Matboli et al. (2017) (39) revealed MAP1LC3A with the high diagnostic power for detection of DKD, as the urinary expression level of MAP1LC3A was significantly lower in DKD than the control group. The mRNA expressions of MAP1LC3B were also down-regulated in the DKD group (40). MAP1LC3 is essential in autophagosome mature and selective engulfment of damaged mitochondria fusing with lysosomes and processing mitophagy (41, 42). Therefore, oxamate could be involved in reprogramming the DKD via the influence of some potential therapeutic pathways, where the DE MAP1LC3A and MAP1LC3B are involved.

Significant Pathways and Protein-Protein Interaction Networks in Human Proximal Tubular Cells

Lysosome (hsa04142) was identified as one of the common significant pathways based on 3,884 DE genes, 3,087 DE genes in the turquoise module and 797 DE genes in the blue module (, , ), where ACP2, AGA, ARSA, ASAH1, ATP6AP1, CD63, CTNS, CTSA, CTSB, CTSC, CTSD, CTSF, CTSH, CTSL, CTSS, CTSZ, DNASE2, FUCA1, FUCA2, GAA, GLA, GNPTG, GUSB, HEXA, HEXB, HGSNAT, HYAL3, LAMP1, LAPTM4A, LAPTM4B, LGMN, LIPA, LITAF, MAN2B1, NAGLU, NAGPA, NAPSA, NEU1, NPC2, PPT1, PSAP, SLC11A2, SMPD1 and TPP1 were enriched in. Lysosome, known as the place for autophagy, has been recognized as a highly dynamic organelle that bidirectionally contact with mitochondria, to sense nutrient stress and control the switch between anabolism and catabolism by regulating mitochondrial biogenesis and autophagy (43, 44).

There is one pathway, p53 signaling pathway (hsa04115), only found in all 3,884 DE gene enrichment (), where AIFM2, BBC3, BCL2L1, BID, CCNE1, CDKN1A, CHEK2, COP1, GADD45A, GADD45B, IGFBP3, PERP, SHISA5, SIAH1, TP53 and TP53I3 were enriched. Chen et al. (45) revealed that MPC-5 cells from the high glucose-induced proliferation-inhibition and apoptosis-promotion via p53 signaling pathway can be protected by silencing the CCNG1 (45). Here, IGFBP3, as one of the key genes, interacted with 6 up-regulated and 12 down-regulated DE genes in the PPI networks (). It is reported that increased oxidative stress from high glucose enhances IGFBP3 expression to induce apoptosis; however, increased IGFBP3 expression by high glucose mediates high-glucose-induced apoptosis in PTECs and induces additional oxidative stress, which may result in amplification of hyperglycemic damage (46–48). The inhibition of LDH in hyperglycemic proximal tubular increased aerobic glycolysis (6), and changes in the extramitochondrial-free NADH/NAD+ ratio signal associated with aerobic glycolysis could control the abundance and activity of p53 by the C-terminal binding protein (CtBP) family of NADH-sensitive transcriptional regulators (49). Previous studies showed increased NADH/NAD+ ratio in response to the oxamate treatment in hyperglycemic rat proximal tubular cells (6). Here, we propose that the NADH-CtBP-p53 pathway may present as one of the hyperglycemia-dependent metabolic reprogramming caused by oxamate.

Conclusions

In summary, our study conducted the genome-wide transcriptome and co-expression analysis for HK-2 cells to identify 3,884 DE genes (adjusted P-value ≤ 0.05), where 1,664 of them were up-regulated and 2,220 of them were down-regulated. In addition, two modules of co-expression patterns were constructed to perform GO terms and pathways. We found that lysosome (hsa04142) can be enriched in DE genes of all, turquoise and blue modules, but p53 signaling pathway (hsa04115) was only detected in all DE genes. IGFBP3 (adjusted P-value = 1.34×10-75 and log2(FC) = 2.64), one of the key genes in p53 signaling pathway (hsa04115), interacted with several DE genes in the PPI networks and could be considered as the candidate biomarker due to its impact in high glucose conditions, after further validations. Together, our results highlight a possibility that anaerobic glycolysis-induced lactate in hyperglycemic HK-2 cells might integrate metabolic signals to mitochondria-lysosome contacts through some unexplored mechanism, which oscillates the coupling of mitochondrial biogenesis and mitophagy, and prolongedly compromises stress resistance.

Data Availability Statement

The raw RNA sequencing data were deposited in the Gene Expression Omnibus (GEO) of National Center for Biotechnology Information (NCBI) with the accession number GSE182138 at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE182138.

Author Contributions

ZW, XW, and GQ conceived and designed the experiments. ZW and JY performed the experiments. XW, DH, and ZW analyzed the data and wrote the manuscript. ZW, DH, DF, XW, and GQ improved the manuscript. All authors read and approved the final manuscript.

Funding

This study was supported by the Overseas Research Project of Scientific and Technological Talents in Health and Family Planning of Henan Province (No. 2015013) and National Natural Science Foundation of China (No. 82100896).

Conflict of Interest

XW was employed by company Konge Larsen ApS.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.