Variation in PTCHD2, CRISP3, NAP1L4, FSCB, and AP3B2 associated with spherical equivalent.

PURPOSE: Ocular refraction is measured in spherical equivalent as the power of the external lens required to focus images on the retina. Myopia (nearsightedness) and hyperopia (farsightedness) are the most common refractive errors, and the leading causes of visual impairment and blindness in the world. The goal of this study is to identify rare and low-frequency variants that influence spherical equivalent.
METHODS: We conducted variant-level and gene-level quantitative trait association analyses for mean spherical equivalent, using data from 1,560 individuals in the Beaver Dam Eye Study. Genotyping was conducted using the Illumina exome array. We analyzed 34,976 single nucleotide variants and 11,571 autosomal genes across the genome, using single-variant tests as well as gene-based tests.
RESULTS: Spherical equivalent was significantly associated with five genes in gene-based analysis: PTCHD2 at 1p36.22 (p = 3.6 × 10(-7)), CRISP3 at 6p12.3 (p = 4.3 × 10(-6)), NAP1L4 at 11p15.5 (p = 3.6 × 10(-6)), FSCB at 14q21.2 (p = 1.5 × 10(-7)), and AP3B2 at 15q25.2 (p = 1.6 × 10(-7)). The variant-based tests identified evidence suggestive of association with two novel variants in linkage disequilibrium (pairwise r(2) = 0.80) in the TCTE1 gene region at 6p21.1 (rs2297336, minor allele frequency (MAF) = 14.1%, β = -0.62 p = 3.7 × 10(-6); rs324146, MAF = 16.9%, β = -0.55, p = 1.4 × 10(-5)). In addition to these novel findings, we successfully replicated a previously reported association with rs634990 near GJD2 at 15q14 (MAF = 47%, β = -0.29, p=1.8 × 10(-3)). We also found evidence of association with spherical equivalent on 2q37.1 in PRSS56 at rs1550094 (MAF = 31%, β = -0.33, p = 1.7 × 10(-3)), a region previously associated with myopia.
CONCLUSIONS: We identified several novel candidate genes that may play a role in the control of spherical equivalent. However, further studies are needed to replicate these findings. In addition, our results contribute to the increasing evidence that variation in the GJD2 and PRSS56 genes influence the development of refractive errors. Identifying that variation in these genes is associated with spherical equivalent may provide further insight into the etiology of myopia and consequent vision loss.

Introduction

Uncorrected refractive errors, the leading cause of visual impairment and blindness worldwide [1], are also associated with several ocular disorders, such as glaucoma, cataract, retinal detachment, and macular degeneration, further increasing the risk of vision loss in later life [2,3]. In the United States, myopia affects 1 in 3 individuals 20 years of age or older, and the prevalence of hyperopia increases with increasing adult age [4]. Although mild refractive errors can be easily managed with non-invasive optical means or refractive surgery, determining the underlying causes of refractive errors may lead to better management of their clinical complications.

The etiology of refractive errors is complex and not fully understood. Many environmental factors have been linked to increased myopia, including high intensity near-work [5], high intelligence [6-8], and low outdoor exposure [5,9-11]. These associations have been inconsistent across studies [9,12,13], and the causative relationship between these factors and spherical equivalent has not been established. In addition to environmental factors, genetic factors play an important role in control of spherical equivalent. The estimated heritability of spherical equivalent in populations of European descent ranges from 0.50 to 0.86 [14-21]. Genome-wide linkage studies have reported evidence of genetic linkage to more than 30 chromosomal regions for myopia and spherical equivalent [22]. Genome-wide association studies (GWASs) targeting more common variants have identified susceptibility loci associated with spherical equivalent or myopia in East Asian populations [23-30] and in populations of European ancestry [30-35]. However, in most cases the biologic basis of these associations have yet to be experimentally validated [36].

Evolutionary theory and empirical data suggest that rare and low-frequency variants could have relatively large effects on the risk of developing complex traits, and collectively rare variants may make significant contributions to disease susceptibility in a given population [37]. However, these variants were poorly characterized by earlier generations of GWAS genotyping arrays [38]. Sequencing of the exome or whole genome directly assays all variants, including rare and low-frequency variants, but it is currently costly to perform on thousands of study subjects. Exome arrays, with enhanced resolution of coding variants, provide a less costly alternative to sequencing to evaluate the role of rare (minor allele frequency (MAF) < 1.0%) and low-frequency (MAF = 1.0 – 5.0%) variants in phenotypes.

