A Genome-Wide Association Study for Susceptibility to Visual Experience-Induced Myopia.

Purpose: The rapid rise in prevalence over recent decades and high heritability of myopia suggest a role for gene-environment (G × E) interactions in myopia susceptibility. Few such G × E interactions have been discovered to date. We aimed to test the hypothesis that genetic analysis of susceptibility to visual experience-induced myopia in an animal model would identify novel G × E interaction loci.
Methods: Chicks aged 7 days (n = 987) were monocularly deprived of form vision for 4 days. A genome-wide association study (GWAS) was carried out in the 20% of chicks most susceptible and least susceptible to form deprivation (n = 380). There were 304,963 genetic markers tested for association with the degree of induced axial elongation in treated versus control eyes (A-scan ultrasonography). A GWAS candidate region was examined in the following three human cohorts: CREAM consortium (n = 44,192), UK Biobank (n = 95,505), and Avon Longitudinal Study of Parents and Children (ALSPAC; n = 4989).
Results: A locus encompassing the genes PIK3CG and PRKAR2B was genome-wide significantly associated with myopia susceptibility in chicks (lead variant rs317386235, P = 9.54e-08). In CREAM and UK Biobank GWAS datasets, PIK3CG and PRKAR2B were enriched for strongly-associated markers (meta-analysis lead variant rs117909394, P = 1.7e-07). In ALSPAC participants, rs117909394 had an age-dependent association with refractive error (-0.22 diopters [D] change over 8 years, P = 5.2e-04) and nearby variant rs17153745 showed evidence of a G × E interaction with time spent reading (effect size -0.23 D, P = 0.022). Conclusions: This work identified the PIK3CG-PRKAR2B locus as a mediator of susceptibility to visually induced myopia in chicks and suggests a role for this locus in conferring susceptibility to myopia in human cohorts.

Refractive errors, such as myopia, hyperopia, and astigmatism, occur when light is not focused accurately on the retinal photoreceptor layer of the nonaccommodated eye. Of these refractive errors, myopia—which is characterized by an axially elongated eye that causes light to focus in front of the retina—is currently a topic of intense research interest, because its prevalence has risen markedly in the past few decades and it is becoming an increasingly frequent cause of untreatable visual impairment and blindness.1 Approximately 40% to 50% of individuals in Western populations, such as the United States and Europe, are myopic,2–4 while the figure is even higher in young-adult populations of highly urbanized areas of east and southeast Asia.5–7 Myopia is associated with an increased risk of glaucoma, subcapsular cataract, retinal detachment, and maculopathy,8 and is now among the leading causes of blindness worldwide.1,9,10

Genetic studies have demonstrated refractive errors to be highly heritable,11,12 yet environmental risk factors, such as insufficient time outdoors13,14 and high levels of near-work/educational attainment,15,16 have also been shown to be important. Reconciling major roles for both genetics and environmental risk factors in myopia is challenging; gene–environment (G × E) interactions offer an attractive explanation for how environmental exposures can exert profound changes on the prevalence of myopia and yet the heritability of refractive errors can still be in the range 50% to 80% from twin and family studies.

Detecting the genetic variants underlying G × E interactions is notoriously difficult.17–21 Thus, while more than 150 independent genetic variants with genome-wide significant “main” (i.e., noninteraction) effects on refractive error have been identified,22,23 the number of variants with known G × E interaction effects is 10-fold less24–28 (and to date, only a handful of these G × E interactions been replicated in independent samples24,27). Tkatchenko et al.28 identified a myopia-predisposing interaction between variants in the promoter of APLP2 and the time children spent reading. Investigation of APLP2 was prompted by an earlier study29 in a primate model of experimentally induced myopia in which APLP2 gene expression was upregulated in the retina of eyes developing myopia. Here, we hypothesized that a genome-wide association study (GWAS) in animal model of myopia would also have the potential to identify candidate genes for myopia in humans, especially genes participating in G × E interactions controlling susceptibility to myopia induced by changes in the visual environment. In prior work, we found genetics explained approximately 50% of the interanimal variation in susceptibility to form-deprivation myopia in outbred chicks.30 Here, we built on these results by carrying out the first GWAS for myopia susceptibility, using the chick form-deprivation model.31

Methods

Experimental Animals

The experimental work was approved by the Animal Subjects Ethics Subcommittee of The Hong Kong Polytechnic University. The care and use of the animals in this study complied with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research.

White Leghorn Gallus gallus–specific pathogen-free eggs obtained from Tin Hang Tech Ltd., China were hatched in batches of approximately 20 per week. Chicks were reared in wire-mesh cages with a suspended infrared heat lamp controlling the temperature at 25°C under a 12/12-hour light/dark diurnal cycle and given access to water and commercial chick starter ad libitum. On day 7 after hatching, chicks were given a general anesthetic and their ocular component dimensions were measured by high-resolution A-scan ultrasonography. They were then monocularly deprived of sharp vision by being fitted with a translucent diffuser over one eye. On day 11 after hatching (i.e., after 4 days of form vision deprivation) the diffuser was removed and the refractive error of each eye was measured using retinoscopy while the chick was awake, followed by A-scan ultrasonography and re-attachment of the diffuser under general anesthesia. On day 12 after hatching, chicks were euthanized by carbon dioxide asphyxiation and a blood sample was collected from each chick for DNA extraction.

Form-Deprivation Treatment

