Visual AbstractAbstract

Background and objectives Chronic kidney disease (CKD) is a global public health issue associated with large economic burdens. CKD contributes to higher risks of cardiovascular complications, kidney failure, and mortality. The incidence and prevalence rates of kidney failure in Taiwan have remained the highest in the world.

Design, setting, participants, & measurements Assessing genetic factors that influence kidney function in specific populations has substantial clinical relevance. We investigated associations of genetic variants with eGFR. The quality control filtering and genotype imputation resulted in 10,008 Taiwan Biobank participants and 6,553,511 variants for final analyses. We examined these loci with in silico replication in individuals of European and African ancestry.

Results Our results revealed one significant locus (4q21.1) and three suggestive significant loci (17q23.2, 22q13.2, and 3q29) for eGFR in the Taiwanese population. In total, four conditional-independent single nucleotide polymorphisms were identified as the most important variants within these regions, including rs55948430 (Coiled-Coil Domain Containing 158), rs1010269 (BCAS3), rs56108505 (MKL1), and rs34796810 (upstream of DLG1). By performing a meta-analysis, we found that the 4q21.1 and 17q23.2 loci were successfully replicated in the European population, whereas only the 17q23.2 locus was replicated in African ancestry. Therefore, these two loci are suggested to be transethnic loci, and the other two eGFR-associated loci (22q13.2 and 3q29) are likely population specific.

Conclusions We identified four susceptibility loci on 4q21.1, 17q23.2, 22q13.2, and 3q29 that associated with kidney-related traits in a Taiwanese population. The 22q13.2 (MKL1) and 3q29 (DLG1) were prioritized as critical candidates. Functional analyses delineated novel pathways related to kidney physiology in Taiwanese and East Asian ancestries.

Introduction

Taiwan has the highest incidence and prevalence of chronic kidney disease (CKD), which is associated with health and economic costs (1,2). The Global Burden of Disease study estimates that 1.2 million people died from kidney failure, an increase of 32% from 2005 to 2015 (3). In 2010, an estimated 2.3–7.1 million people with kidney failure died without access to long-term dialysis (4). CKD is known as a disease multiplier for cardiovascular diseases (2). Hence, targeting CKD progression can considerably reduce cardiovascular mortality and attenuate the progression to kidney failure (5).

The aggregation of CKD in families suggests a genetic component to this complication (6,7). Population-based genome-wide association studies (GWAS) support a role of genetic susceptibility for kidney function decline (8). Genetic determinants of eGFR have been explored in the European general populations (9,10). A meta-GWAS identified 24 novel loci associated with eGFR in European ancestry individuals in addition to 29 known loci (9). Another meta-GWAS revealed the susceptibility loci, such as the MHC region, UNCX, and MPPED2-DCDC5, for eGFR in East Asian populations (11). In addition, two genetic loci on chromosome 17 (Chr. 17; BCAS3) and Chr. 2 (LRP2) were replicated in both Japanese (12) and Korean populations (13). More than 350 novel loci were further identified using a meta-GWAS of over 1 million individuals across ancestry groups from both population-based and hospital-based studies (14,15). These results provided prioritized gene lists for target selection and experimental follow-up.

Evidence shows that the effects of genetic variants vary across study populations and distinct phenotypically defined kidney diseases (8). These findings suggest that the interethnic genetic variability should be considered in the genetic study for kidney-related traits. Previous GWASs were mostly conducted in European ancestries. In this study, we implemented a genome-wide approach to identify susceptibility loci for kidney function in a Taiwanese population.

Materials and MethodsStudy Participants

Taiwan Biobank is a national project that has collected specimens and lifestyle and genetic information in Taiwan. Taiwan Biobank has recruited a community-based cohort since 2005. The participants are volunteers of ages 30–70 years. Written informed consent from participants were provided. This study included a total of 10,686 participants from Taiwan Biobank in 2018. This study was approved by the Institutional Review Board of Taipei Medical University (TMU-JIRB N201802059) and the Taiwan Biobank in Academia Sinica (TWBR10505–05, TWBR10602–02, and TWBR10703–02).

The medical information from Taiwan Biobank mainly consisted of three categories: questionnaires, physical examinations, and blood/urine tests. The presence of hypertension and/or diabetes mellitus was defined as having hypertension and/or diabetes mellitus as a current illness or fulfilling the diagnostic criteria of hypertension (16) and/or diabetes mellitus (17).

