Introduction

Keratoconus is a noninflammatory degenerative disorder that results in bulging and distortion of the corneal surface, leading to irregular astigmatism and progressive myopia. In advanced cases, corneal scarring and even corneal blindness can occur. Keratoconus has an incidence of approximately 1 in 2,000 individuals with a prevalence varying from 8.8 to 2300 per 100,0001, 2. It is a leading indication for corneal transplantation in many countries, especially in Australia, Middle East and Africa3. Management of keratoconus varies from conservative visual correction by spectacles or contact lenses for mild disease, to surgical interventions such as collagen cross-linking, intracorneal rings and keratoplasty for advanced disease. The onset of keratoconus is insidious and the progression is irreversible. Therefore, early diagnosis of keratoconus and its progression is needed. However, the variable risk of keratoconus progression poses a challenge to the personalized management for patients4. Knowing the risk factors for keratoconus would thus be helpful for early detection and monitoring the progression of the disease.

Keratoconus is a multifactorial disease resulting from the interaction of environmental, behavioural and genetic factors. Major environmental and behavioural factors include contact lens wear5 and chronic eye rubbing6. The genetic aetiology is evidenced by the bilaterality, familial aggregation7,8,9, monozygotic twins concordant of the disease10, its association with other genetic diseases such as Down syndrome11 and Leber’s congenital amaurosis12, and the ethnic difference in the prevalence and incidences13. Genetic associations for keratoconus will provide insight into disease mechanisms and help identify biomarkers for early detection of keratoconus onset and monitoring its progression. Of note, about 14% of the patients with keratoconus have a family history9. So far, however, the difference in the genetic basis of familial and sporadic keratoconus is unclear. Since the family history does not affect disease severity, the pooling of all cases in genetic studies is deemed reasonable14.

So far, 6 chromosomal loci have been identified for isolated keratoconus by linkage analysis, namely 2p24 15, 3p14-q13 16, 5q14.3-q21.1 2, 13q32 18, 16q22.3-q23.1 19, and 20q12 20. However, no disease-causing mutation has been identified from these loci. Besides, genome-wide association studies (GWAS) and candidate gene association studies have reported over 150 polymorphisms in more than 60 genes/loci for keratoconus. Among them, 7 genes/loci were identified by GWAS, including the HGF 21, LOX 22, FOXO1 and FNDC3B genes23, and the 3p26, 2q21.3 and 19q13.3 loci24. However, most of these associations were inconsistent across different study cohorts, making the roles of the genes/loci inconclusive.

In this study, we conducted a systematic review and meta-analysis to summarize the genetic association evidence for all variants in genes that were previously reported for keratoconus, and evaluated potential trans-ethnic heterogeneities. We first presented the association results from selected original studies/cohorts in the forest plots and then provided a prioritized list of studies and genes variants for further analysis. For SNPs that have been meta-analyzed in prior studies, our study provides an update of the summary association results by including new studies.

Results

Selection of studies

We retrieved a total of 978 records published between 1950 and 1 June 2016 from MEDLINE, Embase, Web of Science, and HuGENET for review. After removing 339 duplicated records we evaluated 639 citations and selected 36 articles for full-text assessment. Among them, 2 were reviews25, 26 and 32 were molecular genetic studies, including 2 GWAS23, 24 and 30 candidate gene association studies. A total of 64 genes/loci and 156 variants have been identified from the full-text review (Supplementary Table 1; Fig. 1). In the meta-analysis, we excluded 12 of the 36 articles because 7 of them were about gene variants that were not tested in additional independent studies27,28,29,30,31,32,33, 2 reported insufficient genotype data for meta-analysis34, 35, 2 were reviews25, 26, and 1 was an animal study36. We did not receive genotype data after contacting some of the authors34, 35. Finally, 24 studies were included for meta-analysis, involving a total of 53 SNPs in 28 genes/loci (Fig. 1)21,22,23,24, 32, 37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54. Among these 24 studies, 20 were candidate gene studies conducted in different populations, including Whites39, 41, 43, 44, 47, 49,50,51,52, 55, Arabic37, 38, 40, Chinese45, 56, Korean53, 54, Japanese48, Indian46, and Turkish42. The total sample sizes from these candidate gene studies were 3,037 patients with keratoconus and 9,928 controls. The 2 GWAS included 2,333 keratoconus patients and 16,655 controls of Caucasian origin (Table 1)23, 24.

