Assumption-free estimation of the genetic contribution to refractive error across childhood.

PURPOSE: Studies in relatives have generally yielded high heritability estimates for refractive error: twins 75-90%, families 15-70%. However, because related individuals often share a common environment, these estimates are inflated (via misallocation of unique/common environment variance). We calculated a lower-bound heritability estimate for refractive error free from such bias.
METHODS: Between the ages 7 and 15 years, participants in the Avon Longitudinal Study of Parents and Children (ALSPAC) underwent non-cycloplegic autorefraction at regular research clinics. At each age, an estimate of the variance in refractive error explained by single nucleotide polymorphism (SNP) genetic variants was calculated using genome-wide complex trait analysis (GCTA) using high-density genome-wide SNP genotype information (minimum N at each age=3,404).
RESULTS: The variance in refractive error explained by the SNPs ("SNP heritability") was stable over childhood: Across age 7-15 years, SNP heritability averaged 0.28 (SE=0.08, p<0.001). The genetic correlation for refractive error between visits varied from 0.77 to 1.00 (all p<0.001) demonstrating that a common set of SNPs was responsible for the genetic contribution to refractive error across this period of childhood. Simulations suggested lack of cycloplegia during autorefraction led to a small underestimation of SNP heritability (adjusted SNP heritability=0.35; SE=0.09). To put these results in context, the variance in refractive error explained (or predicted) by the time participants spent outdoors was <0.005 and by the time spent reading was <0.01, based on a parental questionnaire completed when the child was aged 8-9 years old.
CONCLUSIONS: Genetic variation captured by common SNPs explained approximately 35% of the variation in refractive error between unrelated subjects. This value sets an upper limit for predicting refractive error using existing SNP genotyping arrays, although higher-density genotyping in larger samples and inclusion of interaction effects is expected to raise this figure toward twin- and family-based heritability estimates. The same SNPs influenced refractive error across much of childhood. Notwithstanding the strong evidence of association between time outdoors and myopia, and time reading and myopia, less than 1% of the variance in myopia at age 15 was explained by crude measures of these two risk factors, indicating that their effects may be limited, at least when averaged over the whole population.

Introduction

Myopia is the most common eye disorder worldwide [1,2] and is a risk factor for several important ocular pathologies, including glaucoma, retinal detachment, and chorioretinal degeneration [3-6]. The etiology of myopia is poorly understood, although genetic and environmental influences are known to be involved [for recent reviews, see 1,2,6–15].

The extent to which a quantitative trait is determined by genetic factors is termed “heritability.” Heritability in the “broad sense” (H2), as calculated in twin studies, refers to the proportion of the phenotypic variation in the trait due to additive, dominant, and gene–gene interaction effects [16] while heritability in the “narrow sense” (h2), as calculated in family studies, refers to the proportion of phenotypic variance explained by additive genetic effects alone [17]. Since genes and environment never operate in isolation, the utility of a heritability estimate may appear questionable [16,18]. However, the narrow sense heritability indicates the potential of genetic testing for identifying individuals at high risk of a specific disease [19]. Refractive error has been shown to be highly heritable [1,2,6-15,20-24]: In a recent review [25], the broad sense heritability ranged from 0.75 to 0.90 and the narrow sense heritability from 0.15 to 0.70.

The availability of techniques to accurately and cost-effectively survey genetic variation across the genome now makes it feasible to quantify the contribution of causal variants tagged by genotyped or imputed single nucleotide polymorphisms (SNPs) to differences in metric traits or disease susceptibility between individuals [26,27]. The inter-subject variance in a trait that can be explained by all the SNPs has been termed SNP heritability (h2SNP). For an infinitely large population in which every polymorphism was genotyped without error, SNP heritability h2SNP would equal heritability in the narrow sense h2. For finite-sized samples in which only common SNPs can be genotyped or imputed accurately, h2SNP therefore corresponds to a lower-bound estimate of h2. Theoretical and empirical work suggests that h2SNP provides a valuable guide to, first, the potential for predicting at-risk individuals using currently available genotyping platforms and, second, the potential for identifying specific disease-associated genetic variants in a genome-wide association study (GWAS) using existing genotyping and imputation methods but applied to larger samples.