One eye was chosen at random and fitted with a diffuser (made from a sheet of 0.8-mm thick polypropylene with a light absorbance of 0.07 log units). The diffuser was affixed to the periorbital skin surrounding the orbit of the treated eye using 3-4 sutures in the 12, 4, 6, and 8 o'clock positions30 while the chick was anaesthetized (intramuscular injection of ketamine 50 mg/kg and xylazine 3.5 mg/kg). After recovery from anesthesia, the treated eye was observed to be able to open freely. The diffuser could be removed for the purpose of taking retinoscopy and ultrasound measurements by cutting the sutures and re-affixed under anesthesia.

Measurements of Axial Length and Refractive Error in Chicks

All measurement procedures were carried out as described previously.30 Streak retinoscopy was performed on both eyes of awake chicks at a working distance of 33 cm, under dim illumination. Spherical refractive error was measured in the horizontal and vertical meridians, and converted to sphero-cylinder format. No correction was made for the small eye artefact of retinoscopy,32 because our interest was in the relative difference in refractive error between treated and control eyes, not their absolute refractive error.

The A-scan ultrasonography system consisted of a 20-MHz transducer of focal length 25 mm fitted with a saline stand-off of 15 mm perfused at a rate of 0.15 mL/min, a Panametrics model 5073PR pulser-receiver (Olympus Corporation, Waltham, MA, USA) and a personal computer fitted with an Acqiris DP-110 data acquisition card (Acqiris, Plan-les-Ouates, Switzerland). Waveforms were sampled at 100 MHz. For each reading 50 traces were averaged, and processed in real time to detect reflections due to the ocular components and convert time differences to distances in millimeters, using custom-written software in Visual Basic (Microsoft Corporation, Redmond, WA, USA).30 The velocity of sound in chick optical media was taken as 1.6078 mm/μs in the lens and 1.5340 mm/μs in the other ocular media.33 Between three and six readings were taken for each eye, and the average of the three highest (longest axial length) readings was calculated. Body weight was measured to the nearest gram prior to anesthesia on days 7 and 11 after hatching.

The treatment-induced change in axial length (ΔAXL) due to form deprivation was calculated as: ΔAXL = ΔAXLT − ΔAXLC, where ΔAXLT and ΔAXLC represent the change in axial length in the treated and control eye, respectively, over the 4-day period of monocular form deprivation.30 Treatment-induced change in anterior chamber depth (ΔACD), lens thickness (ΔLT), and vitreous chamber depth (ΔVCD) were calculated equivalently, for example, ΔVCD = ΔVCDT − ΔVCDC.

Genotyping and Quality Control

DNA extraction from blood samples and PCR-based sexing were carried out as described.34 Because (as described in the Results section) there was an association between ‘batch' (the group of chicks hatched and treated in each week of the study) and ΔAXL in the full sample of chicks, chicks in each weekly batch were ranked in order of ΔAXL and the top and bottom 20% were selected for genotyping. Accordingly, a total of 380 chicks were selected for genotyping, 190 from the high tail and 190 from the low tail of the ΔAXL distribution. This approach of selecting chicks within batches was designed to avoid needing to include a term for batch in the genetic analysis, as the inclusion of such a categoric variable with 48 levels would have reduced statistical power to detect genetic loci (note that the number of chicks per batch ranged from 6 to 28, with a mean of 20 and a median of 22; Supplementary Fig. S3). The 380 DNA samples were assigned to the wells of four 96-well plates, with the plate and well chosen at random in order to avoid “plate effects.” One well of each plate was assigned an “internal duplicate” sample for the purpose of quality control (QC). DNA samples were genotyped on the 600K Affymetrix Axiom Chicken Genotyping Array (Affymetrix, Inc., Santa Clara, CA, USA) by Aros-Eurofins Ltd. Internal QC was carried out by Aros-Eurofins according to the best-practice guidelines of the genotyping array manufacturer, as described in the Axiom Genotyping Solution Data Analysis Guide (P/N 702961 Rev. 3). The sample QC thresholds applied were dish QC (DQC) more than 0.82, call rate more than 97%, percentage of samples on plate passing the above thresholds more than 95%, and average call rate per plate more than 98.5%. The single nucleotide polymorphism (SNP) QC thresholds applied included call rate more than 97% and cluster quality Fisher's linear discriminant (FLD) 3.6 or more. After this internal QC, genotypes for 580,961 genetic markers were released. The average call rate for these markers was 99.5%. The average concordance rate of genotypes for the four intentionally included duplicate samples was also 99.5%.

Additional QC was performed using PLINK v1.9035 to remove markers with a call rate less than 95%, markers with duplicate genomic positions (e.g., tri-allelic variants), and markers lacking annotation information (e.g., markers located on microchromosomes). Markers were not removed based on a test for Hardy Weinberg equilibrium (HWE) because the chicks were expected to be partially inbred (i.e., related) and because alleles associated with a specific phenotype will no longer be in HWE in samples selected from the phenotype extremes.36 Sex inferred from the genotyping array was determined by examining the number of heterozygous genotypes on the Z and W chromosomes. For all 380 samples, sex inferred from the genotyping array was concordant with that determined by the PCR-based assay.

Quanto37 was used for statistical power calculations. At a genome-wide level of statistical significance (α = 0.05/304,963 = 1.6e−07) for an effect size of ΔAXL = 0.055 mm (see below) per copy of the risk allele and a normally distributed trait with mean ± SD of 0.55 ± 0.17 mm, the study had power = 0.19 of detecting a marker with minor allele frequency (MAF) 0.1, power = 0.72 for a marker with MAF = 0.2 and power greater than 0.90 for a marker with MAF greater than 0.3 (i.e., approximately 20%–90% power to detect a marker that explained 2%–5% of the variation in myopia susceptibility). Accordingly, markers with MAF less than 0.1 were removed due to limited statistical power.

