|Home | About | Journals | Submit | Contact Us | Français|
Refractive error is the most common eye disorder worldwide, and a prominent cause of blindness. Myopia affects over 30% of Western populations, and up to 80% of Asians. The CREAM consortium conducted genome-wide meta-analyses including 37,382 individuals from 27 studies of European ancestry, and 8,376 from 5 Asian cohorts. We identified 16 new loci for refractive error in subjects of European ancestry, of which 8 were shared with Asians. Combined analysis revealed 8 additional loci. The new loci include genes with functions in neurotransmission (GRIA4), ion channels (KCNQ5), retinoic acid metabolism (RDH5), extracellular matrix remodeling (LAMA2, BMP2), and eye development (SIX6, PRSS56). We also confirmed previously reported associations with GJD2 and RASGRF1. Risk score analysis using associated SNPs showed a tenfold increased risk of myopia for subjects with the highest genetic load. Our results, accumulated across independent multi-ethnic studies, considerably advance understanding of mechanisms involved in refractive error and myopia.
Refractive error is the most important cause of visual impairment in the world1. Myopia, or nearsightedness, in particular is associated with structural changes of the eye, increasing the risk of severe complications such as macular degeneration, retinal detachment, and glaucoma. The prevalence of myopia has been rising dramatically over the past few decades2, and it is estimated that 2.5 billion people will be affected by myopia within a decade3. Although several genetic loci influencing refractive error have been identified4-10, their contribution to phenotypic variance is small, and many more loci are expected to explain its genetic architecture.
Here the Consortium for Refractive Error and Myopia (CREAM) presents results from the largest international genome-wide meta-analysis on refractive error with data from 32 studies from Europe, the United States, Australia, and Asia. The meta-analysis was performed in three stages: as a first step, we investigated genome-wide association study (GWAS) results of 37,382 individuals from 27 populations of European ancestry (Supplementary Note, Supplementary Table 1) using spherical equivalent as a continuous outcome; as a second step, we aimed to test cross-ethnic transferability of the statistical significant associations from the first stage in 8,376 individuals from 5 Asian cohorts (Supplementary Note, Supplementary Table 1). As a third step, we performed a GWAS meta-analysis on the combined populations (total n = 45,758). Subsequently, we examined the influence of associated alleles on the risk of myopia in a genetic risk score analysis, and lastly, we evaluated gene expression in ocular tissues and explored potential mechanisms by which newly found loci may exert their effect on refractive development.
At step 1, we analyzed ~2.5 million autosomal single nucleotide polymorphisms (SNPs) which were obtained through whole-genome imputation of genotypes to HapMap 2. The inflation factors (γGC) of the test statistics in individual studies contributing to the meta-analysis ranged between 0.992 and 1.050, indicating excellent within-study control of population substructure (Supplementary Table 2). The overall lambda was 1.09, consistent with a polygenic inheritance model for refractive error (Q-Q plot, Supplementary Figure S1). We did not perform a lambda correction as Yang et al. have shown that in this situation substantial genomic inflation can be expected, even in the absence of population structure and technical artifacts11. We identified 309 SNPs that exceeding the conventional genome-wide significance threshold of P=5.0×10−8 in the European ancestry sample. These SNPs were clustered in 18 distinct genomic regions across 14 chromosomes (Figure 1, Table 1). At step 2, we investigated the 18 best associated SNPs in the Asian population: ten showed evidence of association (Table 1). The most significant association in both ancestry groups was at a previously identified locus on chromosome 15q14 in the proximity of the GJD2 gene (SNP rs524952; Pcombined=1.44×10−15)4,12. The locus near the RASGRF1 gene was also replicated in the meta-analysis (SNP rs4778879; Pcombined=4.25×10−11)9, the remaining 16 genome-wide significantly associated loci had not previously been reported in association with refractive error. Those loci that were not significant in the smaller sized Asian population mostly had a similar effect size and direction of effect as in the European ancestry sample. At step 3, we identified 8 additional loci which exceeded genome-wide significance in the combined analysis (Table 2). Regional and forest plots of the associated loci are provided in Supplementary Figures 2 and 3.
Genotype distributions of the risk alleles were evaluated in Rotterdam Studies 1-3 (n=9,307). The clinical utility for the prediction of risk of myopia was evaluated by a weighted genetic risk score analysis based on the aggregate of effects (regression coefficients betas) of individual SNPs derived from the meta-analysis, using the middle risk category as a reference. Risk scores ranged from a mean risk score of 1.88 (95% CI 1.86 - 1.89) in the lowest risk score category to 3.63 (95% CI 3.61-3.65) in the highest risk score category. Having the lowest or the highest genetic risk score was associated with an odds ratio of 0.38 (95% CI 0.18-0.77), and an odds ratio of 10.97 (95% CI 3.727-31.251) of myopia, respectively (Figure 2). The predictive value (area under receiver operating characteristic curve, AUC) of myopia versus hyperopia was 0.67 (95% CI 0.65-0.69), a relatively high value for genetic factors in a complex trait13,14. The genetic variants explain 3.4% of the phenotypic variation in refractive error in the Rotterdam Study.
We examined the expression of genes harboring a genetic association signal by measuring levels of RNA in various eye tissues, and found most of these genes expressed in the eye (Supplementary Table 3). The genes PRSS56, LOC100506035, and SHISA6 were not available in the expression data set; all other genes were expressed in the retina. Subsequently, we assessed the areas where our SNP hits reside for H3K27ac modification marks15, and HaploReg16 annotations for marks of active regulatory elements (Supplementary Table 4, Supplementary Figure 4). We found that many hits contain these elements, and alteration of regulatory function is therefore a suggestive mechanism.
The widely-accepted model for myopia development is that eye growth is triggered by a visually-evoked signaling cascade, which originates from the sensory retina, traverses the retinal pigment epithelium and choroid, and terminates in the sclera, where active extracellular matrix (ECM) remodeling results in a relative elongation of the eye17. Many of the genes in or near the identified loci can be linked to biological processes that drive this cascade. Neurotransmission in the retina is a necessary mechanism for eye growth regulation; the most significantly associated gene GJD2 plays a role herein. This gene forms a gap junction between neuronal cells in the retina, enabling intracellular exchange of small molecules and ions. The other previously-reported gene RASGRF1 is a nuclear exchange factor that promotes GDP/GTP exchange on Ras-family GTPases, and is involved in synaptic transmission of photoreceptor responses18,19. Both GJD2 and RASGRF1 knockout mice show retinal photoreception defects18,20. One of the newly identified genes, GRIA4 (SNP rs11601239; Pcombined=5.92×10−9) also has a potential function in this pathway. This gene is a glutamate-gated ion channel that mediates fast synaptic excitatory neurotransmission21, is present in various retinal cells22, has been shown to be critical for light signaling in the retina23 and emmetropization24. Another gene involved in synaptic transmission is A2BP1 (also known as RBFOX1; SNP rs17648524; Pcombined=5.64×10−10, an RNA-binding splicing regulator which modulates membrane excitability25.
We identified for the first time a number of genes involved in ion transport, channel activity and maintenance of membrane potential. KCNQ5, a potassium channel regulator (SNP rs7744813; Pcombined=4.18×10−9), participates in the transport of K+ from the retina to the choroid, and may contribute to voltage-gated K+ channels in photoreceptors and retinal neurons associated with myopia26,27. CD55 (SNP rs1652333; Pcombined=3.05×10−12) is known to elevate cytosolic calcium ion concentration. Other ion channel genes include CACNA1D, a voltage-sensitive calcium channel regulator; KCNJ2, a regulator of potassium ion transport; CHRNG, a nicotinic cholinergic receptor; and MYO1D, a putative binder of calmodulin, which mediates Ca+ sensitivity to KCNQ5 ion channels.
Retinoic acid is synthesized in the retina and highly expressed in the choroid and has been implicated in eye growth in experimental myopia models28-30. Retinol dehydrogenase 5 (RDH5), a novel refractive error susceptibility gene (SNP rs3138144; Pcombined=4.44×10−12), is involved in the recycling of 11-cis-retinal in the visual cycle31. Mutations in RDH5 cause congenital stationary night blindness (OMIM# 136880), a disease associated with myopia. Other genes involved in retinoic acid metabolism are RORB (RAR-related orphan receptor), and CYP26A1, genes that were significant in the European ancestry studies. Notably, retinoic acid contributes to ECM remodeling by regulating cell differentiation.
ECM remodeling of the sclera is the pathological hallmark in myopia development. LAMA2 (laminin α2, SNP rs12205363; Pcombined=1.79×10−12) is the most prominent gene in this respect. LAMA2 forms a subunit of the heterotrimer laminins which are essential components of basement membranes, stabilizing cellular structures and facilitating cell migration32. The two bone morphogenic genes (BMP2, SNP rs235770; Pcombined=1.57×10−8; BMP3) can also be placed within the ECM architecture. They are members of the transforming growth factor-β (TGFβ) super-family, regulate growth and differentiation of mesenchymal cells, and may orchestrate the organization of other connective tissues than bone such as sclera. Remarkably, BMP2 shows bidirectional expression in retinal pigment epithelium in myopia animal models33.
Genes involved in eye development appeared as a separate entity among the gene functions. SIX6 (SNP rs1254319; Pcombined=1.00×10−8 has been linked to anophthalmia and glaucoma34,35, PRSS56 (protease serine 56, SNP rs1656404; Pcombined=7.86×10−11) to microphthalmia36-38, CHD7 to CHARGE syndrome, a congenital condition with severe eye structural defects, and ZIC2 to brain development including visual perception. For the remaining novel gene associations, a mechanism in the pathogenesis of myopia is not immediately clear. Results from Ingenuity and the Protein Link Evaluator39 (Supplementary Figure 5) visualize the subcellular location of all associated genes, and illustrates their interrelationships. Direct connections between genes were surprisingly infrequent, suggesting molecular disease heterogeneity or functional redundancy in the pathobiological events involved in development of refractive error and myopia.
In summary, we identified 24 new chromosomal loci associated with refractive error through a large-scale meta-analysis of GWAS from international multi-ethnic studies. The significant overlap in genetic loci for refractive error between subjects of European ancestry and Asians provides evidence for shared genetic risk factors between the populations. The tenfold increased risk of myopia for those carrying the highest number of risk alleles depicts the clinical significance of our findings. Further elucidation of the mechanisms by which these loci affect eye growth carries the potential to improve the visual outcome of this common trait.
We performed a meta-analysis on directly genotyped and imputed SNPs from individuals of European ancestry in 27 studies, with a total of 37,382 individuals. Subsequently, we evaluated significant SNPs in 8,376 subjects of Asian origin from 5 different studies, and performed a meta-analysis on all studies combined.
All studies participating in this meta-analysis are part of the Consortium for Refractive Error and Myopia (CREAM). All studies had a population-based design and had a similar protocol for phenotyping (Supplementary Table 1). Eligible participants underwent a complete ophthalmologic examination including a non-dilated measurement of refractive error of both eyes. Exclusion criteria were all conditions that could alter refraction, such as cataract surgery, laser refractive procedures, retinal detachment surgery, keratoconus, or ocular or systemic syndromes. Inclusion criteria were persons aged 25 years and over who had data on refractive error and genotype.
The meta-analysis of step 1 was based on 27 studies of European ancestry: 1958 British Birth Cohort, ALSPAC, ANZRAG, AREDS1a1b, AREDS1c, CROATIA-Korcula, CROATIA-Split, CROATIA-Vis, EGCUT, FECD, TEST/BATS, FITSA, Framingham, GHS 1, GHS 2, KORA, ORCADES, TwinsUK, WESDR, YFS, ERF, DCCT, BMES, RS1, RS2, RS3, and OGP Talana. The second step was formed by 5 Asian studies: Beijing Eye Study, SCES, SIMES, SINDI, and SP2.
General methods, demographics and phenotyping and genotyping methods of the study cohorts can be found in the Supplementary Note and Supplementary Table 1. All studies were performed with the approval of their local Medical Ethics Committee, and written informed consent was obtained from all participants in accordance with the Declaration of Helsinki.
Particulars of genotyping in each cohort, particular platforms used to generate genotyping and methods of imputation can be found in more detail in the Supplementary Table 5. To produce consistent datasets and enable meta-analysis of studies across different genotyping platforms, the cohorts performed genomic imputation on the HapMap Phase 2 available genotypes with MACH40 or IMPUTE41, using the appropriate ancestry as templates.
Each cohort applied stringent quality control procedures prior to the imputation, including minor allele frequency cutoffs, Hardy-Weinberg equilibrium (P > 10−7), genotypic success rate (>95%), Mendelian inconsistencies, exclusion of individuals with more than 5% shared ancestry (exception made for family-based cohorts in which due adjustment for family relationship was made) and removal of all individuals whose ancestry as determined through genetic analysis did not match the prevailing ancestry group of the own cohort. SNPs with low imputation quality were filtered using metrics specific to the imputation method and thresholds used in their previous GWAS analyses. Hence, imputation quality criteria varied slightly among studies, and low-confidence imputed SNPs were omitted in the meta-analysis for individual studies.
Spherical equivalent was calculated according to the standard formula (SE=sphere + ½ cylinder), and the mean of two eyes was used for analysis. When data from only one eye was available, the SE of this eye was used.
Each cohort performed association analyses in which the spherical equivalent (determined as described above) was the dependent variable and genotypes (number of alleles in each of the HapMap2 loci) as the independent variables. Analyses in all cases also adjusted for sex and age at the time of phenotype measurement. In family-based cohorts score-test based association test was used to adjust for within-family relatedness (see Supplementary Note)42,43. Study-specific lambda estimates are shown in Supplementary Table 2.
All study effect estimates were corrected using genomic control and were oriented to the positive strand of the NCBI build 36 reference sequence of the human genome, which was the genomic build on which most available genotyping platforms were based. The coordinates and further annotations of the SNPs were further converted into build 37, the most recent of the available builds at the time of writing.
Meta-analyses used effect size estimations (beta regression coefficients) and standard errors from individual cohorts’ summary statistics. Random-effects were assumed for all the meta-analyses which were performed using GWAMA44. We tested for heterogeneous effects between the two ancestries using METAL45 for Linux. For the purpose of these analyses, we defined significance as equal to or better than the conventional multiple testing genome-wide thresholds of association (P<5.0 × 10−8) for stage 1 and nominally significant probabilities (P<0.05) for stage 2. Manhattan, regional plots and forest plots were made using R (see URLs) and Locuszoom (see URLs)46.
For the Rotterdam Study 1-3, a weighted genetic risk score per individual was calculated using the regression coefficients from the GWAS meta-analysis model for the association of SNPs within the associated 26 loci (Table 1, Table 2; per locus only one SNP was included in the analysis) and the individual allele dosages per genotype to evaluate the relationships between myopia (SE ≤ − 3 D), emmetropia (−1.5 D ≤ SE ≤ 1.5 D) and hyperopia (SE ≥ +3 SD). The weighted risk scores were categorized and mean odds ratios per risk score category were calculated for subjects with myopia versus hyperopia, using the middle risk score category as a reference. Subsequently, the area under the receiver curve (AUC) was calculated for myopia versus emmetropia and myopia versus hyperopia. Lastly, the proportion of variance of spherical equivalent explained by the identified SNPs was calculated. For these analyses, we used SPSS version 20.0.0 (SPSS Inc.).
Independently designed, collected, and reported human ocular tissue array data from two different sources, as well as literature reviews were used to verify evidence of expression of the candidate genes.
Human gene expression data of RPE, photoreceptors and choroid were obtained essentially as described47 and the dataset has been deposited in NCBI’s Gene Expression Omnibus48 (GEO series accession number GSE20191; see URLs). In short, postmortem eye bulbs (retinal pigment epithelium was obtained from six donor eyes, choroid was obtained from three donor eyes and photoreceptors were obtained from three donor eyes), provided by the Corneabank Amsterdam, were rapidly frozen using liquid nitrogen. Donors were between 63 and 78 years old and had no known history of eye pathology. Cryosections were cut from the macula, and histology was used to confirm a normal histological appearance. Retinal pigment epithelium, photoreceptor and choroidal cells were isolated from macular sections using a Laser Microdissection System (PALM). Total RNA was isolated and the mRNA component was amplified, labeled and hybridized to a 44K microarray (Agilent Technologies)49. At least three to six microarrays were performed per tissue. Sample isolation, procedures and expression microarray analysis were carried out according to MIAMI guidelines. As a measure of the level of expression, we sorted all the genes represented on the 44K microarray by increasing their expression, and we calculated the corresponding percentiles (Supplementary Table 3a).
We assessed expression of the associated genes in sclera, cornea and optic nerve tissue in an additional dataset (unpublished data). Adult eyes were obtained from the North Carolina Eye Bank (Winston-Salem, North Carolina). All whole globes were immersed in RNALater (Quiagen, Hilden, Germany) within 6.5 hours of collection, shipped overnight on ice, and dissected on the day of arrival. The retina, choroid and scleral tissues were isolated at the posterior pole using a circular, double embedded technique using round 7 mm and 5 mm biopsy punches. To reduce contamination of retina to the other ocular tissues samples, the second biopsy punch of 5 mm was used in the center of the 7 mm punch after retinal removal. RNA samples (quality control of RNA concentration and 260/280 nm ratios using Nanodrop®) (Invitrogen, Carlsbad, California, USA) were hybridized to whole genome microaray Illumina® HumanHT-12 v4 Expression BeadChips (over 25,000 genes and 48,000 probes) in two batches. The first batch was hybridized to adult RPE, choroid, and sclera RNA samples (n=6). The second batch of newer chips with additional probes was hybridized to adult optic nerve and cornea samples (n=6). The data were exported from Illumina® GenomeStudio and log2 transformed. Sample outliers were determined by principle component analyses using the Hoteling’s T2 test50 (at 95% confidence interval) and removed from further analyses. The data intensity was normalized by Quantile normalization followed by Multichip Averaging51 to reduce chip effects. For each tissue type, the probes with signal intensities below background levels and those with the lowest (5%) signal intensities (detection P<0.10) were excluded. Evidence of expression in the remaining probes was defined by detection P<0.05. Probes with detection P values < 0.10 and > 0.05 required additional tissue expression support from EyeSAGE or literature reports (Supplementary Table 3b).
We used the ‘Integrated Regulation from ENCODE’ track in the UCSC genome browser to look at H3K27ac modificiation marks as a mark of active regulatory elements. Numbers of H3K27ac modification marks were counted between the associated topSNP from a locus and the nearest gene and within (the nearest) gene itself. We also used HaploReg16 annotations to look for other signs of regulatory activity at the site of the associated SNP itself, such as enhancer histone marks, DNAse hypersensitivity sites, binding proteins and motifs changed.
We used two different programs for pathway analsysis; Ingenuity (see URLs), version August 2012, application build 172788, content version 14197757) and the Disease Association Protein-Protein Link Evaluator (DAPPLE)39.
Subcellular localization assignment and functional annotation of myopia associated disease genes as well as molecular pathway analysis was carried out using the Ingenuity knowledge database. The candidate myopia disease genes discovered in this study were entered into the Ingenuity knowledge database (IPA). We used the “IPA toggle subcellular layout” function to show the subcellular location (extracellular, plasma membrane, cytoplasm, nucleus, unknown) of the proteins corresponding to these genes, which yield a first glance which signaling molecules and pathways are involved in myopia. Subsequently, we used the IPA “connect” function to discover potential direct or indirect functional relationships or molecular pathways in between these entries. This yielded surprisingly little hits, which suggest molecular disease heterogeneity and/or functional redundancy in the pathobiological events leading to myopia. Next, we used the IPA “overlay” function to annotate the myopia candidate disease genes with (their involvement in) “functions and diseases”, “canonical pathways” and a range of custom made gene lists from previous studies, including photoreceptor, RPE, and choroidal specific transcripts (partly published52). Lastly, we used the Disease Association Protein-Protein Link Evaluator (DAPPLE)39 to look for physical connections between proteins encoded from disease-genes associated regions.
We gratefully thank the invaluable contributions of all study participants, their relatives and staff at the recruitment centers. Complete funding information and acknowledgements per study can be found in the Supplementary Note.
Author contributions V.J.M.V, P.G.H., R.W., C.J.H., C.C.W.K, A.W.H., D.A.M., T.L.Y., and C.M.v.D. performed analyses and drafted the manuscript. C.C.W.K, D.S., C.J.H., J.E.B.W., S.M.S., C.M.v.D., A.H., D.A.M., S.M., A.P., V.V, C.W., P.N.B., T.Y.W., J.R., T.L.Y., K.O., O.P., S.P.Y., J.G., A.M., M.P., S.K.I., and N.P. jointly conceived the project and supervised the work. J.E.B.W., S.M.S., D.A.M., T.L.Y., C.J.H., C.C.W.K, D.S., J.E.B.W., C.M.v.D., R.W., P.G.H., V.J.M.V., K.O., Y.Y.T., T.Y.W., P.N.B., V.V., N.A., B.A.O., A.H., J.R.V., F.R., A.G.U., N.P., C.M., A.M., T.Z., B.F., J.F.W., Z.V., O.P., A.F.W., C.H., I.R., S.K.I., E.C., J.H.L., R.P.I.J., S.J., M.S., J.J.W., P.M., I.C., J.S.R., P.M.C., C.E.P., G.W.M., A.M., W.A., F.M., M.P., L.C.K., T.D.S., E.Y.D., A.N., O.R., C.C.K., T.M., A.D., R.T.O., Y.Z., J.L., R.L., P.C., V.A.B., W.T.T., E.V., T.A., E.S.T., C.C.W.K., A.M., T.H., R.K., B.E.K.K., J.E.C., K.P.B., L.J.C., C.P.P., D.W.H.H., S.P.Y., J.W., O.P., J.B.J., L.X., H.S.W., S.M.H., A.D.P., M.K., T.L., K.M.M., C.L.S., C.W., N.J.T., D.M.E., B.S.P., J.P.K., G.M., G.H.S., M.K.I., X.Z., C.Y.C., A.W.H., S.M., R.H., J.A.G., and Q.F. were responsible for study-specific data. G.H.S., Q.F., and J.A.G. were involved in the genetic risk score analysis. T.L.Y., A.A.B.B., T.G.M.F.G., and F.H. performed the data expression experiments. A.A.B.B., T.G.M.F.G., A.M., and S.M. were involved in pathway analyses. J.E.B.W., S.M.S., D.A.M., T.L.Y., K.O., T.Y.W., P.N.B., T.G.M.F.G., S.K.I., E.C., J.J.W., A.J.M.H.V., C.C.K., B.E.K.K., S.P.Y., C.W., N.J.T., G.H.S., M.K.I., A.W.H., and J.A.G. critically reviewed the manuscript.
Competing financial interests The authors declare no competing financial interests.
URLs R, http://www.r-project.org; Locuszoom, http://csg.sph.umich.edu/locuszoom; Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20191; Ingenuity, http://www.ingenuity.com