Family-Based Genome-Wide Association Study of South Indian Pedigrees Supports WNT7B as a Central Corneal Thickness Locus.

Purpose: To identify genetic risk factors contributing to central corneal thickness (CCT) in individuals from South India, a population with a high prevalence of ocular disorders.
Methods: One hundred ninety-five individuals from 15 large South Indian pedigrees were genotyped using the Omni2.5 bead array. Family-based association for CCT was conducted using the score test in MERLIN.
Results: Genome-wide association study (GWAS) identified strongest association for single nucleotide polymorphisms (SNPs) in the first intron of WNT7B and CCT (top SNP rs9330813; β = -0.57, 95% confidence interval [CI]: -0.78 to -0.36; P = 1.7 × 10-7). We further investigated rs9330813 in a Latino cohort and four independent European cohorts. A meta-analysis of these data sets demonstrated statistically significant association between rs9330813 and CCT (β = -3.94, 95% CI: -5.23 to -2.66; P = 1.7 × 10-9). WNT7B SNPs located in the same genomic region that includes rs9330813 have previously been associated with CCT in Latinos but with other ocular quantitative traits related to myopia (corneal curvature and axial length) in a Japanese population (rs10453441 and rs200329677). To evaluate the specificity of the observed WNT7B association with CCT in the South Indian families, we completed an ocular phenome-wide association study (PheWAS) for the top WNT7B SNPs using 45 ocular traits measured in these same families including corneal curvature and axial length. The ocular PheWAS results indicate that in the South Indian families WNT7B SNPs are primarily associated with CCT. Conclusions: The results indicate robust evidence for association between WNT7B SNPs and CCT in South Indian pedigrees, and suggest that WNT7B SNPs can have population-specific effects on ocular quantitative traits.

Ocular quantitative traits such as central corneal thickness (CCT), axial length (AXL), and intraocular pressure are heritable intermediate phenotypes (endophenotypes) for common complex eye disorders such as keratoconus, myopia, and glaucoma.[1] CCT is a highly heritable ocular quantitative trait with up to 95% of its phenotypic variance due to genetics.[2] Thin CCT is related to several diseases of the cornea, especially keratoconus[3] and brittle corneal syndrome.[4] Very thin corneas are a hallmark of Ehlers Danlos,[5] and thicker than normal corneas are found in patients with aniridia.[6] Thinner than average CCT can influence development of primary open angle glaucoma[7,8] with more severe disease evident in people with thinner corneas.[9-11]

CCT varies among ethnic populations with individuals of African descent having lower values than European Caucasians and East Asians.[2,12,13] Genome-wide association studies (GWAS) in European Caucasians,[14-16] Asians,[14,17] and Hispanics[18] have identified ZNF469, RXRA-COL5A1, COL8A2, and FOXO1 among others as important loci contributing to CCT. RXRA-COL5A1 and ZNF469 have been associated with CCT in most populations studied while the associations of other loci (COL8A2, FOXO1) may be restricted to specific populations.[19] Recently, WNT7B single nucleotide polymorphisms (SNPs) have been associated with CCT in Latinos,[20] and interestingly some of these same SNPs were associated with AXL and corneal curvature, traits influencing myopic refractive error, in a Japanese population.[21]

Few genetic studies of ocular quantitative traits have been completed in individuals from South India, a population with high prevalence of common ocular conditions, especially cataract and glaucoma.[22-27] In Indian populations, CCT is thinner than the average values for Caucasians[22] suggesting that CCT could be an important factor in the development of CCT-related common ocular disorders in this population. To identify genetic risk loci for CCT in South Indians, we completed a family-based association study using large pedigrees, many with consanguineous matings that are typical for this geographic region. For the top SNPs located in the WNT7B region, we also completed a phenome-wide association study (PheWAS) to examine the range of phenotypes associated with WNT7B SNPs in this South Indian population.

Materials and Methods

Pedigrees and Quantitative Traits

