Myopia, or nearsightedness, has rapidly emerged as a global health concern in the last three decades[@b1]. It is one of the leading causes of visual impairment and is associated with potentially blinding ocular complications including myopic maculopathy, glaucoma, cataract and retinal detachment[@b2]. Evidence from family and twin studies strongly supports the heritability of myopia[@b3]. Estimates for the heritability of the quantitative trait refractive error have been reported to be as high as 90% (ref. [@b4]). On the other hand, the rapid upsurge of myopia in the last few decades in many parts of the world is likely to be a consequence of lifestyle changes such as the increasing educational intensity, in particular in urban East Asia[@b5][@b6], and potentially gene and environment (G × E) interactions.

Major attempts undertaken in genome-wide association studies (GWAS) to elucidate the genetic determination of myopia and refractive error have recently led to the discovery of \>30 distinct susceptibility loci[@b7][@b8]. Nevertheless, collectively these genetic variants are estimated to explain \<12% of phenotypic variance in refractive errror[@b7][@b8]. As myopia is a result of the combination of genetic and environmental factors, interplay between genes and environment may account for a substantial proportion of the phenotypic variance. In recent times, we showed interactions between education and genetic risk score of myopia derived from 26 known GWAS single-nucleotide polymorphisms (SNPs) in the Rotterdam Study[@b9]; the combined effect of genetic predisposition and education on the risk of myopia was substantially greater than anticipated from a simple sum of these two factors. At the gene level, some genes such as *SHISA6-DNAH9* have been shown to interact with education level and exhibit strong genetic effects for myopia among Asians with at least higher secondary education[@b10]. In the current study, we demonstrate that new genetic effects implicated in myopia development could be uncovered by studying interactions between genetic variants and education level.

In the context of the aetiology of refractive errors, education attainment is generally considered a surrogate measure for accumulated near work activity[@b1]. When viewing near objects, the eye generates extra optical power through the process of accommodation to focus the image on the retinal plane, to maintain clear vision[@b11]. There is an accommodative lag (less accommodation produced than needed) in many myopes, resulting in a hyperopic defocus on the retina for near work, which has long been proposed to promote eye growth[@b1][@b12], but whether this occurs before or after the onset of myopia in humans is less clear. The retina has a central role in the mechanism linking such visual input with eye growth and refractive development[@b13]. Several neurotransmitters or molecules have been implicated in this process by animal studies including dopamine, acetylcholine, vasoactive intestinal peptide, GABA (γ-aminobutyric acid) and glucagon[@b14][@b15]. However, an organized framework for the retinal signalling mechanisms underlying refractive error development under various environmental conditions remains to be elucidated.

Factoring in environmental exposures may enhance power for the detection of genes, especially in circumstances where a genetic locus has a differential effect conditional on specific environment exposures[@b16]. Gene--environment-wide interaction studies (GEWIS) using a joint meta-analysis (JMA) approach on SNP main effects and SNP × environment interactions have recently been described[@b17][@b18]. This approach has successfully identified six novel loci associated with fasting insulin and glucose accounting for interactions with body mass index[@b18]. It also led to the identification of two novel loci for pulmonary function that did not emerge from analyses based on the genetic main effects alone[@b19]. The well-documented effects of educational attainment on myopia and refractive error make the proposed interaction an excellent analytical candidate for the GEWIS.

The availability of large-scale GWAS spherical equivalent data sets from the Consortium for Refractive Error And Myopia (CREAM) makes G × E interaction analyses feasible. To identify additional genetic variants for refractive error, we performed GEWIS-based analyses on 40,036 adults of European ancestry from 25 studies and 10,315 adults of Asian ancestry from 9 studies. We identified nine new loci using the JMA approach, where three loci exhibited G × E interaction on refractive error in Asians, including the GABA~C~ receptor subunit ρ1 gene *GABRR1*.

Results
=======

Educational level and its main effects on spherical equivalent
--------------------------------------------------------------