Figure 1
figure 1

Flow diagram of study selection process.

Table 1 Characteristics of eligible studies for the meta-analysis.

Genes reported in keratoconus GWAS

We first meta-analyzed the SNPs that were reported in the four keratoconus GWAS23, 24 and additional independent studies based on the GWAS21, 22, 39,40,41, 47, 49, 55, 56. A total of 27 SNPs in 22 genes/loci were involved. Among them, 16 SNPs in 14 genes/loci showed a summary P value < 0.05 (Table 2). Of note, 3 SNPs in 3 respective genes/loci reached genome-wide significance, including FOXO1 rs2721051 (P = 5.6 × 10−11, I2 = 0), RXRA-COL5A1 rs1536482 (P = 2.5 × 10−9, I2 = 0), and FNDC3B rs4894535 (P = 1.4 × 10−8, I2 = 0) (Table 2 and Fig. 2). The P values for the remaining 13 significantly-associated SNPs ranged from 6.1 × 10−7 (IMMP2L rs757219) to 0.035 (19p12 rs8111998) (Table 2).

Table 2 Allelic associations of gene variations with keratoconus using cohorts from both GWAS and subsequent replication studies.
Figure 2
figure 2

Meta-analysis of the 5 SNPs in 4 genes/loci showed genome-wide significance. Of the 4 genes/loci, 3 were detected in genome-wide association studies, including (A) FOXO1 (rs2721051, P = 5.6 × 10−11, I2 = 0), (B) RXRA-COL5A1 (rs1536482, P = 2.5 × 10−9, I2 = 0) and (C) FNDC3B (rs4894535, P = 1.4 × 10−8, I2 = 0). The (D) COL4A4 (rs2229813, P = 1.3 × 10−12, I2 = 0) gene was identified by candidate gene analysis.

We then performed meta-analysis only using the candidate gene studies, including those based on the GWAS findings or other hypotheses. One SNP from GWAS was significantly associated with keratoconus, i.e., RXRA-COL5A1 rs1536482 (P = 1.5 × 10−5, I2 = 0), while FOXO1 rs2721051 (P = 9.4 × 10−3, I2 = 0), BANP-ZNF469 rs9938149 (P = 0.017, I2 = 27%), COL4A4 (rs2228557, P = 0.020, I2 = 70%) and COL4A3 (c.2685 A > C, P = 0.032, I2 = 0) were nominally significant (Table 3). One SNP, FNDC3B rs4894535, reached a genome-wide significance in the overall population but did not show a significant association in the pooled Chinese and Arabic samples (P = 0.078, I2 = 0; Table 3). The other 4 genes/loci that have been reported in GWAS (i.e., MPDZ-NFIB, COL5A1, LOX and HGF) were also insignificant (P > 0.050; Table 3).

Table 3 Allelic associations of gene variations with keratoconus based on purely candidate gene studies.

Stratification analysis

To reduce the potential impact of trans-ethnical heterogeneity to the overall genetic association, we grouped the study cohorts into Whites and others (including Chinese, Korean, Japanese, Indian and Arabic). The 5 SNPs that were identified from GWAS showed a robust or nominal significance in Whites: FOXO1 rs2721051 (P = 1.5 × 10−9, I2 = 11%), MPDZ-NFIB rs1324183 (P = 1.8 × 10−4, I2 = 49%), BANP-ZNF469 rs9938149 (P = 2.6 × 10−4, I2 = 42%), COL5A1 rs7044529 (P = 9.9 × 10−4, I2 = 12%) and HGF rs3735520 (P = 3.6 × 10−3, I2 = 66%; Supplementary Table 2). Moreover, 2 SNPs in the COL4A4 gene identified by candidate gene studies were strongly associated with keratoconus in Whites, namely rs2229813 (P = 1.3 × 10−12, odds ratio (OR) = 2.38; I2 = 0) and rs2228557 (P = 4.5 × 10−7, OR = 0.54; I2 = 0) (Supplementary Table 2 and Fig. 2). In contrast, SNP rs2229813 showed a nominal association with keratoconus in combined Chinese and Arabic samples (P = 0.047, I2 = 16%). The odds ratio was notably toward an opposite direction (OR = 0.74; Supplementary Table 2). Moreover, most of aforementioned significant SNPs in Whites were not significant in the Chinese and Arabic samples, including FOXO1 rs2721051 (P = 0.31; I2 = 0), BANP-ZNF469 rs9938149 (P = 0.32; I2 = 0), MPDZ-NFIB rs1324183 (P = 0.63; I2 = 82%), and COL5A1 rs7044529 (P = 0.95; I2 = 0) (Supplementary Table 2), indicating ethnic diversities.

