Genetic associations for keratoconus: a systematic review and meta-analysis

Genetic associations for keratoconus could be useful for understanding disease pathogenesis and discovering biomarkers for early detection of the disease. We conducted a systematic review and meta-analysis to summarize all reported genetic associations for the disease. We searched in the MEDLINE, Embase, Web of Science, and HuGENET databases for genetic studies of keratoconus published from 1950 to June 2016. The summary odds ratio and 95% confidence intervals of all polymorphisms were estimated using the random-effect model. Among 639 reports that were retrieved, 24 fulfilled required criteria as eligible studies for meta-analysis, involving a total of 53 polymorphisms in 28 genes/loci. Results of our meta-analysis lead to the prioritization of 8 single-nucleotide polymorphisms (SNPs) in 6 genes/loci for keratoconus in Whites. Of them 5 genes/loci were originally detected in genome-wide association studies, including FOXO1 (rs2721051, P = 5.6 × 10−11), RXRA-COL5A1 (rs1536482, P = 2.5 × 10−9), FNDC3B (rs4894535, P = 1.4 × 10−8), IMMP2L (rs757219, P = 6.1 × 10−7; rs214884, P = 2.3 × 10−5), and BANP-ZNF469 (rs9938149, P = 1.3 × 10−5). The gene COL4A4 (rs2229813, P = 1.3 × 10−12; rs2228557, P = 4.5 × 10−7) was identified in previous candidate gene studies. We also found SNPs in 10 genes/loci that had a summary P value < 0.05. Sensitivity analysis indicated that the results were robust. Replication studies and understanding the roles of these genes in keratoconus are warranted.

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  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 Chinese 32 or Arabs 40 . 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 studies 37,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/loci 22,24,50,51,55 . In 3 studies 22,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 studies 50, 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 cases 50 . 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 stress 57 . For example, Fox proteins regulate the expressions of anti-oxidative enzymes such as superoxide dismutase and thioredoxin reductase 58,59 . Moreover, reduced FOXO1 expression has been reported to induce apoptosis in human trabecular meshwork cells in response to oxidative stress 60 . 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 thinning 61 . We also found association of keratoconus with IMMP2L, a mitochondrial inner membrane protease. Mutation in IMMP2L also accumulates oxidative stress 62 . 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 differentiation 63 , 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 study 65 . 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 heterodimers 66 . Null mice of both RXRA and RXRA/RAR showed abnormal embryonic eye morphologies, including thickening of corneal stroma and absence of anterior chamber 66 . 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 gene 67 . 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 syndrome 68 ; COL4A3 and COL4A4 for Alport syndrome 69 ; HGF for non-syndromic hearing loss 70 ; IMMP2L for Gilles de la Tourette syndrome 71 ; KCND3 for spinocerebellar ataxia 72 ; LOX for thoracic aortic aneurysms and dissections 73 ; MPDZ for leber congenital amaurosis and retinitis pigmentosa 74 ; RAB3GAP1 for Warburg Micro syndrome and Martsolf syndrome 75 ; and ZNF469 for Brittle cornea syndrome 76 . 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 Scientific RepoRts | 7: 4620 | DOI:10.1038/s41598-017-04393-2 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 tables [21][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 otherwise 77 . If only the OR and 95% CI were reported, we estimated the SE following the methods described in our previous papers 77,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 bias 81 . 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 2 82 . 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 estimator 83 , regardless of the Q statistics and I 2 . Of note, to assess the replication results of SNPs identified in the GWAS on keratoconus 23, 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 (P Chi-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.