An important, fundamental difference between SNP-heritability estimation methods and prior methods of estimating heritability is that h2SNP is computed using a sample of essentially unrelated individuals. Thus, h2SNP estimates have the attractive feature of not being biased by the confounding effects of shared family environment. Furthermore, they are not reliant on other assumptions required to obtain comparable estimates in twins [26,27]. As with conventional h2 estimates, the genetic correlation between two traits (i.e., the degree to which the two traits are determined by a shared set of additive genetic variants) can also be calculated using SNP data.

Prior work suggests that at least 40–50% of the between-subject variance in height [27,28], axial length [29], and corneal curvature [29] can be explained by common SNPs. We applied this approach to estimate the SNP heritability of refractive error for children participating in a birth cohort study. We also calculate genetic correlations between refractive error measured at different ages, to infer whether the same or different genetic variants influence susceptibility to refractive error across childhood.

Methods

Subjects and ethical approval

The Avon Longitudinal Study of Parents and Children (ALSPAC) recruited 14,541 pregnant women resident in Avon, England, with expected dates of delivery 1 April 1991 to 31 December 1992. When the oldest children were approximately 7 years of age, an attempt was made to bolster the initial sample with eligible cases who had failed to join the study originally, resulting in an additional 713 children being enrolled. The total number of children who were alive at 1 year of age was 14,701. The enrollment is described in more detail in a cohort profile article [30]. This research adhered to the tenets of the Declaration of Helsinki. Ethical approval for the study was obtained from the ALSPAC Ethics and Law Committee and the Local Research Ethics Committees (LRECs). Separate consent was given for the collection of refractive data from the participant’s own optometrist. The ALSPAC website contains details of all the data available through a fully searchable data dictionary.

Phenotypes

ALSPAC participants were invited to a research clinic when they were aged approximately 7½, 10½, 11½, 12 ½, and 15½ years old. At each visit, refractive error was estimated with non-cycloplegic autorefraction (Canon R50 instrument, Canon USA Inc., Lake Success, NY) as described and validated previously [31-33]. After the outlier readings were removed, the mean spherical equivalent (SE) refractive error was calculated as the autorefraction sphere power plus half of the cylinder power and averaged between the two eyes. For one set of analyses, refractive errors were transformed to normal scores [34] to test whether the non-normal distribution of the trait biased the results. As discussed below, and in the Appendix, autorefraction in the absence of cycloplegia is known to lead to inaccuracy, with the degree of measurement error varying with the age of the child, his or her level of refractive error, the type of autorefractor used, and other subject-specific and random factors [33,35-38].

When the study children were aged 8–9 years, a questionnaire was completed by their mother that included the following questions: (1) “On a weekend day, how much time on average does your child spend each day out of doors in summer?” (2) “On normal days in school holidays, how much time on average does your child spend each day reading books for pleasure?” The response options for these questions were “None at all,” “1 hour,” “1–2 hours,” and “3 or more hours.”

Genetic and statistical analyses

DNA samples were genotyped using Illumina HumanHap550-Quad bead arrays, as described [39,40]. Blood sampling and DNA extraction were carried out as described [30]. Samples of known non-European ancestry, with missingness >3%, minimal or excessive autosomal heterozygosity, cryptic relatedness >10% identity by descent (IBD), or with a sex mismatch were excluded. EIGENSTRAT analysis and multidimensional scaling analysis seeded with individuals from HapMap phase 2 revealed no additional outliers. SNPs with call rate <95%, minor allele frequency (MAF) <1%, or Hardy–Weinberg equilibrium (HWE) p value ≤5×10−7 were excluded. After phasing haplotypes with ShapeIT v2.r644 [41], IMPUTE2 (V2.2.2) [42] was used to impute unobserved marker genotypes for the 8,237 participants whose samples passed the quality control criteria (imputation was performed simultaneously with a sample of 8,196 mothers of ALSPAC participants [43], who were included to improve phasing and imputation accuracy) with the Dec 2013 release of the 1000 Genomes Project haplotypes (Version 1, Phase 3) as the reference set. Data for seven participants were excluded from further analysis after the subjects withdrew their consent to participate. Genomic control lambda (λGC) values computed after GWAS analyses for a range of traits using the ALSPAC children data set showed minimal evidence of population stratification (λGC<1.06).

