Rare variants in optic disc area gene CARD10 enriched in primary open‐angle glaucoma

Abstract Background Genome‐wide association studies (GWAS) have identified association of common alleles with primary open‐angle glaucoma (POAG) and its quantitative endophenotypes near numerous genes. This study aims to determine whether rare pathogenic variants in these disease‐associated genes contribute to POAG. Methods Participants fulfilled strict inclusion criteria of advanced POAG at a young age of diagnosis. Myocilin mutation carriers were excluded using direct sequencing. Whole exome sequencing was performed on 187 glaucoma cases and 103 local screened nonglaucoma controls then joint‐called with exomes of 993 previously sequenced Australian controls. GWAS‐associated genes were assessed for enrichment of rare predicted pathogenic variants in POAG. Significantly enriched genes were compared against Exome Aggregation Consortium (ExAC) public control. Results Eighty‐six GWAS disease or trait‐associated glaucoma genes were captured and sequenced. CARD10 showed enrichment after Bonferroni correction for rare variants in glaucoma cases (OR = 13.2, P = 6.94 × 10−5) with mutations identified in 4.28% of our POAG cohort compared to 0.27% in controls. CARD10 was significantly associated with optic disc parameters in previous GWAS. The whole GWAS gene set showed no enrichment in POAG overall (OR = 1.12, P = 0.51). Conclusion We report here an enrichment of rare predicted pathogenic coding variants within a GWAS‐associated locus in POAG (CARD10). These findings indicate that both common and rare pathogenic coding variants in CARD10 may contribute to POAG pathogenesis.


Introduction
Globally, glaucoma is the most common cause of irreversible and preventable blindness and the second most common cause of all blindness after cataract, accounting for 8% of blindness worldwide in 2010 (Pascolini and Mariotti 2012). The defining feature of glaucoma is an optic neuropathy characterized by optic disc cupping with corresponding visual field loss (Foster et al. 2002). With early identification and intervention, mainly with intraocular pressure (IOP)-lowering therapy, vision loss is usually preventable (Leske et al. 2003). The disease incorporates a spectrum of subtypes ranging from congenital to age-related, normal of raised IOP, open or closed iridocorneal angle and primary or secondary etiology (Kwon et al. 2009). Primary open-angle glaucoma (POAG) is the most common disease subset and the high IOP phenotype is predominant in Caucasian and African ethnicities (Tham et al. 2014). While elevated IOP is the single greatest risk factor for POAG, other associated endophenotypes include optic disc morphology and central corneal thickness (CCT) (Miglior et al. 2007;Jiang et al. 2012;Hollands et al. 2013).
The exact contribution of rare variants to glaucoma is not known. Common tagged single-nucleotide polymorphism (SNP) analysis at the genome level, using SNP microarray technology cannot ascertain the burden of rare variants due to the poor linkage disequilibrium between common and rare variants (Siu et al. 2011). This casecontrol study utilized whole exome sequencing (WES) to investigate the degree of enrichment and proportion of disease burden accounted for by rare pathogenic variants in genes near common variants known to be associated with POAG and related quantitative traits from GWAS as this has not been adequately addressed in the literature.

Materials and Methods
This study adhered to the principles listed by the revised Declaration of Helsinki and the Australian National Health and Medical Research Council (NHMRC) statement of ethical conduct in research involving humans. Ethical approval was obtained from the Southern Adelaide and Flinders University Clinical Research Ethics Committee. All peripheral blood samples were collected for genomic DNA extraction as a part of the Australian and New Zealand Registry of Advanced Glaucoma (ANZRAG) as previously described (Souzeau et al. 2012).