This study adhered to the tenets of the Declaration of Helsinki and has been reviewed and approved by the Institutional Review Boards of Massachusetts Eye and Ear Infirmary and Medical Research Foundation, Sankara Nethralaya, Chennai, India. After obtaining written informed consent, 197 individuals from 15 Indian pedigrees were recruited at Sankara Nethralaya, Chennai, India. CCT was measured by an ultrasonic pachymeter in triplicate and the average value was used. Methods to measure the other traits used in the PheWAS are described in the Supplementary Methods. Collections of samples for replication cohorts are described in the Supplementary Methods.

Genotyping and Quality Control (QC)

Genotyping for the South Indian families was performed at the Ocular Genomics Institute at the Massachusetts Eye and Ear Infirmary using the Illumina HumanOmni2.5-8 Beadchip kit (2,379,855 markers; Illumina, Inc., San Diego, CA, USA). Genotypes were called using GenomeStudio (v2011.1, Illumina, Inc.). The genetic sex of all individuals was consistent with the reported sex. Two samples were removed because genotyping call rates were <99%. The average call rate per sample was >99.8%. QC for 2,352,697 (98.9%) well-clustered SNPs was performed with PLINK (v1.07, provided in the public domain, http://pngu.mgh.harvard.edu/∼purcell/plink/).[28] 25,088 (1.1%) SNPs with call frequency <90% and 881,678 (37.5%) SNPs with minor allele frequency (MAF) <0.01 were removed from the analysis. 164,174 (7.0%) SNPs with Mendelian errors and 58,443 (2.5%) SNPs on chromosome X or Y, or on the mitochondrial chromosome were also excluded. After QC, 1,223,314 SNPs were included in the final analysis. Genotyping for replication cohorts is described in the Supplementary Methods.

Statistical Analysis

The kinship coefficients for pairwise relationships across pedigrees were estimated from the SNP data using the KING software (provided in the public domain, http://people.virginia.edu/∼wc9c/KING/index.html).[29] The heritability for each trait was estimated with restricted maximum likelihood–based linear modeling in the GCTA software (provided in the public domain, http://gcta.freeforums.net/),[30] taking into account all pedigree relationships simultaneously. Inverse-normal transformation of ranks was applied to CCT measurements before analysis. Age and sex were included as covariates in the association tests. The genome-wide association test was performed using the score test in MERLIN (v1.1.2, provided in the public domain, https://csg.sph.umich.edu/abecasis/Merlin/),[31,32] which incorporated genetic relatedness based on the family structure. Because this program applies a restriction on pedigree size, 8 of the 15 pedigrees were split into nonoverlapping fragments of ≤18 bits using the PedSTR program (provided in the public domain, http://mga.bionet.nsc.ru/soft/PedStr/PedStr.tar.gz),[33] which breaks inbreeding loops and identifies subpedigrees having the maximal total relationship between individuals of interest, resulting in a total of 26 effective subpedigrees used in the final analysis. To avoid an excess of false-positive results in regions of strong linkage, the likelihood-ratio test was performed to accurately evaluate the SNPs with suggestive association. The regional SNP association plot was generated using SNAP (provided in the public domain, http://archive.broadinstitute.org/mpg/snap/).[34] The variance in CCT explained by all the SNPs in the Indian population was estimated using GCTA.[30]

Meta-analysis using the inverse-variance weighting method was done using both fixed-effects and random-effects models using Review Manager software (RevMan, version 5.3; Copenhagen: The Nordic Cochrane Centre, The Cochrane Collaboration, 2014). The heterogeneity between data sets was evaluated by heterogeneity index (I2) and Cochran's Q statistic.[35] Heterogeneity among data sets was further examined by evaluating differences in ethnicity (Indians, Latinos, or Europeans), study design (family-based design or population design), imputation quality score, age, and sex in meta-regression models using the R package “metafor.”[36] Forest plots were generated using the R package “metafor.”[36]

PheWAS

Forty-five quantitative traits (including CCT) (Supplementary Table S1) were analyzed for association as described above. Methods for measuring each trait are described in the Supplementary Methods. The average value for each trait for both eyes was used for analysis. Age and sex were included as covariates in the association tests. The association tests were performed using the likelihood-ratio test in MERLIN (v1.1.2).[31,32] The PheWAS plots were generated using the R package ggplot2 (provided in the public domain, https://www.r-project.org/).[37] Phenotypes were grouped along the x-axis by categorization of ocular measures (Biometric traits, Corneal traits, Optic nerve traits, Refractive error traits). Each point in the plot represents the −log10(P) value of a trait measure in association analysis. The lower gray dashed line indicates P = 0.05. The upper black dashed line indicates a single-SNP Bonferroni correction P = 0.001 (0.05/45).

Power Analysis

Power analysis was performed using the Genetic Power Calculator (provided in the public domain, http://pngu.mgh.harvard.edu/∼purcell/gpc/).[38] The total proportion of trait variance was derived from the estimated heritability of these ocular traits in the South Indian pedigrees. For CCT, AXL, and corneal curvature, heritability was 0.54, 0.84, and 0.82, respectively. The quantitative trait locus (QTL) increaser allele frequency was set to the same as the marker allele frequency. Linkage disequilibrium between the QTL and the marker was set at D' = 1.0. The sample size was set as 26 because a total of 26 effective subpedigrees were used in the final analysis. The sibling correlation was set as 0.5. The sibship size was set as 2. An additive effects only (1 df) test was used to calculate the power at the type I error rate of 5 × 10−8 for GWAS or 0.001 (0.05/45 traits) for PheWAS. Power results for all traits are listed in Supplementary Table S1.

Results

Study Sample

One hundred ninety-five individuals from 15 pedigrees (Supplementary Fig. S1) were recruited at Sankara Nethralaya Eye Hospital, Chennai, India for a family-based genetic association study. These pedigrees were unrelated to each other; the maximum kinship coefficient estimated from the SNP data across pedigrees was 0.0344. The pedigree size ranged from 2 to 26 members. Ten of the pedigrees included at least one consanguineous mating. Fifty-eight percent of the subjects were female and 42% male. The average age was 44.9 (±15.0) years and the age ranged from 16 to 85 years. These families were not ascertained on specific eye conditions. CCT was measured by an ultrasonic pachymeter in triplicate for each eye (Supplementary Methods) and the average value for both eyes was used (516.2 [±30.2] μm average; 433–608 μm range; Supplementary Table S1).

Genome-Wide Association Results for CCT

After QC, 1,223,314 SNPs were included in the genome-wide CCT analysis. The results for the family-based association test are shown in Supplementary Figure S2. The genomic inflation factor of 1.05 (QQ plot, Supplementary Fig. S3) suggested that population substructure or other confounding factors were not significant. Six SNPs located on chromosomes 6, 13, 18, and 22 showed suggestive evidence of association with CCT (P < 1.0 × 10−5; Table 1), with the top SNP (rs9330813, P = 1.7 × 10−7, β = −0.57, 95% confidence interval [CI]: −0.78 to −0.36 [A]) located in the first intron of WNT7B on chromosome 22 (Fig. 1). CCT association with rs9330813 was two orders of magnitude greater than any other SNP (Table 1) and accounted for 17% of the phenotypic variance in the South Indian families. WNT7B SNPs have previously only been associated with CCT in a Latino population (Mexican American Glaucoma Genetics Study [MAGGS]),[20] and the top SNP in the Latino study (rs10453441) is 422 bp from rs9330813. rs10453441 is in moderate linkage disequilibrium with rs9330813 in the South Indian dataset (r2 = 0.55) and was nominally associated with CCT in the South Indian pedigrees (P = 5.85 × 10−4, Supplementary Table S2).

Table 1

SNPs With P < 1.0 × 10−5 for Association With CCT in South Indian Pedigrees

Figure 1

Regional SNP association plot for the 22q13 region. A region of 408 kb around the top SNP (rs9330813) is displayed. The degree of linkage disequilibrium (LD) between the top SNP and any SNP tested is indicated by red shading. The recombination rate is displayed by a blue line with scale on the right-hand axis. Characterized genes in the region are represented with a green bar. The P value for rs9330813 (1.71 × 10−7) is shown as a red diamond.

To provide further support for the association of WNT7B with CCT in the South Indian pedigrees, we investigated the association of rs9330813 in the Latino study cohort as well as in four independent European data sets (Fig. 2). In addition, we investigated association of rs10453441 with CCT in an independent Singaporean Indian cohort, and five independent European data sets (Supplementary Fig. S4). The WNT7B SNPs were imputed from previous genotype data for the European cohorts. For both rs9330813 and rs10453441, association with CCT was evident with consistent direction of effects observed in all data sets with the exception of one European cohort for rs10453441 (Supplementary Fig. S4). For both SNPs, strongest association was observed for the South Indian and MAGGS (Latinos) data sets, with smaller effects in European cohorts (Fig. 2; Supplementary Fig. S4). Significant heterogeneity was detected among data sets, and ethnicity, study design, imputation quality score, and age and sex were evaluated for contribution. This analysis suggested that the heterogeneity was most likely caused by imputation quality and study design (meta-regression P = 0.0001 and P = 0.02, respectively). Limiting the meta-analysis to data sets with imputation scores >0.7 for each SNP reduced but did not completely eliminate heterogeneity (Fig. 2; Supplementary Fig. S4). Because of the residual heterogeneity reverse inverse weighted meta-analyses were completed using both fixed and random effects and investigated separately the data sets with imputation scores >0.7 for each SNP. Using the fixed effects model, significant association was observed for CCT and rs9330813 [A] (P = 1.7 × 10−9, β = −3.94, 95% CI: −5.23 to −2.66; Fig. 2), and rs10453441 [G] (P = 2.20 × 10−11, β = −3.11, 95% CI: −4.02 to −2.02; Supplementary Fig. S4). Evidence for association improved when only the data sets with imputation scores >0.7 were included in the meta-analyses: rs9330813[A] (P = 5.0 × 10−12, β = −5.59, 95% CI: −7.17 to −4.00; Fig. 2), and rs10453441 [G] (P = 5.3 × 10−12, β = −3.43, 95% CI: −4.40 to −2.45; Supplementary Fig. S4). Reduced but consistent association was observed using the random effects model for both SNPs: rs9330813 [A] (P = 7.0 × 10−3, β = −8.00, 95% CI: −13.85 to −2.15); rs10453441 [G] (P = 1.0 × 10−4, β = −3.44, 95% CI: −5.21 to −1.68).

Figure 2

Meta-analysis for rs9330813 and CCT. Forest plot showing effect estimates for the South Indian pedigree, as well as for the replication effort. Pooled estimates for β and 95% CI were calculated by fixed-effects, inverse variance weighting meta-analysis. Reduced evidence of association but with similar effects was observed if the meta-analysis was calculated using random effects: P = 7.0 × 10−3, β = −8.00, 95% CI: −13.85 to −3.15. Individual data set results are indicated by black squares and summary values are indicated by black diamonds. ORCADES, Orkney Complex Disease Study; TwinsUK, UK Twin Study.

The top WNT7B SNP, rs9330813 is in strong equilibrium with rs9723267; r2 = 0.96 and 1.0 in the South Indian data set, 1000 Genomes (provided in the public domain, http://csg.sph.umich.edu/abecasis/MACH/download/1000G.2012-03-14.html), and Haploreg v.4.1 (provided in the public domain, http://www.broadinstitute.org/mammals/haploreg/haploreg.php), respectively, that disrupts a Rad21 binding motif and a CTCF (CCCTC-binding factor) binding site, as well as other transcription factor binding sites (RegulomeDB, provided in the public domain, http://regulome.stanford.edu/; Supplementary Fig. S5) suggesting a role in regulation of gene expression. The region of intron 1 that includes the WNT7B SNPs associated with CCT contains multiple DNaseI hypersensitivity sites and features of enhancers as annotated by ENCODE in multiple cell types (Supplementary Fig. S5).

In the South Indian family data set, we also replicated association (P < 0.005) with a number of loci previously associated with CCT including RXRA-COL5A1,[16] ZNF469,[15] GPR15,[13] and GLT8D2,[13] although none of these associations were as significant as those observed for the WNT7B SNPs in this population (Supplementary Table S3). It was estimated that 53.8% of the variance in CCT was explained by all the CCT-associated SNPs in this Indian population.

We also investigated the association of the WNT7B SNPs associated with CCT in this study with primary open angle glaucoma (POAG) in our NEIGHBORHOOD European Caucasian data set of 3853 cases and 33,480 controls.[39] However, similar to other studies,[14] we did not find evidence for association of these SNPs with POAG (P > 0.05).

Recently, SNPs also located in this same region of the first intron of WNT7B have been associated with two other ocular quantitative traits, corneal curvature and AXL, in a GWAS using a Japanese population.[21] The lead SNP in the Japanese study, rs10453441, is the same SNP associated with CCT in the Latino study[20] located 422 bp from rs9330813, the lead SNP in the South Indian pedigrees (Supplementary Fig. S5). To determine if the WNT7B association in our data set was specific for CCT, we performed an age- and sex-adjusted PheWAS (Phenotype-wide association study) using association data for 45 ocular quantitative traits measured in the same families used for the CCT analysis (see Supplementary Table S1 for complete list of traits), including AXL and corneal curvature, the two traits associated with the WNT7B SNP rs10453441 in the Japanese study.[21] For the PheWAS, we investigated the top three WNT7B SNPs (rs9330813, rs9723267, and rs75159625) from our data (Supplementary Table S2) and also the top two SNPs in the Japanese study (rs10453441 and rs200329677). Four of these SNPs are preferentially associated with CCT in the South Indian sample (the remaining SNP, rs200329677, was not significantly associated with CCT or any other trait in this data set) (Fig. 3; Supplementary Fig. S6). In the South Indian data set, the PheWAS data did not support significant association of any WNT7B SNP with any trait other than CCT (P > 0.001) including AXL or corneal curvature as was observed in the Japanese study (Fig. 3; Supplementary Fig. S6) despite having sufficient power (>99.9%) for AXL and corneal curvature to detect the associations previously described (Supplementary Table S1).

Figure 3

PheWAS plot for the top SNP associated with CCT in the South Indian population (rs9330813). The association results for each measured trait (Supplementary Table S1) for this SNP were plotted with the phenotypes (ocular traits) grouped along the x-axis and the –log10(P) value for association analysis on the y-axis. The phenotype group is indicated by the color of the graph point as indicated by the side panel. The lower gray dashed line indicates P = 0.05. The upper black dashed line indicates a single-SNP Bonferroni correction for 45 traits, P = 0.001 (0.05/45). Other traits were not labeled in these figures due to limited space. Categories are grouped according to Supplementary Table S1. IOPg, intraocular pressure measured by Goldman applanation; CRF, corneal resistance factor; K_H, corneal curvature, horizontal; K_V, corneal curvature, vertical; RNFL_VC, retinal nerve fiber layer curvature as measured by the Heidelberg Retina Tomography and analyzed by using Glaucoma Probability Score (GPS).

Discussion

This is the first GWAS for CCT in individuals residing in Southern India, a population at increased risk for blinding ocular disorders.[27,40] In this family-based study that included large consanguineous pedigrees, we identified association of CCT with WNT7B SNPs located in an apparent regulatory region likely to impact gene expression. Pedigrees with consanguineous matings are known to have added power for genetic studies of recessive traits. In this study, we have shown that consanguineous families can also provide genetic insights leading to discovery of loci for quantitative traits. The CCT boxplot for three genotypes of top SNP rs9330813 was consistent with an additive model in this South Indian data set (Supplementary Fig. S7). We estimated that we had at least 82% power to detect the associations between these WNT7B SNPs and CCT in this South Indian data set.

WNT7B codes for a member of the Wnt family of proteins that have critical roles in cell growth, patterning, and differentiation of multiple tissues and organs.[41] The canonical WNT signaling pathway that includes WNT7b (the product of WNT7B) is known to contribute to stem cell proliferation in development.[42] In the eye, WNT7B has been shown to have increased expression in the central cornea and may also be necessary for corneal limbal stem cell development.[43] Interestingly, a rare exonic variant in another WNT family member, WNT10A, has also been associated with CCT in a quantitative trait study of European Caucasians.[44]

The WNT7B SNPs associated with CCT are located in the first intron of the gene in a region with multiple DNaseI hypersensitivity sites and enhancers as annotated by ENCODE. The top SNP is in strong linkage disequilibrium with rs9723267 that impacts Rad21 and CTCF (CCCTC-binding factor) binding sites. Rad21 is one of the subunits of the cohesin complex that together with CTCF associates with active enhancers and promoters forming long-range interactions important for gene regulation.[45] Rad21 and CTCF activity is highest when a general transcription factor (TBP) binding site is also nearby[46] as is the case in the WNT7B region associated with CCT (Supplementary Fig. S5), suggesting that genetic variants in this region could impact gene expression.

In addition to the association between WNT7B and CCT, we also confirmed association with several other loci previously associated with CCT in other populations, in particular ZNF469 and RXRA-COL5A1. Genomic association studies have now been completed for CCT in a variety of ethnic populations including European Caucasians,[13-16] Asians,[17] and Latinos.[18,20] Evidence for association of CCT with ZNF469 and RXRA-COL5A1 has been found in most populations, while other CCT loci such as COL8A2, significantly associated in Asians,[17] may be restricted to specific populations.[19] Our study suggests that WNT7B is an important locus for CCT in the South Indian population.

WNT7B SNPs may also contribute to other ocular phenotypes. In a study conducted in Japanese, SNPs in the same genomic region associated with CCT in our study were associated with AXL and corneal curvature, ocular quantitative traits related to refractive error and myopia.[21] We have previously measured 45 quantitative traits in the collection of Indian pedigrees used for this study including AXL, corneal curvature, and refractive error. This collection of quantitative trait data made it possible to complete an ocular PheWAS for the WNT7B SNPs associated with CCT in our study and the WNT7B SNPs associated with AXL and corneal curvature in the Japanese study. Understanding the range of phenotypic consequences of DNA sequence variants may provide insights into the mechanisms by which a variant or gene leads to disease. The PheWAS approach can test the association of a disease-associated variant with a broad range of phenotypes.[47-49] We found that in the South Indian population, the WNT7B SNPs are specifically associated with CCT and did not show evidence of association with any other traits, including those related to myopia and refractive error. While the Japanese study did not specifically interrogate association with CCT, it appears that the WNT7B SNPs can be associated with additional or different traits in the Japanese population. The opportunity to complete a PheWAS to evaluate the association of the WNT7B SNPs with a broad range of ocular phenotypes was a strength of our study.

Conclusions

In summary, our family-based association analysis using South Indian pedigrees has identified WNT7B as a locus for CCT in this population and an ocular PheWAS conducted in the same data set showed that the WNT7B association is specific for this trait in these South Indian pedigrees. WNT7B is known to be associated with CCT in a Latino population[20] but has not been previously shown to be a CCT locus in Asians or European Caucasians, suggesting that genomic studies in specific ethnic populations can uncover new loci for complex traits that provide additional insights into the underlying genetic architecture of these common conditions.

Click here for additional data file.