The statistical approach used for estimating SNP heritability contrasts fundamentally with that of GWAS gene discovery analyses, in which millions of SNPs are tested one by one, resulting in millions of statistical tests. A SNP-heritability analysis requires only one statistical test to be performed (a linear regression analysis). Specifically, the difference in phenotype between pairs of individuals is used as the dependent variable, and genetic relatedness between the same pairs of individuals is used as the primary explanatory variable (modeled as a random effect using restricted maximum likelihood; REML). Other explanatory variables such as age and sex can be included in the model as fixed effects. PLINK 1.90 [44] and GCTA v1.24.4 [45] were used to compute the pairwise genetic relationship matrix (for individual chromosomes and all chromosomes, respectively). There has been debate in the literature about the advantages and disadvantages of excluding SNPs in linkage disequilibrium when computing the genetic relationship matrix for SNP-heritability analysis [46-48]. In light of this debate, the genetic relationship matrix was calculated as described by Zhu et al. [49]. Thus, SNPs with an imputation reliability score <0.6, MAF <0.01, or HWE test p value <10−6 were excluded, followed by extraction of SNPs present in HapMap3, which left 1,192,778 markers for calculating genetic relatedness in a sample of 8,230 individuals. One subject from each of the 1,028 pairs within the cohort with an estimated genetic relationship >0.025 (the approximate degree of sharing expected on average for second or third cousins) was excluded, leaving 7,202 “unrelated” individuals. SNP-heritability analyses were run either with the raw refractive error readings or after transforming refractive error measurements to normal scores [34], to determine whether the non-normal distribution of the trait influenced the SNP-heritability estimate. Bivariate analyses were run to estimate the genetic correlation between all possible pairs of traits [29]. Univariate and bivariate modeling was performed with GCTA. Sex was included as a covariate; however, for the results reported age was not included since its omission had no effect on the parameter estimates.

Simulations were performed to assess the impact of refractive error measurement error due to the lack of cycloplegia. Briefly, this involved simulating a hypothetical trait such that it had a predetermined SNP-heritability level (ranging from 0.2 to 0.8 in steps of 0.1) and then assessing the SNP heritability before and after measurement error “noise” was added to the trait. The choice of the noise level is discussed in Appendix 1 and Figure 1 and Figure 2. In practice, the genotype data and genetic relatedness information for the 3,404 subjects in the age 15 years analysis were used to simulate phenotypic trait values at a given SNP-heritability level (0.2, 0.3, …, 0.8) by randomly selecting 1,000 SNPs (MAF range 0.01 to 0.50) from the genome as “causal variants,” with the effect size for each SNP drawn from a normal distribution, using LDAK [46]. In R, the standard deviation of the simulated trait was converted to 1.29 D to match the standard deviation of refractive error in 15-year-old ALSPAC participants (Table 1) and a randomly selected measurement error noise value was added, drawn from a normal distribution with a mean zero and a standard deviation of 0.50 D. GCTA was used to calculate the SNP heritability for the trait before and after noise was added. Fifty replicates were simulated at each heritability level.

Figure 1

Comparison of subjective refraction recorded by the subject’s own optometrist and non-cycloplegic autorefraction obtained at the 15-year ALSPAC research clinic. Whether the Avon Longitudinal Study of Parents and Children (ALSPAC) clinic record indicated that the subject was not wearing contact lenses during the autorefraction test is indicated by the symbol color. Note that the single outlying data point for the subjects confirmed as not wearing contact lenses corresponds to a myopic participant recorded as perfectly emmetropic by the optometrist, suggesting an error in data entry or retrieval instead of a highly inaccurate autorefraction measurement.

Figure 2

Difference between non-cycloplegic autorefraction obtained at the 15-year ALSPAC research clinic and subjective refraction recorded by the subject’s own optometrist. Two frequency distributions are plotted according to whether the Avon Longitudinal Study of Parents and Children (ALSPAC) clinic record indicated that the subject was or was not wearing contact lenses during the autorefraction test. Note the reduction in outliers and better conformity to a normal distribution when contact lens wear was excluded.

Table 1

Subject demographics by clinic visit. Values are mean ± standard deviation.

