Multiomics Analyses Identify AKR1A1 as a Biomarker for Diabetic Kidney Disease

Diabetic kidney disease (DKD) is the leading cause of end-stage kidney disease. DKD is a complex disorder associated with numerous genetic loci in genome-wide association studies (GWAS) (1). The pathways that lead to worsening of kidney function in individuals with diabetes with persistent hyperglycemia include increased synthesis of extracellular matrix proteins, enhanced cell proliferation, and dysfunction of endothelial cells, together with tubular atrophy, interstitial fibrosis, and the thickening of the glomerular and tubular basement membranes. These are consequences of inflammation, oxidative stress, and production of advanced glycation end products (2).

A number of proteins involved in those pathways have been identified as biomarkers for DKD, including 17 circulating inflammatory proteins associated with development of end-stage kidney disease in the Joslin and Pima study group (3). Recently, Hirohama et al. (4) identified, via unbiased human kidney tissue proteomics, matrix metalloproteinase 7 (MMP7) as a diagnostic marker of kidney fibrosis and blood MMP7 as a biomarker for future kidney function decline in patients with DKD.

In this study, we applied multiple layers of omics analyses to narrow potential DKD biomarkers by using cell type–specific differential gene expression (single-cell RNA-sequencing [scRNA-seq]) and differentially expressed proteins in human kidney cortex from Mendeley Data. (4) We also incorporated GWAS top hits for estimated glomerular filtration rate (eGFR) from the most recent meta-analysis including Chronic Kidney Disease Genetics (CKDGen), Pan-UK Biobank, Million Veteran Program (MVP), Population Architecture Using Genomics and Epidemiology (PAGE), and Surrogate Markers for Micro- and Macrovascular Hard Endpoints for Innovative Diabetes Tools (SUMMIT) consortia (1) and subsequent cis–expression quantitative trait loci (eQTL) and cis–protein quantitative trait loci (pQTL) identification. We identified the specific aldo-keto reductase family 1 member A1 gene (AKR1A1) as a potential biomarker for DKD. We further investigated this gene by integrating data from our Wake Forest African American Diabetes Heart Study (AA-DHS), in which we are assessing the longitudinal effect of plasma metabolites on DKD progression. This systematic approach identifies cross-linked molecular signatures between metabolomics- and transcriptomics- or proteomics-indicated pathways in DKD and provides deep insights into intrinsic relationships between genetic variants, transcript dysregulation, protein/enzyme deficiency, and metabolic dysfunction in development of DKD.

Figure 1 displays the study design and workflow for biomarker identification and verification. The Institutional Review Board of Wake Forest University School of Medicine approved the local studies, and all patients provided written informed consent. For a detailed reporting of our methods, see Supplementary Methods.

All data generated and analyzed in this study are available from public databases (genomics, transcriptomics, and proteomics) or the corresponding authors (metabolomics and microscopy) upon reasonable request.

Differential gene expression was performed on kidney proximal tubule (PT) cells. A total of 790 differentially expressed genes were identified, with 530 upregulated and 260 downregulated genes in the Kidney Precision Medicine Project cohort (14 DKD cases vs. 20 control participants) in 64,333 kidney cells. Differential protein expression was performed in 33 human kidney cortex biopsy specimens (23 DKD cases vs. 10 participants) from the Mendeley Data. When gene identifiers (IDs) were merged with corresponding protein IDs, 24 gene–protein pairs were found to be in common in terms of ID and direction of differential expression with cutoff false discovery rate (FDR)–adjusted P value of <0.05 and |log₂ fold change [log₂FC]| >0.25 (Supplementary Table 1).

