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.

C hronic kidney disease (CKD) is a global public health problem [1][2][3] , and is associated with an increased risk for cardiovascular disease, all-cause mortality and end-stage renal disease 4,5 . Few new therapies have been developed to prevent or treat CKD over the past two decades 1,6 , 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 CKD [7][8][9][10][11][12][13][14][15] . Subsequent functional investigations point towards clinically relevant novel mechanisms in CKD that were derived from initial GWAS findings 16 , 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 o60 ml min À 1 per 1.73 m 2 ) with B2.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 loci 7-9 ( 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 o1.0 Â 10 À 6 . Of these 48 novel SNPs, 21 were genome-wide significant with P values o5.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 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 o5.0 Â 10 À 8 (Table 1). Of these, 23 fulfilled additional replication criteria (q-value o0.05 in stage 2). Only rs6795744 at the WNT7A locus demonstrated suggestive replication (P value o5.0 Â 10 À 8 , q-value 40.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.
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  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.
Exploratory comparison of the association effect sizes in subjects with and without hypertension based on our previous work 7 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 pressure 17 , myocardial infarction 18 , left ventricular mass 19 , heart failure 20 , fasting glucose 21 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).
( 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-analysis 22 , 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 log 10 Bayes Factor 46, 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 log 10 Bayes Factor 46.
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 o10 À 5 in the discovery metaanalysis, 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 o10 À 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.  Table 12) or closest gene. A complete overview of the genes in each locus is given in the regional association plots ( Supplementary Fig. 4). ySNP 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. zFor 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 I 2 . Q statistic P value 40.05 for all SNPs, except rs11666497 (SIPA1L3, P value ¼ 0.04). 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 eGFRcreaassociated SNPs with P values o10 À 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 transcription 24 , 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 410 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 a ¼ 5.6 Â 10 À 3 ).
Functional characterization of new loci. To prioritize genes for functional studies, we applied gene prioritization algorithms including GRAIL 25 , 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 o0.05 or DEPICT false discovery rate (FDR) o0.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 Effect on log(eGFRcrea) a b (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 hybridization 12 . 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 nonrenal 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 groups 7-13 by identifying 450 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 effort 7 , 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 phenotypes 26 .
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 Human renal glomerular endothelial cells 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 o0.05) at the P value o10 À 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.
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 nutrition 28 . 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 variation 9 ARTICLE o15 ml min per 1.73 m 2 were set to 15, and those 4200 were set to 200 ml min À 1 per 1.73 m 2 . CKD was defined as eGFRcrea o60 ml min À 1 per 1.73 m 2 . Diabetes was defined as fasting glucose Z126 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 B2.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 package 32 in R, before including them into the meta-analysis. Genome-wide meta-analysis was performed with the software METAL 33 , assuming fixed effects and using inverse-variance weighting. The genomic inflation factor l 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 w 2 with 1 degree of freedom, with b and s.e. representing the effect of each SNP on the phenotype and its standard error, respectively 34 . Genomic-control (GC) correction was applied to P values and s.e.'s in case of l41 (first GC correction). SNPs with an average minor allele frequency (MAF) of Z0.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 I 2 statistic.
After removing SNPs with MAF of o0.05 and which were available in o50% of the studies, SNPs with a P value of r10 À 6 were selected and clustered into independent loci through LD pruning based on an r 2 of r0.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 r10 À 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  N ¼ 4,955). Thirteen studies contributed replication data on CKD (N ¼ 33,972; 4,245 CKD cases; studies with o50 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 r5.0 Â 10 À 8 in accordance with previously published guidelines 35 ; (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 o0.05 in the replication stage. Q-values were estimated using the package QVALUE 36 in R. The tuning parameter lambda for the estimation of the overall proportion of true null hypotheses, p0, was estimated using the bootstrap method 37 . 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 QUANTO 38 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).
Associations stratified by diabetes and hypertension status. 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 twosample t-test for correlated data: where b DM and b nonDM represent the SNP effects on log(eGFRcrea) in the two groups, s.e. is the standard error of the estimate and r(.) 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 work 7 . 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 P 53 Þ=varðyÞ 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 i is the estimated effect of the i th SNP on y, y corresponds to the sexand age-adjusted residuals of the logarithm of eGFRcrea and var(SNP i ) ¼ 2 Â MAF SNPi Â (1 À MAF SNPi ) 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 ICBP 17 , myocardial infarction in CARDiOGRAM 18 , left ventricular mass 19 , heart failure 20 , the urinary albumin to creatinine ratio (CKDGen consortium, personal communication) and fasting plasma glucose in MAGIC 21 . 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 (r 2 40.2) within ±1 MB distance were checked for their association with other traits according to the NHGRI GWAS catalog 40 (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 consortium 11 .
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 Z0.01 and imputation quality r 2 Z0.3 were considered for the meta-analysis. After meta-analysis, we removed SNPs with MAF o0.05 and which were available in o50% 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 metaanalysis 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 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms10023 of Trans-ethnic Association studies) software 22 . 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 Z0.01 and imputation quality r 2 Z0.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 Z120,000) were considered after meta-analysis. Genomewide significance was defined as a log 10 Bayes' Factor (log 10 BF) Z6 (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 used 25 . 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 (r 2 40.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 PLINK 42 to identify independently associated SNPs using all SNPs with eGFRcrea association P values o10 À 5 (HapMap release 27 CEU data 43 ; LD r 2 threshold ¼ 0.01; physical kb threshold ¼ 1,000). We then used the DEPICT method 23 to construct associated regions by mapping genes to independently associated SNPs if they overlapped or resided within LD (r 2 40.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 nonoverlapping 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 o1e À 5 (Fig. 4), we estimated their overlap by computing the pairwise Pearson correlation coefficient r between each pair of gene sets followed by discretization into one of three bins; 0.3rro0.5, low overlap; 0.5rro0.7, medium overlap; rZ0.7, high overlap.
DNase I hypersensitivity analysis. The overlap of SNPs associated with eGFRcrea at Po10 À 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 tissues 44 (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 PLINK 42 in windows of 100 kb and maximum r 2 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 r 2 ¼ 0.8 in windows of 100 kb, again using LD relationships in the 1000G EUR panel, restricted to SNPs with MAF 41% 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 o5.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 studies 45 . We therefore investigated the genomic mapping of the known and replicated novel index SNPs, as well as their perfect LD proxies (n ¼ 173, r 2 ¼ 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 Consortium 46 . 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 enrichment 47 . 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 annotation 45 . 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 .
Functional characterization of new loci. Replicated gene regions were prioritized for functional studies using the following criteria: (1) GRAIL identification of a gene in each region of P valueo0.05 or DEPICT, FDR o0.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 mM. 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.