Commingling Analyses of Central Corneal Thickness and Adjusted Intraocular Pressure in an Older Australian Population

Abstract

Purpose.: To test the hypothesis that there is a major genetic determinant for central corneal thickness (CCT) in a large, population-based sample and to test whether such a determinant accounts for results previously ascribed to intraocular pressure (IOP).

Methods.: Measurements of IOP and CCT were collected for 1356 individuals, 49 years of age or older, participating in the Blue Mountains Eye Study (third ascertainment call). Linear regression methods were used to adjust CCT for age. IOP was adjusted for CCT. A commingling analysis was performed with a C++ program, SKUDRIVER, to investigate whether the observed CCT and IOP data, adjusted for covariates, were best described by a one-, two-, or three-distribution model. The fitted models were compared by Akaike information criterion (AIC) values.

Results.: Significant skewness was present in CCT and IOP distributions. The most parsimonious model for age-adjusted CCT was a single-distribution model. For IOP adjusted for CCT the most parsimonious model was a mixture of three distributions with means corresponding to 16.02, 23.36, and 35.59 mm Hg. The proportion of total IOP variance attributable to these distributions was 18.3%. Model parameters were almost identical with those in which IOP was not adjusted for CCT. None of the analyses suggested deviation from Hardy-Weinberg equilibrium.

Conclusions.: The findings revealed no evidence of a single major genetic determinant for CCT. Previous results suggesting the presence of a major locus for IOP do not appear to be confounded by the influence of CCT on IOP.

Primary open-angle glaucoma (POAG) is a progressive optic neuropathy with an established genetic component in its origin. ^1–3 In most 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 more than 25 chromosomal regions have been linked to the disease, only three genes (MYOC, ⁴ OPTN, ⁵ and WDR36 ⁶ ) have been identified. It is likely that these genes 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. 

The etiologic complexity of POAG can be reduced by separately studying quantitative features of the phenotype, such as intraocular pressure (IOP). Elucidating the genetic determinants of this quantitative feature in healthy eyes may improve our understanding of the damage to the optic disc in glaucomatous eyes. Quantitative traits are likely to be more powerful in detecting new genes than the dichotomous POAG trait. ^11,12 They may also have simpler genetic backgrounds, can be studied population-wide, and are less prone to misclassification. 

Commingling analysis as applied to traits of interest in glaucoma has been described in detail elsewhere. ^13,14 In brief, the analysis investigates whether the observed distribution of a quantitative trait is best modeled by a single distribution or by an admixture of multiple distributions. A single distribution suggests the lack of a single major genetic determinant for the trait, whereas multiple distributions suggest the presence of a major genetic determinant, especially if the distributions are in Hardy-Weinberg equilibrium (HWE) for population-ascertained data. 

IOP largely satisfies the criteria to be considered as an endophenotype ¹⁵ for POAG. It is associated with illness in the population, it is heritable, it is present in an individual whether or not disease is present, and it cosegregates, to an extent, with the disease. 

Elevated IOP is a major risk factor for POAG. Although only a small number of patients with OHT develop POAG every year, studies on the conversion from OHT to POAG have demonstrated that the overall risk of developing POAG is approximately five times higher in subjects with IOPs greater than 21 mm Hg than in subjects with lower IOPs ¹⁶ and that the higher the IOP at screening, the greater the risk of POAG. ¹⁷  

Furthermore, IOP is currently the only modifiable risk factor for POAG. The Ocular Hypertension Treatment Study demonstrated that the therapeutic lowering of IOP reduces the risk of conversion from ocular hypertension to glaucoma. ¹⁸  

IOP has been demonstrated to be a heritable trait, ^19–22 and a commingling analysis ¹⁴ has suggested the presence of a major gene accounting for 18% of the variance of IOP in a population. Commingling analysis investigates the strength of evidence that a single gene has a major effect and provides an estimate of the locus-specific heritability. 

However, none of these studies has accounted for the possible influence of central corneal thickness (CCT) on the results. CCT is known to affect tonometric measurements of IOP ^23–35 and has also been proposed as an independent risk factor for POAG. ^36–38 Furthermore, several studies in different ethnic groups have demonstrated the high heritability of CCT. ^39–41  

