key: cord-1021243-jbwk6q4b
authors: Wang, Zhimin; Hao, Dan; Fang, Dong; Yu, Jiating; Wang, Xiao; Qin, Guijun
title: Transcriptome Analysis Reveal Candidate Genes and Pathways Responses to Lactate Dehydrogenase Inhibition (Oxamate) in Hyperglycemic Human Renal Proximal Epithelial Tubular Cells
date: 2022-03-17
journal: Front Endocrinol (Lausanne)
DOI: 10.3389/fendo.2022.785605
sha: 44668e9a7650d477c2aaeb4acea3ce3bf3ab83ba
doc_id: 1021243
cord_uid: jbwk6q4b

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 log(2)(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.

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

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.

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.

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.

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 .

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.

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

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

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

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

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 ( Figure 6A ). 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 ( Figure 6B) . 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 ( Figure 6B ). Based on 3,087 DE genes in the turquoise module, the three most significant GO terms and pathways ( Figures 7A, B) 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 (Figures 8A, B) .

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, 

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.

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.

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.

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.

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.

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. 

CKD in Diabetes: Diabetic Kidney Disease Versus Nondiabetic Kidney Disease

Diabetic Kidney Disease

Mitochondrial Energetics in the Kidney

The Tubular Hypothesis of Nephron Filtration and Diabetic Kidney Disease

Current Challenges and Future Perspectives of Renal Tubular Dysfunction in Diabetic Kidney Disease

Metabolic Consequences of Lactate Dehydrogenase Inhibition by Oxamate in Hyperglycemic Proximal Tubular Cells

The Ugly Duckling of Energy Metabolism

Regulation of Blood Glucose by Hypothalamic Pyruvate Metabolism

WGCNA: An R Package for Weighted Correlation Network Analysis

A Gene-Coexpression Network for Global Discovery of Conserved Genetic Modules

Gene Coexpression Networks for the Analysis of DNA Microarray Data

Weighted Gene Co-Expression Network Analysis Identifies Modules and Functionally Enriched Pathways in the Lactation Process

Integrated Analysis of mRNA and MicroRNA Co-Expressed Network for the Differentiation of Bovine Skeletal Muscle Cells After Polyphenol Resveratrol Treatment

Assessment of the Low Inhibitory Specificity of Oxamate, Aminooxyacetate and Dichloroacetate on Cancer Energy Metabolism

HISAT: A Fast Spliced Aligner With Low Memory Requirements

HTSeq-A Python Framework to Work With High-Throughput Sequencing Data

Transcript Assembly and Quantification by RNA-Seq Reveals Unannotated Transcripts and Isoform Switching During Cell Differentiation

Transcriptomic Changes in Bovine Skeletal Muscle Cells After Resveratrol Treatment

DEGseq: An R Package for Identifying Differentially Expressed Genes From RNA-Seq Data

Differential Expression Analysis for Sequence Count Data

Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data With Deseq2

Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters

Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks

Mitochondrial Dysfunction in Diabetic Kidney Disease

Quantitative Ultrastructure of Human Proximal Tubules and Cortical Interstitium in Chronic Renal Disease (Hydronephrosis)

Oxamate Improves Glycemic Control and Insulin Sensitivity via Inhibition of Tissue Lactate Production in Db/Db Mice

Interleukin 6 Production by Human Proximal Tubular Epithelial Cells In Vitro: Analysis of the Effects of Interleukin-1a (IL-1a) and Other Cytokines

Interleukin 6 and Its Receptor: Ten Years Later

Role of MAP Kinase Pathways in Mediating IL-6 Production in Human Primary Mesangial and Proximal Tubular Cells

Interleukin-6. In: The Curated Reference Collection in Neuroscience and Biobehavioral Psychology

From PDE3B to the Regulation of Energy Homeostasis

Advances in Targeting Cyclic Nucleotide Phosphodiesterases

Phosphodiesterase 3b (PDE3B) Regulates NLRP3 Inflammasome in Adipose Tissue

Regulation of the Expression of the Hepatocellular Sulfate-Oxalate Exchanger SAT-1 (SLC26A1) by Glyoxylate: A Metabolic Link Between Liver and Kidney?

Oxamate-Mediated Inhibition of Lactate Dehydrogenase Induces Protective Autophagy in Gastric Cancer Cells: Involvement of the Akt-mTOR Signaling Pathway

Combination of Metformin, Sodium Oxamate and Doxorubicin Induces Apoptosis and Autophagy in Colorectal Cancer Cells via Downregulation HIF-1a

Beclin 1: A Role in Membrane Dynamics and Beyond

Depletion of Adipocyte Becn1 Leads to Lipodystrophy and Metabolic Dysregulation

Evaluation of Urinary Autophagy Transcripts Expression in Diabetic Kidney Disease

Caffeic Acid Attenuates Diabetic Kidney Disease via Modulation of Autophagy in a High-Fat Diet/Streptozotocin-Induced Diabetic Rat

Atg8 Family LC3/GAB ARAP Proteins Are Crucial for Autophagosome-Lysosome Fusion But Not Autophagosome Formation During PINK1/Parkin Mitophagy and Starvation

The E3 Ubiquitin Ligase NEDD4 is an LC3-Interactive Protein and Regulates Autophagy

Mitochondrial Dynamics at the Interface of Immune Cell Metabolism and Function

Lysosomes as Dynamic Regulators of Cell and Organismal Homeostasis

Silencing CCNG1 Protects MPC-5 Cells From High Glucose-Induced Proliferation-Inhibition and Apoptosis-Promotion via MDM2/p53 Signaling Pathway

Insulin-Like Growth Factor-Binding Protein-3 Mediates High Glucose-Induced Apoptosis by Increasing Oxidative Stress in Proximal Tubular Epithelial Cells

Coordination of Mitophagy and Mitochondrial Biogenesis During Ageing in C. Elegans

Autophagosome Formation and Cargo Sequestration in the Absence of LC3/GABARAPs

P53 Is Regulated by Aerobic Glycolysis in Cancer Cells by the CtBP Family of NADH-Dependent Transcriptional Regulators

The 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.Copyright © 2022 Wang, Hao, Fang, Yu, Wang and Qin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.