All DNA samples were confirmed to have a call rate more than 95% for autosomal markers. However, one chick was excluded due to extreme heterozygosity (beyond 4 SD of the mean level for all 380 samples; potentially indicative of contamination from another DNA sample).

A genetic relatedness matrix (GRM) was created using PLINK v1.90.35 The vast majority of samples had a low level of relatedness. Specifically, of the approximately 72,000 pairwise relationships, there were 60 pairs of chicks with a GRM relatedness coefficient (as defined38) more than 0.1 and 10 pairs with a relatedness coefficient greater than 0.2 (note that half-siblings would have a relatedness coefficient ∼0.25). Only one pair of chicks had a relatedness coefficient more than 0.4 indicative of being first-degree relatives (e.g., full siblings or parent-offspring).

Genome-Wide Association Tests for Myopia Susceptibility in Chicks

Of 580,961 genetic markers, 275,998 were removed for the following reasons: a call rate of 95% or less (10,076), a MAF less than 10% (255,551), duplicate positions (62), located on the W or Z sex chromosome (6696), and no annotation information (3613). This left 304,963 SNP or insertion–deletion (indel) markers for use in the association tests. After exclusion of the single chick sample exhibiting extreme heterozygosity, there were 379 samples for inclusion in the association tests.

Linear regression modeling in the 379 selected chicks demonstrated that ΔAXL was associated with sex, final (posttreatment) body weight (FBW), and the sex × FBW interaction term (note that, by design, ΔAXL was no longer associated with batch in the sample of 379 selected chicks, due to chicks having been selected for extreme ΔAXL within each batch). Genetic-association tests were carried out with GEMMA,39 which uses a linear mixed model to account for relatedness. To account for the use of selective genotyping (i.e., the 379 chicks analyzed were from the phenotype extremes), ΔAXL values were inverse-normal transformed, because the linear mixed model implemented by GEMMA assumes a normally distributed trait. To account for covariates, the inverse-normal transformed ΔAXL was regressed on sex, FBW, and the sex × FBW interaction, and the residuals were analyzed by GEMMA. Inclusion of genotyping plate (with the 4 plates encoded using 3 binary dummy variables) as a covariate in the analysis had minimal effect on the results.

Genetic Association of Candidate Genes With Refractive Error in Humans

