Association of SIX1/SIX6 locus polymorphisms with regional circumpapillary retinal nerve fibre layer thickness: The Nagahama study

SIX1 and SIX6 are glaucoma susceptibility genes. Previous reports indicate that the single nucleotide polymorphism (SNP) rs33912345 in SIX6 is associated with inferior circumpapillary retinal nerve fibre layer (cpRNFL) thickness (cpRNFLT). Although the region of visual field defect in glaucoma patients is directly related to cpRNFL thinning, a detailed sector analysis has not been performed in genetic association studies. In the present study, we evaluated 26 tagging SNPs in the SIX1/SIX6 locus ±50 kb region in a population of 2,306 Japanese subjects with 4- and 32-sector cpRNFLT analysis. While no SNPs showed a significant association with cpRNFLT in the 4-sectored analysis, the finer 32-sector assessment clearly showed a significant association between rs33912345 in the SIX1/SIX6 locus with inferior cpRNFL thinning at 292.5–303.8° (β = −4.55, P = 3.0 × 10−5). Furthermore, the fine-sectored cpRNFLT analysis indicated that SIX1/SIX6 polymorphisms would affect cpRNFL thinning at 281.3–303.8°, which corresponds to parafoveal scotoma in a visual field test of glaucoma patients.

Recent advances in optical coherence tomography (OCT) have allowed ophthalmologists to perform quantitative evaluations of circumpapillary retinal nerve fibre layer thickness (cpRNFLT), which represents glaucomatous optic neuropathy (GON) with high reproducibility and reliability.This morphological testing method generates objective data less influenced by the problems listed above 26 .The association between cpRNFLT and VFD pattern is highly correlative [27][28][29] , suggesting that cpRNFLT would be a better measure to evaluate the genetic contribution to GON endophenotype than visual field testing.
To date, genetic studies on cpRNFLT have shown consistent contributions of the single nucleotide polymorphisms (SNPs) rs33912345 and rs10483727 located within the SIX1/SIX6 loci to cpRNFL thinning in the superior and inferior sectors, but not in the temporal, nasal, and global sectors [30][31][32] .The significance of the SIX1/SIX6 locus in glaucoma was initially discovered by GWAS for VCDR and primary open-angle glaucoma (POAG), and subsequent studies confirmed the association of polymorphisms in this region with glaucoma onset. 9,33,34 Acording to previous GON structure-function correlations [27][28][29] , a finer cpRNFLT analysis with >4 sectors would be better suited to fully evaluate the risk of QOV-threatening VFD patterns.Therefore, the present study examined genetic associations of the SIX1/SIX6 locus with cpRNFLT using 4-and 32-sector analyses.