The search for quantitative trait loci (QTLs) and genes influencing IOP would be greatly aided by studies explicitly addressing the effect of CCT. For this reason, we performed commingling studies on IOP and CCT data from the Australian Blue Mountains Eye Study cohort. We sought to investigate the likelihood that there is a single major genetic determinant of CCT and to determine whether a previous genetic model ascribed to IOP ¹⁴ is robust to adjustment of IOP for CCT, rather than erroneously ascribing results to IOP that are actually due to the effects of CCT, through the surrogate measure of IOP. 

Methods

Study Population

The sample was a part of the Blue Mountains Eye Study (BMES), a population-based survey of vision and common eye diseases in the Blue Mountains region west of Sydney, Australia. The study adhered to the tenets of the Declaration of Helsinki and was approved by the Western Sydney Area Health Service Human Ethics Committee. Written, informed consent was obtained from all participants. The population has been described in detail previously. ⁴² 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. Among the 779 nonparticipants, 501 (11.3%) persons 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. ^43–45 The BMES examined these 3654 patients aged 49+ years during 1992 to 1994 (baseline-known as BMES-1), 2335 (75% of 3111 survivors) at a 5-year examination (during 1997 to 1999 known as BMES-2), and 1935 (75% of 2581 survivors) at a 10-year examination (during 2002 to 2004, known as BMES-3). 

Phenotyping

CCT was measured in 1356 individuals with ultrasonic pachymetry at the 10-year examination (Healey PR, et al. IOVS 2005;46:ARVO E-Abstract 3520). The disparity between the total number of participants in BMES-3 and those with recorded values of CCT arose due to instrumental difficulties at the time of CCT measurement. The cases in which CCT measurements were not performed did not show systematic bias compared with those in which CCT was measured. Three measurements were recorded for corneal thickness of each eye. The left eyes of 6 individuals and right eyes of 10 individuals were excluded from the study, as data for CCT were missing from their phenotypic data. This data deficit reduced the sample sizes to 1350 and 1346 for the left and right eye CCT, respectively. CCT commingling analyses were performed for both eyes separately, both with and without adjustments for covariates. 

IOP measurements were recorded by using a Goldmann applanation tonometer, as described elsewhere. ¹⁴ In their study, Viswanathan et al. ⁴⁶ took a conservative approach: Only one eye of each subject was considered for analysis because of the problems relating to intereye correlation. To enable a comparison of findings between this work and the present study, we selected cases with CCT measurements from BMES-3 from the dataset used for the previous IOP commingling study. Hence, the sample size was reduced to 1305 for IOP commingling analyses with or without adjustment for CCT. 

Statistical Analysis

Before the commingling study, univariate and multivariate linear regression analyses were performed to detect whether any adjustments were needed for the effects of explanatory covariates. Adjustment was made, if necessary, by analyzing the residuals after linear regression. 

The software used to implement the commingling analysis was the C++ program SKUDRIVER, written by one of the authors (ACV), and the program SKUMIX. Both programs are freely available (http://statgen.iop.kcl.ac.uk/skudriver/ provided by the Statistical Genetics Unite, Kings College London, UK). SKUDRIVER takes as input a user-specified range of starting values for each of the following variables (see also Fig. 1): within-genotype variance (V), homozygote mean (U), dominance (D), displacement (T), allele frequency (Q), power transform variables (P and R), and inbreeding coefficient (F). Displacement (T) is defined as the difference between the mean values of the two homozygous distributions. Dominance (D) represents the mean value of the heterozygous distribution relative to the two homozygotes. Thus, the three genotypic means are at U, U+DT, and U+T. Since the input parameters in SKUDRIVER can be specified as either fixed or estimated, the user may constrain the model to a single distribution by fixing the value of T as 0, or may specify a two-distribution model by fixing the value of D as 0 or 1. Q is assigned to be the frequency of the allele associated with the displaced distribution so that in the three-distribution model, under HWE, the proportions of the population within each of the distributions are (1 – Q)², 2Q(1 – Q), and Q ². However, the program also allows deviation from HWE by introducing an inbreeding coefficient F, so that the proportions within the distributions become (1 – Q)² + FQ(1 – Q), 2Q(1 – Q)(1 – F), and Q ² + FQ(1 – Q). 

Figure 1.

 

Histogram of IOP: x-axis, IOP in mm Hg; y-axis, number of persons.

One of the important features of the software is the facility to specify starting values P and R of a power transformation to reduce skewness of the form y = R/P [(x/R + 1) ^P − 1], where 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. ⁴⁷ Significant skewness may be tested with 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). 

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 logL). The best-fitting model was chosen according to the Akaike information criterion (AIC), defined as −2 logL + twice the number of free parameters. ⁴⁸ The AIC penalizes for adding free parameters and thus selects the most parsimonious model that fits the data well. 

