**Purpose.**:
To test the hypothesis that there is a major genetic determinant of vertical disc diameter (VDD) and vertical cup-to-disc ratio (VCDR) in a large, population-based sample.

**Methods.**:
Data were collected from 3654 individuals, 49 years of age or older, participating in the Blue Mountains Eye Study. VDD and VCDR were determined from stereo optic disc photographs. Commingling analyses in SKUDRIVER/SKUMIX were performed in nonglaucomatous eyes to investigate whether the observed VDD and VCDR data were best described by a one-, two-, or three-distribution model.

**Results.**:
VDD data did not show evidence of commingling. After adjustment for the effects of age, VDD and intraocular pressure, the best model for VCDR consisted of a mixture of three distributions in Hardy-Weinberg equilibrium. The proportion of the variance in VCDR explained by this mixing component was 0.58.

**Conclusions.**:
Findings from this study are consistent with the presence of a major gene that accounts for 58% of the variance in VCDR. These results strongly support further efforts to identify the genetic variants responsible for this quantitative trait, which is a key constituent of the phenotype of primary open-angle glaucoma (POAG).

^{ 1–3 }In the majority of cases, POAG is inherited as a complex disease: it is assumed to result from many interactive genetic and environmental factors, none of which individually is necessary or sufficient to cause the disease. Although high-throughput genotyping technologies are becoming increasingly feasible and affordable, and underlying methodologies to unravel complex diseases are developing rapidly, the multifactorial etiology of POAG is still proving a hard nut to crack. More than 25 chromosomal regions have been linked to the disease, but only three genes (

*MYOC*,

^{ 4 }

*OPTN*,

^{ 5 }and

*WDR36*

^{ 6 }) have been identified. These genes most likely contribute to the pathogenesis of POAG in less than 5% of cases in the general population.

^{ 7–10 }Genes accounting for a greater proportion of the known heritable component of POAG thus remain to be identified.

^{ 11–13 }Heritability estimates of VCDR ranged from 0.48 to 0.65.

^{ 12–14 }These estimates were based on family studies and provided information on the additive effects of all involved genes. Although a high heritability itself may already be promising for gene-finding studies, it would be interesting to know whether the heritable component solely involves genes of small effect or also includes one or more genes that have a relatively large effect. The latter would be easier to detect in gene-finding studies, and their presence would therefore even more resolutely support quantitative trait–based strategies.

^{ 15 }a form of model fitting that employs the method of maximum likelihood.

^{ 16 }Commingling analysis investigates the strength of evidence for a single gene of major effect and provides an estimate of the locus-specific heritability, which is the proportion of the total phenotypic variance explained by the effect of the major gene.

^{ 17 }In brief, all permanent noninstitutionalized residents 49 years of age or older were invited to participate. Of the 4433 eligible individuals, 3654 (82.4%) attended baseline eye examinations between 1992 and 1994. Of the 779 nonparticipants, 501 (11.3%) refused, 68 (1.5%) had died, and 210 (4.8%) had moved away from the area. The response rate compares well with the best population-based research in glaucoma.

^{ 18–20 }

^{ 21 }and measurement of intraocular pressure (IOP) with Goldmann applanation tonometry. Visual fields were initially assessed with a 30° suprathreshold screening test (Humphrey 76-point test; Carl Zeiss Meditec, Inc., Dublin, CA). Full-threshold 30-2 visual field tests of each eye were subsequently performed in subjects with suspected glaucoma.

^{ 22 }The vertical disc diameter (VDD) was measured to the nearest 0.01 mm as the longest diameter between the inner limits of the scleral ring in a range between clock hours 11 to 1 and 5 to 7. The optic cup was determined by its contour, with the outer margin taken to be the point where the wall met the plane of the disc surface at the level of the scleral ring. The vertical cup-disc ratio (VCDR) was calculated from the disc and cup measurements. Optic disc measurements were corrected for the magnification effect of the eye-camera system according to spherical equivalent refraction, as described by Bengtsson and Krakau.