Among the 878 independent single nucleotide polymorphisms (SNPs) significantly associated with eGFR (P < 5 × 10⁻⁸) from the GWAS (GWAS Catalog no. GCST90100220) (1), only 10 fell in the 1 MB–centered cis region of the 24 prioritized genes, based on differential gene and protein expression (Supplementary Table 2). Among these 10 SNPs, only rs499600 was a cis-eQTL in the GTEx database for at least one tissue and a cis-pQTL for AKR1A1 in the Atherosclerosis Risk in Communities (ARIC) study plasma samples (http://www.nilanjanchatterjeelab.org/pwas). AKR1A1 was the only gene that survived the eGFR GWAS hit and cis-eQTL/pQTL screen. Additional AKR1A1 cis-SNPs were identified as cis-pQTLs of the protein level in ARIC samples. These SNPs are summarized in Supplementary Table 3. A posterior probability of the causal variant hypothesis (when H4 > 0.8 for colocalization analysis) is displayed in Supplementary Table 3. The linkage disequilibrium plot of these SNPs is shown in Supplementary Fig. 1.

As shown in Supplementary Table 1 and Fig. 2A, AKR1A1 expression was significantly lower only in PT cells, and not other cell types, of patients with DKD compared with those of control participants (FDR = 1.74 × 10⁻⁴³; log₂FC = −0.26). In the Nephrotic Syndrome Study cohort, AKR1A1 expression in kidney tubulointerstitium was positively associated with eGFR in patients with nephrotic syndrome (5) (NephroSeq V5; http://v5.nephroseq.org) in the microarray data (P = 1.11 × 10⁻⁴; r = 0.524) and positively associated with eGFR in patients with focal segmental glomerular sclerosis in the European Renal cDNA Bank (ERCB) RNA-seq data (P = 6.10 × 10⁻⁴; r = 0.744) (NephroSeq V5) (Supplementary Fig. 2).

The Mendeley data set contained 33 human kidney cortical samples: 23 samples from patients with DKD and 10 from healthy participants (https://data.mendeley.com/datasets/83k89shdx5/). Detailed demographics of the study sample have been described previously (4). Differential protein expression results show that AKR1A1 levels were significantly lower in patients with DKD compared with healthy control participants (FDR = 0.029; log₂FC = −0.61; see Supplementary Table 1).

Immunofluorescence on kidney cryosections verified that PT AKR1A1 protein levels in patients with DKD were significantly lower than that in nondiabetic control participants (P = 0.006; see Fig. 2B and C).

AKR1A1 is enriched in PT cells and downregulated only in PT cells in patients with DKD, compared with control participants. We conducted gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis on the basis of cell type–specific differential gene expression. A total of 152 KEGG terms reached the threshold of FDR <0.05, forming 30 clusters. AKR1A1 was involved in six of those pathways, and three pathways were in the top three clusters (Fig. 2D, red, blue, and green clusters) in terms of overall significance ranking. These pathways included chemical carcinogenesis-reactive oxygen species (cluster 1), ascorbate and aldarate metabolism (cluster 3), and pentose and glucuronate interconversions (cluster 3) (Supplementary Tables 4-1 and 4-2).

Figure 3 shows top serum metabolites that predicted development of DKD in the meta-analysis of Wake Forest AA-DHS participants (N = 624 for both discovery and replication sample sets). Between the top 20 metabolites in the two analysis models predictive of DKD in baseline samples, 17 were in common. Among these 17 metabolites, 3,4-dihydroxybutyrate had the highest odds ratio with FDR-adjusted P value of 2.11 × 10⁻⁹ and 5.58 × 10⁻⁹ in models 1 and 2, respectively (Supplementary Table 5).

Research in DKD has spawned numerous databases that are valuable for facilitating the identification of novel biomarkers using multidisciplinary approaches. We applied multiple layers of analyses using publicly available data sets together with Wake Forest AA-DHS metabolomics data to illuminate how a potential DKD biomarker, AKR1A1, was identified. To our knowledge, this is the first time a single gene conferring risk for DKD progression has been investigated through multiple lenses across genomics, single-cell transcriptomics, proteomics, and metabolomics.

As shown in Fig. 1, our strategy was to narrow down genes from incorporation of scRNA-seq–based PT cell–specific differential gene expression and proteomics-based differential protein expression in kidney cortex tissue. The common transcripts and proteins were those that survived cross-validation between transcriptomics and proteomics screens. We chose PT cells for differential gene expression because of their relevance to kidney functional decline or falling eGFR (6–8). Studies by researchers in the Susztak laboratory (9) reveal consistent enrichment of eGFR GWAS signals in PT cells. In addition, PT changes show a strong association with the glomerular filtration rate, and tubule cell dysfunction occurs in parallel with glomerular changes in DKD (10). In advanced DKD, kidney PT cell dysfunction correlates with eGFR decline. Recent genetic studies revealed the role of PT cells, not glomerular cells, in determining glomerular filtration (1). Hyperfiltration and kidney hypertrophy were the earliest changes in DKD (11). PT cells are mainly responsible for renal enlargement and hyperfiltration. PT injury eventually leads to development of tubulointerstitial inflammation and fibrosis in DKD (12). Independent GWAS eGFR hits from a large meta-analysis including five international consortia were applied (1) as additional layers of genomics data to filter findings from single cell transcriptomics and kidney cortex proteomics. We then used the GTEx project (13) (https://gtexportal.org/) and pQTL database from Johns Hopkins University’s ARIC Study (14) (http://nilanjanchatterjeelab.org/pwas/) to connect key molecules that formed a transomics biochemical network. These prioritized cis-eQTL variant-transcript pairs were further matched with plasma proteome cis-pQTLs, because cis-eQTLs and cis-pQTLs are considered important instruments to identify causative genes and proteins causing disease (15,16). Finally, we performed a longitudinal analysis to assess plasma metabolites predictive of DKD progression in our AA-DHS samples. We determined that AKR1A1 deficiency may be the hub in this enzyme-related pathway and contribute to progression of DKD.

As shown in Fig. 4, AKR1A1 is downregulated in PT cells from patients with DKD. Compared with healthy control participants, expression of AKR1A1 protein was also lower in kidney cortex specimens from patients with DKD. This was confirmed by our immunofluorescence imaging of kidney cryosections. A cis-pQTL (rs2229540, A/C) located in exon 3 of AKR1A1 gene was identified to be a functional (posterior probability = 0.99) variant (Asn52Ser) that affects plasma protein levels of AKR1A1 (14). A similar result was obtained from the ProteomeXchange Consortium for this causal variant in liver (17) (Supplementary Table 3). Although SNP rs2229540 was not a cis-eQTL for any tissue in the current GTEx project release (V8). This may be due to the small sample size used by the GTEx project and the low minor allele frequency of 0.05 of this SNP in the general population. A previous GWAS (1) identified SNPs rs499600 and rs659437, both located near AKR1A1, as associated with eGFR and as functional cis-eQTLs for brain and skeletal muscle; and identified cis-pQTLs (P = 1.45 × 10⁻⁵ and 3.19 × 10⁻⁵, respectively; although not robust enough to support a posterior P value for prediction of a functional variant). The fact that rs499600 is not a cis-eQTL in kidney tissue may undermine the importance of the relationship of this gene and the disease phenotype (i.e., eGFR). It is worth noting that rs499600 is located in a cis-regulatory 3′ region of AKR1A1, near an H3K27Ac mark (genome.ucsc.edu), where differential acetylation may result in differential AKR1A1 gene expression. Because we explored differential expression of AKR1A1 solely in PT cells, it will be important to determine cell type–specific eQTL analyses to profile tissues across the human body at single-cell resolution using single nucleus RNA-seq. The GTEx portal is currently working on this. The role of AKR1A1 in the regulation of eGFR was further supported by the Nephrotic Syndrome Study Network (NEPTUNE) cohort (5) (NephroSeq V5), where AKR1A1 expression in the renal interstitium was positively associated with eGFR in patients with nephrotic syndrome in microarray data, as well as eGFR in patients with focal segmental glomerular sclerosis in the ERCB RNA-seq data (NephroSeq V5).

Based on differential gene expression in PT cells, AKR1A1 may be involved in several pathways, including chemical carcinogenesis-reactive oxygen species, ascorbate and aldarate metabolism, and pentose and glucuronate interconversions. This is consistent with our metabolomics data showing that DKD progression can be predicted by 3,4-dihydroxybutyrate, a marker of succinate semialdehyde. AKR1A1 encodes aldo-keto reductase family 1 member A1, also known as alcohol dehydrogenase (NADP⁺) or aldehyde reductase. It catalyzes the NADPH-dependent reduction of a variety of carbonyl-containing compounds to their corresponding alcohols. It has enzymatic activity toward endogenous metabolites with a preference for negatively charged substrates, such as glucuronate and succinate semialdehyde (18). In kidney, AKR1A1 is mainly expressed in the renal cortical PTs, where it may be involved in detoxification processes. When aldehyde reductase (AKR1A1) is insufficient, it causes accumulation of succinate semialdehyde (Fig. 3B). Succinate semialdehyde is considered a reactive carbonyl and may lead to increased oxidative stress, which is viewed as a critical pathogenic factor implicated in the initiation, development, and progression of many kidney diseases, including DKD (19). Knockdown of AKR1A1 in cells increased reactive oxygen species levels (20), suggesting that the enzyme may be involved in redox-regulated activities and that it protects cells against oxidative stress (21). Fujii et al. (22) reported that Akr1a1 knockout mice were more likely to develop DKD, indicating the protective role of AKR1A1 from DKD.

In addition, glucuronate was also a significant predictor for DKD in our metabolomics study. The enzyme (aldehyde reductase) encoded by AKR1A1 is responsible for conversion of d-glucuronate to l-gulonate in the ascorbic acid (vitamin C) synthesis pathway. It reduces d-glucuronate to l-gulonate, which enters the pentose and glucuronate interconversion pathway (Fig. 3B). This pathway is very active in the renal cortex, which is also most abundant in AKR1A1 localized to the PTs (23). SNP-based pathway analysis implicated ascorbate and aldarate metabolism, and pentose and glucuronate interconversion in the pathogenesis of DKD in patients with type 1 diabetes (24). Proteomics and metabolomics studies also identified pentose and glucuronate interconversion as among the most affected pathways in early-stage and advanced DKD (25).

Although we do not have genotype data to link with plasma metabolite (3,4-hydroxybutyrate and glucuronate) levels, the alleles representing decreased AKR1A1 enzyme levels (plasma and liver) and lower eGFR in GWAS studies were featured by β estimates in the same direction (Supplementary Table 3). Thus, it is possible that decreased liver and kidney AKR1A1 enzyme levels work synergistically to worsen kidney function (eGFR). Moreover, association of AKR1A1 SNPs (rs659437 and rs499600) with serum cystatin-C concentration was identified in 363,228 UK Biobank participants (26) (Supplementary Table 3). Cystatin-C is freely filtered by the glomerulus and catabolized by renal proximal tubular cells (27) and has been recognized as a promising early DKD biomarker (28). It is tempting to speculate that AKR1A1, abundantly expressed in PT cells, could prevent DKD by reducing tubulointerstitial damage/fibrosis caused by oxidative stress and inflammation.

One limitation of this study is that various layers of omics studies came from different participants and resources; however, it is extremely difficult for one research group to conduct all layers of omics analyses from a single large cohort. Another limitation is that mixed populations were used in GWAS, scRNA-seq, and kidney tissue proteomics analyses. We have noticed that cis-pQTLs of the plasma proteome did show population deviation between European and African ancestry populations (14). However, AKR1A1 gene expression in tubular interstitium was consistently positively associated with eGFR either in African American people (NEPTUNE cohort) or White people (ERCB). Even though the metabolomic profile might be different, the function of a protein is thought to be universal across different populations.

In summary, by integrating single-cell transcriptomics, proteomics, metabolomics, and genomic architecture with cis-eQTL and cis-pQTL characteristics, we identified AKR1A1, as well as downstream metabolites, as potential biomarkers for DKD progression. The longitudinal analysis on metabolomics provides deep insight into the pathophysiological mechanisms of DKD. We confirmed the findings of previous studies on the importance of ascorbate and aldarate metabolism and pentose and glucuronate interconversion pathways, and detected the involvement of AKR1A1 in the development and progression of DKD at single-cell resolution. Longitudinal studies on plasma AKR1A1 protein as a predictive biomarker for DKD are needed. The role of AKR1A1 in nondiabetic chronic kidney disease warrants further investigation.

Acknowledgments. The authors are thankful to Kidney Precision Medicine Project investigators for sharing the kidney tissue scRNA-seq data; ARIC investigators and Dr. Nilanjan Chatterjee of Johns Hopkins University for sharing the pQTL data; and NEPTUNE investigators for sharing gene expression and kidney phenotype correlation data. The authors are also grateful to Nephroseq for providing a platform for data mining, and to the GTEx project for making the eQTL data publicly available. We are grateful to Dr. Katalin Susztak of the University of Pennsylvania for sharing proteomics data on the Mendeley Database and the generosity of ProteomeXchange Consortium to make proteomics data publicly accessible. Last, we appreciate National Human Genome Research Institute, European Bioinformatics Institute, University of Cambridge, and National Center for Biotechnology Information for making GWAS summary statistics available in the GWAS Catalog. The authors are thankful to the National Science Foundation ACCESS program (project BIO220154) for computational resource support in the public data analyses.

Funding. This work was supported by the National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK; grant DK071891 to B.I.F.) and the National Institute of Aging (grant AG058921 to N.D.P.). The results here are based, in part, on data generated by the Kidney Precision Medicine Project (https://www.kpmp.org; accessed 1 June 2023), which is funded by the NIDDK (grants U01DK133081, U01DK133091, U01DK133092, U01DK133093, U01DK133095, U01DK133097, U01DK114866, U01DK114908, U01DK133090, U01DK133113, U01DK133766, U01DK133768, U01DK114907, U01DK114920, U01DK114923, U01DK114933, U24DK114886, UH3DK114926, UH3DK114861, UH3DK114915, UH3DK114937).

Duality of Interest. B.I.F. has consulted for AstraZeneca Pharmaceuticals and Renalytix; received research funding from AstraZeneca Pharmaceuticals and Renalytix; and holds a U.S. patent with Wake Forest University Health Sciences and related to APOL1 genetic testing. No other potential conflicts of interest relevant to this article were reported.

Author Contributions. L.M. and B.I.F. conceived the project. D.L. performed cell type–specific differential gene expression (scRNA-seq), proteomics, and pathway analyses, and L.L. validated these analyses. F.-C.H. performed metabolomics analysis. M.M. enrolled patients who underwent nephrectomy for kidney sections and subsequent gene expression studies. Y.A.C. performed immunofluorescence staining. L.M. and Y.A.C. performed microscopy imaging. B.I.F. reviewed Kidney Precision Medicine Project pathology slides for pathologic diagnoses. L.M., B.I.F., and D.L. wrote the manuscript. M.M., N.D.P., J.S.P., and D.W.B. revised and edited the manuscript. All authors reviewed and approved the final version for publication. L.M. is the guarantor of this work and, as such, has full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.