Results

Central Corneal Thickness

The sample had a mean of 538.92 μm for left eye CCT with SD 33.72 μm and for right eye 540.14 μm with SD 34.27 μm. Linear regression–based methods found that the CCTs of both eyes were associated with age, and so CCT values were adjusted for this covariate. Stepwise regression was performed to find the best-fitting models of CCT for each eye. Models were chosen on the basis of adjusted R ². Both unadjusted and age-adjusted CCT were standardized before analysis with SKUDRIVER/SKUMIX. The results of the commingling analysis on the unadjusted left eye CCT are shown in Table 1. 

Table 1.

 

Commingling Analysis of Unadjusted CCT Data

Model	V	U	D	T	Q	P	−2 LL + Constant	AIC
1 Distribution, untransformed	0.999	0.000	—	—	—	1	1915.07	1919.07
1 Distribution, transformed	0.999	−0.032	—	—	—	0.297	1911.05	1917.05
2 Distribution, untransformed	0.969	0.000	0	3.341	0.052	1	1911.93	1919.93
2 Distribution, transformed	0.696	−0.050	0	1.162	0.805	−0.102	1909.96	1919.96
3 Distribution, untransformed	0.894	0.000	0.272	3.551	0.049	1	1911.39	1921.39
3 Distribution, transformed	0.696	0.050	1.000	1.159	0.405	−0.103	1909.96	1921.96

  n = 1350. Best-fitting model shown in bold.

In Table 1, when one distribution was specified, significant skewness was removed by the power transformation (χ² = 4.02, P < 0.05). This process is reflected by the findings that the skewness of the untransformed left eye CCT data were 0.186, whereas after the power transformation (with a value of R in the power transformation chosen as 11.0 and a value of P optimized by SKUDRIVER/SKUMIX as 0.297), it was 0.057. The power transform had no significant effects in the two- and three-distribution models. 

In Table 1, the Akaike information criterion (AIC) shows that the transformed one-distribution model is the best-fitting model (minimum AIC value, 1917.05) for unadjusted left eye CCT data. Allowing the inbreeding coefficient (F) to vary (i.e., allowing departure from Hardy-Weinberg proportions) did not significantly improve the fit of any of the non-null models. For left eye CCT, the one-distribution model fit better than any two- or three-distribution model. 

Table 2 shows similar findings for age-adjusted left eye CCT data. Age adjustment of left eye CCT had no remarkable effect on the commingling findings, as the transformed one-distribution model is still the best-fitting model (minimum AIC value, 1917.00) for left eye CCT data. For left eye CCT the one-distribution model fit better than any two- or three-distribution model either adjusted or unadjusted for age. 

Table 2.

 

Commingling Analysis of Age-Adjusted Left Eye CCT Data

Model	V	U	D	T	Q	P	−2 LL + Constant	AIC
1 Distribution, untransformed	0.999	0.000	—	—	—	1	1915.06	1919.06
1 Distribution, transformed	0.999	−0.032	—	—	—	0.302	1911.00	1917.00
2 Distribution, untransformed	0.969	0.000	0	3.582	0.046	1	1911.16	1919.16
2 Distribution, transformed	0.977	−0.021	0	3.440	0.042	0.544	1909.83	1919.83
3 Distribution, untransformed	0.896	0.003	0.255	3.805	0.046	1	1910.63	1920.63
3 Distribution, transformed	0.974	−0.021	0.057	3.459	0.042	0.546	1909.83	1921.83

  n = 1350. Best-fitting model shown in bold.

Similar findings were obtained applying the commingling method to right eye CCT, either unadjusted or adjusted for age. Therefore, it is concluded that the CCT data sets are explained by a single normal distribution better than any other fitted models. This implies that a major locus is unlikely to determine the CCT. 

Intraocular Pressure