Participants
Unrelated Caucasian participants with advanced glaucoma and nonglaucoma Caucasian controls were included for the study. All DNA samples were extracted from peripheral blood samples, using the QIAamp â DNA blood kit (Qiagen, Hilden, Germany) following the manufacturer's protocol. Our case cohort inclusion criteria selected only participants with the most severe disease and the youngest age (mean = 44.4 years, SD = 10.4 years) of onset within the ANZRAG database. For this study, participants were included if there was glaucomatous visual field loss involving at least two of the four central fixation squares with a pattern standard deviation of <0.5% on a reliable Humphrey 24-2 field (Carl Zeiss, Dublin, CA), or a mean deviation of worse than À15 dB in the worst affected eye. For participants without formal visual field testing, their bestcorrected visual acuity had to be worse than 20/200 due to glaucomatous damage. Furthermore glaucomatous damage had to also be present in the less affected eye as measured by a reliable Humphrey 24-2 field, with corresponding neuroretinal rim thinning. With the ANZRAG database participants with the youngest age of diagnosis were selected for this study. Efforts were made to ensure a comparable split of participants with high and normal IOP. High-tension glaucoma (HTG) was defined as having a maximum recorded IOP of >21 mmHg. Individuals with known MYOC mutations were excluded from this study as the genetic cause for their POAG has already been identified.
For the local control cohort, all participants were examined to exclude glaucoma or glaucoma-related phenotypes. Furthermore, in silico analysis included a larger unexamined control cohort from the Australian Osteoporosis Genetics Consortium (Estrada et al. 2012) (AOGC). These controls were female participants with high or low bone mass who were otherwise healthy by self-report.

Data acquisition
Exon capture and enrichment was performed on all glaucoma cases and local controls with the SureSelect Human All Exon V4 enrichment kit (Agilent, Santa Clara, CA, USA) as per manufacturer's protocol and enriched DNA sequenced on a HiSeq2000 (Illumina, San Diego, CA, USA) with paired end 100 bp reads (Macrogen Next Generation Sequencing Services). Local controls served as both technical and phenotypic controls. Raw experimental data were cleaned and joint-called with AOGC exome data that were captured with Nimblegen Human Exome Capture V2 (Roche, Basel, Switzerland) and sequenced on the HiSeq2000 (Illumina, San Diego, CA, USA) (University of Queensland Centre for Clinical Genomics). FASTQ files were aligned to human genome build hg19 with novoalign (version 3.02.08). Duplicate reads were marked with Picard MarkDuplicates (version 1.124). Local indel realignment and base quality recalibration were per-

Data analysis
Using in-house UNIX scripts, only protein coding exonic and splicing site variants were selected for analysis. Genes were included in the list if they were reported by any published GWAS or meta-analysis paper to be the closest to SNPs associated with POAG or disease endophenotype at P < 5 9 10 À8 in any ethnicity. Quality control threshold was set at a Genotype Quality (GQ) score of 20. Focusing on the impact of rare mutations, all variants with known minor allele frequency (MAF) of >1% in dbSNP, 1000 genomes and ExAC were excluded. Pathogenicity filtering removed all synonymous variants and missense variants considered tolerated or only possibly damaging by both SIFT and PolyPhen-2 HVAR, and retained stop site, frameshift and predicted deleterious missense mutations. Raw VCF data from the Exome Aggregation Consortium were annotated with the ANNO-VAR pipeline followed by variant selection using the same in-house scripts to target only rare deleterious mutations in the non-Finnish European subgroup to act as a larger in silico population control. Mutation loads per gene were calculated for the glaucoma cases, local control, AOGC control and ExAC control cohorts by summing the minor allele counts of all qualifying variants in the same gene and dividing by the average number of captured alleles for those variants, thereby adjusting for capture rate. Odds ratios were generated between cases and each control cohort separately, and Fisher's exact test used to calculate P-values, and Bonferroni correction applied for multiple testing of all analyzed genes (n = 86, threshold P = 5.81 9 10 À4 ). Five sets of randomly chosen independent groups of 86 WES-captured genes were selected using a random number generator and analyzed using the same protocol to act as control gene sets. This step further strengthens the positive results by examining the likelihood of false-positive findings.

Validation of variants
All variants within significantly associated genes were independently validated by capillary sequencing (primers listed in Table S1). PCRs were conducted on carrier samples using 40 ng of DNA, 10 lmol of each corresponding forward and reverse primers, 0.5 units of HotStarTaq DNA polymerase (Qiagen, Hilden, Germany), 2 lmol of dNTP and PCR buffer in a 20 lL volume per reaction. Thirty cycles were performed at an annealing temperature of 62°C. Five lL of each PCR product was incubated with 2 lL of Shrimp Alkaline Phosphatase (Affymetrix, Santa Clara, CA, USA) and 0.5 lL of E. coli exonuclease I (New England BioLabs, Ipswich, MA, USA) at 37°C for 60 min. Sequencing was carried out on the fluorescence-based capillary electrophoresis system 3130xl Genetic Analyzer (Life Technologies, Carlsbad, CA, USA) using BigDye Terminator V3.1 (Life Technologies, Carlsbad, CA, USA) according to manufacturer's protocol.