Clinic visit	N	Age (years)	Refractive error (D)	Sex (%male)
Age 7	5103	7.53±0.30	+0.19±0.88	51%
Age 10	4862	10.64±0.25	+0.06±1.10	50%
Age 11	4433	11.74±0.23	−0.04±1.12	49%
Age 12	4438	12.80±0.23	−0.17±1.17	49%
Age 15	3404	15.45±0.31	−0.40±1.29	47%

The variance in refractive error at ages 10, 11, 12, and 15 predicted (or explained) by time spent outdoors was assessed using a general linear model, using the information reported by the mother when the child was aged 8–9 (see the Phenotypes section, above). Sex was included as a covariate. The variance predicted (explained) by time spent reading was gauged using the same method.

Results

At each ALSPAC research clinic, non-cycloplegic autorefraction measurements were obtained for between 3,404 and 5,103 participants who also had genotype data available and were less related to each other than second to third cousins (Table 1). Univariate GCTA modeling suggested the SNP heritability for refractive error was 0.19 to 0.27 (weighted mean h2SNP=0.23, SE≈0.08, p<0.001) when estimated using untransformed trait values and 0.20 to 0.34 (weighted mean h2SNP=0.28, SE≈0.08, p<0.001) after the trait was transformed to a normal distribution (Table 2). These estimates therefore represent the proportion of the inter-subject variance in refractive error that can be explained by common SNPs in linkage disequilibrium (LD) with causal genetic variants that influence refractive error. Accordingly, a lower bound estimate for the heritability of refractive error—free from any bias due to sharing a common home environment—is approximately 25%.

Table 2

SNP-heritability (estimates of variance explained by SNPs) for refractive error, at a range of ages.

Research clinic visit	N	Untransformed trait			Transformed trait
		h2SNP	SE	p-value	h2SNP	SE	p-value
7	5103	0.250	0.067	9.7e-05	0.342	0.068	1.6e-07
10	4862	0.198	0.069	1.5e-03	0.204	0.069	1.1e-03
11	4433	0.194	0.075	4.0e-03	0.260	0.076	2.0e-04
12	4438	0.270	0.076	1.4e-04	0.323	0.076	5.1e-06
15	3404	0.259	0.097	2.9e-03	0.270	0.096	1.8e-03
Mean		0.233	~0.08	<0.001	0.281	~0.08	<0.001

Analyses included sex as a covariate. Estimates are shown for refractive error in Diopters (“untransformed trait”) or after refractive error readings were transformed to normal scores (“transformed trait”). P values correspond to a test of whether the heritability estimate is different from zero. The mean was calculated after weighting by sample size (N), using the formula, Σ(h2N)/ΣN

Between 2,997 and 4,238 unrelated genotyped participants underwent autorefraction at two different ages (Table 3). The genetic correlation between refractive error assessed at each age varied between 0.77 and 1.00 (SE=0.07) for untransformed trait values, and after refractive error readings were transformed into normal distributions. Estimates from closely-spaced clinic visits were generally higher than those from more widely-spaced visits, as would be expected, and genetic correlation estimates for refractive error assessed at or after the age of 10 were all ≥0.94 (Table 3). Thus, at each of the ages examined the same set of genetic variants explained about 25% of the variance in refractive error, instead of one set of genetic variants acting early in childhood and a different set of variants acting later in childhood.

Table 3

Genetic correlations between refractive error assessed at two different ages.

Age	Age 7	Age 10	Age 11	Age 12	Age 15
Age 7	x	0.94 (0.051) n=4238	0.98 (0.072) n=3837	0.91 (0.063) n=3862	0.77 (0.127) n=3007
Age 10	0.87 (0.083) n=4238	x	1.00 (0.035) n=4060	1.00 (0.037) n=4062	1.00 (0.083) n=3146
Age 11	0.89 (0.089) n=3837	1.00 (0.077) n=4060	x	0.99 (0.032) n=3878	0.96 (0.088) n=2997
Age 12	0.78 (0.089) n=3862	1.00 (0.067) n=4062	0.99 (0.058) n=3878	x	1.00 (0.055) n=3127
Age 15	0.77 (0.124) n=3007	1.00 (0.107) n=3146	0.96 (0.099) n=2997	0.94 (0.070) n=3127	X

