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 PGC1a. 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.

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% CO 2 ) 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 Dglucose (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×10 3 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-PGC1a (ab54481), Anti-CASP3 (9665), Anti-BCL2 (ab59348), Anti-BAX (GB11690), Anti-BECN1 (bs-1353R), Anti-MAP1LC3 (14600-1-ap) and Anti-b-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 mg 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) (Supplementary File 1) 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)(20)(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 log 2 (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 log 10 (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 1a (PGC1a), using the re-cultured and retreated cells. The selected genes might be involved in mitophagy and mitochondrial biogenesis for the flexibility of mitochondrial respiration adapting to nutrient stress. Glyceraldehyde-3phosphate 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− DDCt was used to calculate the relative gene expression levels. All the primers for qRT-PCR experiments are shown in the Table 1.

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 b 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 b power parameter (soft threshold) was equal to 8 when the R 2 of the free-scale topology was equal to 0.8 (Supplementary Figure 1).
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).

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 ( Figure 1A) and HG medium ( Figure 1B). It was suggested that oxamate affected the cell viability in a dosedependent manner in either LG or HG culture condition within 24 hours (Figure 1).
Further western blot assay was used to detect the protein expressions of PGC1a, CASP3, BCL2, BAX, BECN1 and MAP1LC3 (Supplementary Figure 2). The results showed that oxamate significantly decreased the protein expression levels of PGC1a 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 (Figures 2A, B). Culturing HK-2 cells in 25 mM HG medium with 40 mM oxamate for 24 hours is suggested to prompt a diabetic state  Table 3. We found that the most up-regulated DE gene was SAT1 (adjusted Pvalue = 5.57×10 -83 and log 2 (FC) = 1.75), while the most downregulated DE gene was DYNC1H1 (adjusted P-value = 3.66×10 -41 and log 2 (FC) = -1.14). In addition, the log 2 (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). Figure 3 displayed the obvious division between up-regulated DE genes and down-regulated DE genes based on log 2 (FC) values. The range of log 2 (FC) values reached to -5 for downregulated DE genes and 5 for up-regulated DE genes ( Figure 3A). 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 ( Figure 3B). 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 PGC1a (P-value < 0.05) showed the lower expressions ( Figure 3C).   Figure 4A). 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 ( Figure 4B). 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 ( Figure 4C). Obviously, DE genes were grouped into 2 modules that had similar co-expressions using the average linkage hierarchical clustering algorithm ( Figures 4B, C), 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 4D).

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 ( Figure 5A). Based on the regulation status of DE genes, 430 up-regulated and 1126 down-regulated significant GO terms ( Figure 5A), and 22 up-regulated and 58 down-regulated significant pathways ( Figure 5B) were achieved respectively. We found that the three most significant GO terms were chromosome segregation (GO:0007059, adjusted P-value =   Figure 6A). 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 ( Figure 6B). Based on 3,087 DE genes in the turquoise module after WGCNA analysis, the three most significant GO terms and pathways ( Figures 7A, B) 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 ( Figure 7B) 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 ( Figures 8A, B), where only RNA transport was still one of the three most significant pathways ( Figure 8B).

Protein-Protein Interaction Networks of Differentially Expressed Genes
The PPI network results were visualized using the top eight upregulated DE genes (SAT1, IGFBP3, AQP3, IFI6, S100A9, AKR1C1, CYP24A1 and AKR1C2) and the top six downregulated 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 (Figure 9). For example,

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 PGC1a 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 A B   (Figures 1, 2); 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 log 2 (FC) = 0.54) was identified as one up-regulated DE gene, whereas PDE3B (adjusted P-value = 1.80×10 -2 and log 2 (FC) = -0.52) as one down-regulated DE gene in our study (Supplementary File 2).
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-a), and lipopolysaccharide (LPS) as the stimuli (27)(28)(29)(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 b 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 log 2 (FC) = 1.75) ( Table 3). 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 A B  been reported that oxamate could induce autophagy via downregulation HIF-1a 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 log 2 (FC) = 0.71) and MAP1LC3B (adjusted Pvalue = 8.47×10 -4 and log 2 (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) ( Figure 3C) 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 ( Figure 6B), 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 ( Figure 9A). 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)(47)(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 upregulated 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 log 2 (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.