Genes in the Gal4 or Gal5 chicken genome build within 100 kb of all genetic variants with P < 1.64e−05 in the chick myopia susceptibility GWAS were selected as candidate genes (where 1.64e−05 corresponded to the ‘suggestive significance threshold', defined as 100× higher than the genome-wide significance threshold of 1.64e−07). This identified eight candidate genes, all of which had human homologues. The genomic coordinates of the human genes were obtained from the University of California, Santa Cruz (UCSC) Genome Browser for genome build GRCh37.3 (hg19).

A gene-based test (MAGMA v1.06)40 was used to assess whether the candidate genes identified in the chick myopia susceptibility GWAS were enriched for genetic markers associated with refractive error in human GWAS studies (a gene-based test was chosen in preference to seeking to find a human genetic variant analogous to the lead variant identified in the chick GWAS, because the presence and function of genetic variants is very unlikely to be conserved between species as evolutionarily distant as human and chicken). MAGMA's gene-based test considers all of the markers within a specified gene locus (here, we considered the genomic interval between the transcription start and stop site of a gene, plus a flanking 20-kb region at the 5′ and 3′ ends) and accounts for the nonindependence of markers in linkage disequilibrium (LD). MAGMA correctly accounts for gene size and for the variable density of genetic variants within genes40 (e.g., some genes may contain many more SNPs than others). Two sets of human GWAS summary statistics were analyzed. First, a meta-analysis of GWAS for spherical equivalent refractive error in n = 44,192 participants of European ancestry aged older than 25 years carried out by the CREAM consortium.23 The contributing studies of the CREAM consortium imputed genotype data to the 1000-Genomes Project phase 3 reference panel. Participants provided informed consent during recruitment into the individual studies.23 Secondly, a GWAS for spherical equivalent refractive error in n = 95,505 UK Biobank participants of European ancestry aged 37 to 73 years who had autorefraction information available and no history of eye disorders, carried out by the UK Eye & Vision (UKEV) consortium.41 Genotype data for UK Biobank participants were imputed to the Haplotype Reference Consortium (HRC). reference panel and a combined 1000 Genomes Project-UK10K reference panel by Bycroft et al.42 The UKEV GWAS41 only analyzed markers on the HRC reference panel43 with MAF 0.05 or more and IMPUTE4 INFO metric greater than 0.9. Ethical approval for the UK Biobank project was obtained from the National Health Service (NHS) National Research Ethics committee (Ref. 11/NW/0382) and all participants provided informed consent. Demographic details for the human replication samples are provided in Supplementary Table S1.

Quantile-quantile (QQ) plots were created for the markers in a gene locus plus a flanking 20-kb region at the 5′ and 3′ ends. QQ plots were also created for genetic markers selected as being in linkage equilibrium in unrelated UK Biobank participants of European ancestry (using the –indep-pairwise 50 5 0.1 command in PLINK 2.0).35 Regional association plots were created using LocusZoom44 for all markers from 100 kb upstream of the PIK3CG gene to 100 kb downstream of the PRKAR2B gene. Fine-mapping was carried out using FINEMAP45 for variants within ± 500 kb of the lead variant rs117909394, with LD between variants computed using PLINK 1.9 for a randomly selected sample of 10,000 unrelated UK Biobank participants of White European ancestry.

Tests for PIK3CG-PRKAR2B Locus Gene-Environment Interactions in Children

Tests for G × E interactions were performed using data from the ALSPAC cohort. The ALSPAC research team recruited 14,541 pregnant women residing in Avon, UK with expected dates of delivery between April 1, 1991 and December 31, 1992. Of the initial 14,541 pregnancies, 13,988 children were alive at 1 year. The original cohort was largely representative of the UK 1991 Census; however, over time there was a trend for greater attrition of families of low socioeconomic position and of non-White ethnic origin. Details of the ALSPAC study cohort have been published.46,47 Ethical approval for the study was obtained from the ALSPAC Ethics and Law Committee and Local Research Ethics Committees. The ALSPAC website contains details of all the data available through a fully searchable data dictionary (in the public domain, http://www.bris.ac.uk/alspac/researchers/data-access/data-dictionary/).

Refractive error was assessed longitudinally in ALSPAC participants using noncycloplegic autorefraction at research clinics attended when the children were approximately 7.5, 10.5, 11.5, 12.5, and 15.5 years of age.12 DNA samples from approximately 10,000 participants were genotyped on the Illumina HumanHap550 quad chip genotyping platform. The time participants spent outdoors each day when they were aged approximately 8-years old was ascertained from a questionnaire completed by the child's mother or carer, with the question, “On a school weekday, how much time on average does your child spend each day out of doors in summer?” The response options were “None at all,” “1 hour or less,” “1 to 2 hours,” or “3 or more hours.” Children were classified as spending a “high” amount of time outdoors if the response was “1 to 2 hours” or “3 or more hours,” and as “low” otherwise. The time participants spent reading at age 8-years old was ascertained from the question on the same questionnaire, “On normal days in school holidays, how much time on average does your child spend each day reading books for pleasure?” Children were classified as spending a “high” amount of time reading if the response was “1 to 2 hours” or “3 or more hours” per day, or a “low” amount of time reading if the response was “None at all” or “1 hour or less.”

There were 7981 children with genome-wide genotype data remaining after excluding those whose data failed QC assessments, related individuals, participants with non-European ancestry, and individuals who withdrew their consent. Of these 7981 individuals, 4989 had refractive error information from at least three visits during the 7.5- to 15.5-year period. Among the latter sample, 2047 and 2308 children were classified as having spent a high versus low amount of time outdoors at age 8 years (904 had missing information about time spent outdoors), while 1624 and 2697 children were classified as spending a high versus low amount of time reading (with 668 having missing information). Demographic characteristics of the ALSPAC sample are provided in Supplementary Table S1. Linear mixed models28 were used to examine refractive trajectories and the influence of interactions between a specific genetic variant (1 of 3 lead SNPs in the PIK3CG-PRKAR2B region) and an environmental exposure (either time spent outdoors or time spent reading). The refractive error of each child at the baseline visit and the linear trajectory of refractive error with age were modeled as individual-level random effects. All models contained fixed-effects polynomial terms for age (for the first, second, third, and fourth order), as well as the predictive variables SNP genotype (coded 0, 1, or 2) and environmental exposure (coded as 0 = low, 1 = high) and terms for a SNP × age interaction and an environmental exposure × age interaction. More complex models including terms for a two-way interaction between SNP genotype and the environmental exposure and a three-way interaction between SNP genotype, the environmental exposure and age were also evaluated. In preliminary tests, model fits were similar with or without a term for sex; therefore, sex was excluded from the final models. Two of three SNPs tested had low MAF (MAF = 0.05 for rs17153745 and MAF = 0.03 for rs117909394), which meant there were fewer than 10 individuals homozygous for the minor allele. Individuals homozygous for the minor allele of rs17153745 and rs117909394 were therefore recoded as heterozygotes to avoid instability in the model fits (i.e., SNP genotype was coded as 1 or 2 rather than 0, 1, or 2 for these SNPs).

Results

Myopia Induced by the Monocular Deprivation of Sharp Vision

A total of 987 chicks, hatched in 48-weekly batches, were monocularly deprived of sharp vision (form deprived) for 4 days to induce myopia in the treated eye. A PCR-based sex test was successful in 959 chicks, and revealed that 52% of chicks were male and 48% female. The posttreatment ocular features of the chicks are presented in Table 1. The average ΔAXL was 0.55 ± 0.17 mm (mean ± SD; P < 2.2e−16). Although male chicks had longer eyes than females both before and after treatment, ΔAXL was similar in the two sexes (P = 0.80; Table 1). Treated eyes were more myopic than control eyes (−4.16 ± 3.02 vs. 6.52 ± 0.96 diopters [D], P < 2.2e−16). Prior to treatment, female chicks were more hyperopic than male chicks (6.62 ± 1.00 vs. 6.42 ± 0.91 D, P < 0.001); however, the level of induced myopia (ΔMSE) was similar in the two sexes (−10.64 ± 3.07 D in males and −10.73 ± 2.97 D in females, P = 0.64; Table 1). The correlation between ΔAXL and ΔMSE was 0.74 (Fig. 1B, P < 2.2e−16).

Table 1

Ocular Component Dimensions After the Form Deprivation Treatment Period

Trait	Eye	Males (n = 501)	Females (n = 458)	All (n = 959)	P Value (M vs. F)
FBW, g		78.0 ± 10.6	77.4 ± 10.5	77.7 ± 10.5	0.42
ACD, mm	Control	1.46 ± 0.04	1.43 ± 0.04	1.45 ± 0.04	<0.001
	Treated	1.57 ± 0.08	1.53 ± 0.08	1.55 ± 0.08	<0.001
	ΔACD (T − C)	0.11 ± 0.07	0.10 ± 0.06	0.10 ± 0.06	0.08
	P value (T vs. C)	<0.001	<0.001	<0.001
LT, mm	Control	2.15 ± 0.05	2.11 ± 0.04	2.13 ± 0.05	<0.001
	Treated	2.15 ± 0.05	2.12 ± 0.05	2.14 ± 0.05	<0.001
	ΔLT (T − C)	0.004 ± 0.04	0.01 ± 0.04	0.01 ± 0.04	0.34
	P value (T vs. C)	<0.001	<0.001	<0.001
VCD, mm	Control	5.49 ± 0.17	5.40 ± 0.17	5.45 ± 0.17	<0.001
	Treated	5.93 ± 0.22	5.84 ± 0.22	5.88 ± 0.22	<0.001
	ΔVCD (T − C)	0.43 ± 0.14	0.44 ± 0.14	0.44 ± 0.14	0.17
	P value (T vs. C)	<0.001	<0.001	<0.001
AXL, mm	Control	9.10 ± 0.19	8.94 ± 0.2	9.02 ± 0.21	<0.001
	Treated	9.65 ± 0.27	9.49 ± 0.27	9.57 ± 0.28	<0.001
	ΔAXL (T − C)	0.54 ± 0.17	0.55 ± 0.17	0.55 ± 0.17	0.80
	P value (T vs. C)	<0.001	<0.001	<0.001
MSE, D	Control	+6.42 ± 0.91	+6.62 ± 1.00	+6.52 ± 0.96	<0.001
	Treated	−4.22 ± 3.08	−4.10 ± 2.95	−4.16 ± 3.02	0.56
	ΔMSE (T − C)	−10.64 ± 3.07	−10.73 ± 2.97	−10.68 ± 3.02	0.64
	P value (T vs. C)	<0.001	<0.001	<0.001

Values are presented as mean ± SD. LT, lens thickness; Δ, Difference between treated eye (T) minus control eye (C).

Figure 1

Frequency distribution of ΔAXL and ΔMSE in the full chick sample (n = 987). Distribution of ΔAXL (A), ΔMSE (C), and relationship between ΔAXL and ΔMSE (B). Chicks selected for having a high susceptibility (n = 190; black) or low susceptibility (n = 190; dark gray) to induced myopia are indicated. Chicks not selected for high or low susceptibility are plotted as light gray symbols.

Regression models were used to identify covariates associated with the primary outcome measure, ΔAXL, in the full sample of 959 chicks. Separate models were fitted for body weight prior to treatment (initial body weight; IBW) and body weight after treatment (FBW), because these traits were highly correlated. ‘Batch'—the variable indicating the group of chicks hatched and treated each week—was modeled as a categorical variable with 48 levels. In the first model, sex, batch, and the interaction between sex and IBW were associated with ΔAXL. Together these covariates explained 3.2% of the variation in ΔAXL. In the second model, sex, batch, FBW, and sex × FBW were associated with ΔAXL, with similar effect sizes to those observed in the first model. The second model explained 5.8% of the variance in ΔAXL.

To reduce costs, not all of the chicks were genotyped; instead only chicks selected as being in the myopia-susceptibility phenotype extremes were genotyped. Chick selection was based on the phenotype ΔAXL, rather than ΔMSE, because we reasoned this could be measured more accurately (active accommodation during retinoscopy meant that refractive error could only be measured to the nearest 0.50 D, which for a typical ΔMSE of −10.00 D corresponded to 5% [0.5× 100/10]. The SD of A-scan AXL measurements was approximately 7 μm, which for a typical ΔAXL of 0.5 mm corresponded to 1.4% [7 × 100/500]. Although not exactly equivalent, these figures of 5% and 1.4% suggested ΔAXL was the more accurately determined trait). It was originally planned that the chicks would be ranked according to ΔAXL, and the top and bottom 20% of the full sample selected for genotyping. However, in view of the association between ΔAXL and batch, instead of selecting chicks from the whole population at once, chicks were selected within each batch of hatchlings separately. Thus, from within each batch, the 20% of chicks with highest ΔAXL and the 20% of chicks with the lowest ΔAXL were selected. A total of 380 chicks were selected; 190 with high susceptibility to form-deprivation myopia and 190 with low susceptibility (Fig. 1). As a result of the within-batch selection scheme, ΔAXL was no longer associated with batch in the sample of selected chicks. The average ΔAXL was 0.31 ± 0.08 mm in the chicks selected for low ΔAXL, and 0.78 ± 0.08 mm in the chicks selected for high ΔAXL (Fig. 1). The mean difference in ΔAXL between chicks in the high and low subsamples was 0.47 mm (95% CI: 0.45–0.48, P < 2.2e−16). The average ΔMSE was −7.14 ± 2.29 D in the low ΔAXL subsample, and −13.55 ± 2.29 D in the high ΔAXL subsample. The difference in ΔMSE between the high and low subsamples was 6.41 D (95% CI: 5.95–6.88, P < 2.2e−16).

One of the 380 genotyped chicks was excluded from the analysis because of an unusual level of genetic heterozygosity. A genome-wide association analysis for ΔAXL that took account of familial relatedness was carried out in the remaining 379 chicks. Analysis was restricted to 304,963 genetic markers on the autosomal chromosomes that were highly polymorphic in the study sample (MAF > 0.1). A single marker passed the Bonferroni-corrected significance threshold of P < 1.64e−07 (Fig. 2; Supplementary Fig. S1). This was marker rs317386235 (P = 9.54e−08), situated on chick chromosome 1 between the PIK3CG and PRKAR2B genes, which was part of a cluster of variants strongly associated with ΔAXL (Supplementary Table S2). A conditional analysis in which rs317386235 genotype was included in the analysis as a covariate reduced the association signal at the PIK3CG-PRKAR2B locus to background levels, suggesting the presence of a single causal variant in the region (lead variant within 100 kb, after conditioning on rs317386235, was rs313346637, P = 0.008). The only other genetic marker with a P value below the suggestive significance threshold of P < 1.64e−05 was rs313633102 on chick chromosome 12 (P = 1.62e−05; Supplementary Table S2). Additionally, exploratory genome-wide association analyses for ΔMSE or for high versus low myopia susceptibility did not yield any markers with P values below the Bonferroni-corrected genome-wide significance threshold, likely because of limited statistical power resulting from the modest sample size.

Figure 2

Manhattan plot from a GWAS for susceptibility to environmentally induced myopia in chicks. Three hundred seventy-nine chicks from the extremes of the phenotype distribution were analyzed. The x-axis indicates genomic position, the y-axis indicates negative log10(P value). The long-dashed horizontal line indicates the Bonferroni-corrected significance threshold (P = 1.64e−07) and the short-dashed horizontal line the suggestive significance threshold (P = 1.64e−05). Markers situated on alternate chromosomes are plotted as light or dark gray symbols.

Genes within 100 kb of the genetic variants with P values below the suggestive significance threshold of 1.64e−05 in the chick myopia susceptibility GWAS were considered as candidate genes (Supplementary Table S2). This yielded the following eight candidate genes, all of which had human homologues: SYPL1, NAMPT, PIK3CG, PRKAR2B, HBP1, COG5, CENPP, and OGN. A gene-based bioinformatics analysis (MAGMA40) was carried out to test whether any of these genes were enriched for markers associated with refractive error in two human GWAS analyses carried out by the CREAM consortium23 and the UK Biobank Eye & Vision (UKEV) consortium.41 The CREAM GWAS provided evidence that PIK3CG (P = 0.004), PRKAR2B (P = 0.024), and COG5 (P = 0.045) were enriched for markers associated with refractive error (Supplementary Table S3). After accounting for multiple testing (8 genes) the PIK3CG gene retained evidence of association (P = 0.033). The UKEV GWAS provided suggestive evidence that PIK3CG was enriched for markers associated with refractive error (uncorrected P = 0.054; Supplementary Table S4). PIK3CG, PRKAR2B, and COG5 are situated close together on chromosome 7q22.3, and hence markers situated in the flanking regions of any one of these genes may have contributed to the evidence supporting the association of the adjacent genes.

Both the chick myopia susceptibility GWAS and the human GWAS gene-based tests provided the strongest support for the involvement of the PIK3CG and PRKAR2B genes, which both have functions related to cAMP signaling.48 Notably, the phosphoinositide 3-kinase (PI3K) signaling axis has been implicated in previous studies of refractive development. Specifically, PRKAR2B gene expression in the retina was found to be upregulated 1.7-fold in chicks wearing +7 D lenses for 24 hours,49 and the PI3K inhibitor Ly294002 was shown to prevent the exaggerated response to negative lens wear in chicks treated with intravitreal insulin.50 Therefore, attention was focused on the PIK3CG and PRKAR2B genes. To visualize the enrichment of GWAS markers associated with PIK3CG and PRKAR2B, QQ plots and regional association plots were created for the CREAM GWAS summary statistics, the UKEV GWAS summary statistics, and a meta-analysis of the two GWAS datasets (Fig. 3; Supplementary Fig. S2). As expected from the gene-based test results, the QQ plots confirmed that there was an excess of markers with low P values compared with what would be expected under the null hypothesis of no association, including when only markers selected as being in linkage equilibrium were considered (note that this selection of markers in linkage equilibrium was done without reference to any phenotypic information). The regional association plots (Supplementary Fig. S2) suggested that markers associated with refractive error in human participants were distributed across a broad region of low recombination rate encompassing PIK3CG, PRKAR2B, and the intergenic region between them. The most strongly associated marker differed in the following three GWAS datasets: the lead markers were intergenic variant rs17153745 in the CREAM GWAS, rs757903 in the promoter region of PIK3CG in the UKEV GWAS, and rs117909394 situated within the coding region of PRKAR2B in the meta-analysis of the CREAM and UKEV GWAS. Table 2 summarizes the level of association with refractive error in each of the three GWAS datasets, for these three markers. Attempts to fine-map the region using a Bayesian approach45 that evaluated the evidence for between 1 and 5 causal variants in the region consistently provided evidence favoring multiple causal variants rather than a single causal variant. However, it was not possible to confidently select a set of potentially causal markers, perhaps due to the low level of recombination across the PIK3CG-PRKAR2B region.

Figure 3

QQ plots of the PIK3CG-PRKAR2B gene region show enrichment of markers with low P values in human GWAS analyses. Data from the UKEV consortium (A, B), the CREAM consortium (C, D), and a meta-analysis of the UKEV and CREAM consortium (E, F) GWAS for refractive error. For the PIK3CG gene region (A, C, E) and the PRKAR2B gene region (B, D, F) the negative log10 P value expected under the null hypothesis of no association is plotted on the x-axis, and the observed negative log10 P value on the y-axis. Black symbols are results for all variants; dark gray symbols are variants preselected as being in linkage equilibrium. The light gray shaded region shows the 95%CI for the distribution expected under the null hypothesis. The black line is the line of unity (x = y). Variants within 20 kb of the transcription start and sites were included in the QQ plots.

Table 2

Features of the Lead Markers in the PIK3CG-PRKAR2B Gene Region Most Strongly Associated With Refractive Error in Human GWAS Analyses

SNP	CHR	POS	EA	NEA	MAF	CREAM			UKEV			Meta-Analysis
						BETA	SE	P	BETA	SE	P	BETA	SE	P
rs757903	7	106500436	C	T	0.26	−0.026	0.017	1.2e−01	−0.040	0.013	3.5e−03	−0.035	0.011	9.9e−04
rs17153745	7	106659112	A	G	0.05	−0.180	0.036	7.8e−07	−0.065	0.028	4.2e−02	−0.108	0.022	1.2e−06
rs117909394	7	106723885	G	A	0.03	−0.212	0.043	8.7e−07	−0.092	0.032	1.1e−02	−0.135	0.026	1.7e−07

From an analysis of markers located 100-kb upstream of the transcription start of PIK3CG to 100-kb downstream of the transcription stop site of PRKAR2B, the markers rs757903, rs17153745, and rs117909394 were the most strongly associated variants in the UKEV GWAS, the CREAM GWAS, and the CREAM + UKEV meta-analysis, respectively. CHR, chromosome; POS, genomic position for genome build GRCh37.3; EA, effect allele; NEA, noneffect allele; BETA, effect size in dioptres per copy of the effect allele; SE, standard error of BETA.

PIK3CG-PRKAR2B Locus Gene–Environment Interactions in Children

Because the aim of the chick GWAS was to identify variants influencing myopia through a change to the visual environment, it was of interest to determine whether variants at the PIK3CG-PRKAR2B locus exerted G × E interaction effects in humans. Hence, G × E analyses were carried out using data from ALSPAC cohort participants whose refractive development had been followed longitudinally from age 7.5- to 15.5-years old and for whom information on time spent reading and time spent outdoors was available.46 To reduce the extent of multiple testing, only the three SNPs most strongly associated with refractive error in the adult GWAS datasets were evaluated (i.e., those described in Table 2).

Baseline models were fitted initially to test whether the three lead SNPs had main effects (i.e., noninteraction effects) associated with refractive error in childhood. None of the SNPs was associated with refractive error at the age of 7 years; however, 2 of 3 SNPs (rs17153745 and rs117909394) were associated with children's refractive error trajectory between ages 7.5 and 15.5 years (Table 3). Specifically, the A allele of rs17153745 was associated with a −0.014 D/y change in refractive error (P = 0.028) and the G allele of rs117909394 was associated with a −0.027 D/y change in refractive error (P = 5.2e−04). For both variants, the allele predisposing children to a more negative (myopic) refractive error matched the allele associated with a more negative refractive error in the CREAM consortium GWAS and in the UKEV consortium GWAS (Table 2). Thus, the ALSPAC cohort analyses provided further independent evidence of association between the PIK3CG-PRKAR2B locus and susceptibility to refractive error.

Table 3

Tests for Main Effects and Age-Dependent Effects of the Lead Variants in the PIK3CG-PRKAR2B Gene Region for Refractive Error in Children Aged 7- to 15-Years Old From the ALSPAC Cohort

SNP	EA	MAF	Main Effect, D			SNP × Age Interaction, D/y
			BETA	SE	P	BETA	SE	P
rs757903	C	0.26	0.014	0.020	4.9e−01	−0.004	0.003	2.1e−01
rs17153745	A	0.05	−0.039	0.045	3.9e−01	−0.014	0.006	2.8e−02
rs117909394	G	0.03	0.087	0.055	1.1e−01	−0.027	0.008	5.2e−04

Models testing for genetic effects of the three SNPs that differed depending on the time a child spent outdoors at age 8-years old (categorized as either high or low) revealed minimal evidence for such interactions (Supplementary Tables S5−S6). Models testing for interactions between the genotype of the lead SNPs and time spent reading, also categorized as either high or low, likewise provided minimal evidence for three-way (SNP × age × time spent reading) interactions (Supplementary Table S7). However, there was suggestive evidence for a SNP × time spent reading interaction at the baseline age of 7.5-years old for one of three SNPs (Fig. 4; Table 4). Specifically, for children classified as spending a high amount of time reading, those carrying one copy of the risk allele of rs17153745 had a refractive error at age 7.5 years approximately −0.23 D more negative than those not carrying any risk alleles (P = 0.022). Because approximately 90% of Europeans actually carry two copies of the myopia-predisposing allele (A) of rs17153745, this G × E interaction effect is more easily understood as suggesting that the approximately 10% of children who carry a copy of the alternate (G) allele of rs17153745 were relatively protected against the more myopic refractive error usually associated with spending a high amount of time reading. However, we caution that the level of evidence supporting this G × E association with time spent reading was modest (P = 0.022) and that this finding could be a false-positive result, especially because statistical power would have been limited by the low MAF of rs17153745 (MAF = 0.05).

Figure 4

Analysis models for detecting G × E interactions in childhood refractive error trajectories. Graphs show the best-fitting models for the two predictor variables, genotype at SNP rs17153745 from the PIK3CG-PRKAR2B locus and time spent reading (classified as high or low). The three models considered were, a model with no G × E effects (A), only a two-way rs17153745 × time reading interaction at age 7.5-years old (B), or both a 2-way rs17153745 × time reading interaction at age 7.5-years old and a three-way rs17153745 × time reading × age interaction (C).

Table 4

Tests for Two-Way Interactions Between Time Spent Reading at Age 8.5-Years Old and SNP Genotype of Lead Variants in the PIK3CG-PRKAR2B Gene Region for Refractive Error in Children Aged 7- to 15-Years Old From the ALSPAC Cohort

SNP (Test Allele)	SNP Effect at Age 7 y, D			Time Reading Effect at Age 7 y, D			SNP × Time Reading Interaction at Age 7 y, D			SNP × Age Interaction, D/y			Time Reading × Age Interaction, D/y
	BETA	SE	P	BETA	SE	P	BETA	SE	P	BETA	SE	P	BETA	SE	P
rs757903 (C)	−0.005	0.028	8.7e−01	−0.109	0.073	1.3e−01	0.045	0.045	3.1e−01	−0.005	0.003	1.2e−01	−0.022	0.004	8.0e−09
rs17153745 (A)	0.063	0.066	3.4e−01	0.398	0.192	3.9e−02	−0.229	0.099	2.2e−02	−0.018	0.007	1.0e−02	−0.023	0.004	5.6e−09
rs117909394 (G)	0.152	0.077	4.9e−02	0.268	0.240	2.6e−01	−0.162	0.123	1.9e−01	−0.032	0.008	1.1e−04	−0.022	0.004	1.7e−08

Discussion

The discovery of the PIK3CG-PRKAR2B region as a susceptibility locus for refractive error supports our hypothesis that genetic studies in animal models can be leveraged to identify genes with roles in human myopia, as previously suggested.28 The GWAS for myopia susceptibility in chicks supported the polygenic nature of this phenotype,30,51 but highlighted the large sample size necessary to distinguish true-positive association signals against the background of false-positives arising from testing hundreds of thousands of genetic markers. Nevertheless, our findings suggest that sample sizes in animal model studies do not need to reach the tens or hundreds of thousands required to detect G × E effects in population-based human studies.17–21

In designing the study, form deprivation was chosen in preference to negative lens wear as the myopia-inducing visual stimulus. Form deprivation has the following advantages: (1) it is an ‘open-loop' stimulus, and hence the level of induced myopia would not be subject to a ‘ceiling effect', (2) there was prior evidence that susceptibility to form-deprivation myopia in chicks is heritable and polygenic,30 whereas such information is lacking for lens-induced myopia, and (3) using high-power negative lenses, or failure to keep the lenses scrupulously clean, risks inadvertently inducing form deprivation.

A strength of this work was that multiple human cohorts were used to assess the evidence for a role of the PIK3CG-PRKAR2B locus in refractive development. Limitations of the study were as follows: (1) not all of the chick population was genotyped in the myopia susceptibility GWAS, which reduced statistical power, future studies with larger sample sizes would be expected to yield a greater number of candidate loci for follow-up; (2) information about the key childhood exposures of time outdoors and time spent reading were not available for the UKEV and CREAM study participants, which meant that G × E interactions had to be tested in the smaller ALSPAC sample rather than in the larger adult samples; (3) autorefraction in ALSPAC participants was carried out without prior cycloplegia, which is known to result in higher measurement error52,53 and therefore would have further reduced statistical power to detect G × E effects; (4) statistical fine-mapping of the PIK3CG-PRKAR2B locus was not successful in prioritizing candidate causal variants, complicating future work to dissect the molecular mechanisms linking genetic variation to myopia susceptibility; (5) in view of the comparatively weak evidence for SNP × time spent reading interactions in the ALSPAC sample, we were not able to demonstrate definitively that PIK3CG-PRKAR2B variants interact with environmental exposures to cause myopia, although there was strong supporting evidence for main effects at the locus in independent human cohorts; and (6) the mechanism by which PIK3CG-PRKAR2B variants influence myopia susceptibility—and indeed, even the tissues in which these genes act to mediate this effect—could not be determined from our genetic investigations. Further studies in animal models will be required to test potential mechanistic hypotheses.

In conclusion, a GWAS for susceptibility to myopia induced by changes in the visual environment of chicks identified a single genomic locus after applying stringent statistical criteria to account for multiple testing. The lead variant, rs317386235 (P = 9.54e−08), was located on chick chromosome 1 between the PIK3CG and PRKAR2B genes. In humans, the PIK3CG-PRKAR2B region was demonstrated to be enriched for variants associated with refractive error in two large GWAS datasets (CREAM: lead variant rs17153745, P = 7.8e−07; UKEV: lead variant rs757903, P = 0.004; CREAM+UKEV meta-analysis: lead variant rs117909394, P = 1.7e−07). In a sample of children with longitudinal measurements of refractive error over an 8-year period, lead variant rs117909394 had an age-dependent association with refractive error (P = 5.2e−04) and variant rs17153745 demonstrated a suggestive G × E interaction with time spent reading at age 7- to 8-years old (P = 0.022). This work further supports studies in animal models of myopia as a route for the discovery of genetic and lifestyle factors responsible for myopia development in humans.

Click here for additional data file.