The standard error of the genetic correlation is shown in brackets. N refers to the sample size. Analyses included sex as a covariate. Estimates are shown for refractive error in Dioptres above the diagonal, and after refractive error was transformed to a normal score below the diagonal.

To assess the impact of measurement error during autorefraction due to the lack of cycloplegia, we simulated cycloplegic and non-cycloplegic autorefraction traits for the 3,404 children included in the age 15 years analysis (Table 1). In each simulation run, 1,000 SNPs were randomly selected as “causal variants” and assigned a phenotypic effect such that the simulated trait had a predetermined h2SNP value. Then the h2SNP was calculated using GCTA before and after phenotypic measurement error noise was added at the level expected to occur because of the lack of cycloplegia (Appendix 1; Figure 1 and Figure 2). As shown in Figure 3, there was a close correspondence between the simulated and observed h2SNP before the measurement noise was added (“cyclo” autorefraction simulation in Figure 3; slope=1.059). However, the observed h2SNP was underestimated after measurement error noise was added to the trait (“non-cyclo” autorefraction simulation in Figure 3; slope=0.851). Applying a conversion factor of 1.244 (1.059/0.851=1.244) to account for the level of underestimation due to the lack of cycloplegia yielded an adjusted weighted mean h2SNP of 0.35 (SE=0.09). Since the measurement error caused by the lack of cycloplegia would be at least as bad at younger ages, and since any non-systematic source of measurement error such as hyperopes accommodating more than myopes would cause h2SNP to be further underestimated, this adjustment is likely to be conservative.

Figure 3

Comparison of SNP-heritability estimation for a simulated refractive error trait before and after the addition of measurement error “noise.” Traits were simulated using a model designed to yield a specified SNP-heritability level (x-axis). The simulated trait was either analyzed directly (“cyclo”) or after the addition of noise designed to mimic the measurement error of non-cycloplegic autorefraction (“non-cyclo”). Measurement error led to the progressive under-estimation of the true SNP heritability.

To help give context to the genetics findings, the variance explained (or predicted) by time spent reading and by time spent outdoors was also estimated (Table 4). Time spent outdoors assessed using a questionnaire completed by the mother when the participants were aged 8–9 years predicted or explained <0.5% of the variance in refractive error at older ages, while the corresponding figure for time spent reading for pleasure assessed at a similar age was <1%.

Table 4

Variance in refractive error explained by time spent reading and time spent outdoors.

Research clinic visit	N	Untransformed trait		Transformed trait
		r2	P-value	r2	P-value
Time spent reading at age 8–9 years-old
10	4616	0.003	7.7e-04	0.003	3.7e-04
11	4203	0.009	6.7e-09	0.007	4.1e-07
12	4209	0.008	1.4e-07	0.007	3.0e-07
15	3298	0.008	2.1e-06	0.005	1.9e-04
Time spent outdoors at age 8–9 years-old
10	4626	0.001	3.7e-01	0.001	9.6e-02
11	4215	0.001	4.4e-01	0.001	2.2e-01
12	4225	0.003	4.6e-03	0.002	1.4e-02
15	3298	0.002	2.1e-02	0.001	3.3e-02

Analyses included sex as a covariate. Estimates are shown for refractive error in Diopters (“untransformed trait”) or after refractive error readings were transformed to normal scores (“transformed trait”).

Discussion

The SNP-heritability estimate of approximately 35% obtained here gives by far the most compelling support to date for a polygenic inheritance component for refractive error. Unlike prior heritability studies, the present estimate is free from assumptions about, and confounded by, shared environment effects, for example, parents and their children, or siblings, sharing the same home environment.

