Molecular Vision 2017; 23:952-962
<http://www.molvis.org/molvis/v23/952>
Received 26 June 2017 | Accepted 13 December 2017 | Published 15 December 2017
Molecular Vision 2017; 23:952-962
<http://www.molvis.org/molvis/v23/952>
Received 26 June 2017 |
Accepted 13 December 2017 |
Published 15 December 2017
Matthew D. Benson,1 Chiea C. Khor,2 Philip J. Gage,3 Ordan J. Lehmann1,4
1Department of Ophthalmology and Visual Sciences, University of Alberta, Edmonton, Canada; 2Genome Institute of Singapore, Singapore; 3Department of Ophthalmology and Visual Science, Department of Cell and Developmental Biology, University of Michigan Medical School, Ann Arbor, MI; 4Department of Medical Genetics, University of Alberta, Edmonton, Canada
Correspondence to: Matthew D. Benson, Department of Ophthalmology and Visual Sciences, University of Alberta, Katz Building 7-030, Edmonton, Canada, T6E 2E1; Phone: 1-780-964-8548; FAX: 1-780-964-8548; email: mbenson@ualberta.ca
Purpose: To evaluate the ability of a targeted genome-wide association study (GWAS) to identify genes associated with central corneal thickness (CCT).
Methods: A targeted GWAS was used to investigate whether ten candidate genes with known roles in corneal development were associated with CCT in two Singaporean populations. The single nucleotide polymorphisms (SNPs) within a 500 kb interval encompassing each candidate were analyzed, and in light of the resulting data, members of the Wnt pathway were subsequently screened using similar methodology.
Results: Variants within the 500 kb interval encompassing three candidate genes, DKK1 (rs1896368, p=1.32×10−3), DKK2 (rs17510449, p=7.34×10−4), and FOXO1 (rs7326616, p=1.56×10−4 and rs4943785, p=1.19×10−3), were statistically significantly associated with CCT in the Singapore Indian population. DKK2 was statistically significantly associated with CCT in a separate Singapore Malaysian population (rs10015200, p=2.26×10−3). Analysis of Wnt signaling pathway genes in each population demonstrated that TCF7L2 (rs3814573, p=1.18×10−3), RYK (rs6763231, p=1.12×10−3 and rs4854785, p=1.11×10−3), and FZD8 (rs640827, p=5.17×10−4) were statistically significantly associated with CCT.
Conclusions: The targeted GWAS identified four genes (DKK1, DKK2, RYK, and FZD8) with novel associations with CCT and confirmed known associations with two genes, FOXO1 and TCF7L2. All six participate in the Wnt pathway, supporting a broader role for Wnt signaling in regulating the thickness of the cornea. In parallel, this study demonstrated that a hypothesis-driven candidate gene approach can identify associations in existing GWAS data sets.
Primary open-angle glaucoma (POAG) is characterized by optic nerve damage and progressive loss of vision, frequently because of elevated intraocular pressure (IOP). POAG represents a leading worldwide cause of irreversible blindness, with central corneal thickness (CCT) a key factor in the pathogenesis. Alterations in CCT affect the accuracy of IOP measurement, with decreased CCT resulting in IOP underestimation while the converse occurs with increased CCT. Such factors impact the diagnosis and treatment of glaucoma and provide an explanation for CCT being a strong predictor for the progression of ocular hypertension to primary open-angle glaucoma [1]. In this context, knowledge of the genetic factors that influence CCT has scientific importance for the understanding of ocular development and clinical relevance to glaucoma.
The cornea is a highly specialized, multilayered tissue with embryologic contributions from multiple origins. During early eye development, neural crest cells migrate between the developing lens and surface ectoderm with transcription factors, including FOXC1 (Gene ID: 2296; OMIM: 601090), LMX1B (Gene ID: 4010; OMIM: 602575), and PITX2 (Gene ID: 5308; OMIM: 601542), orchestrating key aspects of anterior segment development [2]. The corneal epithelium is formed by the differentiation of surface ectoderm, whereas the corneal stroma and the endothelium are derived from neural crest cells. Although altered CCT has been demonstrated in murine models and patients with mutations in Col8a1 [3], Col8a2 [3,4], Foxe3 [5], Pax6 [6], PITX2 [7], FOXC1 [5], and PAX6 (Gene ID: 5080; OMIM: 607108) [6], understanding of CCT genetics is largely derived from population-based association studies where greater than 26 CCT-associated loci have been reported [8-17]. In addition to the relationship between CCT and POAG, several of the genes associated with CCT also have demonstrated associations with keratoconus, a disorder of progressive corneal ectasia [12]. Genes with variants associated with CCT and keratoconus include COL5A1 (Gene ID: 1289; OMIM: 120215), FOXO1 (Gene ID: 2308; OMIM: 136533), ZNF469 (Gene ID: 84627; OMIM: 612078), among others [12]. As the molecularly undefined proportion of CCT is estimated to be up to 90% [12], there is an opportunity to utilize knowledge of corneal development to identify genes that may influence CCT.
One approach for identifying variants, and thus genes, that influence a complex trait such as CCT is a genome-wide association study (GWAS). This hypothesis-free method analyzes single nucleotide polymorphisms (SNPs) across the genome in cases and controls to determine if any are associated with the phenotype being studied. Generally, 1 to 2 million SNPs are used, with a statistical correction applied based on the number of SNPs (107) analyzed. In this example, only associations with a probability of less than 10−7 are considered statistically significant. GWASs have been successful in advancing gene discovery in a range of disorders [18,19]; however, this approach is best suited for genes of major effect. A GWAS is less successful in genetically heterogeneous disorders where an individual gene may account for less than 1% of a phenotype. In such scenarios, the contribution of the Bonferroni correction (due to the 1 million or more SNPs analyzed) and the cohort size (usually fewer than 10,000 patients) frequently make studies underpowered [20]. Because limits exist on the recruitment of ever larger patient cohorts, and the cost of genotyping remains substantial, there is a need for alternative approaches, especially those that take advantage of existing data sets.
Accordingly, the possibility that real associations exist but they are masked by the level of correction needed for a genome-wide study was evaluated [21]. Therefore, instead of conducting a hypothesis-free GWAS, this study investigated whether knowledge of ocular development could be used. Whether candidate genes with anterior segment developmental roles influenced corneal thickness was determined by analyzing regional SNP data from existing GWAS data sets. SNPs in small intervals encompassing each candidate were investigated to determine if they were statistically significantly associated with CCT. The results demonstrated the utility of the approach as a means of generating novel findings that merit further experimental validation and implicate members of the Wnt pathway in the control of corneal thickness at the population level.
GWAS data sets derived from two Singaporean populations of independent ethnicity (Indian and Malay) were used for the analyses. These populations comprised 2,538 individuals from the Singapore Indian Eye Study (SINDI) [22] and 2,542 individuals from the Singapore Malay Eye Study (SiMES) [23]. Ethics approval was obtained from the Singapore Eye Research Institute (SERI) Institutional Review Board (IRB). In the studies, median CCT had been determined with ultrasound pachymetry using five measurements in the right eye of each patient. Population stratification was well controlled for, and there was no genomic inflation [15].
Ten candidate genes were selected based on their roles in the development of the anterior segment: crumbs homolog 1 (CRB1; Gene ID: 23418; OMIM: 604210) [24], dickkopf-related proteins 1 and 2 (DKK1; Gene ID: 22943; OMIM: 605189 [25] and DKK2; Gene ID: 27123; OMIM: 605415[26]), forkhead box protein E3 (FOXE3; Gene ID: 2301; OMIM: 601094) [27], LIM homeobox transcription factor 1-beta (LMX1B) [27], paired box protein 6 (PAX6) [27], paired-like homeodomain transcription factor 2 (PITX2) [27], forkhead box protein O1 (FOXO1) [28], microRNA 184 (miR184; Gene ID: 406960; OMIM: 613146) [29], and forkhead box protein C1 (FOXC1) [27]. To enhance rigor, ten control genes were selected for comparable analysis, and based on their known functions, these genes were anticipated to not contribute to the development of the anterior segment. The genes comprised salivary amylase 1-alpha (AMY1A; Gene ID: 276; OMIM: 104700) [30], cholecystokinin-A receptor (CCKAR; Gene ID: 886; OMIM: 118444) which regulates gallbladder motility [31], five immunoglobulin chain genes (IGHD4–4; Gene ID: 28496 [32], IGHV3–21; Gene ID: 28444[32], IGHV7–81; Gene ID: 28378 [32], IGKV2–24; Gene ID: 28923 [33], and IGKV3D-7; Gene ID: 28877[34]) and three mitochondrial ribosomal proteins (MRPL10; Gene ID: 124995; OMIM: 611825[35], MRPL11; Gene ID: 65003; OMIM: 611826 [36], and MRPL12; Gene ID: 6182; OMIM: 602375[37]).
An interval of 500 kb encompassing each gene was analyzed [38]; an interval that would generally include the promoter and multiple long distance regulatory elements, in addition to the coding region. PLINK (version 1.07; Cambridge, MA), a C/C++ open source GWAS analysis program, was used for evaluating associations [39]. The software analyzed the association of genotypes and CCT, a quantitative trait, by computing a conventional linear regression equation correcting for age and sex. The genotype of each individual was scripted based on the number of variant alleles present where 0 represented wild-type, 1 represented a heterozygote for the variant, and 2 represented a homozygote for the variant. A minor allele frequency of 0.05 was used in the association analysis. Data were presented graphically with a line indicating the threshold for the statistical significance for each interval.
Haploview software (version 4.2; Cambridge, MA) was used to analyze haplotypes, as well as to calculate and visualize linkage disequilibrium (LD) blocks to allow for correction for multiple comparisons [40,41]. To obtain a corrected significance threshold, standard methodology was used to correct for multiple SNP comparisons [42-46]. This involved dividing the conventional significance level (α=0.05) by the sum of the number of LD blocks and inter-block SNPs in the 500 kb interval. This approach appropriately accounts for the effective number of independent SNPs due to non-random variant associations and allele proximity and offers increased statistical power [42] compared to the overly conservative Bonferroni method. As an example of this method, in the Singapore Indian population, the 500 kb interval containing DKK1 had 36 LD blocks and inter-block SNPs; therefore, the corrected significance level, α (corr), is 0.05/36=1.39×10−3.
Variants reaching statistical significance in individual populations were evaluated in a meta-analysis that combined the Singapore India and Malaysian populations (n=5,080). The mean regression coefficient, weighted based on population size, and standard deviation were determined for the linear regression models for the variants, and these were used to calculate 95% confidence intervals. The results were depicted in a forest plot derived from an existing template [47].
A comparison of CCT in murine Dkk2+/− variants and littermate C57BL/6J wild-type mice was performed using murine variants that had been bred and generated as previously described [48]. The study adhered to the Association for Research in Vision and Ophthalmology (ARVO) Statement for Use in Animals in Research. The investigations using mice were approved by the University of Michigan Committee on Use and Care of Animals. Approximately equal distributions of male and female murine eyes were enucleated at 38 to 48 days of age, fixed in 4% paraformaldehyde, and processed into JB-4 plastic. Central corneal sections were collected, visualized using Lee’s stain, with the stained sections examined, and high magnification images taken using a Leica DM6000 (Buffalo Grove, IL) microscope. After the images were imported into MetaMorph software (Sunnyvale, CA), a total of three measurements (left, middle, and right) of the CCT were taken and averaged. Subsequent comparison of the wild-type (n=10) and heterozygote (n=8) mice CCT was performed using a two-tailed t test with the assumption of two-sample equal variance. Animal care and research protocols are in keeping with guidelines published by the Institute for Laboratory Animal Research.
To evaluate the relationship between members of the Wnt signal transduction pathway and CCT, Wnt genes comprising ligands, receptors, and antagonists were sampled for association with CCT [49]. To narrow the focus to those Wnt genes more likely to be associated with CCT, a meta-analysis was performed for the Singapore Malaysian and Indian populations. Those genes with small p values (less than 0.025) were then studied in each population separately as it was deemed that any p value larger than 0.025 in the meta-analysis was unlikely to be statistically significant in each population following LD correction. Association analysis was performed using the PLINK software with the application of conventional linear regression and LD correction as described above for nine genes: receptor tyrosine kinase-like orphan receptor 2 (ROR2; Gene ID: 4920; OMIM: 602337), receptor-like tyrosine kinase (RYK; Gene ID: 6259; OMIM: 600524), wingless-type MMTV integration site family, member 5B (WNT5B; Gene ID: 81029; OMIM: 606361), low density lipoprotein receptor-related protein 5 (LRP5; Gene ID: 4041; OMIM: 603506), transcription factor 7-like 2 (TCF7L2; Gene ID: 6934; OMIM: 602228), frizzled class receptor 8 (FZD8; Gene ID: 8325; OMIM: 606146), dickkopf WNT signaling pathway inhibitor 3 (DKK3; Gene ID: 27122; OMIM: 605416), transcription factor 7-like 1 (TCF7L1; Gene ID: 83439; OMIM: 604652), and casein kinase 1, alpha 1 (CSNK1A1; Gene ID: 1452; OMIM: 600505).
A variant within a 500 kb interval containing DKK1 was statistically significantly associated with CCT following correction by linkage disequilibrium in the Singapore Indian population (rs1896368, β=3.02, p=1.32×10−3; Figure 1A) with a significance level corrected to α (corr)=1.39×10−3. This variant failed to reach statistical significance in the Singapore Malaysian population (p=0.990; Figure 1B).
In addition, two variants within a 500 kb interval and adjacent to the coding region of DKK2 were associated with CCT, one in the Singapore Indian population (rs17510449, β=–3.24, p=7.34×10−4, α (corr)=2.08×10−3; Figure 1C) and the other in the Singapore Malaysian population (rs10015200, β=3.26, p=2.26×10−3, α (corr)=2.5×10−3; Figure 1D). Neither rs17510449 nor rs10015200 reached statistical significance in either population (p=0.536, p=0.857, respectively; Figure 1C,D).
Finally, two variants in the 500 kb interval containing FOXO1 were statistically significantly associated with CCT (rs7326616, β=–7.43, p=1.56×10−4 and rs4943785, β=–3.72, p=1.19×10−3; α (corr)=1.35×10−3) in the Singapore Indian population (Figure 1E), in keeping with a previous study [11]. These variants, rs7326616 and rs4943785, failed to reach statistical significance in the Singapore Malaysian population (p=0.544, p=0.448, respectively; Figure 1F). Only one of the ten control genes had variants within the 500 kb interval that were statistically significantly associated with CCT (Appendix 1). This was detected only in the Singapore Indian population for the mitochondrial ribosomal protein MRPL11, with no association observed for either of its two paralogs, MRPL10 and MRPL12 (Appendix 1). The overall effect of the identified significant variants when the Singapore India and Malaysian populations were combined in a meta-analysis did not reach statistical significance (Appendix 2).
Given the identification of variants in DKK2 that were statistically significantly associated with CCT in the Singapore Indian and Malaysian populations, the CCT of Dkk2 heterozygous mice was compared to wild-type littermates. However, the mean corneal thickness in heterozygote mice did not differ significantly from that of the wild-type mice (Dkk2+/− 116.33 µm (SD=18.710 µm), WT 117.90 µm (SD=10.960 µm), p=0.478; Appendix 3).
As SNPs near DKK1 and DKK2 were associated with CCT and both genes antagonize Wnt signaling, additional members of the Wnt pathway were screened for association with corneal thickness. In the Singapore Indian population, a variant adjacent to the coding region of FZD8 was statistically significantly associated with CCT following correction by linkage disequilibrium (rs640827, β=–3.49, p=5.17×10−4, α (corr)=1.56×10−3; Figure 2A). This variant failed to reach statistical significance in the Singapore Malaysian population (p=4.34×10−2; Figure 2B). The remainder of the Wnt genes (n=8) screened in the Singapore Indian population failed to reach statistical significance. In the Singapore Malaysian population, variants in the 500 kb interval containing RYK (rs6763231, β=4.35, p=1.12×10−3 and rs4854785, β=–3.05, p=1.11×10−3, α (corr)=1.28×10−3) and within the intragenic region of TCF7L2 (rs3814573, β=–3.17, p=1.18×10−3, α (corr)=1.47×10−3), in keeping with a previous study [50], were associated with CCT following statistical correction (Figure 2F,D). The remaining Wnt genes (n=7) evaluated in the Singapore Malaysian population were not associated with the CCT.
Corneal thickness is a highly heritable trait, and several studies have identified variants and genetic loci that are associated with CCT [12,14-16]. However, most of these studies are GWASs and looked for associations at the genome level only. The possibility that genuine associations with CCT are masked by the standard Bonferroni correction used in a GWAS was tested using a focused approach that targeted genes with roles in corneal development. Notably, targeted reanalysis of GWAS data sets has successfully defined the genetic basis of disease ranging from alcohol dependence to acute respiratory distress syndrome (ARDS) [51-54].
From restricted interval analysis, three genes (DKK1 (rs1896368, p=1.32×10−3), DKK2 (rs10015200, p=2.26×10−3; rs17510449, p=7.34×10−4), and FOXO1 (rs7326616, p=1.56×10−4; rs4943785, p=1.19×10−3) out of ten candidate genes were identified as being statistically significantly associated with CCT. Only one control gene, MRPL11, had variants within the 500 kb interval that were significantly associated with CCT; the remaining nine controls did not reach statistical significance. This result agrees with published data that implicated a role in cornea formation for DKK1 [25] and DKK2 [26,48] and described an association between FOXO1 and CCT in a large-scale GWAS [11,12]. Unlike FOXO1, DKK1 and DKK2 represent novel associations with CCT; these associations were not detected in the original GWAS [15]. The association of variants adjacent to MRPL11 with CCT in the Singapore Indian population was unexpected. Because this control gene has no known or putative role in corneal development, and the number of candidate genes exhibiting an association with CCT was substantially higher than in the controls, this finding most likely represents a chance association. Taken together, the present findings suggest that a hypothesis-driven focused analysis of candidate genes can elicit associations that may not be detected with conventional genome-wide analysis.
Variants that were statistically significantly associated with CCT were identified in the Singapore Indian and Malaysian populations; however, none of the identified variants reached statistical significance in either population (Appendix 2). Given that there are inherent ethnic differences in corneal thickness [55,56], it is unsurprising that common variants are associated with CCT in one ethnic group but not in another. Furthermore, the identification of significant variants in one population but not in the other may also reflect an effect of a relatively small cohort size which reduces the statistical power to detect a true association. Such a limitation in the study may be addressed by analyzing additional populations to broaden the sample size and further evaluate this diversity in CCT in a meta-analysis.
Given that DKK2 was statistically significantly associated with CCT, Dkk2 heterozygous mice were used to further examine this potential relationship. Corneal thicknesses of wild-type mice and Dkk2+/− were compared; however, this study failed to detect a difference in mean CCT between these groups (Appendix 3). Several limiting factors may account for this observation. At the technical level, corneal thickness may be affected by the tissue processing needed to create the histological slides. In addition, generating histological sections perpendicular to the cornea is challenging in 2 mm murine globes, and the wide scatter in data points supports this. Substantially larger sample sizes coupled with the use of either optical coherence tomography (OCT) or in vivo approaches (ultrasonic pachymetry) to directly measure corneal thickness would address this. Furthermore, we used mice with a single functional Dkk2 allele to represent SNPs present in patients. It is possible that mice, with a single null mutation, may not accurately recapitulate the effect of a common and much milder variant. Future ex vivo and in vivo studies with Dkk2+/− mice and with Dkk2−/− mice using larger sample sizes will be beneficial to determine whether variants in these genes have a definite causative effect on murine CCT.
DKK1 and DKK2 belong to the Wnt signaling pathway, a signal transduction pathway involved in embryogenesis, growth, and development, where they inhibit canonical Wnt signaling through interactions with Kremen transmembrane proteins and the LRP6 coreceptor [26]. Wnt signaling controls multiple stages of ocular development, including patterning of the optic vesicle, lens formation, and corneal differentiation [57]. Accordingly, nine of 54 members of the Wnt pathway were screened for an association with CCT, and three genes (FZD8, TCF7L2, and RYK) were identified as being statistically significantly associated with CCT (Figure 2A,D,F, respectively). Although FZD8 is expressed in the corneal epithelium [58] and variants in TCF7L2 have been associated with Fuchs’ corneal endothelial dystrophy [59], RYK does not have a reported corneal role. Interestingly, similar to DKK1 and DKK2, TCF7L2, RYK, and FZD8 have all been shown to function in inhibitory roles in canonical Wnt signaling [60-62]. Given this related function, it is plausible that canonical Wnt inhibition may underlie differences in CCT.
CCT is an important prognostic indicator in patients with primary open-angle glaucoma, and several studies have demonstrated more advanced signs of glaucoma with thinner corneas [63-65]. A randomized prospective study of patients with ocular hypertension demonstrated that a thin CCT was a strong predictor for the development of open-angle glaucoma [1]. Knowledge of the genetic factors that influence CCT, including the identification of association of DKK1 and DKK2 with CCT, is thus important not only for enhancing the scientific understanding of ocular development but also for establishing the prognosis and the risk of progression of disease in patients with ocular hypertension and POAG. In this sense, an understanding of the genetic factors of CCT can allow for individual risk stratification, optimizing the management and treatment decisions for patients with glaucoma.
Several potential limitations of the present study should be discussed. First, we used a candidate-driven approach to identify genetic associations with CCT. As a result, this approach would not be helpful without existing knowledge of genes involved in anterior segment development. Thus, this candidate-driven approach complements and does not replace information derived from candidate-independent GWASs. Second, this study examined CCT associations in two Singaporean populations. Further replication of these findings in other data sets with larger sample sizes would be of benefit to enhance reproducibility and allow for generalization to other populations.
In conclusion, variants near six genes (DKK1, DKK2, FOXO1, RYK, TCF7L2, and FZD8) are associated with CCT in the targeted GWAS reanalysis. The fact that the variants identified were not associated with the CCT in either population may reflect genetic diversity between ethnic groups, the effects of relatively small sample sizes, or other factors. An additional limitation of this study, given the nature of association and regression models, is the inability to ascribe causative roles to the identified genes in the CCT. More comprehensive in vivo analyses of vertebrate models (zebrafish or murine) may clarify the relationship of individual genes with CCT. Finally, given that only a subset of Wnt genes were screened, the possibility exists that other members may also be associated with CCT. A more detailed future analysis of the Wnt family may contribute further to the understanding of genetic associations with corneal thickness.
In summary, this study demonstrates the utility of hypothesis-driven mining of existing GWAS data, identifying plausible associations that merit further study. Potentially, such approaches could be used in existing GWAS data sets to elucidate the genetic etiologies of a broad range of ocular phenotypes.
Appendix 2. Singaporean Indian and Malaysian meta-analysis.
Appendix 3. Central corneal thickness in Dkk2 heterozygous mice.
The authors thank Dr. Curtis French and Dr. Sameer Pant for their technical support. Research supported by grants from Alberta Innovates-Health Solutions. The authors affirm that they have no competing interests.