Genetic associations at 53 loci highlight cell types and biological pathways relevant for kidney function.

Reduced glomerular filtration rate defines chronic kidney disease and is associated with cardiovascular and all-cause mortality. We conducted a meta-analysis of genome-wide association studies for estimated glomerular filtration rate (eGFR), combining data across 133,413 individuals with replication in up to 42,166 individuals. We identify 24 new and confirm 29 previously identified loci. Of these 53 loci, 19 associate with eGFR among individuals with diabetes. Using bioinformatics, we show that identified genes at eGFR loci are enriched for expression in kidney tissues and in pathways relevant for kidney development and transmembrane transporter activity, kidney structure, and regulation of glucose metabolism. Chromatin state mapping and DNase I hypersensitivity analyses across adult tissues demonstrate preferential mapping of associated variants to regulatory regions in kidney but not extra-renal tissues. These findings suggest that genetic determinants of eGFR are mediated largely through direct effects within the kidney and highlight important cell types and biological pathways.

Chronic kidney disease (CKD) is a global public health problem123, and is associated with an increased risk for cardiovascular disease, all-cause mortality and end-stage renal disease45. Few new therapies have been developed to prevent or treat CKD over the past two decades16, underscoring the need to identify and understand the underlying mechanisms of CKD.

Prior genome-wide association studies (GWAS) have identified multiple genetic loci associated with CKD and estimated glomerular filtration rate (eGFR), a measure of the kidney's filtration ability that is used to diagnose and stage CKD789101112131415. Subsequent functional investigations point towards clinically relevant novel mechanisms in CKD that were derived from initial GWAS findings16, providing proof of principle that locus discovery through large-scale GWAS efforts can translate to new insights into CKD pathogenesis.

To identify additional genetic variants associated with eGFR and guide future experimental studies of CKD-related mechanisms, we have now performed GWAS meta-analyses in up to 133,413 individuals, more than double the sample size of previous studies. Here we describe multiple novel genomic loci associated with kidney function traits and provide extensive locus characterization and bioinformatics analyses, further delineating the physiologic basis of kidney function.

Results

Stage 1 discovery analysis

We analysed associations of eGFR based on serum creatinine (eGFRcrea), cystatin C (eGFRcys, an additional, complementary biomarker of renal function) and CKD (defined as eGFRcrea <60 ml min−1 per 1.73 m2) with ∼2.5 million autosomal single-nucleotide polymorphisms (SNPs) in up to 133,413 individuals of European ancestry from 49 predominantly population-based studies (Supplementary Table 1). Results from discovery GWAS meta-analysis are publicly available at http://fox.nhlbi.nih.gov/CKDGen/. We performed analyses in each study sample in the overall population and stratified by diabetes status, since genetic susceptibility to CKD may differ in the presence of this strong clinical CKD risk factor. Population stratification did not impact our results as evidenced by low genomic inflation factors in our meta-analyses, which ranged from 1.00 to 1.04 across all our analyses (Supplementary Fig. 1).

In addition to confirming 29 previously identified loci789 (Fig. 1a; Supplementary Table 2), we identified 48 independent novel loci (Supplementary Table 3) where the index SNP, defined as the variant with the lowest P value in the region, had an association P value <1.0 × 10−6. Of these 48 novel SNPs, 21 were genome-wide significant with P values <5.0 × 10−8. Overall, 43 SNPs were identified in association with eGFRcrea (nine in the non-diabetes sample), one with eGFRcys and four with CKD, as reported in Supplementary Table 3. Manhattan plots for CKD, eGFRcys and eGFRcrea in diabetes are shown in Fig. 1b,c and Supplementary Fig. 2, respectively.

Figure 1

Discovery stage genome-wide association analysis.

Manhattan plots for eGFRcrea, CKD and eGFRcys. Previously reported loci are highlighted in light blue (grey labels). (a) Novel loci uncovered for eGFRcrea in the overall and in the non-diabetes groups are highlighted in blue and green, respectively. (b) Results from CKD analysis with highlighted known and novel loci for eGFRcrea. (c) Results from eGFRcys with highlighted known and novel loci for eGFRcrea and known eGFRcys loci.

Stage-2 replication

Novel loci were tested for replication in up to 42,166 additional European ancestry individuals from 15 studies (Supplementary Table 1). Of the 48 novel candidate SNPs submitted to replication, 24 SNPs demonstrated a genome-wide significant combined stage 1 and 2 P value <5.0 × 10−8 (Table 1). Of these, 23 fulfilled additional replication criteria (q-value <0.05 in stage 2). Only rs6795744 at the WNT7A locus demonstrated suggestive replication (P value <5.0 × 10−8, q-value >0.05). Because serum creatinine is used to estimate eGFRcrea, associated genetic loci may be relevant to creatinine production or metabolism rather than kidney function per se. For this reason, we contrasted associations of eGFRcrea versus eGFRcys, the latter estimated from an alternative and creatinine-independent biomarker of GFR (Supplementary Fig. 3; Supplementary Table 4). The majority of loci (22/24) demonstrated consistent effect directions of their association with both eGFRcrea and eGFRcys.

Table 1

The 24 novel SNPs associated with eGFRcrea in European ancestry individuals.

SNP ID*	Chr.	Position (bp)†	Locus name‡	Effect/Non effect allele (EAF)	SNP function§	Stage 1 (discovery)||		Stage 2 (replication)		Combined analysis¶
						Beta	P value	Beta	q-value	Beta	P value#	I2 %**
