Transcriptome-wide association analysis identifies DACH1 as a kidney disease risk gene that contributes to fibrosis

More than 1 in 10 people in the US suffers from chronic kidney disease (CKD). CKD is one of the strongest risk factors for cardiovascular death. The annual mortality of diabetic patients with end-stage renal failure on dialysis can be as high as 15% (1). Currently available therapies slow disease progression, but do not stop or reverse kidney function decline, highlighting the critical need for better mechanistic understanding (2).

Kidney function is a heritable trait; therefore, substantial efforts have been dedicated to mapping genetic variants associated with kidney disease. Genome-wide association studies (GWAS) for kidney function are now available for more than a million participants. These studies have identified several hundred disease susceptibility loci (3–7). GWAS remains a powerful signal discovery method; however, disease-causing genes, cell types, and mechanisms encoded by the genetic signals remain poorly understood.

GWAS studies by themselves are unable to identify disease-causing genes, as close to 90% of GWAS-identified genetic variants are in the noncoding region of the genome. Expression quantitative trait loci (eQTL) are genomic loci that explain variation in gene-expression levels. As GWAS catalogues noncoding variants associated with disease and eQTL identifies the association between variants and gene expression, integration of GWAS and eQTL studies have been powerful in prioritizing potential causal genes from GWAS studies. Several methods have been developed for integrating GWAS and eQTL data sets. Our group has previously generated kidney compartment eQTL data sets and used Bayesian colocalization methods to understand whether genetic signals that control gene expression (in eQTL) and associate with phenotype (in GWAS) originate from the same genetic locus (8–10). Bayesian integration of kidney function–associated GWAS results with compartment-specific eQTL data was able to prioritize 27 genes for estimated glomerular filtration rate (eGFR) GWAS loci (8, 11).

Gene expression imputation followed by a transcriptome-wide association study (TWAS) has been recently proposed as a powerful approach for prioritizing candidate risk genes underlying complex traits. TWAS aims to test whether the expression of a gene mediates the genotype effect on the phenotype development. Different TWAS approaches have been developed. These approaches include Mendelian randomization using summary statistics (summary Mendelian randomization [SMR]) for eQTL and GWAS associations, using the top eQTL signal (12). Newer methods perform a multi-SNP prediction analysis that can more explicitly model linkage disequilibrium (LD) when causal variants are not genotyped. The most popular TWAS methods, such as TWAS FUSION and MetaXcan, use different regression models to “impute” cis-regulatory gene expression into much larger GWAS data sets to achieve this goal (13, 14).

The causal genetic variant identification in GWAS studies is hindered by LD. Most loci contain dozens, sometimes hundreds, of SNPs that pass the genome-wide significance for disease association. Recently, epigenetic information, such as functional annotation of the noncoding genomic regions, has been successfully used to prioritize GWAS variants (15). Functionally important GWAS variants are more likely to be on a cell-type–specific regulatory region; therefore, the intersection of GWAS loci with open chromatin regions identified by transposase insertion (ATAC-Seq) or histone ChIP (ChIP-Seq) can be powerful for causal variant prioritization. Regulatory regions are cell-type specific; therefore, new single-cell analytical methods open new avenues for causal variant prioritization that were missed in prior bulk experiments (15).

Here, we used a multipronged approach to prioritize genes, cell types, and causal variants for the functional annotation of kidney function GWAS (5, 6). First, we performed TWAS to prioritize genes whose expression mediate the genotype effect on phenotype. We investigated human kidney single-cell epigenome annotation by snATAC-Seq to highlight causal variants. Such analyses prioritized a small list of genes, including a cell-fate determination factor, Dachshund homolog 1 (DACH1), for kidney disease development. We used CRISPR-Cas9 gene editing to define the causal role of eGFR GWAS SNPs in regulating the expression of DACH1. Functional analysis using mice with tubule-specific genetic deletion or overexpression and cultured cell models demonstrated that DACH1 loss is associated with a proinflammatory tubule-cell phenotype and cytokine release, leading to macrophage infiltration, fibrosis, and CKD development.