To examine the impact of rare and low-frequency coding variants on spherical equivalent, we conducted an exome array analysis on participants in the Beaver Dam Eye Study (BDES). Additionally, data permitting, we evaluated known regions of association reported in other populations for spherical equivalent and myopia.

Methods

Study participants

The BDES is a population-based cohort established in 1987 to study age-related eye disorders [39]. A total of 5,924 individuals aged from 43 to 84 years were identified from 3,715 households through a private census of Beaver Dam, Wisconsin. Between 1988 and 1990, 4,926 (83.14%) individuals participated in the baseline examination. Familial relationships of 2,783 participants were confirmed and the data was assembled into pedigrees. To improve the sampling efficiency, we genotyped individuals with intraocular pressure or spherical equivalent values at the baseline, resulting in a nested sub-cohort of 2,032 individuals with the full spectrum of spherical equivalent values. The study adhered to the ARVO statement on human subjects and was approved by the institutional review board at the University of Wisconsin. The study complied with the guidelines of the Declaration of Helsinki. Informed consent was obtained from all participants in the study.

Clinical evaluation

All participants in the BDES received a complete eye exam that included standardized evaluation of refraction using the Humphrey 530 refractor (San Leandro, CA) [39]. When the automated refraction yielded visual acuity of 20/40 or worse, a modification of the Early Treatment Diabetic Retinopathy Study (ETDRS) protocol [40] was followed to obtain the best-corrected visual acuity. The mean spherical equivalent (MSE) in diopters for each eye was calculated using the standard formula (MSE = sphere + 0.5 × cylinder). Nuclear lens opacity was determined by grading slit-lamp lens photographs using a standard protocol, resulting in a five-level scale [41]. Eyes without a lens, with an intraocular lens, or with best-corrected visual acuity of 20/200 or worse were excluded. Age, sex, and education were also obtained at the baseline visit.

Genotyping and quality control

DNA from 2,032 BDES participants, 22 blind duplicate samples, and 44 HapMap samples were genotyped on the Illumina Infinium HumanExome-12 v1.1 (Illumina, Inc., San Diego, CA) array at the Genetic Resources Core Facility (GRCF), the Johns Hopkins School of Medicine, the Institute of Genetic Medicine (Baltimore, MD). Of the 2,032 BDES participants, 1,908 (93.9%) were successfully genotyped and had a call rate >98%. After related individuals (n = 126) and individuals with unresolved gender inconsistencies (n = 15), missing spherical equivalent measurements of either eye (n=118), differences in spherical equivalent > ± 4 D between the eyes (n = 15), and missing values in age, sex, education, or nuclear sclerosis (n = 74) were excluded, 1,560 individuals remained in the final analytical cohort.

Of the 242,901 variants genotyped, we excluded the following: sex chromosomal variants (n = 5,465), insertion/deletion variants (n = 134), variants with call rates ≤98% (n = 6,027), monomorphic variants (n = 132,213), and duplicate variants (n = 314). The final data set included 98,748 polymorphic variants. We evaluated Hardy–Weinberg equilibrium (HWE) on these variants and reported the p value estimated from the HWE exact test for significant associations.

To evaluate genetic ancestry, we performed principal component analysis (PCA) using SMARTPCA in EIGENSTRAT (version 4.2) [42]. All BDES participants genotyped were self-reported non-Hispanic white. The first two principal components were included in the analysis to control for potential population stratification, resulting in a genomic inflation factor (λ) of 1.02.

Trait definitions