The eight loci whose smallest P value was observed in the ‘no diabetes' group
rs3850625	1	201,016,296	CACNA1S	A/G (0.12)	Exonic, nonsyn. SNV	0.0080	2.55E−09	0.0071	5.46E−03	0.0083	6.82E−11	0
rs2712184	2	217,682,779	IGFBP5	A/C (0.58)	Intergenic	−0.0049	1.65E−08	−0.0055	2.06E−03	−0.0053	1.33E−10	0
rs9682041	3	170,091,902	SKIL	T/C (0.87)	Intronic	−0.0067	1.36E−07	−0.0046	2.33E−02	−0.0068	2.58E−08	2
rs10513801	3	185,822,353	ETV5	T/G (0.87)	Intronic	0.0070	3.80E−09	0.0046	1.79E−02	0.0072	1.03E−09	0
rs10994860	10	52,645,424	A1CF	T/C (0.19)	UTR5	0.0075	1.00E−11	0.0061	5.46E−03	0.0077	1.07E−12	2
rs163160	11	2,789,955	KCNQ1	A/G (0.82)	Intronic	0.0067	9.02E−09	0.0050	9.89E−03	0.0065	2.26E−09	14
rs164748	16	89,708,292	DPEP1	C/G (0.53)	Intergenic	0.0047	9.92E−09	0.0019	4.19E−02	0.0046	1.95E−08	17
rs8091180	18	77,164,243	NFATC1	A/G (0.56)	Intronic	−0.0054	1.43E−08	−0.0052	5.46E−03	−0.0060	1.28E−09	0
The 16 loci whose smallest P value was observed in the ‘overall' group
rs12136063	1	110,014,170	SYPL2	A/G (0.70)	Intronic	0.0049	2.33E−07	0.0028	2.31E−02	0.0045	4.71E−08	0
rs2802729	1	243,501,763	SDCCAG8	A/C (0.43)	Intronic	−0.0050	7.37E−08	−0.0029	2.05E−02	−0.0046	2.20E−08	9
rs4667594	2	170,008,506	LRP2	A/T (0.53)	Intronic	−0.0045	2.37E−07	−0.0043	5.62E−03	−0.0044	3.52E−08	4
rs6795744‡‡	3	13,906,850	WNT7A	A/G (0.15)	Intronic	0.0071	9.60E−09	0.0019	5.15E−02	0.0060	3.33E−08	18
rs228611	4	103,561,709	NFKB1	A/G (0.47)	Intronic	−0.0055	4.66E−10	−0.0060	8.91E−04	−0.0056	3.58E−12	3
rs7759001	6	27,341,409	ZNF204	A/G (0.76)	ncRNA intronic	−0.0053	2.64E−07	−0.0045	9.10E−03	−0.0051	1.75E−08	0
rs10277115	7	1,285,195	UNCX	A/T (0.23)	Intergenic	0.0095	1.05E−10	0.0079	9.03E−04	0.0090	8.72E−14	0
rs3750082	7	32,919,927	KBTBD2	A/T (0.33)	Intronic	0.0049	2.52E−07	0.0031	1.96E−02	0.0045	3.22E−08	2
rs6459680	7	156,258,568	RNF32	T/G (0.74)	Intergenic	−0.0065	1.96E−10	−0.0019	4.62E−02	−0.0055	1.07E−09	0
rs4014195	11	65,506,822	AP5B1	C/G (0.64)	Intergenic	0.0061	2.19E−11	0.0034	1.42E−02	0.0055	1.10E−11	0
rs10491967	12	3,368,093	TSPAN9	A/G (0.10)	Intronic	−0.0092	3.03E−10	−0.0106	3.93E−04	−0.0095	5.18E−14	0
rs7956634	12	15,321,194	PTPRO	T/C (0.81)	Intronic	−0.0068	2.46E−09	−0.0069	1.51E−03	−0.0068	7.17E−12	0
rs1106766	12	57,809,456	INHBC	T/C (0.22)	Intergenic	0.0062	4.67E−08	0.0058	8.79E−03	0.0061	2.41E−09	9
rs11666497	19	38,464,262	SIPA1L3	T/C (0.18)	Intronic	−0.0064	8.58E−08	−0.0041	1.53E−02	−0.0058	4.25E−08	24
rs6088580	20	33,285,053	TP53INP2	C/G (0.47)	Intergenic	−0.0055	7.17E−10	−0.0027	2.31E−02	−0.0049	1.79E−09	0
rs17216707	20	52,732,362	BCAS1	T/C (0.79)	Intergenic	−0.0084	5.96E−13	−0.0051	6.69E−03	−0.0077	8.83E−15	1

bp, basepairs; Chr, chromosome; EAF, effect allele frequency; eGFRcrea, eGFR based on serum creatinine; GWAS, genome-wide association studies; SNP, single-nucleotide polymorphism; UTR, untranslated region.

*SNPs are grouped by the stratum where the smallest P value in the discovery and combined analysis was observed. In the ‘no diabetes' group, sample size/number of studies were equal to 118,448/45, 36,433/13 and 154,881/58, in the discovery, replication and combined analyses, respectively. In the ‘overall' group, the numbers for the three analyses were equal to 133,413/48, 42,116/14 and 175,579/62, respectively.

†On the basis of RefSeq genes (build 37).

‡Conventional locus name based on relevant genes in the region as identified by bioinformatic investigation (Supplementary Table 12) or closest gene. A complete overview of the genes in each locus is given in the regional association plots (Supplementary Fig. 4).

§SNP function is derived from NCBI RefSeq genes and may not correspond to the named gene.

||Twice genomic-control (GC) corrected P value from discovery GWAS meta-analysis: at the individual study level and after the meta-analysis.

¶For random-effect estimate, see Supplementary Table 4.

#P value of the meta-analysis of the twice GC-corrected discovery meta-analysis results and replication studies.

**Between-study heterogeneity, as assessed by the I2. Q statistic P value >0.05 for all SNPs, except rs11666497 (SIPA1L3, P value=0.04).

‡‡For this SNP, the conditions for replication were not all met (q-value >0.05 in the replication stage).

