Rare, Potentially Pathogenic Variants in ZNF469 Are Not Enriched in Keratoconus in a Large Australian Cohort of European Descent

Abstract

Purpose: The Zinc Finger Protein 469 (ZNF469) gene has been proposed as a candidate gene for keratoconus due to the association of an upstream polymorphism (rs9938149) with the disease in two independent studies, and the role of the gene in the autosomal recessive disease Brittle Cornea Syndrome. Coding variants in ZNF469 have been assessed for association with keratoconus in several small studies, with conflicting results. We assessed rare, potentially pathogenic variants in ZNF469 for enrichment in keratoconus patients in a cohort larger than all previous studies combined.

Methods: ZNF469 was sequenced in 385 Australian keratoconus patients of European descent, 346 population controls, and 230 ethnically matched screened controls by either whole exome sequencing or targeted gene sequencing. The frequency of rare and very rare potentially pathogenic variants was compared between cases and controls using χ² or Fisher's exact tests and further explored using a gene based test (Sequence Kernel Association Test [SKAT]), weighting on the rarity of variants.

Results: A total of 49 rare, including 33 very rare, potentially pathogenic variants were identified across all groups. No enrichment of rare or very rare potentially pathogenic variants in ZNF469 was observed in our cases compared to the control groups following analysis using χ² or Fisher's exact tests. This finding was further supported by the SKAT results, which found no significant difference in the frequency of variants predicted to be damaging between cases and either control group (P = 0.06).

Conclusions: Rare variants in ZNF469 do not contribute to keratoconus susceptibility and do not account for the association at rs9938149.

Keratoconus (OMIM 148300) is a complex disease characterized by the progressive thinning and protrusion of the cornea. This causes high myopia and irregular astigmatism and can severely affect vision and quality of life.¹ Standard treatments such as hard contact lenses are difficult to manage, are often painful for the patient, and corneal transplantation is ultimately required in a large proportion of patients.^2,3 Keratoconus has an incidence of 1 in 50,000^2,4 and usually develops between puberty and early adulthood. The prevalence is reported to be between 17 and 229 per 100,000, and varies greatly between studies, ethnicities, and geographical locations.^2,5–10 In Caucasians, the prevalence is estimated at 54.5 per 100,000.² There is strong evidence for the role of genetic risk factors in keratoconus susceptibility; however, largely the functional variants and genes involved in disease development and progression have not been elucidated.^10–12 Identifying genetic factors underlying keratoconus susceptibility would improve our understanding of the disease pathogenesis and aid the development of novel, nonsurgical treatments, which would ultimately improve the quality of life for patients. 

As keratoconus progresses, the cornea thins and bulges forward forming a cone shape. Therefore, it has been hypothesized that genes involved in central corneal thickness (CCT) may play a role in the etiology of keratoconus. Under this hypothesis, a genome-wide association study (GWAS) by Lu et al.¹³ for CCT assessed genome-wide significant single nucleotide polymorphisms (SNPs) in a large cohort of keratoconus patients.¹³ A SNP upstream of the Zinc Finger Protein 469 (ZNF469) gene, rs9938149, showed a suggestive association with keratoconus; however, the genotype associated with a thinner cornea was associated with decreased keratoconus risk.¹³ This finding was replicated in an independent cohort showing the same direction of association,¹⁴ indicating that the association is likely to be real, if nonintuitive. As ZNF469 is the closest gene to this SNP, it has been hypothesized that genetic variation within the gene may account for the association at rs9938149 as well as contribute to CCT and keratoconus susceptibility. 

The potential role of ZNF469 in keratoconus pathogenesis is further supported by the role of this gene in Brittle Cornea Syndrome type 1 (BCS1, OMIM 229200). BCS1 is a rare, autosomal recessive connective tissue disorder, caused by biallelic loss-of-function variants in ZNF469. A key feature of this syndrome is extremely thin corneas that are prone to spontaneous rupture. This suggests that ZNF469 is important for the structural integrity of the cornea. 