^{ 23 }All photographs were graded by one or both of two trained graders. The chief investigator (PM) adjudicated discrepancies. Interobserver variability was assessed in a masked fashion in a random sample of 100 optic discs and was in the excellent agreement range.

^{ 24 }

*n*= 108). If only one eye of a subject was phakic, this eye rather than its nonphakic fellow was considered for analysis. If both eyes were phakic, one eye was chosen at random for inclusion in the analysis. Eyes with tilted optic discs (

*n*= 78) or with other disc anomalies, such as colobomata (

*n*= 1), disc drusen (

*n*= 1), or optic atrophy (

*n*= 1), were excluded from the dataset. A further 16 eyes were excluded because of high myopia (spherical equivalent greater than −8 D).

^{ 25 }The main analyses in this study were performed on a “normal” population, which excluded patients with glaucoma in either eye (

*n*= 90). The diagnosis of glaucoma was made on the basis of typical glaucomatous visual fields loss on the Humphrey 30-2 test, combined with matching optic disc rim thinning, as described previously.

^{ 24 }For the analyses of VDD, 86 eyes were excluded because no gradable optic disc photographs were available, and for the analyses of VCDR, an additional 5 eyes were excluded because covariate data were incomplete, leaving valid data for 3273 and 3268 subjects for the analyses of VDD and VCDR, respectively.

*q*; genotypic means of

*m*

_{1},

*m*

_{2}, and

*m*

_{3}; and within-genotype variance

*s*

^{2}, a likelihood function

*L*for an individual observation is defined under Hardy-Weinberg equilibrium as

^{ 26 }: where

*x*is the observed trait value of a randomly ascertained member of the population and

*f*(

*x*;

*m*,

*s*

^{2}) is a normal density function with mean

*m*and variance

*s*

^{2}. The overall likelihood of this mixture model is computed as the product of the likelihoods of the individual observations. The maximum likelihood of the model may be compared with that of the null model, which consists of a single normal distribution, for a test of the major gene effect.

^{ 15 }Both programs are available at http://statgen.iop.kcl.ac.uk/skudriver/. The commingling analysis was performed on the adjusted and standardized VDD and VCDR data, and comprised maximum likelihood estimation for each of three models: single-distribution (the null model), two-distribution (i.e., fully dominant or recessive) and three-distribution model.

^{2}, 2Q(1 − Q), and Q

^{2}. However, the program also allows deviation from Hardy-Weinberg equilibrium by introducing an inbreeding coefficient F, so that the proportions within the distributions become (1 – Q)

^{2}+ FQ(1 − Q), 2Q(1 − Q)(1 − F), and Q

^{2}+ FQ(1 − Q).

^{ 27 }

*y*=

*R*/

*P*[(

*x*/

*R*+ 1)

*– 1], where*

^{P}*R*is chosen such that every

*x*/

*R*+ 1 is positive in the sample and

*P*is optimized as part of the maximum-likelihood estimation. This method allows the fit of multiple distributions to be assessed after skewness has been removed, which is important, since skewness in itself may lead to the mistaken conclusion that more than one distribution is present.

^{ 15 }Significant skewness may be tested for by a likelihood ratio test comparing a model in which

*P*is fixed to a value of 1 (untransformed model) with a corresponding model in which

*P*is not constrained (transformed model).

^{ 28 }The SKUMIX program provided a measure of the goodness of fit for the one-, two-, and three-distribution models, both with and without a power transformation, expressed as minus twice the logarithm of the likelihood (−2 log

*L*). Hypothesis testing was achieved by referring the difference in this quantity between two models to a χ

^{2}distribution with degrees of freedom equal to the difference in the number of free parameters. Since multiple comparisons were made, each probability was corrected with a Bonferroni correction, to avoid spuriously significant results.

^{ 29 }The best-fitting model was chosen according to the Akaike Information Criterion (AIC), defined as −2 log

*L*+ twice the number of free parameters.

^{ 30 }The AIC penalizes for adding free parameters and thus selects the most parsimonious model that fits the data well. The Akaike weight (

*w*) was used to assess model selection uncertainty.