Results
Study population.Participants were excluded due to prior intraocular surgery (N = 125), axial length ≥26 mm (N = 258), presence of other ocular disease (N = 64), and cpRNFL image quality (N = 54) as described in the Methods section.A total of 2,306 subjects passed the exclusion criteria and the participant demographics are shown in Table 1.The average age was 57.6 ± 13.6 years (mean ± SD; range, 34-80 years).Mean axial length of the right eye was 23.78 ± 1.04 mm (range, 16.46-25.99mm).Two-thirds of the participants (68.9%) were females; however, this did not affect global cpRNFLT after the adjustment for age and axial length.The mean global cpRN-FLT was 101.6 ± 12.0 μm, and 122.4 ± 20.0 μm, 132.1 ± 20.8 μm, 76.3 ± 13.9 μm, and 75.5 ± 14.3 μm in the superior, inferior, temporal, and nasal sectors, respectively.
Association of SIX1/SIX6 polymorphisms with regional cpRNFLT.Association of the 26 tagging SNPs along the 4 cpRNFLT sectors are shown in Table 2.Only rs12147345 showed a marginal association with cpRNFLT in the temporal region (β = 1.08; 95% confidential interval [CI], 0.26-1.90;P = 0.010).On the other hand, a finer assessment with 32 sectors found 5 additional SIX1/SIX6 SNPs, including rs33912345, to be significantly or marginally associated with a cpRNFLT region, mostly in inferior sectors (Table 3).Among these 6 SNPs, only rs148908311 showed a marginal association with cpRNFLT in the superior region, whereas the other 5 SNPs showed significant or marginal associations with cpRNFLT in the inferior regions.The strongest and only significant association (P < 6.0 × 10 −5 ; 0.05/26 SNPs/32 sectors) was observed between rs33912345 and cpRNFLT in the inferior region at 292.5-303.8°(β = −4.55;95% CI, −2.42-−6.69),P = 3.0 × 10 −5 ).For all SNPs that were significantly or marginally associated with regional cpRNFLTs, we confirmed identical effect directions of the genetic variances to the regional cpRNFLT as to the global cpRNFLT.
We observed that rs10483727-the only SIX1/SIX6 SNP previously associated with glaucoma susceptibility by GWAS-was in close linkage disequilibrium (LD) with rs33912345 (R-sq = 0.99 in the JPT 1000 genomes dataset).Although rs10483727 was not included in the tagging SNPs in the present study, the C risk allele of rs33912345 (corresponding to the T risk allele of rs10483727) resulted in both regional and global cpRNFL thinning.