eGFR GWAS target gene prioritization by TWAS. Here, we annotated eGFR-associated noncoding GWAS signals identified in the large multi-ethnic CKDGen (n = 765,348) and Million Veteran Program (MVP) cohorts (n = 280,722; refs. 5, 6). We used genotype and gene-expression information from microdissected human kidney tubule samples obtained from n = 121 subjects (8). We performed a summary-based TWAS using the FUSION TWAS pipeline (13). This approach uses reference LD and gene-expression panels with GWAS summary statistics to estimate the association between the cis-genetic component of gene expression and kidney function (Figure 1A). The TWAS FUSION pipeline prioritized 51 genes (Bonferroni’s test, P < 0.01) for eGFR in the CKDGen GWAS. Thirty-nine of the 51 prioritized genes also reached significance in the MVP GWAS (Table 1) The full results are shown in Supplemental Tables 1 and 2 (supplemental material available online with this article; https://doi.org/10.1172/JCI141801DS1). As a next step, in our gene-prioritization strategy, we performed SMR analysis followed by the heterogeneity in dependent instrument (HEIDI) test. This analysis prioritized 22 target genes (Bonferroni’s test, P < 0.01; HEIDI, P > 0.05) for eGFR in the CKDGen GWAS. Twenty out of the 22 prioritized genes were also significant in the MVP GWAS ( and Supplemental Tables 3 and 4). Finally, we performed MetaXcan and prioritized 130 genes (Bonferroni’s test, P < 0.01) for eGFR in the CKDGen GWAS, and 114 genes remained significant in the MVP GWAS (Supplemental Tables 5 and 6). Ten genes were identified by all methods, including DAB2, a gene previously identified using Bayesian colocalization and subsequently validated in animal models (ref. 8 and Figure 1G).

Figure 1

TWAS and regional association plots for kidney function and DACH1 expression. (A) Conceptual model for TWAS. (B–D) Regional association plots for DACH1 locus (index SNP; rs626277). (B) eGFR GWAS data from CKDGen study. (C) eQTL data in kidney tubules. (D) eQTL data in kidney glomeruli. Each dot represents 1 SNP. The dots are colored according to their relationship to the index SNP (rs626277). The red dots indicate high correlation (r² > 0.8) (LD) with the index SNP. The left y axis indicates ––log₁₀ (P value). The right y axis indicates recombination rate (cM/Mb). The red dotted horizontal line in the GWAS regional plot indicates the genome-wide significance threshold (5 × 10⁻⁸). Local permutation-based P value was used to define significance for the eQTL (tubule; P = 8.5 × 10⁻⁶, eQTL glom; P = 1.2 × 10⁻⁵). The x axis indicates the genomic location on chromosome 13. The arrow indicates the transcriptional direction for DACH1. (E) Box plot: the x axis represents the SNP (rs626277 genotype A/A, A/C, C/C), and the y axis shows relative DACH1 expression (P = 5.61 × 10^–08) in human kidney tubule samples. (F) Effect sizes of eGFR GWAS (CKDGen) SNPs (y axis) were plotted against eQTL SNPs (x axis). Error bars show the standard errors of SNP effects. Significant SNPs in SMR analysis were plotted. Red triangle indicates the top eQTL signal. (G) Venn diagram of TWAS, SMR, and MetaXcan listing the 10 genes prioritized by all 3 methods.

Table 1

List of genes prioritized by TWAS FUSION integration of eGFR GWAS

Table 2

List of genes prioritized by SMR and HEIDI integration of CKDGen eGFR GWAS

Identification of DACH1 as a kidney disease risk gene. To validate our computational gene prioritization, we focused on the chromosome 13 region, which showed a strong and reproducible association with kidney function (P = 1.96 × 10⁻¹⁹ in the CKDGen GWAS, Figure 1B). We observed that genetic variants at the eGFR GWAS region also showed a strong association with DACH1 expression in microdissected human kidney tubules (P = 1.32 × 10^–8; Figure 1C). Integration of the GWAS and eQTL signals using Bayesian colocalization (10) indicated shared genetic variants for the phenotype (eGFR) and gene expression (DACH1) in kidney tubules (PP4 = 0.973; Figure 1, B and C). In contrast, we did not detect significant eQTL signals in glomeruli (Figure 1D), even though DACH1 was highly expressed in glomeruli (16, 17). The relative DACH1 expression in renal tubules was significantly lower (P = 5.61 × 10^–08) in samples with the eGFR risk SNP allele A at rs626277 (Figure 1E). We further confirmed the rs626277 kidney tubule eQTL effect on DACH1 expression in an independent cis-eQTL study that included patients with nephrotic syndrome (ref. 18; P = 0.0084; Supplemental Figure 1, A and B). Metaanalysis of compartment-based eQTL and 44 GTEx tissues indicated that the eQTL effect of rs626277 on DACH1 was kidney-tubule specific (m = 1.00), potentially explaining the kidney phenotype development (Supplemental Figure 2). Summary data–based Mendelian randomization (12) indicated that direction and effect size consistent effects in eGFR GWAS and eQTLs by analyzing multiple variants (Figure 1F), suggesting a causal relationship.

To narrow the potential disease-causing variants and causal cell types, we analyzed single-cell open chromatin information of human and mouse kidney cells generated by snATAC-Seq (19). We found an open chromatin region, consistent with the promoter around the transcription start site (TSS) of DACH1, in all kidney epithelial cells (Figure 2A). We identified several cell type–specific open chromatin regions, located on intronic regions that were consistent with enhancers, in podocytes and distal convoluted tubules (DCT) in human kidney (Figure 2A). The eGFR GWAS variants did not overlap with the promoter region, but overlapped with the intronic open chromatin region in DCT cells (Figure 2B), suggesting that these SNPs might regulate the expression of DACH1 in DCT cells. On the other hand, we failed to observe a direct overlap between podocyte-specific open chromatin area and eGFR GWAS variants (Figure 2B). In addition to the human kidney, we also analyzed open chromatin data from mouse kidney samples by snATAC-Seq. The cell-type–specific open chromatin regulatory region of Dach1 showed strong conservation (Figure 2C).

Figure 2

Single-cell open chromatin data to narrow causal cell types and variants. (A) Single-cell open chromatin landscape of the human kidney at the DACH1 locus. The top row shows the genome organization, followed by eGFR GWAS-significant SNPs and open chromatin tracks in DCT. Boxes in light purple, blue, and red indicate genomic region deleted in CRISPR-Cas9 experiments. LOH, loop of Henle; CD_PC, principal cells of the collecting duct; CD_IC, intercalated cells of the collecting tubule; PT, proximal tubule. (B) Close-up of the open chromatin area in the vicinity of eGFR GWAS-significant variants. The index SNP from TWAS and SMR was plotted (rs500830 and rs932768, respectively). Note the overlap of risk SNP in DCT (left lower panel) and lack of overlap in podocytes (right lower panel). The box in light red indicates genomic region targeted by CRISPR-Cas9 deletion. (C) Single-cell open chromatin landscape of the Dach1 locus in the mouse kidney. The peak-to-peak coaccessibility loop was inferred by CICERO, and a heuristic cutoff of 0.25 was used. The top row shows the gene tracks, followed by the eGFR GWAS-significant SNP liftover track. (D) CRISPR-Cas9–assisted genomic editing in HEK293 cells stably expressing Cas9. Boxes in light purple, blue, and red indicate the deleted genomic regions. The guide RNAs for positive control, negative control, and risk region were designed. Bar graphs show relative DACH1 expression following guide RNA transfection for negative control, risk region, and positive control. **P < 0.01, 1-way ANOVA and Tukey’s post hoc test.

We reasoned that genetic risk variant directly overlapping with the open chromatin region in DCT could have the highest possibility of regulating DACH1 expression. To experimentally validate our hypothesis, we performed CRISPR-Cas9–assisted genome-editing experiments in HEK293 cells stably expressing Cas9. We designed guide RNAs to delete the intronic open chromatin risk region that included the eGFR GWAS variants (risk region); exon 1 was deleted as a positive control region and a region that lacked open chromatin peaks or significant eGFR GWAS variants was deleted as negative control (Figure 2D). DACH1 expression was markedly lower following the transfection of guide RNA that targeted positive control and risk regions compared with guide RNA that targeted the negative control region (Figure 2D).

In summary, GWAS and eQTL integration with snATAC-Seq indicated that genetic variants associated with kidney function localized to DACH1 regulatory regions in the distal part of the kidney tubule segment. The risk allele correlated with DACH1 levels in tubule samples and genome editing confirmed the role of this region regulating DACH1 levels.

DACH1 is expressed in podocytes and in distal tubule cells in mice and humans. To validate the protein expression and localization of DACH1, we performed double-immunofluorescence staining using mouse and human kidney tissue samples. Consistent with prior reports, we found that DACH1 was highly expressed in glomerular podocytes in mouse and human kidneys (Figure 3 and refs. 16, 17). Double staining with segment-specific markers indicated no detectable overlap with proximal tubule markers such as aquaporin1 (AQP1) and lotus tetragonolobus lectin (LTL) in the mouse kidney. On the other hand, we found that DACH1 was strongly expressed in cadherin-16 (Cdh16) (kidney-specific cadherin–positive [KSP-positive]) distal tubule cells and weakly expressed in the AQP2-positive collecting duct principal cells in the mouse kidney (Figure 3). Similarly, we detected DACH1 expression in calbindin-positive (CAL-positive) distal tubule cells and to some degree in the uromodulin-positive (UMOD-positive) ascending limb of the loop of Henle and AQP2-positive collecting duct principal cells in human kidney samples (Figure 3). Our immunofluorescence analysis indicated that despite DACH1 having an open-promoter area in all epithelial cells, expression of DACH1 was restricted, correlating with the intronic open chromatin area (enhancer), which was present in podocytes and distal kidney cells, such as ascending limb of the loop of Henle, distal tubule, and collecting duct principal cells.

Figure 3

DACH1 expression in mouse and human kidneys. (A–D) Representative immunofluorescence staining of DACH1 (red), DAPI (blue), and tubule-specific markers (green) in healthy adult mouse kidneys. Markers for proximal tubule, LTL (A) and AQP1 (B); distal tubule, Ksp-cadherin (KSP) (C); and collecting tubule, AQP2 (D). The right panels in A–D show higher magnification images of the regions within the white-dotted frames. G, glomerulus. (E–H) Representative immunofluorescence staining for DACH1 (red), DAPI (blue), and tubule-specific markers (green) in healthy human kidney. Markers for proximal tubule, LTL (E); thick ascending loop of Henle, UMOD (F); distal tubule, CAL (G); and collecting tubule, AQP2 (H). Right panels in E–H show higher magnification images of the regions within the white-dotted frames. Scale bars: 20 μm.

Mice with lower Dach1 expression in kidney tubule cells are more susceptible to kidney disease development. Integration of GWAS and eQTL studies suggested a causal role of lower tubule DACH1 levels in kidney disease development. To directly investigate the role of DACH1 in kidney disease, we generated mice with genetically lowered Dach1 expression in kidney tubule cells by crossing Ksp-cre mice with Dach1-floxed mice (20–22). Gene and protein expression analysis confirmed the reduction in Dach1 expression in Ksp^cre/Dach1^fl/WT and Ksp^cre/Dach1^fl/fl mice when compared with littermate WT mice (Figure 4B and Supplemental Figure 3A; see complete unedited blots in the supplemental material). Mice were born at the expected Mendelian ratio and appeared healthy at birth. We did not observe differences in life span at 40 weeks of age or structural changes on periodic acid–Schiff–stained (PAS-stained) kidney sections at 12 weeks of age (Figure 4H).

Figure 4

Mice with tubule-specific Dach1 loss are more susceptible to injury and fibrosis. (A) Experimental design: WT, Ksp^Cre/Dach1^fl/WT (Dach1 HZ), and Ksp^Cre/Dach1^fl/fl (Dach1 KO) mice were injected with FA or NaHCO3 (sham) i.p. Kidneys and serum were collected 7 days after injection. (B) Relative Dach1 mRNA expression in whole kidney tissue samples. (C–E) Profibrotic gene expression in whole kidney tissue samples: Col1a1 (C), Col3a1 (D), and Fn1 (E). Gene expression levels were normalized to Gapdh. (F) Representative Sirius red staining. Scale bar: 20 μm. (G) Quantification (as a percentage of positive area) of Sirius red staining in whole kidney samples. (H) Representative image of PAS-stained kidney sections. Scale bar: 20 μm. (I) BUN levels. Sham-treated group: WT (n = 4), Dach1 HZ (n = 4), Dach1 KO (n = 4); FA-treated group: WT (n = 10), Dach1 HZ (n = 8), Dach1 KO (n = 8). Light gray bars represent WT, red bars represent Dach1 HZ, and green bars represent Dach1 KO group. *P < 0.05; **P < 0.01; ***P < 0.001, 1-way ANOVA and Tukey’s post hoc test.

The functional effect of a quantitative change in gene expression, such as by the eGFR GWAS variant, usually becomes important in injury settings. Therefore, we next analyzed Ksp^creDach1^fl/WT and Ksp^creDach1^fl/fl mice in the folic acid–induced (FA-induced) kidney injury model (Figure 4A). Dach1 levels were lower in kidneys of Ksp^cre/Dach1^fl/WT mice and even lower in kidneys of Ksp^creDach1^fl/fl mice (Figure 4B). Markers of kidney fibrosis, such as collagen1a1 (Col1a1), collagen3a1 (Col3a1), and fibronectin1 (Fn1), were elevated in kidneys of FA-injected Ksp^creDach1^fl/WT mice and were even higher in kidneys of FA-injected Ksp^creDach1^fl/fl mice when compared with FA-injected WT mice (Figure 4, C–E). Sirius red staining to evaluate the degree of fibrosis confirmed the increase in fibrosis in FA-injected Ksp^creDach1^fl/WT and Ksp^creDach1^fl/fl mice (Figure 4, F and G). Histological analysis by PAS-stained kidney sections showed more severe damage in mice with lower tubule Dach1 levels (Figure 4H). Consistently, serum blood urea nitrogen (BUN), a measure of kidney function, was higher in Ksp^creDach1^fl/WT and Ksp^creDach1^fl/fl mice after FA treatment (Figure 4I).

We also analyzed the role of Dach1 in a diabetic kidney injury model induced by streptozotocin (STZ) injection and uninephrectomy (Unx) (Figure 5A). Albuminuria levels, quantified as urine albumin creatinine ratio, gradually increased with time. At 22 weeks of age, the degree of albuminuria in the diabetic Ksp^creDach1^fl/WT and Ksp^creDach1^fl/fl mice was markedly higher than in diabetic WT mice, even though blood glucose levels were similar between the diabetic groups (Figure 5, B and C). Expression of fibrosis-associated genes (Col1a1, Col3a1, and Fn1) was observably higher in kidneys of diabetic Ksp^creDach1^fl/WT and Ksp^creDach1^fl/fl mice (Figure 5, D–F).

Figure 5

Mice with tubule-specific Dach1 loss are more susceptible to diabetic injury. (A) Experimental design: WT, Ksp^Cre/Dach1^fl/WT (Dach1 HZ), and Ksp^Cre/Dach1^fl/fl (Dach1 KO) mice underwent nephrectomy and STZ injection. (B) Urine albumin-to-creatinine ratio at 6, 14, 18, and 22 weeks of age. The x axis indicates age (weeks). Black line, WT; red line, Ksp^Cre/Dach1^fl/WT (Dach1 HZ); green line, Ksp^Cre/Dach1^fl/fl (Dach1 KO). ^aP < 0.05, Dach1 KO vs. WT; ^bP < 0.01, Dach1 KO vs. WT; P < 0.05, Dach1 KO vs. Dach1 HZ; P < 0.05, Dach1 HZ vs. WT. (C) Blood glucose levels of sham-treated and Unx-STZ–treated groups. (D–F) Relative transcript levels of Col1a1 (D), Col3a1 (E), and Fn1 (F) in whole kidney samples of WT, Ksp^Cre/Dach1^fl/WT (Dach1 HZ), and Ksp^Cre/Dach1^fl/fl (Dach1 KO) mice. Gene expression levels were normalized to Gapdh. (B–F) Sham-treated group: WT (n = 3), Dach1 HZ (n = 3), Dach1 KO (n = 3); Unx-STZ treatment group: WT (n = 7), Dach1 HZ (n = 8), Dach1 KO (n = 6). *P < 0.05, **P < 0.01, ***P < 0.001, 1-way ANOVA and Tukey’s post hoc test.

In summary, we generated mice with tubule-specific heterozygous or homozygous deletion of Dach1. We observed increased kidney disease severity in mice with lower tubule Dach1 levels both in the FA and diabetic injury models, confirming the role of Dach1 in kidney disease development.

Elevated Dach1 levels protect kidney tubules from fibrosis development. Next, we generated mice with conditional inducible transgenic expression of Dach1 by placing the Dach1 gene under the tetracycline-inducible promoter (TRE Dach1; Figure 6A). TRE Dach1 mice were mated with Pax8rtTA mice to generate Pax8rtTA/TRE Dach1 double-transgenic mice, in which Dach1 expression is induced in kidney tubules by doxycycline (23). Mice were placed on doxycycline-containing chow at 6 weeks of age. Expression of Dach1 was increased (Figure 6B and Supplemental Figure 3B; see complete unedited blots in the supplemental material) in kidneys of Pax8rtTA/TRE Dach1 mice. Double-transgenic mice showed no baseline abnormalities. We induced kidney disease by FA injection, and mice were sacrificed 7 days following injections. Markers of kidney fibrosis, such as Col1a1, Col3a1, and Fn1, were lower in kidneys of FA-injected Pax8rtTA/TRE Dach1 mice when compared with FA-injected WT mice (Figure 6, C–E). Sirius red staining to quantify the degree of fibrosis confirmed the lower fibrosis in tubule-specific Dach1 transgenic mice following FA injection (Figure 6, F and G). Histological analysis by PAS-stained kidney sections showed less severe damage, and serum BUN levels were lower in FA-injected tubule-specific Dach1 transgenic mice (Figure 6, H and I). In summary, we found that mice with higher Dach1 levels were protected from kidney fibrosis, confirming the role of Dach1 in kidney disease development.

Figure 6

Mice with tubule-specific transgenic expression of Dach1 is protected from fibrosis. (A) Experimental design: Pax8rtTA/TRE-Dach1 (Dach1 OE) mice and control littermates (WT) were placed on doxycycline-containing chow. Mice were injected with FA or NaHCO3 (sham) i.p. Kidneys and serum were collected 7 days after injection. (B) Relative Dach1 mRNA expression in whole kidney tissue samples. (C–E) Profibrotic gene expression in whole kidney tissue samples: Col1a1 (C), Col3a1 (D), Fn1 (E). Gene expression levels were normalized to Gapdh. (F) Representative images of Sirius red staining. Scale bar: 20 μm. (G) Quantitative analysis (as a percentage of positive area) of Sirius red staining of kidney sections. (H) Representative images of PAS-stained kidney sections. Scale bar: 20 μm. (I) BUN levels. Sham-treated group: WT (n = 6), Dach1 OE (n = 6); FA-treated group: WT (n = 6), Dach1 OE (n = 8). Light gray bars represent WT; red bars represent Dach1 OE. *P < 0.05; **P < 0.01; ***P < 0.001, 1-way ANOVA and Tukey’s post hoc test.

Single-cell expression analysis of kidneys of mice with tubule-specific lower Dach1 levels identifies cycling cells. While we did not observe macroscopic phenotypic alterations in tubule-specific Dach1 heterozygous or homozygous mice, we next determined the impact of Dach1 loss on cell type–specific molecular changes by single-cell RNA-Seq (ref. 24 and Figure 7A). We analyzed 18,347 cells from 8-week-old Ksp^creDach1^fl/WT mice and compared with 17,106 cells obtained from WT mice. After quality control (QC), we kept 10,958 cells from Ksp^creDach1^fl/WT mice and 7800 cells from WT mice. Clustering analysis identified all known kidney cell types (Figure 7, B and C, and Supplemental Figure 4, A–C). The list of cell type–specific genes is shown in Supplemental Table 7. With this sample size, we did not identify obvious cell type–specific differentially expressed genes when Ksp^cre/Dach1^fl/WT and WT mice kidneys were compared. We found, however, that the percentage of cells identified as novel was much higher in Ksp^cre/Dach1^fl/WT kidneys (Figure 7B). Indeed, almost 90% of novel cells originated from the Ksp^cre/Dach1^fl/WT (Supplemental Figure 4D). Clusters labeled novel expressed high levels of cell cycle genes, including Mki67, Ccna2, Ccnb2, Top2a, and Stmn1 (Supplemental Table 7).

Figure 7

Single-cell RNA-Seq of WT mice and mice with tubule-specific heterozygous deletion of Dach1. (A) Experimental design: single-cell suspensions were generated from whole kidney tissue samples of WT and Ksp^Cre/Dach1^fl/WT (Dach1 HZ) mice. (B) Dimension reduction (UMAP) clustering. The following clusters were identified: novel cell types (Novel1, Novel2), NK cells, CD8 effector T cells (CD8 effector), B lymphocytes (B1, B2), granulocytes (Granul), macrophages (Macro), monocytes (Mono), B-type intercalated cells (B-IC), A-type intercalated cells (A-IC), principal cells of the collecting tubule (CD PC), connecting tubule (CNT), DCT, ascending loop of Henle (ALOH), descending loop of Henle (DLOH), proximal tubule S1 segment (PT S1), PT S2, PT S3, podocytes (Podo), endothelial cells (Endo). (C) Bubble plots of cell type–specific gene expression. The size of each dot corresponds to the percentage of positive cells, while the color intensity of each dot represents the average gene expression. (D) Relative transcript levels of Stmn1 and Top2a in kidneys of WT and Ksp^Cre/Dach1^fl/WT (Dach1 HZ) mice. Gene expression levels were normalized to Gapdh. (E) Representative immunohistochemistry staining of STMN1 in WT and Ksp^Cre/Dach1^fl/WT (Dach1 HZ) mice. Scale bar: 20 μm. (F) Number of STMN1-positive cells per high power field (HPF) in WT and Ksp^Cre/Dach1^fl/WT (Dach1 HZ) mice. WT (n = 6), Dach1 HZ (n = 6). *P < 0.05, 2-tailed Student’s t test.

To validate the single-cell RNA-Seq analysis, we analyzed an additional 6 WT and 6 Ksp^cre/Dach1^fl/WT mice. We found that the expression of Stmn1 and Top2a (markers of cycling cells) was elevated in Ksp^cre/Dach1^fl/WT mice (Figure 7D). Finally, immunostaining for STMN1 and Ki-67 in kidney tissue confirmed the increase in STMN1- and Ki-67–positive tubule cells in Ksp^cre/Dach1^fl/WT mice (Figure 7, E and F). Overall, our single-cell analysis highlighted more cycling cells in mice with renal tubule–specific heterozygous Dach1 loss.

DACH1 transcriptionally controls cell cycle genes and cell proliferation. To understand the molecular pathways controlled by DACH1, we examined DACH1 ChIP-Seq data (25). The analysis identified 19,460 peaks corresponding to approximately 5000 genes controlled by DACH1 (Figure 8A). To understand the molecular function of genes bound by DACH1, we performed gene ontology analysis (Figure 8B and Supplemental Table 11). The gene ontology analysis identified enrichment for several gene groups. The most prominent enrichment was for genes associated with cell cycle, consistent with the single-cell analysis.

Figure 8

DACH1 transcriptionally controls cell cycle genes in kidney tubule cells. (A) DACH1-ChIP 19,460-binding sites, corresponding to 5111 genes in the K562 cell (human immortalized myelogenous leukemia). Lower panel shows the distribution of peaks compared with TSS distance. (B) Top 10 most significant pathways from the Molecular Signatures Database (MsigDB) using the Genomic Region of Enrichment of Annotations (GREAT) package (see Supplemental Methods). (C) Primary mouse kidney tubule cells were isolated from Dach1^fl/fl mice and infected with adenovirus-GFP (Ade-GFP; control [CTRL]) or adenovirus-Cre-GFP (Ade-Cre-GFP; Dach1 KO). Representative Western blots of DACH1 protein expression in control and Dach1 KO tubule cells. GAPDH was used as loading control. (D) DACH1-ChIP-PCR for Ccnd1 and Ccnd2 in control kidney tubule cells. Three independent experiments were performed. (E) Integrated Genome Browser view of the human CCND2 locus. From top to bottom, the labels indicate the CCND2 locus, DACH1–ChIP-Seq, followed by human kidney snATAC-Seq in DCT, proximal tubule, loop of Henle, principal cell, intercalated cell. The box highlights the overlapping peaks between DACH1–ChIP-Seq and snATAC-Seq. (F) Relative gene expression of Ccna1, Ccnd1, Ccnd2, and Ccne1 in Ade-GFP or Ade-Cre-GFP transfected cells (n = 3). Gapdh was used as internal control. (G) Results of cell proliferation analysis using the MTT assay at 24 and 48 hours after Ade-GFP or Ade-Cre-GFP infection. (n = 3 in triplicate). *P < 0.05, **P < 0.01, 1-way ANOVA and Tukey’s post hoc test.

To validate the direct regulation of genes by DACH1 in kidney tubule cells, we performed DACH1-ChIP-PCR in cultured mouse primary kidney tubule cells. We isolated cells from Dach1^fl/fl mice and infected them with adenoviral GFP or Cre-GFP to generate control and Dach1-deficient cells, respectively (Figure 8C). We found that DACH1 was bound to the regulatory region of key cell cycle regulator genes, including Ccnd1 and Ccnd2, in control tubule cells (Figure 8D), but not in Dach1-deficient tubule cells (Supplemental Figure 5A). By combining DACH1-ChIP and snATAC-Seq information, we confirmed that open-chromatin region in DCT overlapped with DACH1-binding sites at the CCND2 locus (Figure 8E).

To understand whether DACH1 transcriptionally controls cell cycle regulator genes we performed quantitative reverse-transcriptase PCR (QRT-PCR) in control and Dach1-deficient tubule cells. Expression of Ccna1, Ccnd1, Ccnd2, and Ccne1 was higher in Dach1-deficient tubule cells (Figure 8F). These results are consistent with previous reports indicating that DACH1 functions as a transcriptional repressor of cell cycle genes (26–28). Finally, to understand whether changes in expression of genes associated with cell cycle cause functional alterations, we measured cell proliferation in control and Dach1-deficient tubule cells. We found a significant increase in cell proliferation 24 and 48 hours after Dach1 deletion (Figure 8G).

As DACH1 is also expressed in podocytes, we investigated changes in expression of cell cycle regulators in control and DACH1 knockdown human podocytes. We found that the expression of CCND1 (but not other genes) was slightly but significantly higher in DACH1 knockdown podocytes (P < 0.05; Supplemental Figure 6A).

DACH1 has been shown to interact with SMAD4 and TGF-β signaling and promote epithelial-to-mesenchymal transition (EMT) (26). We found no consistent changes in EMT markers, such as Snai1, Twist1, and Acta2, in Dach1-deficient tubule cells by QRT-PCR analysis. Expression of Col1a1 was mildly higher in Dach1-deficient tubule cells (Supplemental Figure 7).

Next, we set to determine whether DACH1 controls cell cycle and proliferation in kidneys of mice with different expression levels of Dach1. We found that expression of cell cycle markers, such as Ccnd2 and Mki67, was higher in kidneys of FA-injected mice (Figure 9, A and B). However, mice with tubule-specific heterozygous and homozygous deletion of Dach1 had even higher expression of Mki67 and Ccnd2 when compared with WT mice (Figure 9, A and B). Consistent with the transcript expressions, the number of Ki-67–positive renal tubule epithelial cells was increased in FA-injected mice with tubule-specific heterozygous and homozygous deletion of Dach1 compared with FA-injected WT mice (Figure 9, C and D). Kidney expression of Ccnd2 and Mki67 was lower in FA-injected tubule-specific Dach1 transgenic mice than in FA-injected WT mice (Figure 9, E and F). These results indicate that DACH1 functions as a transcriptional repressor controlling cell cycle genes in vitro and in vivo.

Figure 9

DACH1 controls cell cycle in vivo. (A and B) Relative transcript expression of Ccnd2 (A) and Mki67 (B) in whole kidney samples of sham-treated and FA-injected WT, Ksp^Cre/Dach1^fl/WT (Dach1 HZ), and Ksp^Cre/Dach1^fl/fl (Dach1 KO) mice. (C and D) Representative Ki-67 immunostaining (C) and quantification of Ki-67–positive cells per high-power field in kidney tissue of sham-treated and FA-injected WT, Ksp^Cre/Dach1^fl/WT (Dach1 HZ), and Ksp^Cre/Dach1^fl/fl (Dach1 KO) mice. Scale bar: 20 μm. (E and F) Relative transcript levels of Ccnd2 (E) and Mki67 (F) in whole kidney samples of sham-treated and FA-injected WT and Pax8rtTA/TRE-Dach1 (Dach1 OE) mice. (A–D) Sham-treated group: WT (n = 4), Dach1 HZ (n = 4), Dach1 KO (n = 4); FA-treated group: WT (n = 10), Dach1 HZ (n = 8), Dach1 KO (n = 8). Light gray bars represent WT, red bars represent Dach1 HZ, and green bars represent Dach1 KO group. (E and F) Sham-treated group: WT (n = 6), Dach1 OE (n = 6); FA-treated group: WT (n = 6), Dach1 OE (n = 8). Light gray bars represent WT,; red bars represent Dach1 OE. *P < 0.05, **P < 0.01, ***P < 0.001, 1-way ANOVA and Tukey’s post hoc test.

DACH1 controls myeloid chemotactic factors and macrophage influx. The CHiP data indicated that genes directly regulated by DACH1 were enriched for proinflammatory cytokines (Figure 8B). To validate that these cytokines are also regulated by DACH1 in kidney tubule cells, we performed DACH1-ChIP-PCR in cultured mouse kidney tubule cells (Figure 8C). We confirmed the direct binding of DACH1 to the promotor region of Ccl2 (monocyte chemotactic protein [Mcp1]), macrophage colony stimulating factor (Csf1), and intercellular adhesion molecule 1 (Icam1) in control tubule cells (Figure 10, A–C), but not in Dach1 knock-out tubule cells (Supplemental Figure 5B). DACH1 binding sites overlapped with open chromatin areas in DCT cells at the CSF1 locus (Figure 10D). Expression of Ccl2, Csf1, and Icam1 was higher in kidney tubule cells in the absence of Dach1 (Figure 10E). CCL2 and CSF1 are well-known myeloid and macrophage chemotactic factors. As seen before, we confirmed that DACH1 acts as an upstream transcriptional repressor. We found that expression of CCL2 was increased in DACH1 knockdown human podocytes; however, we did not observe differences in CSF1 and ICAM1 levels, indicating some shared and cell type–specific DACH1 targets (Supplemental Figure 6B).

Figure 10

DACH1 transcriptionally controls cytokine expression in cultured kidney tubule cells. (A–C) DACH1-ChIP-qPCR for Ccl2 (A), Csf1 (B), Icam1 (C) using primary mouse kidney tubule epithelial cells. The y axis is presented as percentage of input (n = 3). (D) Integrated Genome Browser view of the human CSF1 locus. From top to bottom, the labels indicate CSF1 locus, DACH1-ChIP-seq, followed by human kidney snATAC-Seq (open chromatin region) in DCT, proximal tubule, loop of Henle, principal cells, intercalated cells. The red box highlights the overlapping peaks of DACH1 ChIP-Seq and human kidney snATAC-Seq. (E) The relative expression of Cxcl1, Cxcl2, Ccl2, Ccl5, Tnfa, Csf1, Il1b, Il17c, Il27, and Icam1 in Ade-GFP or Ade-Cre-GFP transfected renal tubule cells from Dach1^fl/fl mice (n = 3). Gapdh was used as internal control. (F) Experimental scheme of Raw 264.7 macrophage chemotaxis assay. The supernatants of Ade-GFP or Ade-Cre-GFP-transfected renal tubule cells from Dach1^fl/fl mice were added to the lower chamber. The number of macrophages (macrophage migration) was analyzed by relative fluorescence units (RFU) after 2 or 24 hours (n = 3 in triplicate). *P < 0.05, **P < 0.01, ***P < 0.001, 1-way ANOVA and Tukey’s post hoc test.

To confirm the functional role of cytokines secreted by tubule cells, we performed macrophage migratory assay. We found that the supernatant obtained from Dach1-deficient kidney tubule cells increased macrophage migration compared with control supernatant (Figure 10F).

To further confirm the role of DACH1 in regulating a proinflammatory tubule phenotype in vivo, we analyzed mice with differing Dach1 gene dosage in kidney tubules. We found that expression of Ccl2 and Csf1 was slightly but significantly (P < 0.05) higher in kidneys of mice with tubule-specific homozygous deletion of Dach1 compared with WT mice at baseline. (Figure 11, A and B). We observed that Ccl2, Csf1, and Icam1 expression levels were markedly increased in kidneys following FA injection and that their expression levels were higher in kidneys of mice with tubule-specific heterozygous and homozygous deletion of Dach1 when compared with WT mice (Figure 11, A–C). On the other hand, mice with tubule-specific transgenic expression of Dach1 showed lower expression of Ccl2, Csf1, and Icam1 in kidneys following FA injection (Figure 11, F–H). As CCL2 and CSF1 are the key macrophage chemotactic factors, we examined macrophage infiltration in kidneys. Macrophage number, as analyzed by F4/80 expression, was slightly higher in kidneys of mice with tubule-specific homozygous deletion of Dach1 at baseline. We observed that the number of macrophages was markedly increased in mice with tubule-specific heterozygous and homozygous deletion of Dach1 following FA treatment (Figure 11, D and E).

Figure 11

DACH1 controls cytokine expression and macrophage infiltration. (A–C) Relative transcript expression of Ccl2 (A), Csf1 (B), and Icam1 (C) in kidneys of sham-treated and FA-treated WT, Ksp^Cre/Dach1^fl/WT (Dach1 HZ), and Ksp^Cre/Dach1^fl/fl (Dach1 KO) mice. (D–E) Representative immunostaining with F4/80 (macrophage marker) (D) and quantification of F4/80-positive area (E) in kidney tissue of sham-treated and FA-treated WT, KspCre/Dach1^fl/WT (Dach1 HZ), and Ksp^Cre/Dach1^fl/fl (Dach1 KO) mice. Scale bar: 20 μm. (F–H) Relative transcript expression of Ccl2 (F), Csf1 (G), and Icam1 (H) in kidneys of sham-treated and FA-treated WT and Pax8rtTA/TRE-Dach1 (Dach1 OE) mice. Gapdh was used as internal control. (A–D) Sham-treated group: WT (n = 10), Dach1 HZ (n = 10), Dach1 KO (n = 10); FA-treated group: WT (n = 10), Dach1 HZ (n = 8), Dach1 KO (n = 8). Light gray bars represent WT, red bars represent Dach1 HZ, and green bars represent Dach1 KO. (E and F) Sham-treated group: WT (n = 6), Dach1 OE (n = 6); FA-treated group; WT (n = 6), Dach1 OE (n = 8). Light gray bars represent WT; red bars represent Dach1 OE. *P < 0.05, **P < 0.01, ***P < 0.001, 1-way ANOVA and Tukey’s post hoc test.

Overall, our results indicate that DACH1 acts as a transcriptional repressor of myeloid chemotactic factors in kidney tubule cells, controlling macrophage infiltration and kidney disease development.

DACH1 levels correlate with fibrosis, proliferation, and inflammation in patient kidneys. In order to understand the impact of altered DACH1 expression in human kidney fibrosis, we performed double-immunofluorescence studies on healthy and CKD human kidney samples. We found that DACH1 expression was lower in the distal tubule segment in kidneys of patients with CKD (Figure 12A).

Figure 12

DACH1 levels correlate with cell cycle, inflammation, and disease severity in human kidney tubule samples. (A) Representative immunofluorescence staining of DACH1 (red), DAPI (blue), and CAL (green). The overlap of DACH1/DAPI/CAL in healthy (control) (upper panel) and CKD (lower panel) human kidney. Right panels show higher-magnification images of the regions within the dotted frames. Scale bar: 20 μm. (B) Gene expression analysis was performed using 95 microdissected human kidney tubules obtained from healthy and diseased samples. (C) Correlation between DACH1 normalized expression levels in human kidney tubules and corresponding eGFR (ml/min/1.73 m²). Correlation coefficient (Cor) = 0.20; P = 0.048. (D) Correlation between COL1A1 and DACH1 normalized expression levels in human kidney tubules. Cor = –0.45; P = 2.8 × 10^–05. (E) Gene Ontology analysis of transcripts correlated with DACH1 levels. The x axis represents the –log₁₀ of the P value. Inflammatory or immune-related pathways are highlighted in red and development or cell proliferation–related pathways are highlighted in blue. (F and G). Correlation between CCND2 (F), CCL2 (G), and DACH1 normalized expression levels in human kidney tubules. Cor = –0.44; P = 5.1 × 10^–05 (F), Cor = –0.21 P = 0.04 (G). Correlation analyses were performed using a linear regression model adjusted for sex, age and race.

We next analyzed gene expression from 95 microdissected human kidney tubule samples, including healthy subjects and subjects with varying degrees of diabetic and hypertensive kidney disease (ref. 29 and Figure 12B). DACH1 expression correlated positively with kidney function (eGFR) and inversely with kidney fibrosis, defined by COL1A1, COL3A1, and FN1 levels, consistent with the protein expression and genetic data (Figure 12, C and D, and Supplemental Figure 8, A and B).

Next, we identified genes whose expression levels correlated with DACH1 levels in microdissected human kidney tubule samples. The expression of 513 genes showed negative correlation with DACH1 levels in microdissected human kidney tubule samples. Functional grouping by ontology analysis indicated that DACH1 expression correlated with genes with immune response, leukocyte activation, and cell proliferation (Figure 12E). We found DACH1 expression levels in human kidney tubule samples strongly correlated with CCND2 and MKI67 expression levels (Figure 12F and Supplemental Figure 8C). Lower DACH1 levels were associated with higher CCL2, CSF1 (Figure 12G and Supplemental Figure 8D), and ICAM1 expression levels (Supplemental Figure 8E). Consistently, expression levels of COL3A1, MKI67, and CCND2 correlated with renal function (Supplemental Figure 8, F–H).

Overall, we observed a conserved decrease in DACH1 expression in patients with CKD. DACH1 expression in human kidney tubule samples strongly correlated with expression of DACH1 target genes, such as CCND2 and CCL2, and kidney function and fibrosis.

In summary, here we performed a multistaged approach to functionally annotating risk variants identified by eGFR GWAS. We performed computational integration of kidney function GWAS and gene expression data via the TWAS method. While prior Bayesian colocalization studies highlighted 24 loci in which gene expression and phenotype signal originate from the same genetic locus (10), TWAS tests whether gene expression mediates the genotype effect on phenotype. Our analysis used TWAS FUSION, SMR, and MetaXcan to integrate GWAS and gene expression data. We leveraged the CKDGen GWAS (5) for discovery and MVP GWAS (6) for validation. Our comprehensive analysis prioritized 39 genes in TWAS FUSION, 20 in SMR, and 114 in MetaXcan. Genes nominated by individual TWAS and colocalization methods could play an important role in disease development, but require functional validation. Here, we decided to focus on DACH1, which was nominated consistently by all TWAS methods and replicated by both GWAS studies. Although the association between SNPs on chromosome 13 and kidney function has been described by prior GWAS (3, 4, 6), the causal variants, responsible gene, cell type, and biological mechanism for kidney disease development have not been identified systematically. Here we addressed these challenges for what we believe is the first time.

Prior studies have implied that rare loss-of-function coding mutations in DACH1 can cause kidney developmental defects. For example, double-homozygous coding missense mutations of DACH1 and BMP4 were found in a patient with bilateral renal cystic dysplasia (30). DACH1 mutations have also been described in branchio-oto-renal (BOR) syndrome, which comprises numerous congenital anomalies characterized by branchial arch deformation (31). During development, DACH1 plays a role in epithelial-mesenchymal interaction patterning and cell-fate determination (27, 32, 33), which explains the kidney development defect in patients with DACH1 loss-of-function mutations. We believe that our work illustrates the convergence of classic monogenic and complex trait diseases by DACH1, such as the kidney disease caused by rare loss-of-function coding mutations and the conferred CKD risk mediated by common noncoding regulatory variants.

Our study evaluated the role of variants in the noncoding (intronic) region on chromosome 13 that are shown to be associated with kidney function in multiple eGFR GWAS (3–6). Our single-cell epigenetic studies–highlighted candidate variants were localized to intronic open chromatin regions (likely enhancers) in DCT cells. These intronic regions are critical for gene expression regulation, so that while all kidney epithelial cells had an open promoter region, the expression of DACH1 strongly correlated with intronic enhancer regions that were present only in podocytes and the distal part of the kidney tubules. We experimentally validated that eGFR risk variant located on the intronic enhancer region in distal tubule cells indeed altered DACH1 expression. DACH1 is highly expressed in podocytes, and the role of DACH1 in podocyte dysfunction has been described earlier (16, 17). We did not observe significant colocalization of the GWAS signals with the glomerular eQTL and podocyte-specific open chromatin areas. Therefore, we believe that the eGFR GWAS signals are related to DACH1 expression in tubules rather than podocytes; however, the role of podocytes cannot be excluded. Our results demonstrated the utility of single-cell epigenetic fine mapping in the detection of potential causal SNPs and disease-responsible cell types.

The effect of heterozygous Dach1 loss in kidney tubule cells was observed on kidney function and fibrosis development in different kidney disease models, such as that induced by FA injection and diabetes. Genetic overexpression of Dach1 in kidney tubules protected from kidney injury development in the FA-induced kidney injury model. The protective effect of Dach1 was partial; at the same time, we only achieved a 2- to 4-fold increase in Dach1 expression at the whole kidney level. We found an increase in the number of cells expressing cell cycle genes in kidney tubules of mice with heterozygous loss of Dach1 (even without injury) by the sensitive scRNA-Seq method. However, using traditional kidney phenotyping methods, such as histological and functional analysis, we failed to identify observable phenotypic changes in mice with heterozygous loss of Dach1. In the future, more sensitive analytical methods will be needed to detect minor molecular alterations in the kidney.

DACH1 has DNA-binding properties and functions as a transcriptional repression and tumor suppressor (27, 32). Prior studies reported the role of DACH1 in cell cycle regulation in podocytes (17), immortalized human renal tubule cells (17), and breast cancer cells (27), supporting our result that DACH1 controlled the expression of cell cycle regulators, such as Ccnd1 and Ccnd2. The number of dividing cells is known to be higher both in acute kidney disease and CKD, and dividing cells are usually less differentiated compared with nondividing cells (34). Increased expression of Ccnd2 and Ki67 indicates that cells entered the cell cycle, but it does not necessarily indicate the completion of cell division. Recent studies highlighted that tubule cells enriched in the G₂/M phase of the cell cycle have proinflammatory properties and play a key role in fibrosis development (35, 36). Our initial results are consistent with an increase in G₂/M arrest (data not shown).

We found that DACH1 also controls the expression of several cytokines, for instance, Ccl2 and Csf1, in renal tubule cells. These cytokines are secreted by injured tubule cells (37, 38). Multiple groups have recently identified the proinflammatory tubule subpopulation in injured kidneys based on kidney single-cell analysis (39). It has not yet been established what drives the proinflammatory transdifferentiation of renal tubule cells. Our data show that DACH1 plays a pivotal repressor role in controlling cytokine expression in renal tubule cells. Consistent with our result, a recent report indicated that DACH1 was associated with a proinflammatory phenotype of high glucose–treated cultured renal tubule cells (40). In mechanistic studies, DACH1 antagonizes FOXM1 promoter occupancy and target gene expression (28). A recent report also highlighted the critical role of FOXM1 in kidney tubule in the context of injury (41). Future studies shall focus on examining the relationship between DACH1 and FOXM1 in kidney tubules. Myeloid chemotactic factors and macrophages play a key role in kidney function decline. Preclinical and clinical studies have persuasively established the role of macrophages, especially CCL2, in kidney disease development (42, 43). DACH1 being a direct regulator of CCL2 can explain the mechanism of DACH1-induced kidney disease and fibrosis development.

In summary, here we show that DACH1 is a kidney disease risk gene. Kidney disease genetic risk variants lower DACH1 expression in kidney tubule cells and increase susceptibility to kidney injury by inducing cycling and cytokine-secreting kidney tubules, causing the influx of inflammatory macrophages and fibrosis development.

TWAS FUSION. We downloaded the FUSION package and the GWAS summary statistics for the MVP and CKDGen studies (5, 6). We generated gene expression weights using our own human kidney tubule RNA-Seq (n = 121) data following the FUSION pipeline (13). Genotypes were imputed to the 1000 Genomes, phase 1, version 3, and restricted to well-imputed (INFO > 0.9) sites. RPKM and log-adjusted gene expression levels were estimated in a generalized linear model controlling for 3 gene expression PCs and rank normalized. We filtered genes that did not exhibit cis-genetic regulation at current sample sizes by keeping only genes with nominally significant estimates of cis-SNP heritability. We refrained from reporting genes in the HLA region due to the complicated LD patterns. To train predictive models, FUSION defines gene expression for samples as a linear function of SNPs (in a 1 Mb region flanking the gene), where the SNP weights, covariates (e.g., sex, age, genotype principal components, genotyping platform, PEER factors), their effects, and random environmental noise are included (44). FUSION estimated weights for expression of a gene in a tissue using multiple penalized linear models (here we used LASSO).

SMR. Summary data–based Mendelian randomization is used to test for potential causal effects of, e.g., gene on complex trait, given a SNP as an instrumental variable, using summary-level data from eGFR GWAS and eQTL studies (12). We performed SMR on eGFR GWAS and eQTLs to estimate the direction of expression changes, using multiple variants in the cis-eQTL region of genes.

MetaXcan. MetaXcan (14) was used to integrate genomic information of cis-eQTL with eGFR GWAS traits, using summary-level data of eGFR GWAS (5, 6) and eQTL study. LD references were estimated based on genotypes of individuals from the 1000 Genome Project (phase 3).

Mice. The KspCre and Pax8rtTA mice were purchased from Jackson Laboratory (23, 45). Dach1^fl/WT or Dach1^fl/fl mice were previously described (22). TRE Dach1 transgenic mice were provided in house. Male mice, 8 to 10 weeks old, and littermates were used for the experiments. Mice were given food and water ad libitum. FA was dissolved in 290 mM NaHCO₃. Mice were injected with FA (250 mg/kg) or NaHCO₃ i.p. Mice were sacrificed 7 days after FA injection. Unx was performed at 6 weeks. STZ was injected at 50 mg/kg i.p. for 5 consecutive days 1 week after Unx. Urine was collected at 14, 18, and 22 weeks of age. Mice were sacrificed at 24 weeks of age. Kidneys were harvested and preserved at –80 °C for RNA and protein analysis or 10% formalin for histology. Serum was collected from the inferior vena cava. The primer sequences used for qPCR are listed in Supplemental Table 8.

Single-cell RNA-Seq of mouse kidney. Single-cell RNA-Seq libraries were generated using the 10× Chromium Single Cell instrument and the 10× Chromium Single Cell 3′ Library Kit according to the previous protocol (24). We generated single-cell suspensions from 2 different mouse kidney samples, WT and KSP^Cre/Dach1^fl/WT (Dach1 HZ). The single-cell sequencing libraries were sequenced on an Illumina HiSeq with a 2 × 150 paired-end kit. The sequencing reads were demultiplexed, aligned to the mouse genome (mm10), and processed to generate gene-cell data matrix using Cell Ranger, version 1.3. We sequenced 35,778 mouse kidney cells. We used the Seurat 3 package to integrate single-cell RNA-Seq data of 2 mouse kidney samples (46). HARMONY(47) was used for batch-effect correction and a cell type–specific gene expression signature was generated using the findmarker function in SEURAT.

Murine primary tubular epithelial cell and human podocyte culture. Three- to four-week-old mice were used to harvest kidney for primary tubular epithelial cell (PTEC) culture. Cells were isolated using 2 mg/ml collagenase I (Worthington Biochemical Product, catalog CLS-1) digestion for 30 minutes at 37°C with gentle stirring. Cells were then filtered through the 100 μm mesh to isolate single cells. Cell suspension was cultured in RPMI 1640 (Corning, catalog 10–040-CM) supplemented with 10% fetal bovine serum (Atlanta Biologicals, catalog S11950), 20 ng/ml EGF (Peprotech, catalog AF-100–15), 1 × ITS (Gibco, Thermo Fisher Scientific, catalog 51500–056), and 1% penicillin-streptomycin (Corning, catalog 30–002-CI) at 5% CO₂ and 37°C. To lower Dach1 expression, cells were infected with Ad5CMV-EGFP (Ad-GFP) or Ad5CMVCre-EGFP (Ad-Cre-EGFP) (University of Iowa Gene Transfer Vector Core, Iowa City, Iowa, USA) at 4 × 10¹⁰ plaque forming units/ml for 24 hours. Infection efficiency was estimated under fluorescence microscope by the presence of GFP-positive cells.

Human podocytes were cultured as described earlier (48). DACH1 knockdown was performed using lentiviral expression vectors (pGFP-C-shLenti) carrying shRNAs targeting DACH1 purchased from Origene (TL313572).

ChIP-qPCR. ChIP assay was performed according to the manufacturer’s protocol (MAGnify Chromatin Immunoprecipitation System, Invitrogen, 492024). Briefly, kidney tubular cells from WT mice were cultured in a 15 cm dish and directly crosslinked with 1% formaldehyde for 10 minutes at room temperature (RT). The reaction was stopped by incubating with glycine for 10 minutes at RT. Cells were lysed with lysis buffer to prepare nuclei. Sonication was performed using Bioruptor UCD300 with high power for 40 cycles (1 cycle; 30seconds on, 30 seconds off) at 4°C to target crosslinked DNA length between 100 and 500 bp. The size of the sonicated DNA was confirmed on 1% agarose gel. Anti-DACH1 antibody was coupled with Dynabeads protein A/G at 4°C for 1 hour. ChIP was performed using 4 μg anti-DACH1 antibody and 4 μg anti-rabbit IgG at 4°C overnight. After reverse crosslinking with proteinase K, DNA was purified with DNA purification magnetic beads and eluted with DNA elution buffer. The sequences of ChIP-qPCR primer were designed using the mouse kidney snATAC-Seq database. The primer sequences for ChIP-qPCR are shown in Supplemental Table 9. The percentage of input was calculated using signals obtained from ChIP and input samples.

CRISPR/Cas9 genomic deletions. HEK293 cells stably expressing Cas9 were a gift of Liling Wang and Lele Song (University of Pennsylvania). Guide RNAs were designed by CRISPOR software and cloned into pLKO.sgRNA plasmid. The plasmids were transfected into Cas9 expressing HEK293 using lipofectamine 3000 at 70% confluent according to the manufacturer’s instructions. Cells were collected 48 hours after transfection, and RNA was extracted using TRIzol. Gene expression was quantified by QRT-PCR. The full sequence of the guide RNAs is included in Supplemental Table 10. Target genomic region deletion was confirmed by Sanger sequencing.

Statistics. Data are expressed as mean ± SEM. Statistical significance was assessed by a 2-tailed Student’s t test for 2-group comparison or 1-way ANOVA for multigroup comparisons followed by Tukey’s post hoc test for subgroup comparisons. P < 0.05 was considered significant.

Accession codes. RNA-Seq data were deposited in the NCBI’s Gene Expression Omnibus database (GEO GSE115098). Mouse kidney single-cell gene expression data were also deposited in the GEO database (GSE107585). In addition, that data can be viewed on the Susztak Lab website (http://susztaklab.com/VisCello/). Mouse kidney snATAC-Seq data were deposited in the GEO database (GSE157079) and can be viewed on the Susztak Lab website (http://susztaklab.com/developing_adult_kidney/igv/). The precomputed human kidney eQTL data are available at the Susztak Lab website (http://www.susztaklab.com/eQTLci/download.php) and in the GEO database (GSE115098). The human kidney single-nuclei ATAC-Seq data are available at the Susztak Lab website (http://www.susztaklab.com/human_kidney/igv/.).