Age, education, and nuclear sclerosis have known effects on spherical equivalent and were significant covariates in our previous linkage analysis [43]. Using linear models, we regressed age (mean centered), education (mean centered), sex, nuclear sclerosis (sum of the nuclear lens opacity score of both eyes), and the first two principal components on the average of spherical equivalent in the right and left eyes. The residuals from this regression were used as the quantitative outcome in our association tests. For replications of loci associated with myopia in prior GWASs, we performed logistic regression analysis that compared myopia to hyperopia and myopia to non-myopia, adjusting for the covariates. We defined individuals with spherical equivalent < −1 D in either eye and < −0.5 D in the opposite eye as myopic and those with spherical equivalent > +1 D in either eye and > + 0.5 D in the opposite eye as hyperopic. Emmetropia was defined as spherical equivalent ranged from −1 D to +1 D in both eyes. The non-myopic group included emmetropes and hyperopes. Individuals who fell outside these categories were excluded from the logistic regression analysis.

Statistical analysis

We applied single-variant and gene-based tests to investigate the association of rare and low-frequency genetic variants with spherical equivalent as a quantitative trait. Single-variant association analysis was performed on variants with MAF ≥1.0% (n = 34,976) under an additive model using PLINK (version 1.07) [44]. To determine the number of independent tests, we tabulated the number of single nucleotide polymorphisms (SNPs) with pairwise linkage disequilibrium (LD) <0.2. There were 23,342 independent variants, which translated to a Bonferroni corrected significance threshold of p<2.1 × 10−6.

Gene-based analysis can be a powerful alternative approach to single-variant analysis in detecting effects of rare variants [45]. Rare variants in our analysis were defined as those with MAF of 1% or less. To examine the impact of rare variants on spherical equivalent, we used the following tests: (1) the burden test, where all variants were collapsed into a gene-based score [46], and (2) the sequence kernel association test (SKAT) [47]. Analysis was restricted to 11,571 genes with at least two rare variants. To optimize the statistical power of SKAT, we applied the Madsen-Browning weighting

which slightly up-weights the rare variants [48]. A p value of <4.3 × 10−6, corresponding to a Bonferroni correction for 11,571 independent gene tests, was used to determine significance.

We evaluated gene regions with GWAS significant associations in published studies or the region spanning two genes if the GWAS association was intergenic [23-34]. For each known GWAS locus, we performed a single-variant association analysis for quantitative refraction and myopia, targeting (1) the exact variant or (2) variants in strong LD (r2 >0.8) with the reported variants. Variants were considered as replicated at p<0.01 if the magnitude and direction of their effects in our study were consistent with those reported in previous studies.

Results

Characteristics of the Beaver Dam Eye Study and the genotyped sub-cohort included in these analyses are presented in Table 1. The baseline mean spherical equivalent in the subset was −0.20 D with a standard deviation of 2.9 D (Appendix 1). Compared to the entire BDES cohort, our subset included a higher proportion (38.6% versus 21.5%, respectively) of myopes and a similar proportion of hyperopes (37.0% versus 37.8%, respectively).

Table 1

The characteristics of the Beaver Dam Eye Study (BDES) participants.

Characteristic	BDES Exome (n=1,560)	BDES Fulla (n=4,452)
Age (years), mean ± SD (range)	59.5±10.6 (43 – 85)	61.1±10.8 (43 – 86)
Female gender, N (%)	884 (56.7)	2,493 (56.0)
Education (years), mean ± SD	12.7±3.1	12.0±2.8
Nuclear Sclerosis b, mean ± SD	4.7±1.6	4.8±1.6
Spherical Equivalent (D)c, mean ± SD (range)	−0.20±2.9 (−13.0 – 10.3)	0.28±2.2 (−13.0 – 10.3)
Myopia, N (%)	602 (38.6)	958 (21.5)
Emmetropia, N (%)	242 (15.5)	1,217 (27.3)
Hyperopia, N (%)	577 (37.0)	1,681 (37.8)

a All BDES individuals with reliable spherical equivalent measurements of both eyes at baseline visit. bThe sum of nuclear lens opacity grading of the right and left eyes at baseline. cThe average of spherical equivalent of the right and left eyes at baseline. Myopia was defined as < −1 D in either eye and < −0.5 D in the opposite eye; hyperopia was defined as > +1 D in either eye and > +0.5 D in the opposite eye; emmetropia was defined as between −1 D and +1 D in both eyes.