Discussion
Our study used a tagging SNP approach to show that polymorphisms in the SIX1/SIX6 region was significantly associated with inferior cpRNFLT and marginally associated with superior cpRNFLT in a community-based Japanese cohort.These results were consistent with previous candidate SNP evaluations on rs33912345 [30][31][32] .In addition, our findings suggest that the 32-sector region-based approach for cpRNFLT enables the detection of underpowered and undermined genetic associations in 4-sector analyses.
The importance of the SIX1/SIX6 locus in glaucoma was initially discovered by a GWAS for VCDR, and subsequent GWAS for POAG confirmed the association between polymorphisms in this locus with glaucoma onset 9,33,34   This is the only locus where an association with cpRNFLT has been established [30][31][32] .In this study, our detailed analysis further specified the cpRNFL region of association, and revealed that SIX1/SIX6 affects cpRNFL thinning at 281.3-303.8°among the inferior region, which may have clinical and biological relevance since retinal nerve fibre layer defects also occur in this region.In addition, according to the Garway-Heath map 27 , RNFL thinning of this region would lead to upper mid-peripheral VFD and so-called early upper nasal step 35 .Because VFD in this region typically results from early-stage glaucoma, SIX1/SIX6 could be associated with initial changes in GON; however, optic fissure closures-known as colobomas-also present in this area.While eyes with apparent colobomas were excluded from our analysis, subclinical cases would likely influence the statistical data.Moreover, since PAX6 mutations are associated with coloboma formation 36 and correlate with SIX6 activation during eye development 37 , SIX1/SIX6 polymorphisms could be involved in the pathophysiology of cpRNFL thinning in the inferior region.Although previous reports have not shown a connection between SIX6 SNPs and coloboma formation 38,39 , further genetic studies with 32-sector cpRNFLT analysis would likely lead to the identification of other genes with key roles in glaucomatous VFD development.
Our findings would suggest that we should evaluate genetic associations to cpRNFLT by dividing it into 32 sectors rather than dividing it into 4 sectors or evaluating cpRNFLT as a whole.Notably, we found a significant association between rs33912345 and inferior region cpRNFLT at 292.5-303.8°(P = 0.025 after Bonferroni correction) in the 32-sectored cpRNFLT analysis, whereas rs33912345 was not significantly associated with the inferior region (225-315°) in the 4-sectored analysis (P = 1.0 after Bonferroni correction).Furthermore, all SNPs with marginal P-values at 292.5-303.8 and 281.3-292.5°showed an equivalent contribution to cpRNFL thinning to that observed with rs33912345.In addition, rs10483727-a RNFLT-susceptible SNP not included in our analysis-showed a significant association with inferior region cpRNFLT at 281.3-303.8°(P = 0.016 and P = 0.016 after Bonferroni correction, respectively) in the 32-sector analysis, but not in the 4-sector analysis (225-315°; P = 1.0 after Bonferroni correction) (Supplementary Tables 1 and 2).
Based on a previous study, cpRNFL thinning at 281.3-303.8°should correspond to mid-peripheral scotoma since central VFD is associated with cpRNFL thinning at 311-40° 27 .VFD is usually classified as central/ mid-peripheral/peripheral scotoma, superior/inferior altitudinal defect, or temporal/nasal hemianopia that shared features of the clinically observed visual field patterns 40 .The altitudinal boundary is separated at 12°    around optic nerve head (ONH) 41 ; thus, cpRNFL should be evaluated at 12° interval or less (at least 30 sectors) to fully evaluate its correspondence to VFD patterns.To our knowledge, this is the first study evaluating the applicability of a region-based approach for GON analysis by examining the genetic contributions of glaucoma susceptibility genes with regional cpRNFLT.Therefore, genetic studies using 32-sectored cpRNFLT might reveal further associations to clinically important glaucoma phenotypes.
To date, three studies have reported significant associations of SIX1/SIX6 polymorphisms to cpRNFL thinning in the upper and lower sectors-but not the nasal and temporal sectors-in 30 POAG cases 32 , 1,243 population controls 31 , and 231 other participants consisting of 20% normal, 44% of suspected glaucoma, and 36% confirmed glaucoma cases 30 .In contrast, our analyses used 2,306 population controls and showed that rs33912345 had the strongest association with cpRNFLT at 292.5-303.8°and a marginal association at 78.8-90.0°.RNFL thinning at 281.3-303.8°would lead to upper mid-peripheral VFDs that could lead to an early upper nasal step 27,35 .However, further genetic studies on visual field testing are necessary to confirm whether the SIX1/SIX6 locus would be an appropriate locus to determine the genetic factors underlying GON in the upper visual fields.
There are several limitations to this study.First, the VFD data was not obtained in our cohort.As the clinical importance of VFD has been widely accepted in GON studies, further confirmation analyses on VFD are needed.Nevertheless, we believe that cpRNFLT is an objective value with high reliability and repeatability and can facilitate the identification of hidden genetic associations for VFD.Second, only data from the right eye of subjects were analysed due to the time constraints placed on OCT acquisition; however, because both eyes are equally affected by genotype, analyses on the left eyes should yield similar findings.Third, the number of participants analysed was rather small compared to the large sample size of the Nagahama cohort.This is mainly because only subjects with genome-wide SNP data and an axial length <26 mm were included in the study.A larger sample analysis might further elucidate its associations to other regions.Fourth, population-based study methods are best evaluated with disease-free subjects.Despite the known association between SIX1/SIX6 SNPs and glaucoma development, these patients were not specifically excluded from the study population since we did not perform visual field testing or slit lamp biomicroscopy required for this diagnosis.Thus, a subsequent study of only healthy subjects will be necessary to confirm the clinical impact of our findings.Lastly, an optimization of sector number may be beneficial in future cpRNFLT studies and it is possible that wider sectors would be sufficient to yield the same result, whereas a finer sector analysis might be able to find other associations.
In conclusion, we confirmed that rs33912345 and rs10483727-the only known cpRNFLT susceptibility SNPs-showed the strongest association with cpRNFLT of those within the SIX1/SIX6 locus.Notably, only the 32-sector cpRNFLT analysis was capable of detecting the significant associations these SNPs with inferior cpRNFL thinning at 292.5-303.8°and 281.3-303.8°,respectively, as the results of 4-sector cpRNFLT analysis were insignificant.Collectively, this suggests that fine regional association analyses are a more effective strategy to assess glaucoma endophenotype in genomic studies and may facilitate the identification of novel genetic associations in disease pathogenesis.