In this study, we were not able to evaluate the potential difference in the genetic basis of familial and sporadic cases as the data from familial cases were limited.

Assessment of potential biases and sensitivity analysis

For quality assessment every study was awarded a star for each of the items, i.e., case definition, ethnicity, and ascertainment of genotype (Supplementary Table 3) according to the Newcastle Ottawa Scale (NOS) system. All the 24 studies were awarded 5 or more stars out of a maximum of 8. Regarding Hardy-Weinberg Equilibrium (HWE), the control groups in 3 study cohorts showed deviation from HWE when tested for FOXO1 (rs2721051)56, COL4A3 (rs10178458 and rs55703767)52, COL4A4 (rs2229813, rs2228555, and rs2229814)52, and VSX1 (rs12480307)45. Therefore, in the sensitivity analysis we first excluded all the cohorts with HWE deviation and recalculated the summary ORs for the 7 SNPs in 4 genes. The associations were not altered (Supplementary Table 4). Furthermore, we omitted each study one at a time sequentially and recalculated the summary outcomes. The significance or insignificance of the summary outcomes was not altered in the sensitivity analysis (data not shown). We did not detect significant small study effects (e.g. publication bias) according to the shapes of funnel plots (Supplementary Figure 1) and the P values from the Egger’s tests, except for COL4A3 (rs55703767), LOX (rs2956540) and VSX1 (rs6138482) in the subgroup analysis by ethnicity (Supplementary Table 1).

Discussion

In this study, we meta-analyzed a total of 53 SNPs in 28 genes/loci for their genetic associations with keratoconus. We identified 8 SNPs in 6 genes/loci that were associated with keratoconus, i.e., FOXO1 rs2721051, FNDC3B rs4894535 and BANP-ZNF469 rs9938149 for the overall combined cohorts, and RXRA-COL5A1 rs1536482, IMMP2L rs757219 and rs214884, and COL4A4 rs2229813 and rs2228557 for Whites. Also, we found nominally significant associations in another 10 genes/loci, including KCND3, RAB3GAP1, UBXD2, MPDZ-NFIB, COL5A1, LOX, HGF, COL4A3, 13q33.3, and 19p12. In contrast, SNPs in 10 genes/loci that were reportedly associated with keratoconus were insignificant in our meta-analysis, including BHLHB2, BIRC8, IL1A, IL1B, KIF26B, LRRN1, PPP3CA, VSX1, 12p13.3 and 3q26.2.

Among the 6 significant genes/loci for keratoconus, 5 were originally identified by GWAS, including FOXO1, FNDC3B, BANP-ZNF469, RXRA-COL5A1, and IMMP2L. In our meta-analysis involving data from the GWAS and independent replication studies, 3 genes/loci (i.e., FOXO1, FNDC3B, BANP-ZNF469) showed consistent effects with low heterogeneity across different study cohorts. Three of them, FOXO1 rs2721051, FNDC3B rs4894535 and BANP-ZNF469 rs9938149, have been tested in both Whites and Asian populations. However, none of them showed a significant association in Chinese32 or Arabs40. Of note, FOXO1 rs2721051 was rare in Chinese with a minor allele frequency of 0.1%32. The lack of significant association in Asians could be due to the small sample size. In this meta-analysis, we also identified a SNP rs2229813 in the COL4A4 gene that showed a summary P value of genome-wide significant in Whites (P = 1.3 × 10−12; OR = 2.38). This gene was identified only in candidate gene studies37, 43, 45, 50, 52. Interestingly the summary P value in the pooled non-Caucasian samples was nominally significant (P = 0.047), but the OR was toward the opposite direction (OR = 0.74). This may suggest trans-ethnic diversities in the genetic components of keratoconus. In the COL4A4 gene, another SNP rs2228557, which was proposed in candidate gene studies, showed a significant summary P value (P = 4.5 × 10−7) in Whites, suggesting COL4A4 could be a genuine susceptibility gene for keratoconus in Whites. However, rs2228557 has only been tested in a Chinese population showing an insignificant association with an opposite OR (1.09)45. Therefore, whether COL4A4 is a biomarker with differential effects on keratoconus among different ethnic groups has yet to be confirmed. Interestingly, these 2 SNPs (i.e., rs2229813 and rs2228557) have not been reported in the published GWAS papers. In GWAS, only SNPs with P values surpassing a certain threshold would have been subjected to replication. Therefore, it would be intriguing to check the COL4A4 SNPs in the GWAS data and assess their association with keratoconus.