The baseline characteristics of 50,351 participants from 34 studies in our meta-analysis are shown in [Table 1](#t1){ref-type="table"}. A total of 40,036 participants were of European descent and 10,315 were of Asian descent; the age of the participants ranged from 20 to 99 years. We grouped individuals into two educational categories: a higher education group that included individuals who completed higher secondary or university education and a lower education group comprising those with lower secondary education or below (see Methods). Among Europeans, the proportions of participants in the higher education group ranged from 16.5% (FITSA[@b20] and OGP Talana[@b21]) to 94.4% (AREDS[@b22]) with an average of 50.7% ([Supplementary Table 1](#S1){ref-type="supplementary-material"}). In Asians, the proportions of individuals in the higher education group ranged from 6.7% (SiMES[@b23]) to 75.9% (Nagahama[@b24]) with an average of 30.0%. Across all studies, individuals in the higher education group had a spherical equivalent refractive error that was on average 0.59 dioptres (D) more myopic, or less hyperopic, compared with those in the lower education group (*β*=−0.59 D; 95% confidence interval (CI): −0.64 to −0.55 D). High education level was associated with a twofold more myopic spherical equivalent in individuals of Asian as compared with European ancestry (Asians: *β*=−1.09 D, 95% CI: −1.20 to −0.98 D; Europeans: *β*=−0.49 D, 95% CI: −0.54 to −0.44 D; [Fig. 1](#f1){ref-type="fig"}). Among Asian studies, we also observed heterogeneity of education effects for refractive error. The education effects on spherical equivalent in Singapore Chinese were significantly larger than that in other Asian studies (Singapore Chinese: *β*=−1.75 D, 95% CI: −1.92 to −1.58 D; other Asian cohorts: *β*=−0.60 D, 95% CI: −0.75 to −0.46 D).

GEWIS in Europeans
------------------

After stringent quality control (QC) filtering, ∼6 million SNPs in each study were eligible for the genome-wide JMA test ([Supplementary Table 2](#S1){ref-type="supplementary-material"}). The JMA for SNP main effects and SNP × education interactions in 40,036 European Ancestry individuals showed an association with spherical equivalent at 12 previously implicated loci ([Fig. 2a](#f2){ref-type="fig"}, [Supplementary Table 3](#S1){ref-type="supplementary-material"} and [Supplementary Fig 1](#S1){ref-type="supplementary-material"}). We also identified four previously unreported loci associated with spherical equivalent achieving genome-wide significance (*P*~JMA~\<5.0 × 10^−8^; *P*~het~≥0.086; [Table 2](#t2){ref-type="table"}): *FAM150B-ACP1*, *LINC00340*, *FBN1* and *DIS3L-MAP2K1*. The significant association for JMA testing at these loci in Europeans was primarily driven by SNP effects in both the lower and higher education strata (4.40 × 10^−8^≤*P*~main~≤1.35 × 10^−6^ and 7.61 × 10^−11^≤*P*~main~≤1.75 × 10^−6^, respectively). SNP × education interaction was not significant (*P*~int~≥0.208). The estimated effect sizes of SNP effects on spherical equivalent were highly similar across education strata.

GEWIS in Asians
---------------

The JMA for spherical equivalent in 10,315 individuals from the Asians cohorts identified genome-wide significant association for three genes: *AREG*, *GABRR1* and *PDE10A* (*P*~JMA~\<5.0 × 10^−8^; [Table 3](#t3){ref-type="table"} and [Fig. 2b](#f2){ref-type="fig"}). SNP × education interaction effects associated with spherical equivalent were observed at all three loci, with genetic effects significantly larger within participants who had a higher level of education compared with those with a lower education level: *AREG* (rs12511037, *β*~int~=−0.89±0.14 D, *P*~int~=6.87 × 10^−11^), *GABRR1* (rs13215566, *β*~int~=−0.56±0.14 D, *P*~int~=8.48 × 10^−5^) and *PDE10A* (rs12206610, *β*~int~=−0.72±0.13 D, *P*~int~=2.32 × 10^−8^). The genotype and phenotype associations were highly significant in the higher education stratum (main genetic effects, 1.97 × 10^−10^≤*P*~main~≤8.16 × 10^−8^) but were considerably weaker in the lower education stratum (0.008≤*P*~main~≤0.243). There was no evidence of inter-study heterogeneity at index SNPs within *AREG*, *GABRR1* or *PDE10A* loci (*Q*-test: *P*~het~≥0.122).

*GABRR1* and *PDE10A* index SNPs were not associated with spherical equivalent in European samples, for either the JMA test, SNP main effect or SNP × education interaction ([Table 3](#t3){ref-type="table"}). *AREG* SNP rs12511037 was excluded in the meta-analysis of European studies after QC filtering; hence, a proxy SNP, rs1246413, in linkage disequilibrium (LD) with rs12511037 in Asians (*r*^2^=0.97) was tested but not associated with spherical equivalent (*P*~JMA~=0.527; *P*~int~=0.176). The meta-regression including study-level characteristics as covariates in the model confirmed the heterogeneity between populations of European and Asian ancestry (*GABRR1*: *P*=0.006; *PDE10A*: *P*=0.0419; [Supplementary Table 4](#S1){ref-type="supplementary-material"}). For *PDE10A*, besides ethnicity, average spherical equivalent of each study also explained the inter-study heterogeneity for the interaction effects (*P*=0.025).

We examined whether the underlying assumption of G × E independence held at these three G × E interaction loci. We performed a meta-analysis of logistic regression analysis for education level on *AREG* SNP rs12511037, *GABRR1* SNP rs13215566 and *PDE10A* SNP rs12296610, adjusting for age, gender and population stratification in Asian cohorts (*n*=10,315). Our analysis did not reveal any significant associations between these loci and education level (*P*≥0.102, *P*~het~≥0.170; [Supplementary Table 5](#S1){ref-type="supplementary-material"}). Furthermore, the three loci were not associated with educational attainment in a large meta-analysis of GWAS recently conducted in European cohorts[@b25]. Thus, our G × E results are unlikely to be biased due to dependence between gene and education.

We also evaluated the association for spherical equivalent in Asian cohorts for four loci identified from European populations. Two of them showed significant associations in the JMA test: *FAM150B-ACP1* (*P*~JMA~=0.031) and *DIS3L-MAP2K1* (*P*~JMA~=0.0042; [Table 2](#t2){ref-type="table"}). The SNP effect sizes in lower and higher education strata in Asians were similar at *FAM150B-ACP1*. The signal at the *DIS3L-MAP2K1* locus was mainly driven by SNP × education interaction in Asians (*P*~int~=7.95 × 10^−4^), whereas the interaction effect was not statistically significant in Europeans (*P*~int~=0.208).

Combined GEWIS of all cohorts
-----------------------------

We subsequently conducted a JMA in the combined data including both the European and Asian participants of all 34 studies. This analysis revealed two additional SNPs: *ARID2-SNAT1* (*P*~JMA~=4.38 × 10^−8^) and *SLC14A2* (*P*~JMA~=2.54 × 10^−8^). Both loci showed suggestive association with spherical equivalent in European cohorts, with the same direction of effect and similar effect sizes in Asian cohorts ([Table 2](#t2){ref-type="table"}). We also detected genome-wide significant associations with spherical equivalent for 17 known loci[@b8] identified in our previous CREAM GWAS ([Supplementary Table 3](#S1){ref-type="supplementary-material"}). The regional plots of the identified novel loci are presented in [Supplementary Fig 2](#S1){ref-type="supplementary-material"}.

Gene and education interactions for GWAS known loci
---------------------------------------------------

We also evaluated the interactions between education and previously reported genetic association with spherical equivalent at 39 loci identified from recent two large GWAS studies[@b7][@b8]. Two SNP × education interactions were nominally significant ([Supplementary Table 6](#S1){ref-type="supplementary-material"}): *TJP2* in Europeans (rs11145488, *P*~int~=6.91 × 10^−3^) and *SHISA6-DNAH9* in Asians (rs2969180, *P*~int~=4.02 × 10^−3^). In general, the index SNPs tested at 39 loci had larger SNP × education interaction effect on spherical equivalent in Asians versus Europeans (meta-regression *P* for fold changes\<0.001; [Supplementary Fig. 3](#S1){ref-type="supplementary-material"}). For 20 SNPs with the same direction of the interaction effect, the magnitudes of interaction effects were fourfold larger on average in Asians than in Europeans (*P*=0.003).

Gene and near work interactions for three identified loci
---------------------------------------------------------

High education levels may reflect an estimator for the greater accumulative effect of near work[@b26][@b27]. We thus examined whether there was evidence for SNP × near work interactions associated with spherical equivalent at the three loci (*AREG*, *GABRR1* and *PDE10A*) in paediatric cohorts (SCORM[@b28], Guangzhou Twins[@b29] and ALSPAC[@b30]; combined *n*=5,835; [Supplementary Table 7](#S1){ref-type="supplementary-material"}). Tentative support for a SNP × near work interaction was observed for *PDE10A* (rs12206610, *P*~int~=0.032, *P*~het~=0.658), with the stronger genetic effect in children spending more hours on reading, writing or compute use. Weaker support for an interaction was noted at *GABRR1* (rs13215566, *P*~int~=0.309, *P*~het~=0.655), although the direction of meta-analysed interaction effect was largely consistent across paediatric studies with that observed in adults. We did not observe the interaction at *AREG* (rs12511037, *P*~int~=0.795, *P*~het~=0.062).

Gene expression in human tissues
--------------------------------

Using the Ocular Tissue Database[@b31], we examined the expression of the associated genes in 20 normal human donor eyes. The majority of genes identified were expressed in human retina, sclera, choroid or retinal pigment epithelium (RPE) ([Supplementary Table 8](#S1){ref-type="supplementary-material"}). Among these genes, *GABRR1*, *ACP1* and *SNAT1* had the highest expression in the retina. The Probe Logarithmic Intensity Error-normalized messenger RNA expression levels in the retina ranged from 121.66 to 236.69. Of note, *MAP2K1* was widely expressed in the retina, sclera and choroid/RPE.

Discussion
==========

This study represents the most comprehensive genome-wide scan of gene and education interactions to date, for refractive error. Here we identified novel genetic loci associated with refractive error by testing the joint contribution of SNP main effects and SNP × education effects in large multi-ethnic populations. Three loci (*AREG*, *GABRR1* and *PDE10A*) showed strong interactions with education in populations of Asian descent, with larger genetic effects within participants who had a higher level of education compared with those with a lower education level; no interactions achieved statistical significance in Europeans for top JMA associations or known myopic loci. Apart from confirming known associations at 17 previous published loci, we identified six new loci (*FAM150B-ACP1*, *LINC00340*, *FBN1*, *DIS3L-MAP2K1*, *ARID2-SNAT1* and *SLC14A2*) significantly associated with spherical equivalent using the combined multi-racial cohort.

A recent meta-analysis of GWAS in multi-ethnic populations comprises 32 studies (*n*=45,756) from CREAM and a large GWAS in Europeans (*n*=45,771) have reported a total of 39 genetic loci associated with refractive phenotypes[@b7][@b8]. The current genome-wide meta-analysis included ∼5,000 more subjects than the previous GWAS of main effects. We identified nine additional novel loci using the JMA approach. These loci can be placed within the biological context of the visually evoked signalling cascade that begins in the retina and mediates sclera remodelling[@b32]. The newly identified genes are involved in retinal neurotransmission (*GABRR1* and *SNAT1*), extracellular matrix remodelling (*FBN1*, *MAP2K1* and *AREG*), circadian rhythm (*PDE10A*) and platelet-derived growth factor receptor signalling (*ACP1*) ([Supplementary Table 9](#S1){ref-type="supplementary-material"}). Network analysis revealed that most of the novel genes may tend to be co-expressed and co-localized with the known myopia susceptibility genes through multiple biological networks such as *LAMA2*, *GJD2*, *RASGRF1*, *BMP3*, *RDH5*, *ZMAT4*, *RBFOX1*, *RDH5* and so on ([Supplementary Fig. 4](#S1){ref-type="supplementary-material"}). Our data, in line with previous findings, substantiate the assumption of heterogeneity in the molecular mechanisms involved in refractive error and myopia.

Among the novel loci, *GABRR1* on chromosome 6q15 (53 kb), encoding GABA~C~ receptor subunit ρ1, is an interesting functional candidate suggestive of a role in myopia development. Modulation of synaptic plasticity via GABA-mediated inhibition would be well placed to alter the 'gain\' of the visually guided feedback system controlling refractive development[@b33]. The lead SNP rs13215566 in *GABRR1*, together with seven SNPs within the LD block (*r*^2^≥0.8), are intronic, potentially affecting regulatory motifs (such as zfp128 and gcm1) that may influence transcriptional regulation ([Supplementary Table 10](#S1){ref-type="supplementary-material"}). The variant rs13215029, in LD (*r*^2^=1) with rs13215566, is associated with *cis*-acting expression of *GABRR1* (*P*=2.3 × 10^−4^) in skin tissues ([Supplementary Table 9](#S1){ref-type="supplementary-material"})[@b34]. A recent pharmacological study provided evidence that retinal dopaminergic and GABAergic neurotransmitters play a substantial role in the modulation of refractive development in form-deprivation myopia[@b35][@b36][@b37]. Stone *et al*.[@b35] have reported that antagonists to GABA-A, -B and -C receptors inhibited form-deprivation myopia in chicks, with greatest effect in the equatorial dimension. GABA receptors, expressed in bipolar and ganglion neuron cells, also interact with dopamine pathways in the retina[@b36]. A recent proteomics study determined that levels of GABA transporter-1 were significantly reduced in myopic murine retina after atropine treatment, implying that GABA signalling is involved in the anti-myopic effects of atropine[@b37]. Altogether, these data suggest that GABA~C~ receptor ρ1 may regulate the development of myopia through functional feedback from RPE to neuron cells in the retina. Further studies are needed to investigate the effect of genetic deletion of *GABRR1* on refractive eye development and the role of the GABAergic pathway in myopia development using gene knockout mice. Therefore, our result in humans is in line with animal experiments, supporting the notion that the GABAergic neurotransmitter signalling pathway in the retina could be a potential factor in the progression of myopia.

SNP rs10889855 on chromosome 6 is an intronic variant within the *ARID2* gene (AT-rich interactive domain 2) and 500 kb downstream of *SNAT1* (solute carrier family 38, member, alias *SLC38A1*). SNAT1 is a transporter of glutamine, a precursor of GABA[@b38]. It is also highly expressed in human retina. In our previous meta-analysis in CREAM[@b8], we identified variants in another glutamate receptor gene *GRIA4* (encoding glutamate receptor, ionotropic); altogether, current evidence supports the notion that retinal neurotransmitters GABA and glutamine may be involved in the refractive development.

The strongest association signal for gene and environment interactions was from rs12511037, located 14 kb downstream the *AREG* gene (amphiregulin). AREG is a ligand of the epidermal growth factor receptor promoting the growth of normal epithelial cells, which is critical for cell differentiation and proliferation such as regrowth of the wounded cornea[@b39]. A link has been found between the muscarinic acetylcholine receptors, dominant in myopia progression, and the epidermal growth factor receptor in muscarinic system[@b40][@b41].

Another novel association, rs16949788 on chromosome 15, derives from a region that spans *DIS3L* and *MAP2K1*. *MAP2K1* encodes mitogen-activated protein kinase 1, which binds to muscarinic receptors during proliferation[@b42] and inhibits the proliferation of human scleral fibroblasts exposed to all-*trans* retinoic acid[@b43]. All-*trans* retinoic acid is a modulator of ocular growth, inhibiting the proliferation of human scleral fibroblasts[@b44].

*FBN1* (Fibrillin 1), a member of the fibrillin family, encodes a large extracellular matrix glycoprotein. Mutations in *FBN1* cause Marfan\'s syndrome, a disorder of connective tissue affecting the ocular, skeletal and cardiovascular systems[@b45]. As a candidate gene for myopia, attempts to study its association with myopia previously produced inconclusive results[@b46][@b47], probably owing, in part, to underpowered studies with insufficient sample sizes. Using data from a large multi-ethnic population, our results support the role of *FBN1* in myopia development.

The genome-wide significant SNPs from the JMA approach did not exhibit any interactions with education in Europeans, in contrast to the significant interactive effect among Asians. In particular, the G × E interactions at *AREG*, *GABRR1* and *PDE10A* were only evident in Asian populations. There are a number of possible reasons for the observed differences. First, the variation of LD patterns and joint effects of genetic variants might affect the transferability of G × E signals across populations. Similar LD patterns were seen at *GABRR1* and *PDE10A* regions across populations whereas a long stretch of LD flanking AREG was present only in Asian populations ([Supplementary Fig. 5](#S1){ref-type="supplementary-material"}). As the true causal variants transferrable across populations are probably not implicated in our study, the identified novel myopia risk loci provide a much-needed starting point for follow-up and functional downstream analyses. Second, we used education level in adults as a surrogate measure of the underlying risk factors for influencing refractive development. Ideally, near work intensity would be measured prospectively in children before the onset of myopia. Studies included in our analyses do not have additional data on relevant childhood exposures. Thus, the best available surrogate measure of cumulative near work exposure in adult cohorts is educational level. We nevertheless believe that education is a highly reliable proxy for the relevant exposures underlying refractive development and is universally associated with refractive error in our study. Third, the differences of G × E interactions in Asian versus Europeans may reflect quantitative differences in near work intensity during childhood. For example, 6- and 7- year-old children in England and Australia reported less near work activity outside of school (1.0--2.3 h per day)[@b48][@b49] compared with children in Singapore and China (2.7--3.5 h per day)[@b50][@b51]. A similar trend was observed in older children[@b48][@b52][@b53][@b54][@b55]. We thus speculate that the total exposure to near work activity may be greater in East Asians compared with European-derived populations with the same levels of education; hence, G × E interaction estimates would tend to be inflated in Asian populations compared with European groups. Fourth, other environmental factors such as outdoor activities could also interact with genes. East Asian children tend to have less exposure to outdoor activities compared with their European peers[@b56]. However, the majority of adult cohorts did not report time outdoors and thus could not be accounted for in the current study. Finally, the population mean of refractive error is less myopic in Europeans (0.10 D) versus Asians (−0.60 D). Of note, for the previously known myopia loci, the magnitudes of interaction effects were fourfold larger on average in Asians than in Europeans ([Supplementary Fig. 2](#S1){ref-type="supplementary-material"}). The impact of G × E interactions may be seen at certain severity levels of myopia.

The risk alleles of rs12511037 in *AREG*, rs1321556 in *GABRR1* and rs12206610 in *PDE10A* had no or weak influence on myopic shift in the lower education group compared with the higher education group. This suggests that the hereditary predisposition to myopia could be latent for the risk allele carriers, if they are less exposed to the myopiagenic environment associated with high-level education. A lack of strong SNP × near work associations in children might be due to the inadequate statistical power in paediatric cohorts of relatively small sample sizes, or the possibility that environmental risk exposures other than near work might underlie the SNP × education interaction seen in the adult Asian samples.

In summary, we identified nine novel loci associated with refractive error in a large multi-ethnic cohort study by GEWIS approach. Our data provide evidence that specific genetic variants interact with education, to influence refractive development, and further support a role for GABA neurotransmitter signalling in myopia development. These findings provide promising candidate genes for follow-up work and may lead to new genetic targets for therapeutic interventions on myopia.

Methods
=======

Study participants
------------------

Thirty-four studies from members of CREAM, comprising 40,036 individuals of European ancestry from 25 studies and 10,315 individuals of Asian ancestry from 9 studies, were made available for this analysis ([Table 1](#t1){ref-type="table"} and [Supplementary Note 1](#S1){ref-type="supplementary-material"}). Individuals aged \<20 years were excluded and so were those who had undergone cataract surgery, laser or other intra-ocular procedures that could alter refraction. Many of these studies were also included in the previous CREAM GWAS on spherical equivalent[@b8]. All studies adhered to the tenets of the Declaration of Helsinki and were approved by their local research ethics committees. The exact names of the Institutional Research Board committees can be found under [Supplementary Note 2](#S1){ref-type="supplementary-material"}. All participants provided a signed consent form before the start of the study.

Phenotyping and education levels
--------------------------------

All participants underwent ophthalmological examinations ([Supplementary Table 1](#S1){ref-type="supplementary-material"}). Non-cycloplegic refraction was measured by autorefraction and/or subjective refraction. Spherical equivalent was calculated as the sphere power plus half of the cylinder power for each eye. The mean spherical equivalent of the right and left eyes was used as a quantitative outcome. When data from only one eye was available, the spherical equivalent of that eye was used. For each study, the participants reported the highest level of education achieved or the years of schooling through a self-reported questionnaire, or in an interview.

We dichotomized education for all participants in each study. The higher education group consisted of those who had achieved the highest educational level of A-levels, high school (higher secondary education), vocational training (for example, diploma), university degree or those with ≥12 years spent in formal education (beginning from first grade). Those who had achieved the highest educational level of O-level, middle school (lower secondary education) or those with \<12 years of formal education were classified into the lower education group. If both number of formal study years and education levels were available in the cohort, we classified participants based on years of formal education. For the four cohorts of relatively young European participants (BATS, DCCT, RAINE and WESDR; total *n*=2,349), almost all of them had completed 12 or more years of schooling. We thus chose to categorize individuals with tertiary or university education as the higher education group in these studies. Sensitivity analysis excluding these four cohorts did not appreciably change our meta-analysis results.

Genotyping and imputation
-------------------------

Detailed information on the genotyping platforms and QC procedures for each study is provided in [Supplementary Table 2](#S1){ref-type="supplementary-material"} and [Supplementary Note 1](#S1){ref-type="supplementary-material"}. Each study applied stringent QC filters for GWAS. In general, duplicate DNA samples, individuals with low call rate (\<95%), gender mismatch or ethnic outliers were excluded. SNPs were excluded if low genotyping call rate (\>5% missingness), monomorphic SNPs, with minor allele frequency (MAF) \<1% or in Hardy--Weinberg disequilibrium (*P*\<10^−6^). After QC filtering, the array genotypes of each study were imputed using the 1000 Genomes Project data as reference panels (build 37, phase 1 release, March 2012) with the software Minimac[@b57] or IMPUTE2 (ref. [@b58]). Approximately six million SNPs that passed imputation quality thresholds (MACH: *r*^2^\>0.5 or IMPUTE info score \>0.5) and with MAF ≥5% were eligible for the meta-analysis ([Supplementary Table 2](#S1){ref-type="supplementary-material"}).

Statistical models
------------------

For each study, a linear regression model for each genotyped or imputed SNP was constructed with the mean spherical equivalent as the outcome. We assumed an additive genetic model where the number of risk alleles is an ordinal variable (0, 1 and 2) for directly genotyped SNPs or a continuous variable of allele dosage probability ranging from 0 to 2 for imputed SNPs. The primary analytic model included SNP, education and SNP × education interaction term, as well as age and sex as covariates. Additional adjustments for the top principal components of genomic marker variations were performed in individual studies when applicable (that is, when there was evidence of population stratification).

We used the following additive genetic model to test for a joint effect of SNP (*β*~SNP~) and SNP × education interaction (*β*~SNP\ ×\ education~) on mean spherical equivalent:

where *Y* is the mean spherical equivalent and education is a dichotomous variable (0=lower education group and 1=higher education group); *cov* is a set of covariates such as age, sex and first top five principal components when applicable. For family-based studies, the kinship matrix was estimated empirically from the SNP data and included as a random effect in the generalized mixed model[@b59]. To test an effect of SNP × education interaction, we assessed *β*~SNP\ ×\ education~ from [equation (1)](#eq1){ref-type="disp-formula"}.

The linear regression analyses in each study were conducted with Quickes or ProbABEL for the unrelated samples and MixABEL for family-based data (see URLs). The command 'robust\' was used in the above software to calculate the robust ('sandwich\', Huber-White) s.e. of *β*~SNP~ and *β*~SNP\ ×\ education~, and error covariance of *β*, to correct the potential inflation of false positive rate for the interaction *P*-value[@b60].

In addition, each study also tested the main effect of education on spherical equivalent by adjusting for age and gender using the linear regression model:

where the definition of each variable is the same as in [equation (1)](#eq1){ref-type="disp-formula"}. We performed meta-analysis of the education effects on mean spherical equivalent in Europeans, Asians (Singapore Chinese versus others) and combined data using a fixed-effect model with inverse-variance weighting (R package 'meta\').

GEWIS join meta-analyses
------------------------

We adopted the JMA approach[@b17][@b61], to simultaneously test both SNP main effects and SNP × education interactions for spherical equivalent with a fixed-effect model, using SNP and SNP × education regression coefficients (*β*~SNP~ and *β*~SNP\ ×\ education~, respectively) and a *β*\'s covariance matrix from each study. A Wald\'s statistic, following a *χ*^2^-distribution with two degrees of freedom, was used to test the joint significance of the *β*~SNP~ and *β*~SNP\ ×\ education~. The JMA was performed with METAL[@b62], as previously described by Manning *et al*.[@b61]. A Cochran\'s *Q*-test was used to assess heterogeneity of the *β*-coefficients across studies for the SNP and interaction effects. To test for interaction between the SNP and education, we conducted a secondary meta-analysis of the SNP × education interaction effects for spherical equivalent (*β*~SNP\ ×\ education~, one degree of freedom), with a fixed-effects model using inverse-variance weighting in METAL; this is a traditional meta-analysis to investigate SNP × education interactions *per se*. Effects and s.e. of the SNP effect on spherical equivalent in the lower education group (*β*~SNP~) and higher education group (*β*~SNP~+*β*~SNP\ ×\ education~) were derived from the JMA output[@b61]. We used the *P*-value of 5 × 10^−8^ as a significant threshold for JMA test. For the SNP and SNP × education effects for the identified top loci, the *P*-value threshold for significance was set at 0.0055=0.05/9 (9 index SNPs underlying analyses).

We performed a meta-regression to explore sources of heterogeneity in our meta-analysis for three loci showing G × E interactions (R package 'metafor\'). Meta-regression included the following study-specific variables as covariates: study sample size, proportion of individuals in the higher education group, average spherical equivalent, education main effects on spherical equivalent (higher education level versus lower), ethnicity (Asian versus European), study design (independent samples versus family-based studies), study year and average age. Meta-regression was also conducted to test the fold changes of the interaction *β*-coefficients in Asians versus Europeans for the 39 known myopia loci.

The study-specific genomic control inflation factors *λ*~gc~ for the joint test for SNP and interaction terms ranged from 1.009 to 1.125 with an average of 1.019 ([Supplementary Table 2](#S1){ref-type="supplementary-material"}), calculated by the ratio of the observed median *χ*^2^ divided by the expected median of the 2df *χ*^2^-distribution (1.382). Genomic control correction was applied to each individual study[@b63]. For studies of small sample sizes (*n*\<500) with *λ*~gc~ \>1.05, we further, before starting the meta-analysis, excluded SNPs showing significant joint *P*-value \<1 × 10^−5^ but neither the SNP nor SNP × education effects supported such an association. Quantile--quantile plots showed only modest inflation of the test statistics in the JMA test (Europeans: *λ*~gc~=1.081; Asians: *λ*~gc~=1.053; Combined: *λ*~gc~=1.092; [Supplementary Fig. 1](#S1){ref-type="supplementary-material"}), similar to previous GEWIS studies with comparable sample sizes[@b18][@b19]. We excluded a small number of markers in the meta-analysis with *P*~het~\<0.0001. The *λ*~gc~ for the SNP × education interaction term in the individual studies ranged from 1.01 to 1.08, indicating little evidence of test statistic inflation on SNP × education effect for each study.

Annotation of genetic variants and gene expression in humans
------------------------------------------------------------

The coordinates and variant identifiers are reported on the NCBI B37 (hg19) genome build and annotated using UCSC Genome Browser[@b64]. We identified variants within each of the LD blocks (*r*^2^≥0.8) in European and Asian populations of the 1000 Genomes Project (100 kb flanking the top SNP at each locus), to apply functional annotations of transcription regulation using HaploReg[@b65] and Encyclopedia of DNA Elements[@b66] data. We also generated funcational association and co-expression network using GeneMANIA[@b67], to determine whether the disease-related genes identified in this study and previous GWAS[@b7][@b8] are functionally connected.

To assess gene expression in human tissues, we examined the Ocular Tissue Database and the EyeSAGE database[@b31][@b68]. The estimated gene and exome-level abundances are available online (see URL). Normalization of gene expression used the Probe Logarithmic Intensity Error method with genomic control-background correction[@b31]. Relationships between genotype and *cis* regulation of gene expression levels ([Supplementary Table 9](#S1){ref-type="supplementary-material"}) were assessed using expression quantitative trait locus associations in multiple human tissues from UK samples[@b34], as well as gene expression profiles obtained from GTExPortal database[@b69].

Additional information
======================

**How to cite this article:** Fan, Q. *et al*. Meta-analysis of gene--environment-wide association scans accounting for education level identifies additional loci for refractive error. *Nat. Commun.* 7:11008 doi: 10.1038/ncomms11008 (2016).

Supplementary Material {#S1}
======================

###### Supplementary Information

Supplementary Figures 1-5, Supplementary Tables 1-10, Supplementary Notes 1-2 and Supplementary References

We gratefully thank all the participants or volunteers who participated in the studies. A full list of funding information and acknowledgements by study can be found in the [Supplementary Note 2](#S1){ref-type="supplementary-material"}.

URLs

Quickest: <http://toby.freeshell.org/software/quicktest.shtml>

ProbABEL: <http://www.genabel.org/packages/ProbABEL>

MixABLE: <http://www.genabel.org/packages/MixABEL>

METAL: [www.sph.umich.edu/csg/abecasis/metal/](http://www.sph.umich.edu/csg/abecasis/metal/)

SCAN: <http://www.scandb.org/newinterface/about.html>

GeneCards: <http://www.genecards.org/>

LocusZoom: <http://csg.sph.umich.edu/locuszoom>

Ocular tissue database: <https://genome.uiowa.edu/otdb>

R: <http://www.rproject.org>

HaploReg v2: <http://www.broadinstitute.org/mammals/haploreg/haploreg.php>

GTExPortal database: <http://www.gtexportal.org/home>

GeneMANIA: <http://www.genemania.org>

**Author contributions** Q.F., V.J.M.V., R.W., C.-Y.C, C.J.H., C.C.W.K. and S.-M.S. conceived the project. Q.F., V.J.M.V., R.W., V.A.B., P.G.H., J.A.G., R.H., V.V., A.P.K., S.M.H., T.L., Y.L., T.H., J.X., C.D.,L.P., J.W.,P.G., C.V., M.M., A.W.H., X.G., J.M., R.P.I., C.L.S., C.B., A.D.P., J.E.B.-W., Y.-Y.T., X.W., X.W., E.T.-H.B., C.E.H.H. and W.L. performed analyses. J.E.H., K.M.W., K.Y., O.P., H.C., I.R., Z.V., J.F.W., P.K.J, G.M., B.S.P., D.M.E., T.-H.S.-A., A.M., A.C.-G., M.B., O.R., M.K., I.S., T.Z., T.M., J.S.R., C.G., M.P., E.M.v.L., N.A., A.G.U., F.R., A.H., J.R.V., Y.X.W., M.K.I., C.S., P.G., V.T., L.Z., R.W.B., R.S., E.-S.T., E.V., E.M., C.-C.K., C.H., R.N.L., P.J.F., B.E.K.K., R.K., H.-S.W., P.M., A.M., T.A., T.L.Y., M.H., O.P., C.M.v.D., J.J.W., C.W., J.B.J.,Y.-Y.T., D.A.M., K.O., N.Y., N.P., T.-Y.W, P.N.B., D.S., C.-Y.C. and C.J.H. were responsible for collecting clinic data and performing genotyping in each study. Q.F., V.J.M.V., R.W., C.C.W.K. and S.-M.S. drafted the paper. P.G.H., V.A.B., J.A.G., V.V., R.H., A.P.K., A.D.P., S.M.H., T.Y., J.J.W., T.-Y.W., P.N.B., C.J.H. and C.-Y.C. critically reviewed the manuscript.

![Forest plot of education main effects on spherical equivalent across studies.\
The *β*-coefficient represents the differences of dioptres in refractive error comparing individuals in higher education group versus lower education group in Europeans (*n*=40,036), Asians (*n*=10,315) and all studies (*n*=50,351). The studies are sorted by effect size of education on spherical equivalent within Europeans and Asians studies, respectively.](ncomms11008-f1){#f1}

![Manhattan plots of −log~10~(*P*) for the JMA on SNP main effects and SNP × education effects on spherical equivalent in (a) European ancestry populations and (b) Asian populations.\
The horizontal red line indicates the genome-wide significance level of *P*~JMA~\<5 × 10^−8^. The horizontal blue line indicates the suggestive significance level of *P*~JMA~\<1 × 10^−5^. Novel loci reaching genome-wide significance are labelled in red and known loci are in grey.](ncomms11008-f2){#f2}

###### Characteristics of study participants.

  Study                        *N*          Study year        Age (s.d.)     Age range  Male (%)   Spherical equivalent
  -------------------------- ------- ------------------------ ------------- ----------- ---------- ----------------------
  *Europeans (*n=*40,036)*                                                                         
   ALIENOR                     509          2006--2008        79.2 (4.1)      73--93    43.2       0.98 (1.98)
   ALSPAC                     1,865         1999--2000        45.9 (4.5)      32--59    0          −0.76 (2.16)
   AREDS                      1,842            1992           68.1 (4.7)      55--81    41.0       0.54 (2.15)
   BATS                        383          1992--2013        24.8 (7.8)      20--67    41.3       −0.67 (1.58)
   BMES                       1,896         1992--2009        66.8 (8.9)      49--94    43.8       0.58 (1.94)
   CROATIA-Korcula             807          2007--2008        56.2 (13.3)     25--94    34.9       −0.13 (1.59)
   CROATIA-Split               787          2008--2009        51.9 (13.0)     25--80    38.6       −1.27 (1.59)
   DCCT                       1,057         1982--1993        35.4 (5.8)      25--49    54.1       −1.47 (1.80)
   EGCUT                       904          2002--2013        56 (17.0)       25--99    38.8       0.33 (3.36)
   EPIC                       1,083         2004--2011        68.8 (7.5)      50--88    43.8       0.34 (2.27)
   ERF                        2,604         2002--2005        48.9 (14.4)     25--87    45.0       0.12 (2.03)
   FES                        2,479   1973--1975/ 1989--1991  54.8 (9.3)      28--84    55.3       0.27 (2.37)
   FITSA                       188          2000--2001        68.5 (3.3)      63--76    0          1.44 (2.08)
   GHS1                       3,178         2007--2008        55.3 (10.9)     35--74    50.4       −0.38 (2.47)
   GHS2                       1,354            2008           54.6 (10.8)     36--74    49.6       −0.39 (2.51)
   KORA                       2326          2004--2006        55.1 (11.8)     35--84    49.4       −0.26 (2.18)
   OGP Talana                  456             2002           52.6 (16.3)     25--89    57.3       −0.20 (0.24)
   ORCADES                    1,124            2009           56.5 (13.2)     29--92    39.1       0.10 (2.07)
   RAINE                       348          2010--2012        20.4 (0.34)     20--22    49.1       0.03 (1.29)
   RS1                        5,702         1991--1993        68.7 (8.7)      55--99    41.0       0.83 (2.55)
   RS2                        2,021         2000--2002        64.3 (7.9)      55--95    46.0       0.48 (2.51)
   RS3                        2,918         2006--2009        56.9 (6.6)      45--86    44.0       −0.28 (2.60)
   TwinsUK                    2,154         1998--2010        53.8 (11.4)     25--84    8.4        −0.96 (2.78)
   WESDR                       561          1979--2007        31.7 (7.0)      25--65    50.3       −1.65 (2.07)
  ** **YFS                    1,490            2011           41.9 (5.0)      34--49    44.6       −1.09 (2.16)
                                                                                                    
  *Asians (*n=*10,315)*                                                                            
   BES                         589          2006--2011        62.1 (8.5)      50--90    34.0       −0.06 (1.86)
   Nagahama                    723          2008--2010        49.2 (15.2)     30--74    33.6       −1.93 (2.46)
   SCES-610K                  1,710         2009--2011        57.5 (7.0)      44--84    51.6       −0.72 (2.69)
   SCES-OmniE                  543          2011--2012        59.3 (8.9)      46--83    51.2       −0.89 (2.74)
   SiMES                      2,256         2004--2006        46.8 (10.2)     40--80    49.1       −0.03 (1.81)
   SINDI                      2,088         2007--2009        55.8 (8.8)      43--84    51.5       0.04 (2.07)
   SP2-1M                      811          1992--1998        46.8 (10.2)     25--80    62.3       −1.80 (2.84)
   SP2-610                     854          1992--1998        48.4(11.3)      25--82    19.6       −1.44 (2.89)
   STARS                       741          2007--2009        38.5 (5.2)      26--58    52.4       −2.80 (2.85)

ALIENOR, antioxydants, lipids essentiels, nutrition et maladies oculaiRes; ALSPAC, avon longitudinal study of parents and children; AREDS, age-related eye disease study; BATS, Brisbane adolescent twins study; BMES, blue mountains eye study; DCCT, diabetes control and complications trial; EGCUT, estonian genome center of the university of Tartu; EPIC, EPIC-Norfolk eye study; ERF, erasmus rucphen family study; FES, Framingham eye study; FITSA, finnish twin study on aging; GHS, Gutenberg health study; KORA, cooperative health research in the region of Augsburg; OGP Talana, ogliastra genetic park, talana study; ORCADES, orkney complex disease study; RAINE, RAINE eye health study; RS, Rotterdam study; TwinsUK, Twins UK study; WESDR, Wisconsin epidemiologic study of diabetic retinopathy; YFS, young finns study; BES, Beijing eye study; SCES, Singapore Chinese eye study; SiMES, Singapore Malay eye study; SINDI, Singapore Indian eye study; SP2, Singapore prospective study program; STARS, strabismus, amblyopia and refractive error study in preschool singaporean children. s.d., standard deviation; age in years; spherical equivalent in dioptres.

Details of each study cohort are described in [Supplementary Note 1](#S1){ref-type="supplementary-material"}.

###### Six genetic loci associated with spherical equivalent from the JMA in the European populations and combined analyses.

  SNP (Chr:BP)               Gene                   Allele      FREQ     Subgroup            Europeans (*n*=40,036)   Asians (*n*=10,315)   All (*n*=50,351)                                                          
  -------------------------- ---------------------- -------- ----------- ------------------ ------------------------ --------------------- ------------------ ------- --------------- ------- ------- --------------- -------
  rs60843830 (2:286756)      *FAM150B-ACP1*         C/G       0.66/0.74  JMA                                             3.71 × 10^−8^           0.086                     0.031       0.980           1.27 × 10^−9^   0.395
                                                                         Lower education             −0.11               4.73 × 10^−8^                         −0.09       0.010               −0.10   1.65 × 10^−9^      
                                                                         Higher education            −0.09               1.75 × 10^−6^                         −0.06       0.509               −0.09   9.83 × 10^−7^      
  rs10946507 (6:22100367)    *LINC00340 (6p22.3)*   A/G       0.47/0.16  JMA                                             3.07 × 10^−8^           0.213                     0.433       0.396           2.24 × 10^−8^   0.249
                                                                         Lower education             −0.08               7.08 × 10^−7^                         −0.04       0.313               −0.08   6.13 × 10^−7^      
                                                                         Higher education            −0.09               1.19 × 10^−8^                         −0.08       0.450               −0.09   1.20 × 10^−8^      
  rs8023401 (15:48703823)    *FBN1*                 G/A       0.87/0.95  JMA                                             1.66 × 10^−9^           0.180                     0.572       0.979           2.85 × 10^−9^   0.495
                                                                         Lower education             −0.15               4.40 × 10^−8^                         −0.06       0.304               −0.13   8.17 × 10^−8^      
                                                                         Higher education            −0.16              7.61 × 10^−11^                         −0.03       0.828               −0.14   2.02 × 10^−9^      
  rs16949788 (15:66590037)   *DIS3L-MAP2K1*         T/C       0.91/0.94  JMA                                             1.34 × 10^−8^           0.721                    0.0042       0.219           2.19 × 10^−8^   0.245
                                                                         Lower education             −0.15               1.35 × 10^−6^                         0.21        0.103               −0.13   4.88 × 10^−6^      
                                                                         Higher education            −0.17               1.89 × 10^−9^                         −0.59       0.014               −0.16   3.90 × 10^−9^      
  rs10880855 (12:46144855)   *ARID2-SNAT1*          T/C       0.51/0.43  JMA                                             7.83 × 10^−7^           0.790                     0.019       0.779           4.38 × 10^−8^   0.867
                                                                         Lower education             −0.09               1.26 × 10^−7^                         −0.06       0.067               −0.09   8.42 × 10^−9^      
                                                                         Higher education            −0.07               1.60 × 10^−5^                         −0.16       0.033               −0.07   3.55 × 10^−6^      
  rs10853531 (18:42824449)   *SLC14A2*              G/A       0.80/0.83  JMA                                             7.82 × 10^−6^           0.052                    0.0023       0.812           2.54 × 10^−8^   0.111
                                                                         Lower education             −0.11               1.27 × 10^−6^                         −0.15   9.01 × 10^−4^           −0.11   3.38 × 10^−9^      
                                                                         Higher education            −0.08               2.12 × 10^−6^                         −0.11       0.288               −0.09   7.14 × 10^−6^      

*β*, *β*-coefficient corresponds to the effect in spherical equivalent (dioptres) for 1 additional copy of the risk allele in the higher or lower education group. FREQ, allele frequency of the risk allele in European/Asian cohorts; JMA, joint meta-analysis on SNP effect and SNP × education interaction effect on spherical equivalent; *P*~het~, *P*-value for the test of heterogeneity at each SNP; SNP, single-nucleotide polymorphism. Allele, risk allele/other allele.

###### Three genetic loci associated with spherical equivalent with a significant SNP × education interaction in Asians and results in European populations.

  SNP (Chr:BP)                                             Gene     Allele     FREQ     Subgroup                    Asians (*n*=10,315)                    Europeans (*n*=40,036)      
  ----------------------------------------------------- ---------- -------- ----------- ------------------ ------- --------------------- ------- -------- ------------------------ -------
  rs12511037[\*](#t3-fn3){ref-type="fn"} (4:75334864)     *AREG*     C/T     0.91/0.95  Lower education     0.07           0.243                  −0.05            0.323               
                                                                                        Higher education    −0.70     1.97 × 10^−10^              −0.03            0.579               
                                                                                        SNP × education     −0.89     6.87 × 10^−11^      0.704    0.02            0.176            0.284
                                                                                        JMA                           5.55 × 10^−10^      0.405                    0.527            0.186
  rs13215566 (6:89918638)                                *GABRR1*    C/G     0.94/0.84  Lower education     −0.13          0.030                  −0.03            0.258               
                                                                                        Higher education    −0.68      1.46 × 10^−8^              −0.01            0.817               
                                                                                        SNP × education     −0.56      8.48 × 10^−5^      0.134   −0.02            0.459            0.457
                                                                                        JMA                            3.81 × 10^−8^      0.122                    0.502            0.630
  rs12206610 6:166016800                                 *PDE10A*    C/T     0.90/0.87  Lower education     0.16           0.008                   0.01            0.759               
                                                                                        Higher education    −0.59      8.16 × 10^−8^               0.01            0.810               
                                                                                        SNP × education     −0.72      2.32 × 10^−8^      0.920   −0.002           0.421            0.111
                                                                                        JMA                            9.21 × 10^−9^      0.902                    0.954            0.305

*β* (higher education/lower education), *β*-coefficient corresponds to the effect in spherical equivalent (dioptres) for 1 additional copy of the effect allele in the higher/lower education group; *β* (SNP × education), *β*-coefficient corresponds to the difference in spherical equivalent (dioptres) for 1 additional copy of the effect allele in the higher versus lower education group; FREQ, allele frequency of the effect allele in Asian/European cohorts; JMA, joint meta-analysis on SNP effect and SNP × education interaction effect on spherical equivalent; LD, linkage disequilibrium; *P*~het~, *P*-value for the test of heterogeneity; SNP, single-nucleotide polymorphism.

*β* and *P*-values for SNP × education interaction were calculated by the meta-analysis of conducting a 1df Wald\'s test of single interaction parameter. Allele is listed as effect allele/other allele.

^\*^SNP rs12511037 was not present in European studies after quality control. Here we present the results of a proxy SNP rs1246413 (*T*/*G*, frequency of risk allele *T*=0.95) in LD with rs12511037 (*r*^2^=0.97).

[^1]: These authors contributed equally to this work.

[^2]: These authors jointly supervised this work.

[^3]: A full list of consortium members appears at the end of the paper.