As an appealing candidate, coding variants in ZNF469 in keratoconus have since been investigated. Two of these studies reported an association of potentially pathogenic variants in ZNF469 with keratoconus, while two showed no association with disease.^15–18 Lechner et al.¹⁷ assessed 112 cases and 96 local controls in a combined United Kingdom and Swiss cohort while Vincent et al.¹⁸ examined a largely Polynesian cohort involving 43 cases and 46 controls from New Zealand. Both of these studies concluded that rare potentially pathogenic variants in ZNF469 were enriched in keratoconus. Conversely, work by Karolak and colleagues¹⁶ did not identify any difference in the frequency of nonsynonymous variants between 42 cases, 49 individuals with high myopia, and 268 controls without eye disease in a Polish cohort. Furthermore, a study by Davidson et al.¹⁵ showed that uncommon and rare potentially pathogenic variants in 11 families with keratoconus (six of Middle Eastern heritage and five of mixed ethnicities) did not segregate with disease. This study also examined four heterozygous carriers of nonsense or frame-shift variants in ZNF469, and found no evidence of keratoconus.¹⁵ 

The potential involvement of ZNF469 in keratoconus remains unclear and further assessment of this candidate gene is required. It is likely that the contrasting findings from previous studies are due to their small sample sizes. To address this, our study examined the frequency of rare potentially pathogenic variants in ZNF469 in 385 keratoconus cases compared to ethnically matched controls. By assessing more cases than all previous studies combined, the present study allows for a more comprehensive analysis of the potential contribution of rare variants in ZNF469 to keratoconus risk. 

Methods

All investigations adhered to the principles of the Declaration of Helsinki and were approved by the Southern Adelaide Clinical Human Research Committee and the Human Research Ethics committee Tasmania. All participants gave written informed consent. 

Case Cohort

Australian keratoconus patients of European descent were referred by treating optometrists and ophthalmologists and recruited through the Flinders Eye Clinic (Adelaide, Australia) and optometry clinics in Adelaide and Melbourne. Additionally, members of Keratoconus Australia were recruited from across Australia via mail. Clinical examinations were performed by an experienced ophthalmologist. Individuals were diagnosed with keratoconus if they had videokeratographic features of keratoconus or any of the following clinical signs: conical corneal protrusion, central or paracentral stromal thinning or other distinctive features such as Fleischer's ring, Vogt's striae, epithelial or subepithelial scarring, or oil droplet sign and/or scissoring of the retinoscopic reflex. In addition, individuals with a history of corneal transplantation for keratoconus were also classified as cases. Peripheral blood was donated by patients and DNA was extracted using the QiaAmp DNA Maxi Kit (Qiagen, Hilden, Germany). 

Control Cohorts

Two Australian control cohorts of European descent were used in this study. The first, the “population controls”, were an ethnically matched cohort from the Anglo-Australasian Osteoporosis Genetics Consortium (n = 346). These individuals are all female with moderately high or low bone mineral density measurements (1.5<|BMD|<4.0), but were not examined for eye disease. The second control group, the “screened controls”, were a cohort of Australians of European descent involved in other genetic studies who were examined for eye disease and found to be unaffected by experienced ophthalmologists at the Flinders Eye Clinic (n = 230). While these individuals had no clinical evidence of keratoconus, the majority of these individuals had advanced glaucoma (n = 195). The remaining individuals either had no evidence of eye disease (n = 22) or were unaffected individuals from families with congenital cataract (n = 9) and nanophthalmos (n = 4). 

Whole Exome Sequencing

Whole exome sequencing (WES) data were available for 99 keratoconus cases, the screened control cohort (n = 230) and the population control cohort (n = 346). WES for the keratoconus patients and the screened control cohort was conducted by Macrogen, Inc., using the SureSelect Human All Exon V4 enrichment kit (Agilent Technologies, Inc., Santa Clara, CA, USA) with paired end sequencing on a HiSeq 2000 (Illumina, San Diego, CA, USA). Raw reads were aligned to the hg19 reference genome with BWA-MEM using the Churchill Pipeline¹⁹ and variants were joint-called with SAMtools and BCFtools (versions 1.3.1).²⁰ The population control cohort was sequenced at the University of Queensland Centre for Clinical Genomics using the TruSeq Exome Enrichment (Illumina) on an Illumina HiSeq 2000. Raw reads were aligned to the reference genome (hg19) using novoalign (version .02.08), and GATK²¹ (version 3.2-2) was used for variant calling and quality score calibration according to GATK's “Best Practices Guidelines”.^22,23 To assess coverage across ZNF469 in the WES data (cases and the two control groups), the mean read depth for each variant position was determined using VCFtools²⁴ and graphed using custom R²⁵ scripts using the R packages ggplot2²⁶ and cowplot.²⁷ For each cohort, regions with a mean depth of <10 were excluded from analyses. These regions were different for the two control cohorts. Variants were only included in analyses if a sequencing depth of ≥10 and a quality score ≥20 was obtained. 