^{ 31 }It represents the probability that the model is the best among the whole set of models.

*P*< 0.001), VDD (B = 0.003,

*P*< 0.001) and IOP (B = 0.004,

*P*< 0.001). The standardized residuals of this multivariate regression model were used for further analysis. The distribution of the standardized and adjusted VCDR data had a skewness of −0.068 and a kurtosis of −0.097.

Characteristic | |
---|---|

Total number | 3273 |

Age (y), mean ± SD | 65.5 ± 9.4 |

Age 49–59 y, n (%) | 964 (29.5) |

Age 60–69 y, n (%) | 1227 (37.5) |

Age 70–79 y, n (%) | 832 (25.4) |

Age 80+ y, n (%) | 250 (7.6) |

Male, n (%) | 1438 (43.9) |

White, n (%) | 3248 (99.3) |

VDD (mm), mean ± SD | 1.51 ± 0.17 |

VCDR, mean ± SD | 0.43 ± 0.14 |

Intraocular pressure (mm Hg), mean ± SD | 16.0 ± 2.7 |

^{2}

_{1}= 9.16,

*P*= 0.027 after Bonferroni correction). This is reflected by the findings that the skewness of the untransformed data was 0.191, whereas after the power transformation (with

*R*fixed as 11.0 and

*P*optimized by SKUMIX/SKUDRIVER as 0.35) the skewness was 0.002. The power transform did not significantly improve the fit of the data when two or three distributions were specified. For the untransformed data, the two-distribution model fitted the data significantly better than the one-distribution model (χ

^{2}

_{2}= 10.92,

*P*= 0.047 after Bonferroni correction), but the three-distribution model did not fit better than the two-distribution model. Considering the transformed models only, neither dataset provided evidence of commingling. Allowing the inbreeding coefficient (

*F*) to vary (i.e., allowing departure from Hardy-Weinberg equilibrium) did not significantly improve the fit of any of the commingled models. According to the AIC, the most conservative and parsimonious model that fitted the data well was the one-distribution transformed model. The Akaike weights show that this model was only 1.12 (0.29/0.26) times more likely than the two-distribution untransformed model to be the best, indicating a considerable degree of model-selection uncertainty.

Model | −2 × Log Likelihood + Constant | χ^{2} (df) Compared with Untransformed* | χ^{2} (df) Compared with 1 Distribution† | χ^{2} (df) Compared with 2 Distributions‡ | Akaike Information Criterion§ | Akaike Weight‖ |
---|---|---|---|---|---|---|

1 Distribution, untrans | 4643.69 | 4647.69 | 0.01 | |||

1 Distribution, trans | 4634.52 | 9.16 (1)¶ | 4640.52 | 0.29 | ||

2 Distribution, untrans | 4632.77 | 10.92 (2)¶ | 4640.77 | 0.26 | ||

2 Distribution, trans | 4631.11 | 1.66 (1) | 3.42 (2) | 4641.11 | 0.22 | |

3 Distribution, untrans | 4631.87 | 11.82 (3) | 0.90 (1) | 4641.87 | 0.15 | |

3 Distribution, trans | 4631.08 | 0.79 (1) | 3.44 (3) | 0.02 (1) | 4643.08 | 0.08 |

*P*−0.57, power transform variable

*R*(fixed) 11.0, inbreeding coefficient (fixed) 0. These parameters gave rise to the distributions shown in Figure 3. When the back-transformed, unstandardized regression residuals were considered, the middle distribution (which represents individuals carrying 1 copy of the rare allele) had a mean of −0.15, and the leftmost distribution (which represents individuals carrying two copies of the rare allele) had a mean of −0.29. As the total variance was 1.05, the residual variance of 0.44 implied that the variance due to the commingling was 0.61 and the locus-specific heritability was 0.58.

Model | −2 × Log Likelihood + Constant | χ^{2} (df) Compared with Untransformed* | χ^{2} (df) Compared with 1 Distribution† | χ^{2} (df) Compared with 2 Distributions‡ | Akaike Information Criterion§ | Akaike Weight‖ |
---|---|---|---|---|---|---|