Association plots of the 24 newly identified genomic regions that contain a replicated or suggestive index SNP appear in Supplementary Fig. 4. The odds ratio for CKD for each of the novel loci ranged from 0.93 to 1.06 (Supplementary Table 4). As evidenced by the relatively small effect sizes, the proportion of phenotypic variance of eGFRcrea explained by all new and known loci was 3.22%: 0.81% for the newly uncovered loci and 2.41% for the already known loci.

Associations stratified by diabetes and hypertension status

The effects of the 53 known and novel loci in individuals with (stage 1+stage 2 n=16,477) and without (stage 1+stage 2 n=154,881) diabetes were highly correlated (correlation coefficient: 0.80; 95% confidence interval: 0.67, 0.88; Supplementary Fig. 5) and of similar magnitude (Fig. 2; Supplementary Table 5), suggesting that identification of genetic loci in the overall population may also provide insights into loci with potential importance among individuals with diabetes. The previously identified UMOD locus showed genome-wide significant association with eGFRcrea among those with diabetes (Supplementary Fig. 2; rs12917707, P value=2.5 × 10−8), and six loci (NFKB1, UNCX, TSPAN9, AP5B1, SIPA1L3 and PTPRO) had nominally significant associations with eGFRcrea among those with diabetes. Of the previously identified loci, 13 demonstrated nominal associations among those with diabetes, for a total of 19 loci associated with eGFRcrea in diabetes.

Figure 2

Association eGFRcrea loci in subjects with and without diabetes.

Novel (a) and known (b) loci were considered. Displayed are effects and their 95% confidence intervals on ln(eGFRcrea). Results are sorted by increasing effects in the diabetes group. The majority of loci demonstrated similar effect sizes in the diabetes as compared with non-diabetes strata. SNP-specific information and detailed sample sizes are reported in Supplementary Table 5.

Exploratory comparison of the association effect sizes in subjects with and without hypertension based on our previous work7 showed that novel and known loci are also similarly associated with eGFRcrea among individuals with and without hypertension (Supplementary Fig. 6).

Tests for SNP associations with related phenotypes

We tested for overlap with traits that are known to be associated with kidney function in the epidemiologic literature by investigating SNP associations with systolic and diastolic blood pressure17, myocardial infarction18, left ventricular mass19, heart failure20, fasting glucose21 and urinary albumin excretion (CKDGen Consortium, personal communication). We observed little association of the 24 novel SNPs with other kidney function-related traits, with only 2 out of 165 tests reaching the Bonferroni significance level of 0.0003 (see Methods and Supplementary Table 6).

To investigate whether additional traits are associated with the 24 new eGFR loci, we queried the NHGRI GWAS catalog (www.genome.gov). Overall, nine loci were previously identified in association with other traits at a P value of 5.0 × 10−8 or lower (Supplementary Table 7), including body mass index (ETV5) and serum urate (INHBC, A1CF and AP5B1).

Trans-ethnic analyses

To assess the generalizability of our findings across ethnicities, we evaluated the association of the 24 newly identified loci with eGFRcrea in 16,840 participants of 12 African ancestry population studies (Supplementary Table 8) and in up to 42,296 Asians from the AGEN consortium11 (Supplementary Table 9). Seven SNPs achieved nominal direction-consistent significance (P<0.05) in AGEN, and one SNP was significant in the African ancestry meta-analysis (Supplementary Table 9). Random-effect meta-analysis showed that 12 loci (SDCCAG8, LRP2, IGFBP5, SKIL, UNCX, KBTBD2, A1CF, KCNQ1, AP5B1, PTPRO, TP53INP2 and BCAS1) had fully consistent effect direction across the three ethnic groups (Supplementary Fig. 7), suggesting that our findings can likely be generalized beyond the European ancestry group.

To identify additional potentially associated variants and more formally evaluate trans-ethnic heterogeneity of the loci identified through meta-analysis in European ancestry populations, we performed a trans-ethnic meta-analysis22, combining the 12 African ancestry studies with the 48 European Ancestry studies used in the discovery analysis of eGFRcrea. Of the 24 new loci uncovered for eGFRcrea, 15 were also genome-wide significant in the trans-ethnic meta-analysis (defined as log10 Bayes Factor >6, Supplementary Table 10), indicating that for most of these loci, there is little to no allelic effect heterogeneity across the two ethnic groups. No additional loci were significantly associated with log10 Bayes Factor >6.

Bioinformatic and functional characterization of new loci

We used several techniques to prioritize and characterize genes underlying the identified associations, uncover connections between associated regions, detect relevant tissues and assign functional annotations to associated variants. These included expression quantitative trait loci (eQTL) analyses, pathway analyses, DNAse I hypersensitivity site (DHS) mapping, chromatin mapping, manual curation of genes in each region and zebrafish knockdown.

eQTL analysis

We performed eQTL analysis using publically available eQTL databases (see Methods). These analyses connected novel SNPs to transcript abundance of SYPL2, SDCCAG8, MANBA, KBTBD2, PTPRO and SPATA33 (C16orf55), thereby supporting these as potential candidate genes in the respective associated regions (Supplementary Table 11).

Pathway analyses

We used a novel method, Data-driven Expression Prioritized Integration for Complex Traits (DEPICT)23, to prioritize genes at associated loci, to test whether genes at associated loci are highly expressed in specific tissues or cell types and to test whether specific biological pathways and gene sets are enriched for genes in associated loci. On the basis of all SNPs with eGFRcrea association P values <10−5 in the discovery meta-analysis, representing 124 independent regions, we identified at least one significantly prioritized gene in 49 regions, including in 9 of the 24 novel genome-wide significant regions (Supplementary Table 12). Five tissue and cell type annotations were enriched for expression of genes from the associated regions, including the kidney and urinary tract, as well as hepatocytes and adrenal glands and cortex (Fig. 3a; Supplementary Table 13). Nineteen reconstituted gene sets showed enrichment of genes mapping into the associated regions at a permutation P value <10−5 (Supplementary Table 14; Fig. 4), highlighting processes related to renal development, kidney transmembrane transporter activity, kidney and urogenital system morphology, regulation of glucose metabolism, as well as specific protein complexes important in renal development.