Although we were not able to evaluate the potential difference in the genetic basis of familial and sporadic cases, we found 2 familial cohorts being tested for different genes/loci22, 24, 50, 51, 55. In 3 studies22, 24, 55, the authors tested the associations of a few genes/loci (e.g. LOX and COL5A1) with keratoconus in a familial cohort using a generalized estimating equation accounting for familial correlations. Some of the significant SNPs identified in unrelated cases also showed significant association with keratoconus in the familial cohort. In another 2 studies50, 51, the authors reported a mutation, “627 + 23 G > A”, in VSX1 that was segregated in cases in several families. However, the mutation did not show significant association with keratoconus in the analysis using all the cases50. The results from the 2 cohorts indicated that the genetic association profiles of sporadic and familial keratoconus could be different.

Results of the present meta-analysis have led to a list of genes and loci associated with keratoconus that can be considered for functional investigations. Further biological investigation on these genes may throw light on new disease mechanisms for keratoconus. For example, FOXO1, RXRA and FNDC3B are the 3 genes that showed genome-wide significant association with keratoconus. FOXO1 is a member of the Forkhead Box (Fox) transcription factor family. Proteins from this family contain a conserved forkhead domain, which is a 110 amino acid DNA-binding domain. Fox proteins are known to be important regulators of the cellular oxidative stress57. For example, Fox proteins regulate the expressions of anti-oxidative enzymes such as superoxide dismutase and thioredoxin reductase58, 59. Moreover, reduced FOXO1 expression has been reported to induce apoptosis in human trabecular meshwork cells in response to oxidative stress60. It has been shown that increased oxidative damage to trabecular meshwork cells results in elevation of intraocular pressure and changing the anterior chamber angle, which would lead to corneal thinning61. We also found association of keratoconus with IMMP2L, a mitochondrial inner membrane protease. Mutation in IMMP2L also accumulates oxidative stress62. Therefore, FOXO1 and IMMP2L might regulate the oxidative stress in the anterior chamber, which affects the intraocular pressure and the corneal thickness. FOXO1 has also been linked to adipocyte differentiation63, which is affected by the gene FNDC3B 64. In this study, FNDC3B is another keratoconus associated gene. The link between adipogenesis and keratoconus is currently unclear. However, FNDC3B was associated with elevated intraocular pressure in a GWAS study65. Hence, FNDC3B may influence the intraocular pressure, the anterior chamber angle and the corneal thickness. Another keratoconus gene is RXRA, which encodes a nuclear retinoic acid receptor protein. There are two classes of nuclear retinoic acid receptors: RXR and RAR, which bind to each other to form RXR/RAR heterodimers66. Null mice of both RXRA and RXRA/RAR showed abnormal embryonic eye morphologies, including thickening of corneal stroma and absence of anterior chamber66. These results suggest a potential role of RXRA and retinoic acid signaling in the ocular development. However, the link among retinoic acid signalling, ocular development, and the abnormal corneal in keratoconus remains to be explored.

It is noteworthy that all of the identified SNPs in the 16 genes/loci are located in intronic, inter-genic, or in 3′- or 5′-untranslated regions. One SNP in HGF, rs3735520 (c.−1652C > T), was reported to modulate the severity of interstitial lung disease in patients with systemic sclerosis by altering the transcriptional efficiency of the HGF gene67. Whether they are in linkage disequilibrium with other coding variants in the relevant genes remained to be elucidated by sequencing analyses.