We tested 34,976 variants for their individual effect on spherical equivalent (Figure 1). No variant reached significance (p<2.1 × 10−6). However, suggestive associations were identified with two variants on chromosome 6p21.1 in the TCTE1 gene (OMIM 186975) region (Appendix 2). A missense SNP, rs2297336 (MAF = 14.1%), causes a phenylalanine to serine substitution (p.Phe261Ser). Each copy of the minor allele G of rs2297336 was associated with a 0.62 D (95% CI = [–0.88, −0.36], p = 3.7 × 10−6) decrease in spherical equivalent. Functional prediction programs predict this variant to be “deleterious” (SIFT [49]) and “probably damaging” (PolyPhen [50]) indicating this variant is likely to have a functional impact. This variant is within a transcription factor (TF) binding site [51]. In addition, a missense SNP rs324146 (p.Pro35Leu, MAF = 16.9%) in LD with rs2297336 (r2 = 0.80) in this study sample also had suggestive evidence of association. Each additional copy of the minor allele G of rs324146 was associated with a 0.55 D lower spherical equivalent (95% CI = [–0.79, −0.30], p = 1.4 × 10−5).

Figure 1

Manhattan plot for the single-variant association results. A total of 34,976 variants with minor allele frequency (MAF) ≥1.0% were tested for single-variant associations. The significance threshold across the exome was set to p<2.1×10−6 (red line). Suggestive associations were discovered on two variants in the TCTE1 gene region at 6p21.1: rs2297336 (p = 3.7 ×10−6) and rs324146 (p = 1.4 ×10−5).

A total of 11,571 gene regions were tested in the gene-based analysis. In these analyses, five gene regions showed significant (p<4.3 × 10−6) evidence of association (Table 2, Appendix 3) with spherical equivalent. The most significant association was with the FSCB gene (OMIM 611779) on chromosome 14q21.2 (p = 1.5 × 10−7). Six protein-altering (one stop-gain and five missense) variants contributed to the gene-based finding.

Table 2

Significant gene-based association results (p<4.3 × 10−6).

Gene	Locationa	Best Test	P value	No. of variants	Amino Acid Alteration (N)b
PTCHD2	1p36.22	SKAT	3.6 × 10−7	4	p.Arg260Pro(2), p.Arg542His(16), p.Pro1034=(1), p.Val1178Ile(2)
CRISP3	6p12.3	SKAT	4.3 × 10−6	3	p.Pro197Thr(3), p.Gly141Ala(2), p.Arg103Gln(1)
NAP1L4	11p15.5	SKAT	3.6 × 10−6	2	p.Ala369Val(1), p.Thr111Ile(2)
FSCB	14q21.2	SKAT	1.5 × 10−7	6	p.Val714Ile(2), p.Ala392Thr(8), p.Ala310Val(1), p.Ala281Val(19), p.Lys251Arg(1), p.Gln74Ter(8)
AP3B2	15q25.2	Burden	1.6 × 10−7	3	p.Ser946Thr(3), p.Lys880Arg(1), p.Lys262Thr(2)

aGene locations reported in GRCh37. bAmino acid alteration reported in dbSNP. n=no. of minor allele copies in the analyzed cohort.

We identified four additional gene regions with significant association with spherical equivalent: the AP3B2 gene (OMIM 602166) on chromosome 15q25.2 (p = 1.6 × 10−7), the PTCHD2 gene (OMIM 611251) on chromosome 1p36.22 (p = 3.6 × 10−7), the NAP1L4 gene (OMIM 601651) on chromosome 11p15.5 (p = 3.6 × 10−6), and the CRISP3 gene on chromosome 6p12.3 (p = 4.3 × 10−6). Although the SKAT and the burden test supported a significant association between rare variants in the AP3B2 gene and quantitative refraction, the p value was more significant in the burden test, possibly suggesting that the variants included in the test influence spherical equivalent in the same direction and with similar magnitudes. Although these findings are intriguing, due to the low frequency of rare variants in these four genes (Table 2, Appendix 3), the observed associations were driven by a small number of observations and should be interpreted with caution.