Figure 3

Bioinformatic analysis of eGFR-associated SNPs.

Connection of eGFR-associated SNPs to gene expression and variant function across a variety of tissues, pathways and regulatory marks was considered. (a) The DEPICT method shows that implicated eGFR-associated genes are highly expressed in particular tissues, including kidney and urinary tract. Shown are permutation test P values (see Methods). (b) Enrichment of eGFRcrea-associated SNPs in DHS according to discovery P value threshold. SNPs from the eGFR discovery genome-wide scan meeting a series of P value thresholds in the range 10−4–10−16 preferentially map to DHSs, when compared with a set of control SNPs, in 6 of 123 cell types. Represented are main effects odds ratios from a logistic mixed effect model. Cell types indicated with coloured lines had nominally significant enrichment (* indicate P values <0.05) at the P value <10−16 threshold and/or were derived from renal tissues (H7esDiffa2d: H7 embryonic stem cells, differentiated 2 days with BMP4, activin A and bFGF; Hae, amniotic epithelial cells; Hrce, renal cortical epithelial cells; Hre, renal epithelial cells; Hrgec, renal glomerular endothelial cells; Rptec, renal proximal tubule epithelial cells; Saec, small airway epithelial cells). (c) ENCODE/Chromatin ChIP-seq mapping: known and replicated novel eGFRcrea-associated SNPs and their perfect proxies were annotated based on genomic location using chromatin annotation maps from different tissues including adult kidney epithelial cells. P values from Fishers' exact tests for 2 × 2 tables are reported (significance level=5.6 × 10−3, see Methods). There is significant enrichment of variants mapping to enhancer regions specifically in kidney but not other non-renal tissues.

Figure 4

Gene set overlap analysis.

The 19 reconstituted gene sets with P value<10−5 were considered. Their overlap was estimated by computing the pairwise Pearson correlation coefficient ρ between each pair of gene sets followed by discretization into one of three bins: 0.3≤ρ<0.5, low overlap; 0.5≤ρ<0.7, medium overlap; ρ⩾0.7, high overlap. Overlap is shown by edges between gene set nodes and edges representing overlap corresponding to ρ<0.3 are not shown. The network was drawn with Cytoscape48.

DNase I hypersensitivity and H3K4m3 chromatin mark analyses

To evaluate whether eGFRcrea-associated SNPs map into gene regulatory regions and to thereby gain insight into their potential function, we evaluated the overlap of independent eGFRcrea-associated SNPs with P values <10−4 (or their proxies) with DHSs using publicly available data from the Epigenomics Roadmap Project and ENCODE for 123 cell types (see Methods). DHSs mark accessible chromatin regions where transcription may occur. Compared with a set of control SNPs (see Methods), eGFRcrea-associated SNPs were significantly more likely to map to DHS in six specific tissues or cell types (Fig. 3b), including adult human renal cortical epithelial cells, adult renal proximal tubule epithelial cells, H7 embryonic stem cells (differentiated 2 days), adult human renal epithelial cells, adult small airway epithelial cells and amniotic epithelial cells. No significant enrichment was observed for adult renal glomerular endothelial cells, the only other kidney tissue evaluated.

Next, we analysed the overlap of the same set of SNPs with H3K4me3 chromatin marks, promoter-specific histone modifications associated with active transcription24, in order to gather more information about cell-type specific regulatory potential of eGFRcrea-associated SNPs. Comparing 33 available adult-derived cell types, we found that eGFRcrea-associated SNPs showed the most significant overlap with H3K4me3 peaks in adult kidney (P value=0.0029), followed by liver (P value=0.0117), and rectal mucosa (P value=0.0445). Taken together, these findings are suggestive of cell-type-specific regulatory roles for eGFR loci, with greatest specificity for the kidney.

Chromatin annotation maps