Results
A total of 187 unrelated participants with advanced POAG and 103 screened nonglaucoma controls were sequenced with whole exome capture and enrichment. Combined with the large unscreened cohort of AOGC controls (n = 993), the total number of controls was 1096. More participants had high-tension glaucoma (HTG, n = 122) than normal-tension glaucoma (NTG, n = 65). To date, 101 genes have been reported by GWAS studies as being near SNPs that are statistically associated with POAG, central corneal thickness, IOP or optic disc morphology. Of these 101 genes, 86 had qualifying (MAF <0.01 and predicted pathogenic) variants captured in glaucoma and control cohorts in this study (Table S2).
The sequencing strategy utilized for cases and local controls in this study yielded good read quality, coverage, and depth. Mean percentage of mappable reads was high at 99.4% with an average total on-target reads per sample of 4.12 9 10 9 and an average depth of 73 reads per target base. Average coverage at >109 depth was 97.9% of all targeted exonic regions. The average numbers of coding SNPs and indels per participant was 19,605 and 465, respectively. AOGC controls had an average depth of 24 reads per target base and 109 coverage of 75.1%. After filtering for rare predicted pathogenic variants only, there was an average of 159 qualifying variants per participant. A total of 1159 variants were within the 86 included genes in case, local or AOGC control cohorts. There was no significant enrichment (OR = 1.12, P = 0.51) in the total carrier rate of qualifying variants in this gene set in glaucoma cases (64.2%, 120 individuals) compared with local and AOGC controls (61.5%, 648 individuals) although there was a trend towards increased mutation load in POAG cases (Table S3).
Four genes (CARD10, CWC27, RERE, and USP37) were nominally enriched (uncorrected P < 0.05) in the glaucoma cohort (Table 1). Only CARD10 (caspase recruitment domain containing protein 10) remained statistically significant after Bonferroni correction for 86 tested genes (P = 6.94 9 10 À5 , corrected P = 5.97 9 10 À3 ) with an odds ratio of 13.2 (3.5-50.2). Eight out of 187 POAG cases (4.28%) carried a rare predicted pathogenic variant in CARD10 confirmed by Sanger sequencing ( Figure S1). Comparatively, qualifying variants in CARD10 were carried by 0.27% (3/1096) of the jointcalled controls. A total of five nonsynonymous variants in this gene were identified in the POAG cohort with all but one absent in the controls (Table 2). These variants were predicted to be pathogenic by SIFT or PolyPhen2 as per the filtering criteria. Additionally, the variants were located in highly conserved regions, all having GERP (Genomic Evolutionary Rate Profiling) scores of >2 (mean GERP score = 4.8). When tested against qualifying variants in the unscreened public domain ExAC non-Finnish European cohort (n = 33370) filtered with our pipeline, CARD10 remained significantly enriched in the glaucoma cohort, albeit with a lower odds ratio (OR = 3.3, P = 3.9 9 10 À3 ). CARD10 variant carriers had younger age at diagnosis (34.9 vs. 44.8 years, P = 0.36) and higher IOP (27.4 vs. 24.8, P = 0.11) than non-CARD10 cases, although the difference was not statistically significant due to the small number of CARD10 carriers. All eight CARD10 qualifying variants were successfully validated in the carrier cases using capillary sequencing ( Figure S1) and submitted to ClinVar database (http://www.ncbi.nlm.nih.gov/clinvar/) with accessions SCV000266585-SCV000266589. Subanalyses of the HTG and NTG cohorts with local and AOGC controls showed significant association of CARD10 variants with HTG subgroup (OR = 15.2, corrected P = 0.01). Five control gene sets of 86 randomly selected genes were examined and none contained any gene that was significantly associated with POAG after Bonferroni correction (Tables S4-S8).

Discussion
Genome-wide association studies have been examined for common variants linked to POAG, NTG and various associated endophenotypes; IOP, vertical cup-to-disc ratio (VCDR), neuroretinal rim area, optic disc area and central corneal thickness. Poor linkage disequilibrium exists between common and rare variants; this has been demonstrated experimentally (Siu et al. 2011). For this reason, previous GWAS designs have reduced power to detect signals at single-nucleotide variants (SNVs) with MAF of <1%. This study employed WES to patch the "gap" left by GWAS and to explore whether rare variants (MAF <1%) within 86 putative POAG risk loci were also enriched within a severe POAG cohort. Although the majority of rare variants from the 86 POAG risk loci do Table 1. Genes at genome-wide significant associated loci that reached nominal significance for enrichment of rare variants in the glaucoma cohort. Mutational load represents the total allele frequency burden of qualifying variants in each gene. Odds ratios are adjusted for capture quality described in the Methods section. P-values are corrected using Bonferroni multiple testing correction.

Odds ratio P-value
Corrected P-value CARD10 2.14 0.14 13.19 6.94 9 10 À5 5.97 9 10 À3 CWC27 1.34 0.14 6.89 9.47 9 10 À3 0.81 RERE 1.60 0.23 5.49 0.0114 0.98 USP37 1.87 0.18 7.32 5.11 9 10 À3 0.44 not show significant enrichment in cases, the plausibility that GWAS candidates may contribute to disease via a common disease rare variant hypothesis has been demonstrated by our study. A previous candidate gene study of SIX6 reported an enrichment of rare variants in POAG in that gene (Carnes et al. 2014). The most prevalent variant, with a carrier frequency of 1.6% in cases, identified by Carnes et al. (2014) (rs146737847:G>A) was found in 1.6% of our POAG cases, but also in 1.9% of local screened controls and 1.6% of AOGC controls with an overall OR of 0.98. The other three rare variants found in that study at a carrier frequency of 0.4% were not detected in our cases or controls. In summary, our study found no enrichment of rare variants in SIX6 in our POAG cases versus controls.

Genome-wide association of CARD10
SNPs near CARD10 were first found to be significantly associated with optic disc area on GWAS in Singaporean Asians, a result replicated in a Dutch Caucasian cohort . Optic disc area is relevant in POAG as a large optic disc is correlated with increased susceptibility to glaucoma in Caucasian populations (Healey and Mitchell 1999). Meta-analysis of optic disc morphology further implicated common variants in CARD10 in VCDR but not POAG after adjustment for optic disc area . This result has since been replicated in a separate case-control study in an Indian population (Philomenadin et al. 2015). Cup-to-disc ratio is fundamental in POAG as it is often used as a diagnostic criterion due to its strong positive correlation with disease incidence (Miglior et al. 2007;Hollands et al. 2013). The SNPs from the original GWAS linking CARD10 to optic disc area were situated between 3262 bp and 7204 bp upstream of CARD10 ). These SNPs are in strong linkage disequilibrium with all common SNPs up to and including rs9610775, which is a missense coding variant, p.R289Q in the CARD10 gene with a population MAF of 13% (Fig. 1). The signal identified in the optic disc area GWAS may possibly relate to CARD10 coding mutations. The SNPs found to be significantly associated with VCDR on meta-analysis were 260 kb upstream of CARD10 with no linkage disequilibrium to the gene ). More likely, these variants are involved in gene expression and regulation. Therefore, the published GWAS data suggests that both coding variants in CARD10 and its expression level may contribute to optic disc pathology.
GWAS was able to highlight the association between regulatory and common coding SNPs in CARD10 with crucial optic disc parameters for the development of POAG Springelkamp et al. 2014;Philomenadin et al. 2015). This study complements the GWAS findings by implicating rare coding CARD10 variants that likely disrupt gene function, and are associated with increased POAG risk within our cohort of advanced glaucoma participants. Both common and rare variants concurrently contribute to disease risk via differing mechanisms affecting the same gene -CARD10. Although more HTG participants harbored rare CARD10 variants than NTG participants, the presence of NTG CARD10 variant carriers suggests that IOP elevation is not required for the potential of CARD10 to cause glaucoma. No definitive inference can be made from this discrepancy due to the small sample size of the subcohorts. A plausible hypothesis would be that the increased susceptibility to retinal ganglion cell apoptosis in the CARD10 mutation carrying individuals is compounded by elevated IOP to cause glaucoma at an earlier age. The participant inclusion criteria of this study selected exclusively for advanced glaucoma phenotype and young age of diagnosis resulting in an enrichment of CARD10 variant carriers that also maintained elevated IOP.
Primary open-angle glaucoma is a disease of enhanced retinal ganglion cell apoptosis (Foster et al. 2002). The NFjB pathway is profoundly involved in regulation of cellular apoptosis. Traditionally, it was only thought to be a promoter of cell survival and proliferation via downstream transcription of anti-apoptotic proteins. As such, overexpression of NFjB leads to enhanced growth of cancerous cells (Escarcega et al. 2007). However more recent evidence in the central nervous system suggests NFjB may also play a proapoptotic role depending on the nature of the noxious stimuli (Kaltschmidt et al. 2005). In certain situations such as neuronal ischemia and Parkinson's disease, upregulation of NFjB led to neuronal death via p53 signaling. Evidence in tumor cells indicates that CARD10 promotes the antiapoptotic effect of NFjB signaling. Therefore loss of function or downregulation of CARD10 leads to increased apoptosis, especially in mice neural crest cells (Grabiner et al. 2007). Conversely, proapoptotic activity of NFjB signalling has been recognized in the rat retina. NMDA (N-methyl-d- aspartate)-induced expression of NFjB p65 led to retinal ganglion cell death which was ameliorated with knockdown of p65 by antisense oligonucleotides (Kitaoka et al. 2007). Retinal ganglion cell damage from NFjB is likely due to glial cell activation and interleukin (IL)-1b secretion (Kitaoka et al. 2007). Given the intimate and complex relationship between CARD10, NFjB and apoptosis, pathogenic coding mutations in CARD10 are likely to affect apoptosis signaling. The coding mutations identified within the glaucoma cohort in this study have not been previously characterized, however they may augment retinal ganglion cell death via up or down-regulation of NFjB. This distinction is important to ascertain in order to translate these findings toward therapeutic strategies in glaucoma management.