What is the relationship between SNP heritability (h2SNP), broad sense heritability (H2), and narrow sense heritability (h2)? Simulation studies and GWAS-GCTA analyses with large sample sizes for traits such as height [27,50] have helped clarify this issue. h2SNP gives a lower-bound estimate of h2. Thus, taking a figure of 54% (95% confidence interval [CI] 36 to 72%) for the narrow sense heritability of refractive error (calculated for a sample of 890 ALSPAC participants who had refraction information available for a sibling [51]) the gap between our h2SNP estimate of 35% and the h2 figure of 54% is likely due, first, to upward bias in the estimation of h2 from shared family environment effects and, second, to two sources of downward bias in h2SNP: namely, the limited ability of genetic relatedness indices to capture the additive genetic effects of rare variants, and statistical “noise” associated with a modest sample size. For height, studies in larger samples and employing higher-density genotyping platforms have shown that h2SNP steadily increases toward h2 [50]. The gap between h2 and H2 is generally attributed to non-additive sources of genetic variation, which are included in H2 but absent from h2. Unlike h2SNP, estimates of h2 and H2 are likely to be inflated due to misallocation of unique or common environment variance. Although our results do not allow us to address the extent of this upward bias (only giving a lower bound estimate for h2 and H2), the large gap between the narrow sense heritability and the broad sense heritability estimates for refractive error (e.g., 15–70% and 75–90%, respectively [25]) implies that dominance and/or interaction effects are likely to be important in refractive development. As with other types of heritability estimate, it is not strictly appropriate to compare SNP-heritability estimates across populations and/or across generations, since sources of variance may not be constant in different settings.

Our h2SNP estimate of about 35% provides an approximate upper limit for the capacity to predict the future refractive error of premyopic children using existing genotype platforms and imputation methods alone. Technical advances and larger sample sizes would be expected to raise h2SNP and improve the ability to predict children at risk of myopia. Inclusion of gene–gene and/or gene–environment interaction effects would be expected to further improve predictive capacity. Notably, recent studies employing Bayesian or mixed models, similar in spirit to GCTA, have demonstrated the potential of genetic prediction from genome-wide common SNPs, which means that identifying individual risk SNPs through GWAS is not a prerequisite for effectively predicting disease [52].

Crudely collected information on two key environmental risk factors for myopia, time spent outdoors and time reading for pleasure, had limited capacity to predict or explain future refractive error, consistent with previous reports [53,54]. However, our finding that these exposures explained <1% of the variance in refractive error likely underestimates the true effects for several reasons. First, the parental questionnaire used to collect the information had a categorical scoring scale rather than a continuous scale, which would have weakened statistical power to predict fine-scale differences. Second, exposure at only one particular age was examined, which may not have been representative of the true exposure. Third, for some participants the level of exposure may have been reported inaccurately. Finally, there may be nonlinear relationships between the time children engage in these activities and refractive error development, such as threshold effects, which would not have been taken into account when we assessed the variance they explained.

Our study had two main limitations. First, the assessment of refractive error at ALSPAC research clinics was done without cycloplegia, which is known to produce a subject-specific measurement error that can either under- or overestimate true refractive error. Lack of cycloplegia will therefore have reduced phenotypic accuracy, especially at younger ages and in hyperopes, and made our SNP-heritability estimates overly conservative (Appendix 1). Note that any systematic “over minus” component of the measurement error due to lack of cycloplegia during autorefraction would have had no impact, since heritability is a measure of explained variance, i.e., between-subject differences, whereas a systematic effect relates to the mean phenotype in the population. Second, the latest age at which ALSPAC participants’ refractive error was assessed was 15-years-old, which is before the age at which myopia typically stabilizes [55]. Therefore, further work will be needed to discover whether the set of genetic variants found here to influence refractive development over the 7- to 15-year-old age range continues to exert an effect into adulthood. Such information will be valuable for genetic prediction and for understanding the extent to which genetically-determined causal mechanisms are shared between early-onset and later-onset myopia.

Conclusion

The SNP heritability for refractive error was estimated to be approximately 35%. This value—calculated using only information about refractive error and an index of genetic similarity between each pair of individuals in the sample—provides the strongest evidence to date that refractive error has an important genetic component. In comparison, the two main environmental risk factors identified to date for myopia—time spent reading and time spent outdoors—together explained or predicted about 1% of the variance in refractive error. Nevertheless, the SNP-heritability figure is lower than that for axial length and corneal curvature (h2SNP about 40–50%) calculated in the same population [29] suggesting that environmental effects exert more influence in determining refractive error than they do for the individual ocular components. The same set of genetic variants was found to influence refractive error at the various ages examined, rather than different sets of genetic variants acting at specific ages.