The exome array contained either the exact variant or at least one variant in strong LD for eight previously reported GWAS loci: SIX6 (OMIM 606326) [30], GJD2 (OMIM 607058) [30,31,34], RASGRF1 (OMIM 606600) [30,32,34], CNDP2 (OMIM 169800) [30], PDE11A (OMIM 604961) [34], PRSS56 (OMIM 613858) [30,34], 4q25 [24], and MIPEP (OMIM 602241) [26]. Variants of GJD2, RASGRF1, and PRSS56 had been previously associated with spherical equivalent and myopia. In our analysis, we had evidence of replication for two of these eight regions (Table 3). We replicated the association of rs634990 located 37.2 kb upstream of the GJD2 gene on chromosome 15q14 with spherical equivalent (Appendix 4). In the Beaver Dam Eye study, the minor allele C of rs634990 was associated with a lower spherical equivalent (MAF = 47%, β = –0.29, 95% CI = [–0.48, −0.11], p = 1.8 × 10−3). In the comparison between myopic and hyperopic individuals, the odds ratio (OR) of having myopia was 1.57 greater (95% confidence interval (CI) = [1.15, 2.15], p = 4.5 × 10−3) for those carrying one copy of the C allele of rs634990 and 1.85 greater (95% CI = [1.27, 2.70], p = 1.4 × 10−3) for the homozygous C carriers. As extreme phenotype samples have been shown to increase power in genetic association studies [52,53], when emmetropic individuals were included in the control group, the association of rs634990 with myopia became smaller and marginally significant (heterozygous carriers: OR = 1.36 [1.03, 1.78], p = 0.03; homozygous carriers: OR = 1.48 [1.07, 2.04], p = 0.02). In addition, rs1550094 on chromosome 2q37.1 in the PRSS56 gene region was significantly associated with spherical equivalent (Appendix 4). Each copy of the minor allele G of rs1550094 on average reduced spherical equivalent by 0.33 D (MAF = 31%, 95% CI = [–0.53, −0.12], p = 1.7 × 10−3). Although rs1550094 was associated with myopia in the analysis of myopia versus hyperopia, this association was not statistically significant for the homozygous carriers (OR = 1.28 [0.80, 2.05], p = 0.30), only for the heterozygous carriers (OR = 1.44, 95% CI = [1.09, 1.91], p = 0.01). A similar trend was observed in the comparison of the myopia versus non-myopia (heterozygous carriers: OR = 1.41 [1.10, 1.79], p = 5.8 × 10−3; homozygous carriers: OR = 1.33 [0.88, 2.00], p = 0.17).

Table 3

Association results for genetic variants in published GWAS loci for spherical equivalent or myopia.

Loci	Gene	GWAS variantsb	GWAS effect size	Tested variant (r2)c	MAF	β (SE)	P value
Spherical equivalent
14q23.1	SIX6 [30]	rs1254319	−0.09	rs1254319 (1.00)	29% [A]	−0.17 (0.10)	>0.01
15q14	GJD2 [31]	rs634990	−0.23	rs634990 (1.00)	47% [C]	−0.29 (0.09)	1.8 × 10−3
15q24.2	RASGRF1 [32]	rs939658	−0.15	rs939658 (1.00)	50% [A]	−0.04 (0.10)	>0.01
18q22.3	CNDP2 [30]	rs12971120	0.11	rs2278161 (1.00)	22% [C]	0.14 (0.12)	>0.01

aLocations are reported in GRCh37 coordinates. bThe index variant reported in published GWAS. cExact variant or top variant in strong LD with the reported GWAS variant (r2 >0.80 in 1000 Genome Pilot 1 CEU population). dHazard ratio. eAllelic odds ratio. fHomozygous/Heterozygous odds ratio.

Discussion

In this study, we performed an exome array analysis on 1,560 BDES participants to identify rare and low-frequency variants associated with spherical equivalent. We discovered novel associations with rare variants in five gene regions, PTCHD2, CRISP3, NAP1L4, FSCB, and AP3B2, and suggestive associations with common variants in TCTE1. In addition, we replicated the previously reported association of GJD2 with spherical equivalent and myopia and discovered an association of spherical equivalent on PRSS56, a locus that had been previously identified for myopia.