In addition to assessing individual regulatory marks separately, we annotated the known and replicated novel SNPs, as well as their perfect proxies in a complementary approach. Chromatin annotation maps were generated integrating >10 epigenetic marks from cells derived from adult human kidney tissue and a variety of non-renal tissues from the ENCODE project (see Methods). The proportion of variants to which a function could be assigned was significantly higher when using chromatin annotation maps from renal tissue compared with using maps that investigated the same epigenetic marks in other non-renal tissues (Fig. 3c), again indicating that eGFRcrea associated SNPs are, or tag, kidney-specific regulatory variants. The difference between kidney and non-renal tissues was particularly evident for marks that define enhancers: the proportion of SNPs mapping to weak and strong enhancer regions in the kidney tissue was higher than in all non-kidney tissues (Fishers' exact test P values from 3.1 × 10−3 to 7.9 × 10−6, multiple testing threshold α=5.6 × 10−3).

Functional characterization of new loci

To prioritize genes for functional studies, we applied gene prioritization algorithms including GRAIL25, DEPICT and manual curation of selected genes in each region (Supplementary Table 12). For each region, gene selection criteria were as follows: (1) either GRAIL P value <0.05 or DEPICT false discovery rate (FDR) <0.05; (2) the effect of a given allele on eGFRcrea and on eGFRcys was direction-consistent and their ratio was between 0.2 and 5 (to ensure relative homogeneity of the beta coefficients); (3) nearest gene if the signal was located in a region containing a single gene. Using this approach, NFKB1, DPEP1, TSPAN9, NFATC1, WNT7A, PTPRO, SYPL2, UNCX, KBTBD2, SKIL and A1CF were prioritized as likely genes underlying effects at the new loci (Supplementary Table 12).

We investigated the role of these genes during vertebrate kidney development by examining the functional consequences of gene knockdown in zebrafish embryos utilizing antisense morpholino oligonucleotide (MO) technology. After knockdown, we assessed the expression of established renal markers pax2a (global kidney), nephrin (podocytes) and slc20a1a (proximal tubule) at 48 hours post fertilization by in situ hybridization12. In all cases, morphant embryos did not display significant gene expression defects compared with controls (Supplementary Table 15).

Discussion

We identified 24 new loci in association with eGFR and confirmed 29 previously identified loci. A variety of complementary analytic, bioinformatic and functional approaches indicate enrichment of implicated gene products in kidney and urinary tract tissues. A greater proportion of the lead SNPs or their perfect proxies map into gene regulatory regions, specifically enhancers, in adult renal tissues compared with non-renal tissues. In addition to the importance in the adult kidney, our results indicate a role for kidney function variants during development.

We extend our previous findings, as well as those from other groups78910111213 by identifying >50 genomic loci for kidney function, many of which were not previously known to be connected to kidney function and disease. Using a discovery data set that is nearly double in size to our prior effort7, we are now able to robustly link associated SNPs to kidney-specific gene regulatory function. Our work further exemplifies the continued value of increasing the sample size of GWAS meta-analyses to uncover additional loci and gain novel insights into the mechanisms underlying common phenotypes26.

There are several messages from our work. First, many of the genetic variants associated with eGFR appear to affect processes specifically within the kidney. The kidney is a highly vascular and metabolically active organ that receives 20% of all cardiac output, contains an extensive endothelium-lined capillary network, and is sensitive to ischaemic and toxic injury. As a result, hypertension, cardiovascular diseases and diabetes each affect renal hemodynamics and contribute to kidney injury. However, many of the eGFR-associated SNPs in our GWAS could be assigned gene regulatory function specifically in the kidney and its epithelial cells, but not in human glomerular endothelial cells or the general vasculature. In addition, variants associated with eGFR were not associated with vascular traits, such as blood pressure or myocardial infarction. Taken together, these findings suggest that genetic determinants of eGFR may be mediated largely through direct effects within the kidney.

Second, despite the specificity related to renal processes, we also identified several SNPs that are associated with eGFR in diabetes, and our pathway analyses uncovered gene sets associated with glucose transporter activity and abnormal glucose homeostasis. Uncovering bona fide genetic loci for diabetic CKD has been difficult. We have now identified a total of 19 SNPs that demonstrate at least nominal association with eGFR in diabetes. The diabetes population is at particularly high risk of CKD, and identifying kidney injury pathways may help develop new treatments for diabetic CKD.

Finally, even though CKD is primarily a disease of the elderly, our pathway enrichment analyses highlight developmental processes relevant to the kidney and the urogenital tract. Kidney disease has been long thought to have developmental origins, in part related to early programming (Barker hypothesis)27, low birth weight, nephron endowment and early growth and early-life nutrition28. Our pathway enrichment analyses suggest that developmental pathways such as placental morphology, kidney weight and embryo size, as well as protein complexes of importance in renal development may in part contribute to the developmental origins of CKD.

A limitation of our work is that causal variants and precise molecular mechanisms underlying the observed associations were not identified and will require additional experimental follow-up projects. Our attempt to gain insights into potentially causal genes through knockdown in zebrafish did not yield any clear CKD candidate gene, although the absence of a zebrafish phenotype upon gene knockdown does not mean that the gene cannot be the one underlying the observed association signal in humans. Finally, our conclusions that eGFRcrea-associated SNPs regulate the expression of nearby genes specifically in kidney epithelial cells are based on DHSs, H3K4me3 chromatin marks and chromatin annotation maps. Since these analyses rely mostly on variant positions, additional functional investigation such as luciferase assay that assess transcriptional activity more directly are likely to gain additional insights into the variants' mechanism of action.

The kidney specificity for loci we identified may have important translational implications, particularly since our DHS and chromatin annotation analyses suggest that at least a set of gene regulatory mechanisms is important in the adult kidney. Kidney-specific pathways are important for the development of novel therapies to prevent and treat CKD and its progression with minimal risk of toxicity to other organs. Finally, the biologic insights provided by these new loci may help elucidate novel mechanisms and pathways implicated not only in CKD but also of kidney function in the physiological range.

In conclusion, we have confirmed 29 genomic loci and identified 24 new loci in association with kidney function that together highlight target organ-specific regulatory mechanisms related to kidney function.

Methods

Overview

This was a collaborative meta-analysis with a distributive data model. Briefly, an analysis plan was created and circulated to all participating studies. Studies then uploaded study-specific data centrally; files were cleaned, and a specific meta-analysis for each trait was performed. Details regarding each step are provided below. All participants in all discovery and replication studies provided informed consent. Each study had its research protocol approved by the local ethics committee.

Phenotype definitions

Serum creatinine was measured in each discovery and replication study as described in Supplementary Tables 16 and 17, and statistically calibrated to the US nationally representative National Health and Nutrition Examination Study data in all studies to account for between-laboratory variation92930. eGFRcrea was estimated using the four-variable Modification of Diet in Renal Disease Study Equation. Cystatin C, an alternative biomarker for kidney function, was measured in a sub-set of participating studies. eGFRcys was estimated as 76.7 × (serum cystatin C)−1.19 (ref. 31). eGFRcrea and eGFRcys values <15 ml min per 1.73 m2 were set to 15, and those >200 were set to 200 ml min−1 per 1.73 m2. CKD was defined as eGFRcrea <60 ml  min−1 per 1.73 m2.

Diabetes was defined as fasting glucose ⩾126 mg dl−1, pharmacologic treatment for diabetes or by self-report. In all studies, diabetes and kidney function were assessed at the same point in time.

Genotypes

Genotyping was conducted in each study as specified in Supplementary Tables 18 and 19. After applying appropriate quality filters, 45 studies used markers of highest quality to impute ∼2.5 million SNPs, based on European-ancestry haplotype reference samples (HapMap II CEU). Four studies based their imputation on the 1000 Genomes Project data. Imputed genotypes were coded as the estimated number of copies of a specified allele (allelic dosage).

Genome-wide association analysis

By following a centralized analysis plan, each study regressed sex- and age-adjusted residuals of the logarithm of eGFRcrea or eGFRcys on SNP dosage levels. Logistic regression of CKD status was performed on SNP dosage levels adjusting for sex and age. For all traits, adjustment for appropriate study-specific features, including study site and genetic principal components was included in the regression and family-based studies appropriately accounted for relatedness.

Stage 1 discovery meta-analysis

GWAS of eGFRcrea were contributed by 48 studies (total sample size, N=133,413); 45 studies contributed GWAS data for the non-diabetes subgroup (N=118,448) and 39 for the diabetes group (N=11,522). GWAS of CKD were comprised by 43 studies, for a total sample size of 117,165, including 12,385 CKD cases. GWAS of eGFRcys were comprised by 16 studies for a total sample size of 32,834. All GWAS files underwent quality control using the GWAtoolbox package32 in R, before including them into the meta-analysis. Genome-wide meta-analysis was performed with the software METAL33, assuming fixed effects and using inverse-variance weighting. The genomic inflation factor λ was estimated for each study as the ratio between the median of all observed test statistics (b/s.e.)2 and the expected median of a χ2 with 1 degree of freedom, with b and s.e. representing the effect of each SNP on the phenotype and its standard error, respectively34. Genomic-control (GC) correction was applied to P values and s.e.'s in case of λ>1 (first GC correction). SNPs with an average minor allele frequency (MAF) of ⩾0.01 were used for the meta-analysis. To limit the possibility of false positives, after the meta-analysis, a second GC correction on the aggregated results was applied. Between-study heterogeneity was assessed through the I2 statistic.

After removing SNPs with MAF of <0.05 and which were available in <50% of the studies, SNPs with a P value of ≤10−6 were selected and clustered into independent loci through LD pruning based on an r2 of ≤0.2 within a window of ±1 MB to each side of the index SNP. After removing loci containing variants that have been previously replicated at a P value of 5.0 × 10−8 (refs 7, 8), the SNP with the lowest P value within each locus was selected for replication (‘index SNP'). If a SNP had an association P value of ≤10−6 with more than one trait, the trait where the SNP had the lowest P value was selected as discovery trait/stratum. Altogether, this resulted in 48 SNPs: 34 from eGFRcrea, 9 from eGFRcrea among those without diabetes, 4 from CKD and 1 from eGFRcys.

Stage 2 replication analysis

In silico replication analysis for any of the studied traits was carried out using eight independent studies whose genotyping platforms are provided in Supplementary Table 19. De novo genotyping was performed in seven additional studies (N=22,850 individuals) of European ancestry (Supplementary Table 20), including the Bus Santé, ESTHER, KORA-F3 (subset of F3 without GWAS), KORA-F4 (subset of F4 without GWAS), Ogliastra Genetic Park (OGP, without Talana whose GWAS was included in the discovery analysis), SAPHIR and SKIPOGH studies (Supplementary Table 20). Summarizing all in silico and de novo replication studies (Supplementary Table 1), replication data for eGFRcrea were contributed by 14 studies (total sample size=42,166), which also contributed eGFRcrea results from non-diabetes (13 studies, N=36,433) and diabetes samples (13 studies, N=4,955). Thirteen studies contributed replication data on CKD (N=33,972; 4,245 CKD cases; studies with <50 CKD cases were excluded) and five on eGFRcys (N=14,930).

Association between eGFRcrea, CKD and eGFRcys and each of the 48 SNPs in the replication studies was assessed using the same analysis protocol detailed for the discovery studies above. Quality control of the replication files was performed with the same software as described above.

We performed a combined fixed-effect meta-analysis of the double-GC corrected results from the discovery meta-analysis and the replication studies, based on inverse-variance weighting. The total sample size in the combined analysis of eGFRcrea was 175,579 subjects (154,881 in the non-diabetes stratum and 16,477 in the diabetes stratum; the sum of these two sample sizes is smaller than the sample size of the overall analysis because some studies did not contribute both strata), 151,137 samples for CKD (16,630 CKD cases) and 47,764 for eGFRcys. Three criteria were used to ensure validity of novel loci declared as significant: (1) P value from the combined meta-analysis ≤5.0 × 10−8 in accordance with previously published guidelines35; (2) direction-consistent associations of the beta coefficients in stage 1 and stage 2 (one-sided P values were estimated to test for consistent effect direction with the discovery stage); (3) q-value <0.05 in the replication stage. Q-values were estimated using the package QVALUE36 in R. The tuning parameter lambda for the estimation of the overall proportion of true null hypotheses, π0, was estimated using the bootstrap method37. When the third criterion was not satisfied, the locus was declared ‘suggestive'.

Power analysis

With the sample size achieved in the combined analysis of stage 1 and stage 2 data, the power to assess replication at the canonical genome-wide significance level of 5.0 × 10−8 was estimated with the software QUANTO38 version 1.2.4, assuming the same MAF and effect size observed in the discovery sample. Power to replicate associations ranged from 87 to 100% for eGFRcrea associated SNPs (median=98%), from 72 to 96% for the CKD-associated SNPs, and was equal to 59% for the eGFRcys-associated SNP (Supplementary Table 3).

For all the 24 novel and 29 known SNPs, the difference between the SNP effect on eGFRcrea in the diabetes versus the non-diabetes groups was assessed by means of a two-sample t-test for correlated data at a significance level of 0.05. We used the following two-sample t-test for correlated data:

where bDM and bnonDM represent the SNP effects on log(eGFRcrea) in the two groups, s.e. is the standard error of the estimate and ρ(.) indicates the correlation between effects in the two groups, which was estimated as 0.044 by sampling 100,000 independent SNPs from our DM and nonDM GWAS, after removing known and novel loci associated with eGFRcrea. For a large sample size, as in our case, t follows a standard normal distribution.

A similar analysis was performed to compare results in subjects with and without hypertension, based on results from our previous work7. The correlation between the two strata was of 0.01.

Proportion of phenotypic variance explained

The percent of phenotypic variance explained by novel and known loci was estimated as , where is the coefficient of determination for each of the 53 individual SNPs associated with eGFRcrea uncovered to date (24 novel and 29 known ones), b is the estimated effect of the ith SNP on y, y corresponds to the sex- and age-adjusted residuals of the logarithm of eGFRcrea and var(SNPi)=2 × MAFSNPi × (1−MAFSNPi)39. Var(y) was estimated in the ARIC study and all loci were assumed to have independent effects on the phenotype.

Test for SNP associations with related traits

We performed evaluations of SNP association with results generated from consortia investigating other traits. Specifically, we evaluated systolic and diastolic blood pressure in ICBP17, myocardial infarction in CARDiOGRAM18, left ventricular mass19, heart failure20, the urinary albumin to creatinine ratio (CKDGen consortium, personal communication) and fasting plasma glucose in MAGIC21. In total, we performed 165 tests, corresponding to 7 traits tested for association against each of the 24 novel SNP, with the exception of myocardial infarction for which results from 3 SNPs were not available (Supplementary Table 6). Significance was evaluated at the Bonferroni corrected level of 0.05/165=0.0003.

Lookup of replicated loci in the NHGRI GWAS catalog

All replicated SNPs, as well as SNPs in LD (r2>0.2) within ±1 MB distance were checked for their association with other traits according to the NHGRI GWAS catalog40 (accessed April 14, 2014).

SNP assessments in other ethnic groups

We performed cross-ethnicity SNP evaluations in participants of African ancestry from a meta-analysis of African ancestry individuals and from participants of Asian descent from the AGEN consortium11.

African ancestry meta-analysis

We performed fixed-effect meta-analysis of the genome-wide association data from 12 African ancestry studies (Supplementary Table 8) with imputation to HapMap reference panel, based on inverse-variance weighting using METAL. Only SNPs with MAF ⩾0.01 and imputation quality r2⩾0.3 were considered for the meta-analysis. After meta-analysis, we removed SNPs with MAF <0.05 and which were available in <50% of the studies. Statistical significance was assessed at the standard threshold of 5.0 × 10−8. Genomic control correction was applied at both the individual study level before meta-analysis and after the meta-analysis.

Transethnic meta-analysis

We performed a trans-ethnic meta-analysis of GWAS data from cohorts of different ethnic backgrounds using MANTRA (Meta-Analysis of Trans-ethnic Association studies) software22. We combined the 48 European ancestry studies that contributed eGFRcrea, which were included in stage 1 discovery meta-analysis, and the 12 African ancestry studies mentioned above for a total sample size of 150,253 samples. We limited our analysis to biallelic SNPs with MAF ⩾0.01 and imputation quality r2⩾0.3. Relatedness between the 60 studies was estimated using default settings from up to 5.9 million SNPs. Only SNPs that were present in more than 25 European ancestry studies and 6 African ancestry studies (total sample size ⩾120,000) were considered after meta-analysis. Genome-wide significance was defined as a log10 Bayes' Factor (log10BF) ⩾6 (ref. 41).

Gene Relationships Across Implicated Loci (GRAIL)

To prioritize the gene(s) most likely to give rise to association signals in a given region, the software GRAIL was used25. The index SNP of all previously known kidney function associated regions, as well as the novel SNPs identified here was used as input, using the CEU HapMap (hg18 assembly) and the functional datasource text_2009_03, established before the publication of kidney function-related GWAS. Results from GRAIL were used to prioritize genes for follow-up functional work.

Expression quantitative trait loci analysis

We identified alias rsIDs and proxies (r2>0.8) for our index SNPs using SNAP software across 4 HapMap builds. SNP rsIDs and aliases were searched for primary SNPs and LD proxies against a collected database of expression SNP (eSNP) results. The collected eSNP results met criteria for statistical thresholds for association with gene transcript levels in their respective original analyses (for references see Supplementary Table 11). Correlation of selected eSNPs to the best eSNPs per transcript per expression quantitative trait loci (eQTL) data set were assessed by pairwise LD. All results are reported in Supplementary Table 11.

DEPICT analysis

In this work, we first used PLINK42 to identify independently associated SNPs using all SNPs with eGFRcrea association P values <10−5 (HapMap release 27 CEU data43; LD r2 threshold=0.01; physical kb threshold=1,000). We then used the DEPICT method23 to construct associated regions by mapping genes to independently associated SNPs if they overlapped or resided within LD (r2>0.5) distance of a given associated SNP. After merging overlapping regions and discarding regions that mapped within the major histocompatibility complex locus (chromosome 6, base pairs 20,000,000–40,000,000), 124 non-overlapping regions remained that covered a total of 320 genes. Finally, we ran the DEPICT software program on the 124 regions to prioritize genes that may represent promising candidates for experimental follow up studies, identify reconstituted gene sets that are enriched in genes from associated regions and therefore may provide insight into general kidney function biology, and identify tissue and cell-type annotations in which genes from associated regions are highly expressed. Specifically, for each tissue, the DEPICT method performs a t-test comparing the tissue-specific expression of eGFRcrea-associated genes and all other genes. Next, for each tissue, empirical enrichment P values are computed by repeatedly sampling random sets of loci (matched to the actual eGFRcrea loci by gene density) from the entire genome to estimate the empirical mean and s.d. of the enrichment statistic's null distribution. To visualize the nineteen reconstituted gene sets with P value <1e−5 (Fig. 4), we estimated their overlap by computing the pairwise Pearson correlation coefficient ρ between each pair of gene sets followed by discretization into one of three bins; 0.3≤ρ<0.5, low overlap; 0.5≤ρ<0.7, medium overlap; ρ⩾0.7, high overlap.

DNase I hypersensitivity analysis

The overlap of SNPs associated with eGFRcrea at P<10−4 with DHSs was examined using publically available data from the Epigenomics Roadmap Project and ENCODE. In all, DHS mappings were available for 123 mostly adult cells and tissues44 (downloaded from http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwDnase/). The analysis here pertains to DHS's defined as ‘broad' peaks, which were available as experimental replicates (typically duplicates) for the majority of cells and tissues.

SNPs from our stage 1 eGFRcrea GWAS meta-analysis were first clumped in PLINK42 in windows of 100 kb and maximum r2 of 0.1 using LD relationships from the 1,000 Genomes EUR panel (phase I, v3, 3/14/2012 haplotypes) using a series of P value thresholds (10−4, 10−6, 10−8, ... and 10−16). LD proxies of the index SNPs from the clumping procedure were then identified by LD tagging in PLINK with r2=0.8 in windows of 100 kb, again using LD relationships in the 1000G EUR panel, restricted to SNPs with MAF >1% and also present in the HapMap2 CEU population. A reference set of control SNPs was constructed using the same clumping and tagging procedures applied to NHGRI GWAS catalog SNPs (available at http://www.genome.gov/gwastudies/, accessed 13 March 2013) with discovery P values <5.0 × 10−8 in European populations. In total, there were 1,204 such reference SNPs after LD pruning. A small number of reference SNPs or their proxies overlapping with the eGFRcrea SNPs or their proxies were excluded. For each cell-type and P value threshold, the enrichment of eGFR SNPs (or their LD proxies) mapping to DHSs relative to the GWAS catalog reference SNPs (or their LD proxies) was expressed as an odds ratio from logistic mixed effect models that treated the replicate peak determinations as random effects (lme4 package in R). Significance for enrichment odds ratio was derived from the significance of beta coefficients for the main effects in the mixed models.

Interrogation of human kidney chromatin annotation maps

Different chromatin modification patterns can be used to generate tissue-specific chromatin-state annotation maps. These can serve as a valuable resource to discover regulatory regions and study their cell-type-specific distributions and activities, which may help with the interpretation especially of intergenic variants identified in association studies45. We therefore investigated the genomic mapping of the known and replicated novel index SNPs, as well as their perfect LD proxies (n=173, r2=1 for proxies) using a variety of resources, including chromatin maps generated from human kidney tissue cells (HKC-E cells). Chromatin immune-precipitation sequencing (ChIP-seq) data from human kidney samples were generated by NIH Roadmap Epigenomics Mapping Consortium46. Briefly, proximal tubule cells derived from an adult human kidney were collected and cross-linked with 1% formaldehyde. Subsequently, ChIP-seq was conducted using whole-cell extract from adult kidney tissue as the input (GSM621638) and assessing the following chromatin marks: H3K36me3 (GSM621634), H3K4me1 (GSM670025), H3K4me3 (GSM621648), H3K9ac (GSM772811) and H3K9me3 (GSM621651). The MACS version 1.4.1 (model-based analysis of ChIP-Seq) peak-finding algorithm was used to identify regions of ChIP-Seq enrichment47. A FDR threshold of enrichment of 0.01 was used for all data sets. The resulting genomic coordinates in bed format were further used in ChromHMM v1.06 for chromatin annotation45. For comparison, the same genomic coordinates were investigated in chromatin annotation maps of renal tissue, as well as across nine different cell lines from the ENCODE Project: umbilical vein endothelial cells (HUVEC), mammary epithelial cells (HMEC), normal epidermal keratinocytes (NHEK), B-lymphoblastoid cells (GM12878), erythrocytic leukemia cells (K562), normal lung fibroblasts (NHLF), skeletal muscle myoblasts (HSMM), embryonic stem cells (H1 ES) and hepatocellular carcinoma cells (HepG2). We tested whether the proportion of SNPs pointing to either strong or weak enhancers in the human kidney tissue cells was different from that of the other nine tissues by means of a Fishers' exact test for 2 × 2 tables, contrasting each of the nine cell lines listed above against the reference kidney cell line, at a Bonferroni-corrected significance level of 0.05/9=5.6 × 10−3.

Replicated gene regions were prioritized for functional studies using the following criteria: (1) GRAIL identification of a gene in each region of P value<0.05 or DEPICT, FDR <0.05); (2) an eGFRcrea to eGFRcys ratio between 0.2 and 5 with direction consistency between the beta coefficients; (3) nearest gene if the signal was located in a gene-poor region. The list of genes selected for functional work can be found in Supplementary Table 12. This same prioritization scheme was also used to assign locus names. Morpholino knockdowns were performed in zebrafish.

Zebrafish (strain Tübingen, TU) were maintained according to established Harvard Medical School Institutional Animal Care and Use Committee protocols (protocol # 04626). Male and female fish were mated (age 6–12 months) for embryo production. Embryos were injected at the one-cell stage with MOs (GeneTools) designed to block either the ATG start site or an exon–intron splice site of the target gene (Supplementary Table 21). In cases where human loci are duplicated in zebrafish, both orthologues were knocked down simultaneously by combination MO injection. MOs were injected in escalating doses at concentrations up to 250 μM. Embryos were fixed in 4% paraformaldehyde at 48 h post fertilization for in situ hybridization using published methods (http://zfin.org/ZFIN/Methods/ThisseProtocol.html). Gene expression was visualized using established renal markers pax2a (global kidney), nephrin (podocytes) and slc20a1a (proximal tubule). The number of morphant embryos displaying abnormal gene expression was compared with control embryos by means of a Fisher's exact test.

Additional information

How to cite this article: Pattaro, C. et al. Genetic associations at 53 loci highlight cell types and biological pathways relevant for kidney function. Nat. Commun. 7:10023 doi: 10.1038/ncomms10023 (2016).