Genotyping, Quality Control Filtering, and Genotype Imputation

Genotyping experiments were conducted by Academia Sinica (Taiwan) on an Affymetrix Taiwan Biobank chip that contained nearly 653,000 single nucleotide polymorphism (SNP) sites. A flow diagram of the study is shown in Figure 1. SNPs with a genotype call rate of <0.95 or a minor allele frequency of <0.05 were excluded. In addition, SNPs with a Hardy–Weinberg equilibrium exact test P<10−5 in the control subset were also removed. SNPs that passed the quality control filtering were further subjected to haplotype phasing using SHAPEIT2 (18) and genotype imputation using IMPUTE2 version 2.3.1 (19) on the basis of a reference panel from the Taiwan Biobank and 1000 Genome Project phase 3 East Asian (1000GP3-EAS) population. After genotype imputation, we filtered out variants with a minor allele frequency of <0.01, and the subsequent removal of indels and triallelic variants left 6,553,511 SNPs for the final analyses.

Figure 1.

Flow diagram of the study design, quality control, and imputation steps for selection and analysis of Taiwan Biobank participants. EAS, East Asian; eQTL, expression quantitative trait locus; DM, diabetes mellitus; GCTA-COJO, genome-wide complex trait analysis-conditional and joint association analysis; GWA, genome-wide association; HTN, hypertension; MAF, minor allele frequency; METAL, a tool for meta-analysis genomewide association scans; PC, principal component; QC, quality control; SNP, single nucleotide polymorphism; TWB, Taiwan Biobank.

For sample quality control, participants with related individuals (pihat >0.1875), a high missing genotype rate of >0.01, an elevated heterozygosity rate that deviated from the mean ± three SDs, or inconsistent data on sex were all excluded, as were 15 participants for whom data for serum creatinine were unavailable, resulting in 10,008 participants in the final analyses.

Definition of Phenotypes

The eGFR was estimated using the four-variable Modification of Diet in Renal Disease study equation (20). Values of eGFR of <15 ml/min per 1.73 m2 were set to 15 ml/min per 1.73 m2 and those of >200 ml/min per 1.73 m2 were set to 200 ml/min per 1.73 m2 to avoid undue influence from outliers, and further, they were natural log transformed.

Replication Study