Although the mechanisms underlying the significant associations of the 16 identified genes/loci with keratoconus are largely unknown, it might be useful for understanding their pathogenic effects by referring to disease pathways identified for other conditions that share the same genes/loci. Eleven genes have been implicated in other diseases, including: COL5A1 for Ehlers-Danlos syndrome68; COL4A3 and COL4A4 for Alport syndrome69; HGF for non-syndromic hearing loss70; IMMP2L for Gilles de la Tourette syndrome71; KCND3 for spinocerebellar ataxia72; LOX for thoracic aortic aneurysms and dissections73; MPDZ for leber congenital amaurosis and retinitis pigmentosa74; RAB3GAP1 for Warburg Micro syndrome and Martsolf syndrome75; and ZNF469 for Brittle cornea syndrome76. The other 6 of the 16 identified genes, namely FOXO1, RXRA, FNDC3B, BANP, UBXD2, and NFIB of the MPDZ-NFIB locus, have not been directly linked to other human diseases.

In this study, we have identified and evaluated the genetic associations for keratoconus by conducting thorough assessments of the existing evidence. We have taken multiple measures to control for potential biases, including subgroup analysis, sensitivity analysis, and Egger’s test. However, this study has some limitations. First, our results could be more applicable to Whites, therefore most of the significant findings should be replicated in other populations with sufficient statistical power, such as the Asian populations. Second, the sample sizes in most of the candidate gene studies were small, especially in Asian populations. We observed lack of associations of almost all SNPs when summarizing the data from Asian cohorts. Therefore, larger cohorts are needed for further validation. Third, although we employed funnel plots and Egger’s tests to identify publication bias, there could still be remaining publication bias due to the reduced power when testing small number of studies in a meta-analysis. Moreover, the COL4A4 variants might not reach the genome-wide significance in the reported GWAS. The non-availability of the data for these variants could be a potential source of publication bias.

In conclusion, we have prioritized 8 SNPs in 6 genes/loci as significant genetic markers for keratoconus in Whites, including FOXO1 rs2721051, RXRA-COL5A1 rs1536482, FNDC3B rs4894535, IMMP2L rs757219 and rs214884, and BANP-ZNF469 rs9938149, and COL4A4 rs2229813 and rs2228557. We also identified 10 genes/loci with suggestive evidence of association with keratoconus. This study has thus provided an up-to-date list of candidate genetic markers for further investigations of their biological roles in keratoconus. More studies are warranted to confirm the reported genetic associations in different populations.

Methods

Searching methods for identifying studies

We searched for relevant records in the MEDLINE, Embase, Web of Science, and HuGENET databases via the Ovid platform. We used the Boolean logic to connect the free terms and controlled vocabularies (i.e. Medical Subject Heading terms): (“polymorphism(s)” OR “mutation” OR “genotype(s)” OR “genetic(s)” OR “gene(s)” OR “allele(s)” OR “DNA”) AND (“keratoconus”) (Supplementary Table 5). We also manually scanned the reference lists of the potentially eligible research articles, reviews and meta-analyses from the initial screening to ensure inclusion of all relevant publications. We did not use language filter in the literature search. The last search was performed on June 1, 2016.

Eligibility criteria

We set the following criteria for eligible studies for meta-analysis: (1) original case-control studies that evaluated the association of gene polymorphisms with keratoconus; (2) the study subjects were unrelated and recruited from explicitly defined populations; and (3) allele or genotype counts or frequencies in both case and control groups were reported or calculable; or odds ratio and 95% confidence intervals (CI) and/or standard error (SE) were reported. We excluded animal studies, case reports, reviews, abstracts, conference proceedings, and editorials.

Study selection, data collection and risk of bias assessment