1 distribution, untrans | 4635.09 | 4639.09 | 0.10 | |||

1 distribution, trans | 4633.91 | 1.18 (1) | 4639.91 | 0.07 | ||

2 distribution, untrans | 4630.03 | 5.06 (2) | 4638.03 | 0.17 | ||

2 distribution, trans | 4627.84 | 2.19 (1) | 6.07 (2) | 4637.84 | 0.19 | |

3 distribution, untrans | 4627.55 | 7.54 (3) | 2.48 (1) | 4637.55 | 0.22 | |

3 distribution, trans | 4625.33 | 2.23 (1) | 8.58 (3) | 2.51 (1) | 4637.33 | 0.25 |

^{ 12–14 }Hardy-Weinberg proportions are respected, and environmental factors of major effect are unknown, it is likely that a gene of major effect explains this commingling. Second, the exclusion of patients with glaucoma is a possible source of selection bias. The exclusion was necessary because we do not know whether the same processes are responsible for VCDR in healthy and glaucomatous eyes. When patients with glaucoma were included, commingling analysis of VCDR resulted in a three-distribution model with dominance −0.59, displacement 2.05, and allele frequency 0.17. The rightmost distribution in this model (

*n*= 93; SD = 0.75; mean = 2.05, corresponding to 0.27 unstandardized residuals) was very similar to the distribution of the included glaucoma population (

*n*= 83; SD = 1.05; mean = 2.04, corresponding to 0.27 unstandardized residuals). This result may indicate that patients with glaucoma form a separate distribution and that the SKUMIX program is able to correctly disentangle this admixture. However, a major genetic origin of the commingling in this heterogeneous population could be disputed. The results for VDD did not change after the glaucoma cases were included.

^{ 11–13 }Our result may be explained by the conservative design of the SKUMIX program, which implements the commingling analysis. This design, which “had to guard against claiming separate distributions where none exist”

^{ 15 }has been tested experimentally.

^{ 32 }When both commingling and segregation analyses were applied to simulated pedigree data in which a major locus was segregating, more than 20% of the samples provided evidence of segregation of a single locus but not of commingling. Our results for VDD therefore do not preclude a genetic determinant of major effect. This finding is also suggested by the Akaike weights, which show a considerable model selection uncertainty and provide some support for a two-distribution model (fully dominant or recessive gene) as well. Another possible explanation of the discrepancy with previous heritability studies is that the latter assessed the additive effects of all involved genes. A collection of several loci with small effects rather than a single major gene determining VDD may lead to high heritability estimates without evidence of commingling.

^{ 27,33 }

^{ 12–14 }Our estimate of 58% for its locus-specific heritability suggested that an important part of this additive genetic variance may be attributable to the effect of a single locus. Moreover, our study provided a model elucidating the allele frequencies, dominance, and displacement associated with this locus.

^{ 34–37 }However, one might question whether a large VCDR in these studies was an actual risk factor or rather an early sign of POAG, and consequently, whether the smaller VCDR associated with the major locus in our study would actually reduce POAG risk. Identifying the gene and exploring its function could shed light on this issue.

^{ 38,39 }If a dichotomous trait–based association analysis were to be considered, it would be desirable to compare individuals having at least one copy of the rare allele with individuals having no copies. The former group would consist of individuals with VCDRs smaller than the lower extreme of the rightmost distribution, say those with VCDR values less than three residual standard deviations from the mean of this distribution; that is, standardized VCDR regression residuals of less than −0.07 − 3 · 0.66 = −2.05. This translates to VCDRs that are at least 0.26 smaller than would be expected based on age, VDD, and IOP. The second group would be those with VCDRs greater than the upper extreme of the middle distribution; that is, those with standardized VCDR regression residuals greater than −1.32 + 3 · 0.66 = 0.66, corresponding to VCDR values of at least 0.08 greater than predicted from the covariates.