The rare variants in the identified genes may not explain much of the phenotypic variation in quantitative refraction in the general population due to their low frequencies. However, the novel loci associated with spherical equivalent may not only be of great influence in a small number of families that harbor these rare variants but also provide valuable insights into novel genes and pathways involved in the etiology of refractive errors. The AP3B2 gene is highly expressed in adult RPE cells [54,55]. The AP3B2 protein is a subunit of the neuron-specific AP3 complex involved in neuronal membrane trafficking via mediating the direct transport of lysosomal proteins from the trans-Golgi network to the lysosomes and lysosome-related organelles [56,57]. A recent study found binding activity between the AP3B2 protein and a membrane scaffold protein, 4.1G, whose deficiency can lead to mislocalization of photoreceptor terminals and visual acuity impairments in mice [58]. The AP3B2 gene has been implicated in the lysosome pathway [59-62]. At the significance level of 0.01, we also found evidence of association with two other genes in the lysosome pathway in the gene-based analysis: SLC17A5 (OMIM 604322) on chromosome 6q13 (p = 2.0 × 10−3) and MFSD8 (OMIM 611124) on chromosome 4q28.2 (p = 2.7 × 10−3; Appendix 5). These findings may suggest the potential involvement of one or more genes of the lysosome pathway in the development of refractive errors.

The gene PTCHD2 on chromosome 1p36.22 is located in the myopia locus 14 (MYP14, OMIM 610320), a known linkage region for myopia and spherical equivalent [63,64]. PTCHD2 is also overlapped with a linkage region for the Volkmann type of cataract (CTRCT8, OMIM 115655) and a locus for primary infantile glaucoma (GLC3B, OMIM 600975), suggesting its potential role in the development of several ocular conditions. This gene is highly expressed in retinal ganglion cells [54]. The PTCHD2 protein interacts with several other proteins in the Hedgehog signaling pathway and may be involved in regulation of lipid transport and cholesterol homeostasis [65,66].

The genes FSCB, NAP1L4, and CRISP3 are all expressed in the retina [54] but do not have a clear mechanism by which they may influence spherical equivalent. The product of the FSCB gene is colocalized with a calcium-binding protein, CABYR, in the mouse testis, and was postulated to be involved in the later stages of fibrous sheath assembly [67]. Alterations of the FSCB gene have also been reported in patients with osteosarcoma [68] and in families with cerebral palsy [69]. The NAP1L4 gene encodes a member of the nucleosome assembly protein (NAP) family. This protein is a highly conserved histone chaperone protein that regulates the histone H2A-H2B heterodimer and mediates nucleosome formation [70,71]. This gene is located near the imprinted gene domain of 11p15.5, an important tumor-suppressor gene region that has been linked to growth-related disorders such as Beckwith-Wiedemann syndrome and various human cancers [72]. The CRISP3 gene is widely distributed in tissues with an exocrine function [73] and is suggested to play an role in the pathogenesis of Sjögren’s syndrome [74,75] and prostate cancer [76-79].

To further infer the potential mechanisms of PTCHD2, CRISP3, NAP1L4, FSCB, and AP3B2 gene underlying refractive errors, we identified genes that were predicted with strong evidence (confidence score >0.70) to be functional partners of these five genes from the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database [80]. Due to the limits of the exon array, not all pathway members were included among the 11,571 genes in our gene-based analysis. For those tested, we found no association at the significance level of 0.01 (Appendix 6), suggesting that these five genes may contribute to the development of refractive error either independently or through other functional partners that have yet to be identified. As refractive errors often develop due to incorrect morphology (i.e., incorrect axial length), further studies are necessary to understand how variation in these genes may influence ocular growth and development.