Two reviewers (S.S.R. and S.T.U.M.) independently screened all the titles and abstracts of identified studies. Disagreements were resolved via discussions with a senior investigator (L.J.C.). After identifying potentially eligible articles, the 4 reviewers (S.S.R., S.T.U.M., X.T.Y., and L.M.) extracted the data separately and cross-validated the data. Consensus was reached via thorough discussion among all the reviewers. In this study, we used ‘Whites’ to represent individuals/populations whose ancestral origins are in the continent of Europe. We designed a customized datasheet for data extraction, which included the first author, year of publication, country of study, ethnicity, definition of case and control, sample sizes in the case and control groups, genes/loci, polymorphisms, allelic and genotypic counts and frequencies, ORs and 95% CIs or SEs of the polymorphisms and corresponding genetic models, and results of the test for HWE in the control group. First, we extracted all the polymorphisms and genes/loci reported in the potentially eligible studies searchable by the end of our search date. For GWAS, we extracted all the variants that were shown to be tested in replication cohorts in the result section and supplementary tables21,22,23,24. For candidate gene study, we extracted all the reported variants. No significance threshold for the genetic association has been applied during the data extraction. We also checked for potentially duplicated cohorts among the studies via comparing research groups and description of study populations. In the studies that had reported 2 or more independent cohorts, we extracted the data of each cohort separately. Second, we selected those polymorphisms that could be meta-analyzed. Third, the missing allele/genotype counts were calculated using the allele frequencies and sample sizes, assuming no deviation from HWE unless reported otherwise77. If only the OR and 95% CI were reported, we estimated the SE following the methods described in our previous papers77, 78. We attempted to contact the authors for additional information if necessary. If the HWE result was not reported, we tested it using the extracted data in the control group by the Chi-square test. Moreover, we used the NOS system (accessed via http://www.ohri.ca/programs/clinical_epidemiology/oxford.asp) to evaluate the quality of the case-control studies (Supplementary Appendix 1)79, 80. We assigned one star to a study if it met one requirement in the NOS from 3 dimensions (i.e., selection, comparability and exposure). The maximum number of stars that a study could obtain was 8. A study of <5 stars in overall or earned no star in any one of the items (i.e., case definition, ethnicity, or ascertainment of genotype) was considered as of suboptimal quality and having high risk in introducing bias81.

Data analysis

We conducted meta-analysis for the SNPs that had been reported in 2 or more study cohorts from at least 2 separated reports. The genetic association was assessed using the allelic (a vs. A) model, where “a” and “A” represented the associated and the reference alleles, respectively. We evaluated the inter-cohort heterogeneity using the I 282. An I 2 value of lower than 25%, between 25% and 50%, and greater than 50% indicated low, moderate, and high heterogeneity, respectively. However, to obtain more conservative results we calculated the summary OR and 95% CI for each SNP only using the random-effect model, in which the weighted effect of a SNP was estimated by inverse variance (IV) and τ2 from the DerSimonian-Laird estimator83, regardless of the Q statistics and I 2. Of note, to assess the replication results of SNPs identified in the GWAS on keratoconus23, 24, we first combined the data from both the GWAS and replication studies, and then removed the data from the initial GWAS. Subgroup analysis by ethnicity was then conducted in Whites and Asian populations (i.e., populations of Asian ancestries including 2 or more ethnic groups from Arabic, Chinese, Korean, Japanese, or Indian populations). We adopted the funnel plots and Egger’s test to assess potential biases (e.g. publication bias)84, 85. A P value of <0.05 in the Egger’s test indicated statistically significant bias. We also conducted the sensitivity analysis to confirm the associations by sequentially omitting each of the included studies one at a time and recalculated the summary outcomes. We then omitted the studies that deviated from HWE (PChi-squre ≤ 0.05), or of suboptimal quality. A finding is more likely to be true when the result is stable in the sensitivity analysis.

Customized analytical scripts based on the metafor package in the R software (v3.0.0, http://cran.r-project.org/) were generated for the meta-analysis.

As a strategy to account for multiple testing, we used the Bonferroni corrected alpha as the cut-off value for confirming the genetic associations. To calculate the adjusted alpha value, we divided 0.05 by the number of SNPs tested (N = 53) and also by the maximum number of different tests a SNP could be done (N = 7). The adjusted significant threshold was therefore 1.35 × 10−4. The P values > 1.35 × 10−4 and ≤ 0.05 were considered nominally significant. We consider a P value < 5 × 10−8 as genome-wide significance.