Targeted Gene Screen in Additional Cases

ZNF469 was screened in a total 341 cases as part of a candidate gene screen using the HaloPlex Target Enrichment System (Agilent Technologies, Inc.) with a custom designed probe panel on pooled DNA samples (Supplementary Table S1). A brief summary of the design is presented in Supplementary Figure S1. DNA pools containing equimolar DNA samples from eight keratoconus patients were prepared as published previously.²⁸ Each DNA pool was indexed with a unique indexing primer cassette, allowing for multiplexed sequencing. Sequencing was conducted in batches of 11 DNA pools on the MiSeq platform (Illumina) using a MiSeq V2 Reagent kit (300 cycles) with paired-end reads. SureCall (Agilent Technologies, Inc.) was used for sequencing analysis using standard trimmer parameters, the BWA-MEM algorithm to align reads and the SNPPET SNP Caller (part of SureCall) to call variants. Variants were called if a minimum read depth of 10 and quality score of 20 was reached. As variants were called from pooled DNA samples, it was expected that if a single alternate allele was present it would be observed on 6.25% of the reads mapping to that position. To account for this, the minimum allele frequency for heterozygous variants was set to 0.035 and the threshold for indel calling was ≥ 0.04. Suspected artefacts and real variants were selected for validation by direct sequencing in all individuals included in the DNA pool, to identify thresholds for the inclusion of variants. In addition, 55 cases included in this candidate gene screen were included in the WES; providing additional cross-validation between the two sequencing strategies. However, variants from these samples were only counted in the analyses once. An in-depth description of the validation of variant calls from the pooled gene screen and the development of thresholds for variant inclusion is available in the Supplementary Materials. 

Variant Validation by Direct Sequencing

To validate variants called in the case data set, primers were designed using Primer3Plus.²⁹ DNA samples were amplified with either GoTaq (Promega, Madison, WI, USA), HotStar Taq Kit (Qiagen), or MyTaq HS Mix (Bioline, London, UK) and purified using Agencourt AMPure XP (Beckman Coulter, Brea, CA, USA) according to the manufacturer's instructions. Purified amplicons were sequenced using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) on either an ABI 310 Genetic Analyzer or an ABI 3500 (Applied Biosystems). DNA sequences were aligned to the human reference genome (hg19), and the chromatograms were manually inspected at the position of each variant using Sequencher 4.10.1 (Gene Codes Corporation, Ann Arbor, MI, USA). 

Variant Filtering Strategies

While it has recently been suggested that ZNF469 is a single-exon gene,³⁰ the capture methods in this study were all designed under the assumption that ZNF469 has two exons with a short 84bp intron. Therefore, only exonic variants were included in analyses. To ensure high quality sequence data and robust comparisons, only regions with sufficient coverage across the relevant data sets were included in each comparison. ANNOVAR³¹ was used to annotate variants with the minor allele frequency (MAF) of the non-Finnish European population in the Exome Aggregation Consortium database³² (ExAC NFE), as well as the predicted pathogenicity of the variants using Sorting Tolerant from Intolerant (SIFT)³³ and the HumDIV algorithm from Polymorphism Phenotyping v2 (PolyPhen2).³⁴ To make our results comparable to previous studies, two filtering strategies were used to define rare and very rare potentially pathogenic variants. Filtering strategy 1 included nonsynonymous variants with a MAF < 0.01 (rare) with a damaging or probably/possibly damaging prediction by SIFT and/or PolyPhen2. As SIFT and PolyPhen2 do not assess insertions/deletions (indels), all indels with a MAF <1% were included in filtering strategy 1. Filtering strategy 2 was identical to filtering strategy 1, except variants were required to have a MAF < 0.001 (very rare). 