Our results contribute to an increasing body of evidence that variants in the GJD2 gene are associated with spherical equivalent. Multiple SNPs in vicinity of the GJD2 gene on chromosome 15q14 have been reported to be associated with refractive error and myopia. rs634990 was first discovered as a susceptibility locus for spherical equivalent in a Dutch population-based study (MAF = 47%, β = –0.25, p = 1.03 × 10−7) [31]. This association was later replicated by the Blue Mountain Eye Study (BMES) [81] and the Consortium for Refractive Error and Myopia (CREAM) [82], as well as in this analysis of the BDES cohort. This variant, rs634990, is also associated with myopia. When we compared the myopes to hyperopes, as was reported by Verhoeven et al. (OR = 1.88 [1.64, 2.16], p<0.001 for homozygous carriers; OR = 1.33 [1.19. 1.49], p<0.001 for the heterozygous carriers) [82], our findings were consistent with this previous study, confirming the association of rs634990 with myopia and indicated a protective effect of its C allele on hyperopia. We did not have data on rs524952, a locus that is strongly correlated with rs634990 (r2 = 0.90) and associated with spherical equivalent [30,81,82] and myopia [34,83] in European populations. Associations of GJD2 variants with refractive errors and myopia have also been replicated in studies conducted in Asian populations [84,85]. GJD2 encodes a member of the connexin family of proteins (hCx36) that form hexameric intercellular gap junction channels. Expression of the GJD2 gene has been found in the inner and outer plexiform layers of the retina in mice, rats, and humans [86]. Although the function of GJD2 in the human retina remains unknown, animal studies suggest an essential role for gap junction channels containing Cx36 (counterpart protein of hCx36 in mice) in the transmission of visual signals through photoreceptor cells, particularly the rod-cone pathways [87]. Disruption of Cx36 expression in mice can cause visual transmission deficits [88]. However, further investigation is necessary to understand how variation in the GJD2 gene influences spherical equivalent.

Our results support the potential involvement of PRSS56 in spherical equivalent. Associations with spherical equivalent or myopia have been found for two SNPs (rs1656404 and rs1550094) in the PRSS56 gene region. We did not have data on rs1656404, a locus that was first reported for spherical equivalent (MAF = 21.0%, β = –0.151, p = 2.38 × 10−9) in a multiancestry cohort [30] and later for myopia [83]. However, rs1550094 was discovered for myopia development (i.e., age of onset) in the GWAS of 23andMe participants of European ancestry (HR = 1.09 [1.07–1.11], MAF = 30.5%, p = 1.3 × 10−15) [34] and is associated with spherical equivalent and myopia in the present study. Due to the low correlation (r2<0.50) between these two SNPs, rs1550094 and rs1656404 were independently associated with a decreased spherical equivalent. Experimental evidence also supports a potential role for PRSS56 in controlling spherical equivalent. The PRSS56 gene encodes a trypsin-like serine peptidase expressed in the human retina, cornea, sclera, and optic nerve [89]. Mutations of the PRSS56 gene cause autosomal recessive posterior microphthalmos and familial nanophthalmos, two rare ocular disorders characterized by abnormally small eyes and severe hyperopia due to reduced axial length [89,90]. In animal studies, mutations of the Mfrp gene caused a progressive accumulation of Prss56 throughout postnatal development and together led to shortening of photoreceptor outer segments and retinal degeneration [91].

We employed an extreme value sampling strategy to maximize the phenotypic differences among individuals. This could lead to the oversampling of variants that are more common among individuals with extreme trait values and thus increase the statistical power of our study. In this study, we directly examined coding variants for association with spherical equivalent. We applied single-variant analysis and gene-based analysis to optimize the power to detect associations for low-frequency and rare variants. However, the exome array was designed using exome sequencing data from ethnically diverse populations. Given that our study population is ethnically homogenous, more than 50% of the genotyped variants were monomorphic, and exclusion of these variants resulted in the mean coverage of six variants in each gene region. This low coverage can be particularly influential when variants in a gene region are mostly singletons or doubletons, resulting in an inflated association driven by a limited number of extreme observations. Therefore, we are cautious with the interpretation of such genes in our gene-based analysis. In addition, variants in sequencing data have to reach certain frequencies to be selected into the exome array. Thus, there is a lack of coverage in the exome array of extremely rare variants, which prohibited us from examining these variants for association with quantitative refraction.

Additional studies focused on assessing the role of low-frequency coding variants across diverse populations are necessary to further validate the reported findings. Such studies would complement the existing GWAS that highlight common variations, and together, these approaches may begin to elucidate the causal alleles and the underlying mechanisms for refractive errors. A better understating of the genetics of spherical equivalent, in combination with the interaction of environmental exposure, may be useful in determining individual risks for myopia and hyperopia and lead to novel or improved measures of prevention, management, and treatment.