Strengths and limitations
A strength of this study is the enrichment of the disease cohort with well-selected extreme disease phenotypes, and it is the largest study to examine whole exomes of a well matched case-control cohort of this nature in the glaucoma field. ANZRAG is the largest collection advanced glaucoma patients worldwide. Our controls are hypernormal and have been individually examined to ensure the absence of glaucoma. Finally, all significantly enriched variants were confirmed via capillary sequencing thereby eliminating the possibility of false-positive calls.
This study has some limitations which were largely technical. Whole exome sequencing has enabled the rapid examination of large numbers of genes at a reasonable cost. Its advantage over microarray technology for genetic studies lies in its ability to capture rare coding variants. However this also introduces the challenge of variable capture. AOGC samples were sequenced on a different platform with a different capture probe set to the case cohort. Joint-calling of the two datasets was used to limit the amount of artifact generated in silico. However, captureand sequencing platform-specific differences remained. Some difference in capture rate at specific locations was evident and translated to incomplete capture of some variants in the AOGC cohort. Adjustment for the capture rate was performed in the calculation of odds ratios and Pvalues. ExAC samples are more heterogeneous in data acquisition compared with this study cohort. Inability to joint-call our data with ExAC may have overestimated the actual mutation burden in this public domain control cohort. Hence the odds ratio between case and controls falls considerably when ExAC data is included and may represent an underestimation. Alternatively, incomplete coverage of AOGC controls may have contributed to an overestimation of the discovery odds ratio. Nevertheless, a statistically significant difference remained between POAG and controls in CARD10 mutational burden even with comparison with ExAC controls.

Conclusion
Overall, the majority of genes near common variants associated with glaucoma and its endophenotypes found in GWAS did not show enrichment of rare variants in POAG. Rare predicted pathogenic variants in CARD10, a gene located near SNPs associated with optic disc area and VCDR in genome-wide association studies, were enriched in a cohort of advanced POAG. Further research is required to elucidate the exact influence of mutations in CARD10 on NFjB signaling, apoptosis, and the development of POAG.

Supporting Information
Additional Supporting Information may be found online in the supporting information tab for this article: Figure S1. Chromatograms of all eight CARD10 mutation carriers in the POAG cohort. Table S1. PCR primers for validation of CARD10 (NM_ 014550.3) variants by direct sequencing. Table S2. List of genes found on GWAS to be associated with POAG or endophenotype. Table S3. Summary of disease burden in genes near GWAS associated SNPs. Table S4. 86 randomly selected genes (Set 1 of 5). Table S5. 86 randomly selected genes (Set 2 of 5). Table S6. 86 randomly selected genes (Set 3 of 5). Table S7. 86 randomly selected genes (Set 4 of 5). Table S8. 86 randomly selected genes (Set 5 of 5).