Statistical Analyses

For both filtering strategies, the number of alternate alleles was compared between cases and the two control data sets separately using χ² or Fisher's exact tests, where appropriate. The odds ratio (OR) and 95% confidence interval (95% CI) were calculated. In addition, the Sequence Kernel Association Test (SKAT)³⁵ was used to assess if potentially pathogenic variants were enriched in the case WES data compared to the two control cohorts. Using the two control cohorts as a separate comparison, SKATBinary was run using the quantile adjusted moment matching (QA) method and the default weight parameter, Beta (1,25). This weighting applies strong weights to rare variants, nonzero weighting to uncommon variants (MAF 0.01–0.05), and almost zero weights to common variants.³⁵ 

Genetic Power Calculations

Power calculations were conducted using the case-control for discrete traits module of the genetic power calculator (http://zzz.bwh.harvard.edu/gpc/).³⁶ These calculations assumed an additive module where the prevalence of keratoconus was 0.00067 (1 in 1500) and D' prime was equal to 1. The high-risk allele frequency was set for each filtering strategy separately, using the frequency of variants identified in the screened controls. 

Variant Visualization

Variants included in filtering strategy 1 were plotted as a bar plot for each group using custom R²⁵ scripts and the R package ggplot2.²⁶ Each variant was plotted along the x-axis according to its genomic position, and the frequency of the variant in the study group was indicated by the height of the bar. Indels were plotted according to the position of the first affected base. To indicate the conservation of each variant position, the bars were colored using a color gradient corresponding to the 100-way vertebrate PhastCons^37,38 score as available from the University of California Santa Cruz (UCSC) genome browser.³⁹ The plots for each group were aligned to a schematic of ZNF469 using the R package cowplot,²⁷ to allow for comparison of the variants across the gene. The position of the zinc finger motifs was plotted onto the gene schematic based on the positions obtained from UniProt⁴⁰ (entry Q96JG9). 

Results

A total of 385 cases, 346 population controls, and 230 screened controls were included in analyses. Demographic details are provided in Table 1. 

Table 1

 

Demographics of the Keratoconus Cases and Control Groups at the Time of Examination

 

As we applied a stringent coverage filter to all sequence data sets, regions containing the majority of the first exon, the beginning of exon two, and the regions encoding the zinc finger domains were excluded from analyses that compared the cases to the population control cohort (Fig. 1). Therefore, comparisons using the population controls included a total of 6700 bases across ZNF469 with high quality sequence data, which corresponds to 56.5% of the coding regions. In contrast, a total of 11,448 bases (96.5% of the coding regions) obtained sufficient coverage across the gene in both the cases and screened controls. For comparisons using the screened controls, only two small regions (85 bp and 279 bp) either side of the intron were excluded from analysis due to poor coverage (Fig. 1). By removing these regions from analysis, the single region that did not meet the minimum depth threshold in the pooled gene screen data was also excluded. A summary of the coverage metrics for the gene screen data set is presented in Supplementary Table S2. 

Figure 1

 

Sequencing coverage in the WES data sets. (A) Coverage across ZNF469 in each data set, based on the mean read depth at variant positions. The horizontal dashed line indicates the minimum depth threshold accepted (10 reads). Regions where the mean read depth in the population controls was below this threshold are shaded light gray and were excluded from analyses in which these controls were included. Regions with insufficient coverage in all cohorts are shaded darker gray. (B) A schematic of ZNF469 where dark gray boxes indicate the position of the zinc finger motifs and the small intron is indicated by the horizontal line.

Across all cohorts, 49 variants (46 single nucleotide variants and three small deletions) fulfilled the criteria for rare potentially pathogenic variants in filtering strategy 1 (Table 2). Of these, 19 variants were observed in cases only, while 17 were unique to the control cohorts and the remaining 13 were in both cases and controls. Eleven variants had been observed in previous studies in either cases or controls, with an additional two variants that are not identical, but occur at the same amino acid as a previously reported variant. A total of 33 rare potentially pathogenic variants were observed 73 times in cases and 23 variants were called 48 times in the screened controls, which was not significantly different (P = 0.66; Table 3). When considering the comparison between the cases and the population controls, 17 variants were identified 31 times in cases and 10 variants were observed 22 times in the controls, which was not significantly different between groups (P = 0.47; Table 3). 

Table 2

 

Variants Included in Analysis Under Filtering Strategy 1 (Rare Potentially Pathogenic Variants) Including the Position, Nucleotide and Protein Changes, Population Frequency, the 100-Way Vertebrate PhastCons Score, Cohort Frequency (Freq) and Allele Count (AC)

 

Table 3

 

Association Analyses Using χ² or Fisher's Exact Test Under Each Filtering Strategy

 

Filtering strategy 2 identified 33 very rare potentially pathogenic variants across all cohorts (Table 2). Sixteen of these variants were unique to keratoconus patients and two, p.(C1693F) and p.(P2548L), were identified in two cases each. Fifteen variants were observed in the control cohorts only, including p.(P3372L), which was identified in two individuals. Two variants, p.(P626_G628del) and p.(E3781K), were identified in both a case and a control. All other variants were observed in a single individual. Two variants, p.(S2242Y) and p.(P3372L), were previously observed in the study by Lechner et al.¹⁷ and p.(E935G) is located at the same amino acid as a variant identified in the study by Davidson and colleagues.¹⁵ In total, 18 variants were identified 20 times in cases and 12 variants were observed 13 times in the screened controls. For the population control comparison, 10 variants were observed in 12 cases and five variants were identified in the controls. Similar to the results for rare variants, very rare variants were not enriched in cases compared to either the screened controls (P = 0.96) or the population controls (P = 0.15; Table 3). Furthermore, the SKAT analyses demonstrated no significant enrichment of variants predicted to be damaging (P = 0.06; Table 4). 

Table 4

 

The Results of the SKAT Analyses

 

For the power calculations, the high-risk allele frequency was assumed to be the frequency of variants identified by the filtering strategies in the 230 screened controls. This was determined to be 0.104 (48/460) for strategy 1 and 0.028 (13/460) for strategy 2. Using these values, our study had 80% power to detect a relative risk of 1.5 and 2.0, respectively. 

For each cohort, the allele frequencies were plotted against the variant position and mapped to a schematic of the gene (Fig. 2). The first exon showed a similar pattern of variation between the cases and screened controls (Fig. 2). Due to insufficient coverage of this region in the population controls, variation in this group could not be assessed. Of the 13 variants in cases in exon 1, two variants, p.(S242I) and p.(R503W), were located at highly conserved nucleotides (PhastCons score >0.99) and p.(E679K) was at a relatively conserved position with a PhastCons score of 0.655. The p.(R503W) variant was observed in a single case (frequency of 0.0013), whereas both p.(S242I) and p.(E679K) were identified in one case and one screened control (frequency of 0.0022). 

Figure 2

 

Summary of identified variants. (A) Schematic of ZNF469 with zinc finger motifs shaded dark gray, (B–D) bar plots indicating the position and alternate allele frequency (AAF) of rare potentially pathogenic variants in 784 case alleles, 460 screened control alleles, and 692 population control alleles, respectively. The bars are colored according to the 100-way vertebrate PhastCons Score for the corresponding position, where dark blue is a score of 0 and red is a score of 1. Light gray shading on the plots for the control cohorts indicate regions with insufficient coverage for variant calling that were excluded from analysis when using these data. The gene schematic and all graphs are aligned vertically to share the same x-axis to allow for comparison.

Few variants were located in the proximal region of exon 2, with only one variant identified in cases. In contrast, a cluster of variants was observed in cases in the distal half of exon 2 (Fig. 2). Specifically, variants were observed in the region spanning just proximal of the first zinc finger motif to midway between the second and third zinc finger motif. Noticeably fewer variants were observed in the control cohorts within this region (Fig. 2). Three variants identified in cases within this cluster, p.(D2902Y), p.(E3026A) and p.(E3137G), were highly conserved (PhastCons scores > 0.96). The p.(D2902Y) variant was present at a similar frequency in cases (0.0038) and the population controls (0.0029), while p.(E3026A) and p.(E3137G) were present in one case and absent in controls. Two additional variants, p.(E3781K) and p.(E3898K), were identified at highly conserved residues at the very distal end of exon two (PhastCons score = 1). These variants were both identified a single case and neither of these variants were observed in the population controls despite good coverage in this region; however, the p.(E3781K) variant was observed in a one screened control. 

Across ZNF469, three variants—p.(A566V), p.(P665L), and p.(R3426Q)—were identified at a frequency >0.01 in both the cases and the screened control cohort, despite a MAF < 0.01 in the public ExAC NFE database. The frequency of these variants could not be assessed in the population controls due to insufficient coverage; however, all three variants were observed in cases at a similar or lower frequency than the screened controls. 

Discussion

In our cohort of Australians of European descent, the overall frequency of rare potentially pathogenic variants in ZNF469 was not different between keratoconus cases and two independent control cohorts using two filtering strategies. Despite a much larger cohort than the earlier studies, our findings were unable to replicate previous reports that show an association of potentially pathogenic variants in ZNF469 and keratoconus development with relative risks of up to 12.^17,18 Our study had 80% power to detect a relative risk of 1.5 for rare potentially pathogenic variants, and 2.0 for very rare potentially pathogenic variants, commensurate with expected effect sizes for rare variants in a complex disease. In our study, we show that the reported large effect sizes do not exist when the whole gene is considered in both cases and controls and therefore if this gene is contributing to keratoconus risk, then effect sizes are likely to be quite small. The main finding of this paper is that the reported large effects of rare variants are likely spurious and are brought about by biased reporting of variants detected in fully sequenced cases not being present in controls who are not fully sequenced or due to the exclusion of variants that are observed in control individuals. Our work shows that controls have just as many variants meeting the potentially pathogenic criteria as keratoconus patients. 

Due to the size of our study, we were also able to map the location of the rare potentially pathogenic variants identified in all three cohorts and show that they span the whole gene, with particular aggregation in the first exon and the distal half of the second exon. Only eight variants observed in cases were located at highly conserved nucleotides. Four of these variants were identified at similar frequencies in both cases and controls and four were only observed in a single case. While it is possible that these rare variants may contribute to keratoconus susceptibility in these few cases, on the whole, the evidence indicates that rare, potentially pathogenic variants in ZNF469 do not make a substantial contribution to keratoconus risk. This finding is consistent with the work by Davidson and colleagues,¹⁵ which showed that uncommon variants (MAF <0.025) did not segregate with disease in families with keratoconus and therefore, at least in isolation, do not contribute to keratoconus susceptibility. Furthermore, the results of the Polish study¹⁶ indicate that potentially pathogenic variants are not enriched in keratoconus, and that ZNF469 is highly allelic in the general population. 

The first reports to assess coding variants in ZNF469 were published by Lechner et al.¹⁷ and Vincent et al.¹⁸ Subsequent studies, including this one, have used similar filtering strategies to allow for direct comparison. Our filtering strategy 1 was based on the criteria used by Vincent et al.,¹⁸ while filtering strategy 2 was based on the method used by Lechner et al.¹⁷ with two key changes. Firstly, Lechner et al.¹⁷ removed any variants from analysis that were present in both their cases and controls. We did not do this as keratoconus is a complex disease and therefore it is likely (and expected) that unaffected individuals will carry risk associated variants without ever developing disease. Secondly, Lechner et al.¹⁷ only used SIFT predictions to classify variants as potentially pathogenic. SIFT uses protein conservation to calculate pathogenicity by comparing a query sequence to similar sequences with similar function.³³ As ZNF469 is a highly variable gene with low conservation in lower mammals and vertbrates,⁴¹ SIFT may misclassify deleterious variants in regions of low conservation. In contrast, PolyPhen2 uses the properties of the substituted amino acids and the proximity to functional domains or structural features, as well as protein conservation to predict pathogenicity.³⁴ As the structure and function of the ZNF469 protein remains largely unknown,⁴¹ PolyPhen2 is likely to better assess regions with poor conservation, particularly the regions that flank the zinc finger domains. Therefore, our study used both SIFT and PolyPhen2 to better capture the pathogenicity of nonsynonymous variants identified in ZNF469. Our study used these robust and complementary methods to replicate the analytical strategies of the previous studies; however, our findings do not support any enrichment of rare potentially pathogenic variants in ZNF469 in keratoconus cases. 

We assessed rare potentially pathogenic variants in ZNF469 in the largest cohort of keratoconus patients to date by combining WES data and sequencing data from a targeted gene screen using pooled DNA samples. Following extensive validation experiments, we demonstrated a high level of consistency of variant calls for individuals sequenced by both methods, validating the utility of pooling DNA samples to maximize cost-effectiveness. Furthermore, these validation experiments were used to develop stringent thresholds that were then applied to all data to ensure only high quality variants were included in analyses. While variants were not assessed for validation in the control cohorts, the inclusion of additional control variants that were filtered out due to these stringent thresholds would only strengthen our findings of no association. Furthermore, to minimize bias due to poor coverage of ZNF469 in the population controls, a strict coverage threshold was applied such that only regions with sufficient coverage were assessed when comparing to this control group. While this may have resulted in the exclusion of important variants, the vast majority of these regions had sufficient coverage in the screened controls and therefore were still assessed. 

As suggested by Davidson et al.,¹⁵ it is likely that variation in ZNF469 is underrepresented in public databases as a result of the poor coverage of the gene by the older WES capture techniques. According to the ExAC Browser, the mean coverage for ZNF469 (ENSG00000225614) is 7.3 reads and the proportion of individuals with at least 10× coverage is less than 20%. Our experience with the poor coverage of ZNF469 in the population control data set, as well as the occurrence of three variants that were identified in our cases and screened controls at a frequency of > 0.01, despite being annotated with a MAF < 0.01 in the ExAC NFE database, supports this. As one might expect, these variants could not be assessed in the population controls due to insufficient coverage. All three of these variants were located at relatively nonconserved nucleotides (PhastCons scores <0.4). Two of these variants, p.(P665L) and p.(R3426Q), were observed at similar frequencies in both our cases and controls and therefore were hypothesized to be benign polymorphisms. These variants were similarly classified in the study by Lechner and colleagues, but were not identified in other reports.^15–18 The third variant, p.(A566V), was identified at more than twice the frequency in the screened controls (0.028) than the cases (0.012). Vincent et al.¹⁸ reported this variant in one Indian and two Caucasian keratoconus cases, while Davidson et al.¹⁵ identified the variant in two cases and two unaffected individuals from two separate consanguineous families of Middle Eastern origin. In addition, the work by Lechner et al.¹⁷ excluded the variant from analysis due to a MAF > 0.01 in their control cohort. As our screened control cohort largely consists of advanced glaucoma patients, it is possible that these variants may be involved in glaucoma susceptibility. However, in the original GWAS that identified the association between the SNP upstream of ZNF469 (rs9938149) and CCT as well as keratoconus, glaucoma cases were also assessed and no association was identified.¹³ Therefore, it is more likely that these variants are benign, uncommon polymorphisms. 

Conclusions

Our study included more keratoconus cases than the total number previously studied and demonstrated no significant difference in the overall frequency of potentially pathogenic variants in ZNF469 compared to controls. In addition, our findings confirm the highly variable nature of ZNF469 and the poor capture of the gene by previous generation WES capture methods, which is likely to have resulted in underestimates of alternate allele frequencies in public databases. Overall, this study indicates that rare potentially pathogenic variants are unlikely to contribute to keratoconus pathogenesis and do not account for the association at the distant upstream SNP, rs9938149. 

Acknowledgments

The authors thank the staff and participants of the Anglo-Australasian Osteoporosis Genetics Consortium, Australian and New Zealand Registry of Advanced Glaucoma and other genetic eye disease studies at Flinders University. 

Supported by the National Health and Medical Research Council (NHMRC) of Australia Centre for Research Excellence (GNT1023911) and Project grant (GNT1104700). KPB is supported by a Senior Research Fellowship from the NHMRC and JEC by an NHMRC Practitioner Fellowship. SEML is supported by the Australian Government Research Training Program Scholarship and the Pennicott Foundation. 

Disclosure: S.E.M. Lucas, None; T. Zhou, None; N.B. Blackburn, None; R.A. Mills, None; J. Ellis, None; P. Leo, None; E. Souzeau, None; B. Ridge, None; J.C. Charlesworth, None; M.A. Brown, None; R. Lindsay, None; J.E. Craig, None; K.P. Burdon, None