Methods
Ethical considerations.Written informed consent was obtained from all participants.Study procedures adhered to the tenets of the Declaration of Helsinki and were approved by the ethics committee of Kyoto University Graduate School of Medicine and the Nagahama Municipal Review Board.

Study participants.
The study population consisted of healthy Japanese volunteers enrolled in the Nagahama Prospective Cohort for Comprehensive Human Bioscience (the Nagahama Study).Participants were recruited between 2008 and 2010 from the general population of Nagahama City, a rural city of 125,000 inhabitants located in central Japan.Community residents from 30-74 years of age, living independently and without physical impairment or dysfunction were eligible.Of the 9,804 included participants, nine withdrew consent to participate, and 26 were excluded because genetic analysis showed an ethnic background other than Japanese.Participants were offered a follow-up assessment 5 years after the baseline evaluation, and 8,294 of the original 9,769 cohort members participated (84.9%).
In the present study, we used a dataset of the follow-up measurement.Study subjects consisted of 2,807 individuals with genome-wide SNP genotyping, axial length, phakic status, and OCT data available by April 2016.Other exclusion criteria included prior intraocular surgery (except for cataract surgery), high myopia (axial length ≥26 mm), and presence of other ocular diseases affecting retinal nerve fibre layer thickness based on fundus photography-such as optic atrophy, anterior ischemic optic neuropathy, optic disc coloboma, retinal vein occlusion, proliferative or severe non-proliferative diabetic retinopathy, retinitis pigmentosa, and other optic neuropathies.Subjects outside of the Japanese ethnic cluster or with poor quality cpRNFL images, which could result from cataracts or small pupils, were also excluded.Glaucomatous status did not serve as exclusion criteria since visual field information was unavailable at the time of analysis.Ultimately, a total of 2,306 subjects with a ≤ 0.9 sample call rate and estimated relatedness (PI-HAT) >0.35 were included in the study population All subjects were assessed by standardized ophthalmic examination, including an objective determination of the refractive error and corneal curvature (Autorefractor ARK-530; Nidek, Gamagori, Aichi, Japan), fundus imaging (CR-DG10; Canon, Tokyo, Japan), and axial length measurements by partial coherence interferometry (IOL Master; Carl Zeiss Meditec, Inc., Dublin, California, USA).The cpRNFL in the right eye was imaged by spectral-domain optical coherence tomography (SD-OCT) (RS-3000 advanced; Nidek, Gamagori, Aichi, Japan).
Circumpapillary retinal nerve fibre layer thickness.The RS-3000 advanced OCT (Nidek) was used to obtain circular B-scans 11.5° in diameter (3.45 mm in the Gullstrand's eye) centred on the optic disc, i.e., a circumpapillary scan.Each B-scan was obtained by averaging 50 images in "Regular mode" to reduce speckle noise.The cpRNFL thickness was defined as the distance between the inner border of the internal limiting membrane (ILM) and the outer border of the RNFL in B-scan images, measured automatically with built-in software, and then manually corrected for all images.We excluded eyes with extensive peripapillary atrophy (affecting cpRNFL scans), RNFL schisis, peripapillary epiretinal membrane, or thickened posterior vitreous membrane that affected segmentation, and low quality or other RNFL segmentation errors.
RNFL thickness was measured at 1,024 points along the 360° OCT circle scan, which were subsequently sectioned into 4 or 32 sectors (Fig. 1).Data for each sector were averaged and the associations between the SIX1/SIX6 SNPs and mean cpRNFLT in each region were analysed.
Genotyping and imputation.DNA samples were prepared and genotyped as described previously 42 .Briefly, 3,712 baseline samples were genotyped using at least one of the three genotyping platforms, HumanHap610K Quad Arrays, HumanOmni2.5 M Arrays, or HumanExome Arrays (Illumina, Inc., San Diego, CA).To ensure high-quality genotype data, a series of quality control (QC) filters, including sample success rate (>95%), individual call rate (>99%), minor allele frequency (MAF) cut-off (>0.01),Hardy-Weinberg equilibrium p-values (>1 × 10 −6 ), and estimated relatedness (PI-HAT < 0.35) were applied to the data from each platform.In addition, seven ancestry outliers were identified by principal component analysis with the HapMap Phase 2 release 28 with the Japanese in Tokyo, Japan (JPT) reference dataset using EIGENSTRAT ver.2.0.QC in PLINK 39,43 (ver.1.07;available at http://pngu.mgh.harvard.edu/~purcell/plink/).As a result, 3,267 baseline and 2,807 follow-up samples passed the QC filters.SNP genotype imputation was performed for the Japanese samples using MaCH 44 .Genotypes of 89 JPT samples from the 1000 Genomes Project (May 2011 release) were used as reference sequences.Imputed SNPs with an R-squared value less than 0.5 were excluded from the following association analyses.
Our dataset contained 488 SNPs within ±50 kb of the SIX1/SIX6 locus (chr14: 60925938-61166155; NCBI build 37).The Tagger program in Haploview 45 was used to identify 26 tagging-SNPs encompassing the 242 known SNPs with an MAF > 0.05 (mean R-sq = 0.962).Since rs33912345 was previously reported as a cpRNFLT-susceptible SNP, it was selected as a tagging SNP for a positive control.The location of these 26 SNPs within the SIX1/SIX6 locus and a linkage disequilibrium plot of this genetic region are shown in Fig. 2. Statistical analysis.Linear regression analyses were performed to determine the associations between regional cpRNFLT and the 26 SNP genotypes assuming additive regression models for a per-minor-allele with an adjustment for age, sex, and axial length.We evaluated regional cpRNFLT associations using 4 divided sectors and finer 32 sectors with Bonferroni corrections.P values < 6.0 × 10 −5 (0.05/26 SNPs/32 sectors) and <4.8 × 10 −4 (0.05/26 SNPs/4 sectors) were considered statistically significant for cpRNFL in the 32-sector and 4-sectors analyses, respectively, whereas P-values < 0.0016 (0.05/32 sectors) and <0.0125 (0.05/4 sectors) were considered marginally significant, respectively.

Figure 1 .
Figure 1.Illustration of the 32 sectors along the RS-3000 advanced (Nidek) optical coherence tomographic (OCT) circular image used for circumpapillary retinal nerve fibre layer thicknesses (cpRNFLT) analysis.(A,B) Circular B-scans 11.5° in diameter (3.45 mm in the Gullstrand's eye) centred on the optic disc were obtained from the right eye of each participant.(C) cpRNFLT was measured at 1,024 points along the 360° OCT circle scan and divided into 4 or 32 sectors.

Figure 2 .
Figure 2. Construction of the SIX1/SIX6 region and tagging SNPs.(A) Location of the 26 tagging SNPs within the SIX1/SIX6 locus are shown relative to NCBI build 37 of the human genome.(B) A linkage disequilibrium map of the SIX1/SIX6 locus ±50 kb region encompassing the 242 SNPs examined in our analysis was plotted with Haploview 4.2 software.A total of 16 haplotype blocks could be distinguished and 10 SNPs were not included in any of the blocks.

Table 1 .
Demographics of study participants.Data represent the mean ± SD.IOP, intraocular pressure, cpRNFLT, circumpapillary retinal nerve fibre layer thickness.† Data were acquired from the right eye.

Table 2 .
Associations of a SIX1/SIX6 locus and circumpapillary retinal nerve fibre layer thickness including 4 sectors.CHR, chromosome; SNP, single nucleotide polymorphism; BP, base pair.*Positions and alleles are given relative to the positive strand of NCBI build 37 of the human genome.† Linear regression analyses were applied assuming additive effect of the per minor allele variant, adjusted for age and sex.Significant (P < 4.8 × 10 −4 ) or suggestive (P < 0.0125) associations are shown in bold.