An in silico replication study was conducted using two independent European (10) and African (9) ancestry datasets from the CKDGen database (http://ckdgen.imbi.uni-freiburg.de) for loci that satisfied P<10−6 in the discovery phase. To combine the findings in the discovery and replication phases, we performed a GWAS meta-analysis with public summary statistics.

Statistical Analyses

A linear regression model was fitted to the natural logarithm of eGFR assuming an additive genetic model using SNPTEST (21). Covariates of age, sex, hypertension, diabetes mellitus, and the first five principal components were included for adjustment. Genome-wide significance levels were set to P<5×10−8, and genome-wide suggestive levels were set to P<10−6 in the analyses. Lead variants were defined as variants that reached the smallest P values in each genetic locus. The independent association signal in the locus of interest was identified with a conditional analysis using GCTA (22). A potential secondary association signal within a given genomic region within ±1 Mb was reported if the variant with the smallest conditional P value was genome-wide significant (P<5×10−8) after conditioning on the lead variant from the discovery phase.

Genomic heritability was defined as the proportion of phenotypic variance attributable to genetic factors. In order to estimate the phenotypic variance explained by autosomal variations, we used a restricted maximum likelihood analysis implemented in GCTA (23). A GWAS meta-analysis was performed on the basis of inverse variance weighting using METAL (24). Successful replication was designated as P<5×10−8, with a consistent directional effect among the discovery and replication datasets. Manhattan and quantile-quantile plots were generated using the qqman package implemented in R. Regional association plots were constructed using LocusZoom (25).

Functional Annotation

Genetic variants were annotated using ANNOVAR (26). LDlink was used to investigate allele frequencies and linkage disequilibrium (LD) patterns for a genomic region of interest across population groups (27). The genomic features and potential function of the variant were predicted using HaploReg version 4.1 (28) and RegulomeDB (29). With those tools, we inputted a query SNP and explored chromatin states, conservation, and regulatory motif alterations within sets of variants with r2≥0.6 (LD was calculated from a 1000 Genomes phase 1 population of Asian ancestry).

Expression Quantitative Trait Locus

Publicly available expression quantitative trait locus association datasets from the Genotype Tissue Expression Project version 7 (30), NephQTL (31), Ko et al. (32), and the human kidney expression quantitative trait locus atlas (33) were explored to identify overlap between gene expression and identified eGFR association results.

Replicated Loci in the National Human Genome Research Institute Genome-Wide Association Study Catalog

Associations of all identified SNPs and SNPs in LD (r2>0.2) within ±1 Mb of other traits were confirmed according to the National Human Genome Research Institute (NHGRI) GWAS Catalog (34).

ResultsBaseline Characteristics of Participants

In total, data from 10,008 postfiltering Taiwan Biobank participants were enrolled in this study (Figure 1). Baseline characteristics of the participants are shown in Table 1. The prevalence of CKD was 2%, which was determined for CKD stages 3–5 (eGFR <60 ml/min per 1.73 m2). The principal component analysis indicated no substructure in the population (Supplemental Figure 1).

Table 1.

Baseline characteristics of the study participants in the Taiwan Biobank

Discovery Phase of the Genome-Wide Association Study–Identified Loci Associated with Kidney Function

We conducted a quantitative trait locus analysis with eGFR as the continuous dependent variable and adjusted for seven covariates. A quantile-quantile plot showed that there was no severe genomic inflation of the association results (λGC=1.03) (Supplemental Figure 2). The Manhattan plot revealed one locus reaching genome-wide significance and three loci reaching suggestive significance (Figure 2). All SNPs with a minor allele frequency of >0.01 explained 12% of the heritability for eGFR on the basis of the GCTA approach. Regional association plots demonstrated zoomed-in proportions of Chr. 4q21.1, Chr. 17q23.2, Chr. 22q13.2, and Chr. 3q29 (Figure 3). Each lead signal was supported by additional SNPs in high LD with marginally significant P values associated with eGFR. In total, 54 SNPs attained suggestive significance. After filtering out SNPs located within the HLA region, 49 variants achieved suggestive significance in the discovery phase for replication.

Figure 2.

Manhattan plots of eGFR in a Taiwanese population revealed one genome-wide and three suggestive significant loci. The upper red line corresponds to a genome-wide significance threshold of 5×10−8, while the lower blue line corresponds to a suggestive association threshold of 10−6.

Figure 3.

One genome-wide and three suggestive significant loci for eGFR in a Taiwanese population. Regional association plots of the loci of interest before (left panel) and after (right panel) conditioning on the lead variant (the purple diamond). Colors represent the degree of linkage disequilibrium (LD; r2) between each variant and the lead variant. LD (r2) was calculated on the basis of the Genome build hg19 version of the 1000 Genomes Project (Asian ancestry, November 2014). Chr, chromosome.

We further identified one locus on 4q21.1 that was significantly associated with eGFR P=9.98×10−11. The locus was led by the intronic variant rs55948430 in the Coiled-Coil Domain Containing 158 (CCDC158) gene. Additionally, 19 and 20 variants within the same locus were located across SHROOM3, FAM47E-STBD1, and STBD1 that demonstrated genome-wide associations and suggestive associations, respectively. We applied a conditional analysis on the lead SNP (rs55948430) in the 4q21.1 locus, however, no additional signal was detected in the region.

Additionally, we identified three loci (17q23.2, 22q13.2, and 3q29) that were suggestively associated with eGFR, where each of the lead SNPs was conditioned. We affirmed each of the lead SNPs as a single independent signal within each region. Lead SNPs on Chr. 17 (rs1010269) and on Chr. 22 (rs56108505) were found to be located on the introns of BCAS3 and MKL1, respectively. The conditionally independent lead SNP on Chr. 3 (rs34796810) resided in the intergenic region near DLG1.

Our results revealed four conditionally independent SNPs that associated with eGFR in the Taiwanese population (Table 2). Three of them were within introns of genes, whereas the other one resided in an intergenic region. The effect alleles of rs55948430[T] and rs1010269[A] were associated with lower eGFR, with effect sizes of −0.032 and −0.018, respectively. The effect alleles of rs56108505[C] and rs34796810[G], on the other hand, were associated with the higher levels of eGFR (β=0.019 and β=0.020, respectively).

Table 2.

Four conditionally independent single nucleotide polymorphisms associated with eGFR were identified in the discovery cohort and the replication (European)/combined stage

Replication of the Detected Variants in Two Independent Datasets

In the European population study, 44 SNPs were available for in silico replication. The replication phase among Europeans resulted in 4q21.1 with 35 variants including rs55948430 in the CCDC158 gene, and 17q23.2 with five variants including rs1010269 in the BCAS3 gene. Both variants reached the genome-wide significant threshold after the meta-analysis (Table 2). We successfully replicated the variants on two loci for kidney function. After the combined analysis, SHROOM3 lead variant rs17319721 had a lower meta-analysis P value (βmeta=−0.012, Pmeta=1.26×10−36) over the discovery P of 5.07×10−8 for the 4q21.1 locus. On the 17q23.2 locus, the previously reported lead SNP rs9895661 had a lower meta-analysis P value (βmeta=−0.013, Pmeta=1.47×10−27) over discovery P of 1.87×10−7. In contrast, the other two suggestive significant loci (22q13.2 and 3q29) identified in the discovery phase were seen only in the Taiwan Biobank cohort.

In the replication analysis among those of African ancestry, 15 SNPs in two loci with genotype data were available for a meta-analysis. Because the lead SNP rs1010269 was not in the African ancestry dataset, we therefore used rs11657044 (r2=0.94, estimated by our GWAS data) as a substitute SNP. After combined analysis, three variants in LD (r2>0.6, estimated by our GWAS data) with lead variant rs1010269 on the 17q23.2 locus reached the genome-wide significance threshold (Table 3). On the other hand, 4q21.1 for eGFR was replicated successfully in both Taiwanese and European populations but not in the African population.

Table 3.

Results of association analyses in the replication (African ancestry) and combined studies

Genome-Wide Association Study Catalog Query

According to the GWAS Catalog, we found multiple traits in genome-wide significant associations with our identified SNPs and variants (Supplemental Table 1). In support of our findings, those SNPs were mostly reported to associate with kidney-related traits, such as serum creatinine, eGFR on the basis of creatinine, BUN levels, the incidence of CKD, and urinary albumin excretion. In addition, the 4q21.1 locus has been reported to associate with magnesium level and red blood cell–related phenotypes (hematocrit, red blood cell counts, and the hemoglobin concentration). The Chr. 17q23.2 locus associated with gout, urolithiasis, serum uric acid levels, and calcium levels; these genes, at least, potentially contribute to the regulation of kidney function. Importantly, consistent with the GWAS from BioBank Japan Project–based GWAS (35) and the recent transethnic studies (14), BCAS3 locus was identified for eGFR and creatinine quantitative trait susceptibility.

Kidney Expression Quantitative Trait Locus Analysis

According to the NephQTL dataset, the lead SNP (rs55948430) of the 4q21.1 locus was found to affect CCDC158 and STBD1 gene expressions in the kidney tubulointerstitium. The allele rs55948430[T] is correlated with lower CCDC158 expression levels in the kidney tubulointerstitium. On the same locus, the European ancestry lead variant (rs17319721) was likely to affect CCDC158, STBD1, and SHROOM3 gene expressions. On the 22q13.2 locus, patients carrying the C allele of the lead SNP (rs56108505) tended to have higher E1A Binding Protein P300 expression levels in kidney tubulointerstitial samples. Finally, the top signal of rs34796810 on 3q29 also influenced Rubicon Autophagy Regulator and Forty-Two-Three Domain Containing 1 gene expressions in glomeruli and the tubulointerstitium, respectively (Supplemental Table 2).

Functional Annotation

On 4q21.1, variants in high LD with the lead variant (rs55948430), including rs1828538 (r2=0.96 as estimated by 1000GP3-EAS), rs56049812 (r2=0.92 as estimated by 1000GP3-EAS), and rs17319721 (r2=0.62 as estimated by 1000GP3-EAS), likely affect protein binding. Among potential functional variants, the 5′-untranslated region variant in FAM47E-STBD1 (rs56049812) was found to be enriched in active transcription start sites, enhancers (H3K4me1 and H3K27ac), and promoters (H3K4me3 and H3K9ac) of histone marks in multiple tissues (Roadmap Epigenomics Consortium 2015).

The BCAS3 intronic variants on 17q23.2, such as rs11657044 (r2=0.89 as estimated by 1000GP3-EAS) and rs9895661 (r2=0.61 as estimated by 1000GP3-EAS), were both enriched in enhancers of histone marks. Finally, on 22q13.2, variants in complete LD with the lead variant (rs56108505), including rs73169089 and rs73167096 (r2=1.0 as estimated by 1000GP3-EAS), were predicted to affect protein binding as well. Moreover, enrichment of enhancers of chromatin marks was also found among these variants.

Discussion

This study revealed four genetic loci that associated with kidney functions in a Taiwanese population, including 4q21.1 (CDCC158, SHROOM3, FAM47E, and FAM47E-STBD1), 17q23.2 (BCAS3), 22q13.2 (MKL1), and 3q29 (DLG1). A novel lead SNP on the 4q21.1 locus was identified in this study. The conditional lead signal of rs55948430 is located in the intronic region of the CCDC158 gene, which encodes CCDC158. CCDC158 variants were reported in a recent transethnic GWAS to associate with kidney function–related traits, such as cystatin C levels, creatinine levels, and eGFR (36). Indeed, CCDC158 expresses in the kidney; however, the biologic function of CCDC158 is unclear. The present lead signal of rs55948430 is in LD (r2=0.73 as estimated by our GWAS data) with the previously reported eGFR variant (rs17319721) of SHROOM3 (Supplemental Figure 3) (9). SHROOM3 encodes a PDZ domain–containing protein that binds F actin and regulates its subcellular distribution. The risk allele of rs17319721 was found to correlate with higher SHROOM3 expression levels in patients undergoing transplantation. Increased SHROOM3 expression facilitated TGF-β1 signaling and contributed to interstitial fibrosis of the transplant and dysfunction of the allograft (37). Introgression of the wild-type Shroom3 gene into a rat with a fawn-hooded hypertensive genetic background restored glomerular function. The wild-type Shroom3 allele rescued glomerular defects induced by knockdown of endogenous Shroom3 in zebrafish (38). Furthermore, SHROOM3 protein was reported to maintain the podocyte architecture through modulating the actin cytoskeleton (39). Consistent with previous studies, the SNP in SHROOM3 was found to affect CKD susceptibility, urine albumin-creatinine ratio, and albuminuria (40,41).

The second genetic locus was found on 17q23.2 (BCAS3). 17q23.2 locus was reported to associate with eGFR in multiple GWASs, including of European and non-European populations (9⇓⇓⇓⇓⇓–15). rs1010269 was in high LD (r2=0.94 as estimated by our GWAS data) with the most significant independent SNP (rs11657044) that associated with eGFR in a European meta-GWAS (9). Stanzick et al. (15) prioritized the BCAS3 variant to the high relevance for the protein (CADD score ≥15). The findings indicated strong susceptibility of rs1010269 and rs11657044 variants. The lead variant was also found to correlate with hemoglobin concentrations (42), kidney stones (43), and gout (44). In addition, BCAS3 gene encodes breast carcinoma–amplified sequence 3 or the microtubule-associated cell migration factor, which is associated with breast cancer. Other common variants of BCAS3 were also reported to associate with gout formation, serum urate levels, and cardiovascular diseases (45). The phenotype-variant correlations highlight the pleiotropic effect of the BCAS3 gene.

The third eGFR-susceptible locus was found on 22q13.2 (MKL1). MKL1, known as myocardin-related transcription factor A (MRTFA), encodes MRTFA. This protein interacts with the transcription factor myocardin, a key regulator of smooth muscle cell differentiation. Intronic variants of MKL1 were reported to associate with serum creatinine and GFR in a Japanese population (35) and transethnic meta-GWAS (14). Moreover, MRTFA was reported to contribute to fibroblast activation and kidney fibrosis (46). Indeed, MRTFA deficiency attenuates AKI in mice with diabetic nephropathy (47).

The final locus was found on 3q29 (upstream of DLG1). The closest gene, DLG1, encodes discs-large MAGUK scaffold protein 1, which plays an important role in epithelial cell polarity and neuronal synaptic function. Functional studies showed that a DLG1 deletion causes severe misalignment of the mouse ureteric smooth muscle cells, resulting in impaired urinary transport and hydronephrosis (48). Okada et al. (11) identified DLG1 variants that affect kidney function in the discovery phase.

By assessing the allelic frequency of the identified lead SNP, we found the diversity among different ethnic groups (Supplemental Table 3). Importantly, the Taiwanese population-specific variants on 22q13.2 exhibited higher-effect allelic frequencies in the eastern population than the western population (0.22–0.25 versus 0.01–0.09). This allelic distribution might explain why those two loci were only identified in the Taiwanese population. On the other hand, the 4q21.1 locus showed consistent associations across ancestry groups, despite variability in its allelic distribution. Interestingly, a previous epidemiologic study indicated that the prevalence of CKD (stages 1–5) (49) was higher with higher eGFR risk allele frequency. For example, Europeans with the highest prevalence of CKD (18%) have the highest risk allele frequency of 0.43, whereas those with African ancestry with the lowest prevalence of CKD (9%) have the lowest risk allele frequency of 0.03.

The strengths of this study are the densely imputed genetic data and the stringent quality control process, leading to the identification for Taiwanese-specific kidney-related loci. Despite these strengths, there are some limitations. First, our study indicated an underestimation of heritability as compared with family study (44%) (7) that might result from the complexity of GFR phenotypes. For multifactorial traits, such as the GFR, most SNPs only contribute to small effect sizes; thus, a larger sample size is required. Second, rare variants and structural variations were not captured by these analyses. Third, adjusting for hypertension and diabetes mellitus may introduce potential bias due to the opening of a backdoor path through unmeasured confounders of the relation between hypertension/diabetes mellitus and eGFR (50). Finally, functional studies for verifying the pathogenicity of variants are still needed. In short, we identified four susceptibility loci associated with kidney-related traits in a Taiwanese population. Our findings provide new insights into the molecular mechanisms of kidney functions.

Disclosures

W.-C. Chang, Y.-C. Chen, W.-H. Chou, C.-C. Kao, H.S.-C. Wong, and M.-Y. Wu report employment with Taipei Medical University. M.-S. Wu reports employment with Taipei Medical University and honoraria from Astellas, AZ, Baxter, and Roche. The remaining author has nothing to disclose.

Funding

This work was supported by Ministry of Science and Technology, Taiwan: MOST109-2314-B-038-131; MOST109-2628-B-038-012; National Health Research Institutes, Taiwan: NHRI-EX111-10926HT; Taipei Medical University, Taiwan (12310-106079;12310-10739) for Yusuke Nakamura Chair Professorship.

Acknowledgments

We thank Dr. Che-Mai Chang (Taipei Medical University, Taiwan) and Dr. Chang, Jan-Gowth (China Medical University, Taiwan) for their helpful comments during the revision.

Author Contributions

W.-C. Chang, M.-S. Wu, and M.-Y. Wu conceptualized the study; C.-H. Chao, Y.-C. Chen, C.-C. Kao, W.-H. Chou, and H.S.-C. Wong were responsible for data curation; W.-C. Chang, C.-C. Kao, and M.-S. Wu were responsible for investigation; Y.-C. Chen was responsible for formal analysis; C.-H. Chao and Y.-C. Chen were responsible for checking the methodology; W.-C. Chang and M.-Y. Wu were responsible for project administration; M.-S. Wu, W.-C. Chang, and M.-Y. Wu were responsible for resources; Y.-C. Chen and W.-H. Chou were responsible for software and IRB; C.-H. Chao and Y.-C. Chen were responsible for validation; Y.-C. Chen and H.S.-C. Wong were responsible for visualization; W.-C. Chang, C.-C. Kao, M.-S. Wu, and M.-Y. Wu were responsible for funding acquisition; M.-S. Wu and W.-C. Chang provided supervision; W.-C. Chang, Y.-C. Chen, C.-C. Kao, H.S.-C. Wong, and M.-S. Wu wrote the original draft; W.-C. Chang, M.-Y.Wu, and W.-H. Chou reviewed and edited the manuscript; and Y.-C. Chen, W.-C. Chang, M.-Y. Wu, C.-C. Kao, and W.-H. Chou were involved in revision.

Received February 18, 2022.Accepted September 16, 2022.Copyright © 2022 by the American Society of Nephrology