The mean value of IOP in the sample was 16.02 mm Hg with SD 2.64. The skewness for this data set was 0.711. Regression-based methods showed that IOP was associated with CCT (Pearson's r = 0.11, P < 0.01). Histograms of IOP and CCT for this sample are shown in Figures 1 and 2, respectively. 

Figure 2.

 

Histogram of CCT: x-axis, CCT in micrometers; y-axis, number of persons.

The commingling analysis was performed for the IOP dataset both unadjusted and adjusted for CCT. 

As shown in Tables 3 and 4, AIC values suggest that a mixture of three distributions is the best-fitting model for both the unadjusted and adjusted IOP data sets. When the one-distribution model was specified, substantial skewness was removed by the power transformation (χ² = 27.86, P < 0.000). However, a power transformation had a negligible effect in the two- and three-distribution models. In the latter model, the value of the P parameter in the power transformation was optimized to 1.139. Allowing the inbreeding coefficient (F) to vary did not significantly improve the fit of any of the non-null models, which suggests that HWE was maintained. 

Table 3.

 

Commingling Analysis of Unadjusted IOP Data

Model	V	U	D	T	Q	P	−2 LL + Constant	AIC
1 Distribution, untransformed	0.999	0.00025	—	—	—	1	1851.03	1855.03
1 Distribution, transformed	0.969	−0.06750	—	—	—	−0.529	1823.17	1829.17
2 Distribution, untransformed	0.832	0.00026	0	3.415	0.121	1	1809.93	1817.93
2 Distribution, transformed	0.847	−0.0207	0	3.275	0.112	0.529	1808.21	1818.21
3 Distribution, untransformed	0.788	−0.00087	0.377	7.237	0.013	1	1804.56	1814.56
3 Distribution, transformed	0.781	0.00533	0.369	7.514	0.013	1.139	1804.45	1816.45

  n = 1305. Best-fitting model shown in bold.

Table 4.

 

Commingling Analysis of CCT–Adjusted IOP Data

Model	V	U	D	T	Q	P	−2 LL + Constant	AIC
1 Distribution, untransformed	0.999	0.00000	—	—	—	1	1851.21	1855.21
1 Distribution, transformed	0.966	−0.06970	—	—	—	−0.576	1820.95	1826.95
2 Distribution, untransformed	0.841	−0.00000	0	3.711	0.108	1	1807.25	1815.25
2 Distribution, transformed	0.851	−0.02327	0	3.455	0.104	0.477	1804.98	1814.98
3 Distribution, untransformed	0.791	−0.00116	0.376	7.406	0.012	1	1802.26	1812.26
3 Distribution, transformed	0.614	−0.12159	−0.691	2.662	0.060	−1.536	1800.99	1812.99

  n = 1305. Best-fitting model shown in bold.

The best-fitting model was the three-distribution one (minimum AIC value, 1814.56, in Table 3) in the light of AIC values for unadjusted IOP data without transformation (P = 1). The parameters of this model were: residual variance, 0.788; homozygote mean, –0.00087; dominance, 0.377; displacement, 7.237; allele frequency, 0.013; power transform variable P (fixed), 1.0; and power transform variable R (fixed), 11.0. As the total variance was 0.976 and residual variance 0.788, the variance due to genotypic means was 0.188. Thus, the percentage of variation explained by the major effect was 18.8%. The estimate of Q as 0.013 corresponds to probabilities of 97.52%, 2.47%, and 0.01% of an individual in the population having genotype AA, Aa, or aa, respectively. On the original, unstandardized IOP scale the three component means correspond to 16.01, 23.20, and 35.12 mm Hg, respectively. 

In the light of the AIC criterion, the best-fitting model for IOP data adjusted for CCT was the three-distribution model (minimum AIC value, 1812.26, in Table 4). The parameters of this model were residual variance, 0.791; homozygous mean, −0.00116; dominance, 0.376; displacement, 7.406; allele frequency, 0.012; power transform variable P (fixed), 1.0; and power transform variable R (fixed), 11.0. As the total variance was 0.974 and residual variance 0.791, the variance due to genotypic means (i.e., mixing component) was 0.183, which implies that the 18.3% variation of IOP adjusted for CCT can be explained by the major effect. The estimate of Q as 0.012 corresponds to 97.67%, 2.31%, and 0.02% probabilities that an individual in the population had genotype AA, Aa, and aa, respectively. On the original IOP scale, for a person with a CCT of 540 μm, three component means correspond to 16.02, 23.36, and 35.59 mm Hg, respectively. 

Allowing the inbreeding coefficient (F) to vary did not significantly improve the fit any of the non-null models, which suggests that HWE was maintained. So, for both the IOP data sets (either unadjusted or adjusted for CCT) the three-distribution model in HWE without transformation was the best model to describe the IOP data. The values of the best parameters were almost identical between adjusted and unadjusted datasets. 

Discussion

The most parsimonious model for the CCT data set from the Blue Mountains Eye Study was the one-distribution model, which implies that there is no single major genetic determinant of CCT. A possible objection to this conclusion is that the analysis included an insufficient number of individuals to support a negative result. However, previous simulation work suggests that the current analysis is adequately powered ⁴⁹ and negative results have been postulated with smaller samples. ⁵⁰  

It is not unusual for a trait, even though highly heritable like CCT, to be found to be influenced by numerous genetic determinants of relatively small effect. A recent example is height. ⁵¹ Findings such as these are useful to elucidate the biology underlying the trait, but are of less use in risk stratification or the search for therapeutic targets. Our results suggest that the heritable component of CCT is determined by multiple genes of relatively small effect and various effect sizes. One of the limitations of commingling analysis is that, although genetic effect size may be estimated (via the model parameter for genetic displacement, T), if a single major gene is implicated, effect sizes cannot be estimated if the results do not suggest a single major gene. 

The best-fitting model for the IOP data set was a mixture of three distributions, which would be compatible with the presence of a major gene that determines IOP. This confirms findings in previous work. ¹⁴ Although IOP was associated with CCT, adjustment of IOP for CCT did not provide any notable changes in the model in comparison to the model for unadjusted IOP. 

There was no evidence of violation of HWE throughout the analyses. The Hardy-Weinberg principle rests on certain assumptions: a diploid organism, sexual reproduction, nonoverlapping generations, bi-allelic genes, identical allele frequencies between males and females, nonassortative mating, large (theoretically, infinite) population size, negligible migration and mutation, and lack of selection pressure on the alleles under consideration. If any of these assumptions is seriously incorrect, departure from HWE may be observed. This does not seem to be the case in the Blue Mountains Eye Study cohort. Previous commingling analyses on IOP ¹⁴ and optic disc morphology ¹³ do not suggest departure from HWE. Furthermore, a recent study ⁵² of the complement factor H (CFH) single-nucleotide polymorphism rs1061170 (Y402H) in this cohort identified CC, CT, and TT genotypes in 13.6%, 46.7%, and 39.7% of the population, respectively. These proportions are in close HWE (χ² P = 0.94 for T allele frequency 0.63). 

The adherence to HWE in modeling studies of both IOP and optic disc morphology (traits which are unlikely to have an identical genetic basis) and in a CFH polymorphism (unrelated to glaucoma endophenotypes) suggests that population substructure is unlikely to be a confounding influence in this cohort. This probability is supported by the largely Caucasian ethnicity of the cohort and by a recent study of the mitochondrial DNA of the cohort ⁵³ which found that 94.2% of samples could be assigned to one of the major European haplogroups. 

A limitation of commingling analysis is that it uses models applied to trait distributions to infer the presence of possible genetic effects, rather than using real genotypic data: it is thus a precursor to formal genotyping studies. The findings of the present study closely resemble those of a previous study and are biologically plausible. ¹⁴ A quantitative trait locus (QTL), which accounts for 18.8% of the IOP trait variance, would be amenable to detection by QTL linkage methods. ⁵⁴ Evidence of commingling does suggest the possibility that a single locus has a major effect on the trait, and commingling analysis can provide guidance in the choice of initial parameter estimates for segregation analysis. Furthermore, these findings justify the effort and expense of a high-density genome-wide association study to identify the putative QTL associated with IOP in this cohort. This study is currently under way as part of the Wellcome Trust Case Control Consortium. Details are available at https://www.wtccc.org.uk/ccc2/projects/ccc2_gc.shtml. 

In conclusion, the present study found no evidence of a single major genetic determinant of CCT distributions under the mixed genetic model in the Blue Mountains Eye Study data set, whereas findings for IOP suggesting the a mixture of three distributions and the inference of a single major gene are robust to adjustment of IOP for CCT. 

Footnotes

 Supported by Special Trustees of Moorfields Eye Hospital.

Footnotes

 Disclosure: M.L. Rahman, None; C. Bunce, None; P.R. Healey, None; P. Mitchell, None; P.C. Sham, None; P. McGuffin, None; A.C. Viswanathan, None