Genetic Correlations Among Corneal Biophysical Parameters and Anthropometric Traits

Abstract

Purpose: The genetic architecture of corneal dysfunction remains poorly understood. Epidemiological and clinical evidence suggests a relationship between corneal structural features and anthropometric measures. We used global and local genetic similarity analysis to identify genomic features that may underlie structural corneal dysfunction.

Methods: We assembled genome-wide association study summary statistics for corneal features (central corneal thickness, corneal hysteresis [CH], corneal resistance factor [CRF], and the 3 mm index of keratometry) and anthropometric traits (body mass index, weight, and height) in Europeans. We calculated global genetic correlations (r_g) between traits using linkage disequilibrium (LD) score regression and local genetic covariance using ρ-HESS, which partitions the genome and performs regression with LD regions. Finally, we identified genes located within regions of significant genetic covariance and analyzed patterns of tissue expression and pathway enrichment.

Results: Global LD score regression revealed significant negative correlations between height and both CH (r_g = −0.12; P = 2.0 × 10⁻⁷) and CRF (r_g = −0.11; P = 6.9 × 10⁻⁷). Local analysis revealed 68 genomic regions exhibiting significant local genetic covariance between CRF and height, containing 2874 unique genes. Pathway analysis of genes in regions with significant local r_g revealed enrichment among signaling pathways with known keratoconus associations, including cadherin and Wnt signaling, as well as enrichment of genes modulated by copper and zinc ions.

Conclusions: Corneal biophysical parameters and height share a common genomic architecture, which may facilitate identification of disease-associated genes and therapies for corneal ectasias.

Translational Relevance: Local genetic covariance analysis enables the identification of associated genes and therapeutic targets for corneal ectatic disease.

Introduction

The causes and consequences of corneal biomechanical dysfunction are complex. Although the biophysical features underlying corneal shape contribute directly to ocular diseases such as keratoconus, they also demonstrate associations with systemic disorders, particularly those of the connective tissue.¹ For instance, Marfan syndrome, a systemic disorder of impaired microfibril formation classically resulting in elongated bone structures, is associated with reduced corneal curvature.^2–4 Ehlers-Danlos syndrome and Williams syndrome, both of which result from defective collagen and elastin synthesis and are also associated with taller stature, contribute to thinner and steeper corneal contours.^5–7 In contrast, Turner syndrome and Down syndrome, both canonically associated with lower body heights, are associated with specific corneal phenotypes such as changes in central corneal thickness (CCT), compared to healthy controls.^8–10 More broadly, recent population-scale evidence suggests an association between reduced corneal refractive power, indicative of corneal flatness, and increased body height.¹¹ Such conditions collectively demonstrate a varied genetic landscape encompassing both autosomal and sex-chromosomal defects, suggesting the existence of an underlying genomic relationship between anthropometric features and corneal biophysical traits.¹² 

Currently, the specific genomic drivers of corneal biomechanical dysfunction remain poorly understood, despite evidence that genetic polymorphisms explain a large proportion of such phenotypic variation.^13–15 Keratoconus itself is believed to share a genomic underpinning with systemic connective tissue disorders. One genome-wide association study (GWAS) identified significant keratoconus associations for 12 novel loci, implicating both collagen matrix production and cell differentiation pathways.¹⁵ Another GWAS in 2021 identified 32 genetic variants associated with corneal curvature, most notably Wnt and fibroblast growth factor signaling.¹⁶ In 2022, a GWAS relying on biophysical endophenotypes as surrogates of keratoconus identified 150 novel loci and similarly implicated pathways involved in connective tissue synthesis.¹⁷ Furthermore, corneal curvature, anterior chamber depth, and axial length are themselves highly heritable,¹⁸ with corneal curvature and axial length demonstrating significant genetic correlation.¹⁹ 

Although association studies have identified some genomic variants associated with keratoconus, traditional GWAS are limited in their ability to resolve common, low-effect variants that likely underlie a large proportion of corneal ectasia heritability.^13,20,21 Such variants may be studied indirectly through genetic correlation analysis, which compares correlations between measured variant effects for distinct phenotypes.²² Such analysis can be performed globally, comparing correlation among all measured variants to assess the genomic similarity of given phenotypes, or locally, measuring region-specific correlations to implicate specific patterns of variation shared in both conditions. To clarify the shared genomic underpinnings of corneal ectasia and systemic connective-tissue development, we performed a genetic correlation analysis of corneal biophysical features and anthropometric features. 

Methods

Ethical Statement

Because summary-level genetic data used in this study are publicly available and deidentified, institutional review board review for this project was not required. Bulk corneal gene expression data were made available from a previous study, which was approved by the institutional review board at Augusta University.²³ 

GWAS Datasets

We obtained summary statistics for publicly available GWAS of keratoconus-related and anthropometric traits, including CCT, corneal hysteresis (CH), corneal resistance factor (CRF), 3 mm index of keratometry result, body mass index (BMI), weight, and height. All GWAS studies considered in the analysis are summarized in Table 1, and the summary data, as well as details on variant calling and preprocessing, are available at the respective citations, as well as via the Neal Lab's data repository for these GWAS studies, available at http://www.nealelab.is/uk-biobank/. 

Table 1.

 

Summary of Traits Analyzed in Genome-Wide Association Studies

 

Both corneal biophysical parameters and anthropometric traits showed significant heritability explained in their respective GWAS data, with explained heritability proportions ranging from 0.14 (BMI) to 0.46 (height). The 3 mm index of keratometry result was an outlier, with an explained heritability of 0.0049. 

Global Genetic Correlations

Genetic correlations occur in the presence of shared single nucleotide polymorphism (SNP) heritability. SNP heritability is the proportion of variability in a trait that is due to genetic contributions, and it can be estimated by classical twin studies or from regressing test statistics from SNPs in a GWAS across linkage disequilibrium (LD) blocks.²⁴ As implemented in the LDSC package (https://github.com/bulik/ldsc),^24,25 we calculated GWAS-explained heritability using all SNPs included in the GWAS. We calculated genome-wide genetic correlation as the covariance of two traits divided by the product of the traits’ heritabilities. Because LD blocks vary with ancestry, we confined our analysis to GWAS with participants of European ancestry. We used genome-wide LD score regression as implemented in the LDSC package. The Bonferroni correction (α = 0.05, n = 21) was applied to account for multiple testing. Power calculations for all pairwise analyses are provided in Supplementary Table S1. 

Local Genetic Correlations

We used ρ-HESS (https://huwenboshi.github.io/hess/)²⁶ to estimate the shared heritability of two traits across partitions of the genome. This method applies a fixed effect model to assess the covariance to the combined effects of SNPs on the two traits. In local genetic covariance analysis, the genome is partitioned into approximately independent LD regions with an average width of 1.6 Mb. This results in 1703 regions across all chromosomes, excluding the sex chromosomes. These regions were the same as the ones implemented by default in ρ-HESS as we had no rationale to prioritize different genome partitions. We conducted local genetic covariance analyses for CRF and height given their high global genetic correlation. We also stratified CRF by sex and performed global genetic correlation analysis with the other corneal and anthropometric traits. 

We prioritized individual genes for downstream analysis by the significance of local covariance between CRF and height. CRF was selected as the endophenotype of choice for this analysis because it had the largest set of significant SNPs of the biophysical traits investigated.¹⁷ Specifically, we identified all unique genes in the GRCh37 reference assembly located within regions exhibiting significant local r_g after Bonferroni correction (α = 0.05, n = 1703). Genome mapping was performed using PyEnsembl (release 97). A full list of the genes is provided as supplementary data. 

Tissue Expression Analysis

We performed single-cell expression analysis of genes in regions of high local covariance using an atlas of gene expression in 19 cell types of the anterior chamber of healthy human eyes.²⁷ These cell types included corneal epithelium, B-cells, macrophages, natural killer cells, mast cells, Schwalbe cells, fibroblasts, collector channel cells, vascular endothelium, Schlemm canal cells, juxtacanalicular tissue cells, beam cells (types A and B), ciliary muscle cells, pericytes, melanocytes, myelinating and non-myelinating Schwann cells, and neurons. In addition to single-cell analysis of healthy tissue, we also compared the expression of genes in regions of high local covariance using bulk corneal tissue from 10 keratoconus patients and 8 healthy controls.²³ Comparisons were performed using Mann-Whitney U tests, with a two-sided significance level of 0.05. 

Functional Gene Annotation

We compared the overlap of gene candidates in regions of significant local covariance with genes associated with either keratoconus or body height in the existing literature. Literature-based gene-phenotype associations were obtained from the DisGeNET ontology, which integrates clinical associations, experimental evidence, and computational predictions to produce a confidence score for gene-phenotype associations.²⁸ 

Functional pathway enrichment among gene candidates was measured using gene set enrichment analysis.²⁹ Analyses were performed in WebGestalt, using PANTHER pathway gene sets to evaluate pathway involvement and DrugBank drug-target gene sets for drug associations, in both cases using the default running parameters.^30–32 

Results

Global Genetic Correlations Between Traits

Both corneal biophysical parameters and anthropometric traits showed highly significant correlations among themselves, as expected. CH was highly positively correlated with CCT (r_g = 0.64; standard error (SE) = 0.06; P = 3.0e-24), as was CH with CRF (r_g = 0.88; SE = 0.01; P < 1e-100) and CCT with CRF (r_g = 0.71; SE = 0.06; P = 8.7e-36). Furthermore, there were extreme positive correlations between weight and both height (r_g = 0.43; SE = 0.02; P < 1e-100) and BMI (r_g = 0.81; SE = 0.02; P < 1e-100), as well as a weaker negative correlation between height and BMI (r_g = −0.09; SE = 0.02; P = 6.3e-6), consistent with previous evidence.²⁵ 

Genetic correlation analysis revealed significant inverse global associations between corneal biophysical features and anthropometric traits (Table 2). Specifically, we observed significant negative genetic correlations between CRF and height (r_g = −0.11; SE = 0.02; P = 6.9e-7), as well as between CH and height (r_g = −0.12; SE = 0.02; P = 2.0e-7). No other pairwise correlations reached significance. The sex-stratified CRF and height genetic correlations (Supplementary Table S2) were similar (females only: r_g = −0.11; SE = 0.02; P = 6.0e-6; males only: r_g = −0.11; SE = 0.03; P = 6.0e-4). 

Table 2.

 

Genetic Correlations Among Keratoconus-Related and Anthropometric Traits^*

 

Local Genetic Correlations for CRF and Height

We next searched for regions of significant local genetic covariance between height and corneal biophysical parameters, which represent drivers of the global genetic correlation between the two phenotypes. Given the high correlation between CRF and CH (r_g = 0.88) and the higher heritability for CRF, we restricted our analysis to the local drivers of the CRF-height genetic association. After adjusting for multiple hypothesis testing, 68 genomic regions showed significant local genetic covariance (Table 3; Supplementary Table S3). The regions had an average length of 1.7 Mb (standard deviation [SD] = 0.75 Mb) and contained an average of 2930 (SD = 1160) unique SNPs and 43 (SD = 38) unique genes per region. Although the global correlation was negative between CRF and height, we observed mixed directionality of covariance at the local level, with 51 significant negative and 17 significant positive covariances. 

Table 3.

 

Ten Regions of Strongest Local Genetic Covariance Between Corneal Resistance Factor and Height

 

Tissue Expression of Correlated Genes

We identified 2874 unique genes in regions of significant local covariance. To determine whether these genes could represent shared genetic drivers between corneal dysfunction and height, we evaluated tissue expression patterns among genes in highly correlated regions (Fig. 1A) using an ocular anterior chamber single-cell gene expression dataset. Among the 19 cell types investigated, genes in regions of high covariance showed the greatest mean expression in corneal endothelial (Schwalbe) cells and the lowest in B lymphocytes. Mean expression in Schwalbe corneal endothelial cells of genes in regions with high covariance was significantly higher than in all other cell types (P ≤ 0.035 in all cases). Furthermore, expression in corneal epithelial cells was also high, being significantly greater (P ≤ 0.001 in all cases) than that in all but four other cell types (Schwalbe, Type B beam, juxtacanalicular tissue, and neuronal cells). We finally compared expression levels of genes in regions of high local covariance in corneal tissue between patients with keratoconus and healthy controls, observing that genes in regions of high genetic covariance showed 13% higher median expression levels in corneas with keratoconus than in healthy corneas (P = 0.003; Figs. 1B, 1C). 

Figure 1.

 

Expression in human ocular tissue of genes in regions of significant genetic covariance. (A) Distribution of single-cell expression levels in ocular cell types for genes in regions of significant covariance reveals high expression in cornea-related cell types, including Schwalbe cells and corneal epithelial cells. BeamCella, type A beam cells; BeamCellb, type B beam cells; CollectorChnIAqVein, collector channel cells; CribiformJCT, juxtacanalicular tissue cells; NKT, natural killer cells; ScEndo, Schlemm canal cells; SchwannCell-my, myelinating Schwann cells; SchwannCell-nmy, nonmyelinating Schwann cells. (B, C) Overall expression levels of genes in regions of high covariance were significantly greater in bulk corneal tissue from keratoconus patients than in corneal tissue from controls.

Functional Pathway Enrichment Among Correlated Genes

We also evaluated whether genes in regions of high local genetic covariance included genes with known associations with both keratoconus and body height in a literature-based ontology of gene-phenotype associations.²⁸ Among the 2874 genes contained in locally correlated regions, 18 had known associations with keratoconus, whereas 153 had known associations with body height, and three genes, NOX4, FNDC3B, and ADAMTS17, had known associations with both. 

We next measured functional pathway enrichment among the list of regions of high covariance using gene set enrichment analysis. Regions with high covariance were significantly enriched for signaling pathways involved in cell-cell adhesion and connective tissue maturation, including cadherin signaling, Wnt signaling, and Rho-GTPase–mediated cytoskeletal regulation (Fig. 2A). We also analyzed enrichment among drug-target gene sets for existing pharmacologic agents, with the most enriched compounds including the estrogen receptor modulator tamoxifen, as well as several metal ions involved in extracellular matrix composition, including copper and zinc (Fig. 2B). 

Figure 2.

 

Enrichment of functional gene sets in regions of significant genetic covariance. (A) Gene set enrichment analysis of PANTHER pathways among genes in regions of high covariance reveals enrichment of pathways related to cell-cell adhesion and connective tissue development. (B) Analysis of drug target gene sets shows enrichment of genes modulated by zinc, copper, and calcium ions in correlated regions, as well as enrichment of targets for tamoxifen.

Discussion

Structural corneal dysfunction is associated with significant visual morbidity in corneal ectatic diseases like keratoconus. While corneal biophysical properties are highly heritable, their genomic basis remains poorly understood.³³ Our findings demonstrate that structural attributes of the human cornea share a common genomic architecture with general anthropometric traits. Both CH and CRF showed significant global genetic correlations with height. We further identified 68 genomic regions with significant levels of local genetic covariance between CRF and height. These regions contained 2874 unique genes, which among them demonstrated enrichment of several functional pathways involved in corneal development, including Wnt/cadherin signaling and metal ion homeostasis. Furthermore, expression of these genes was high in corneal endothelial and epithelial cell types, which have been implicated in a causative pathway for keratoconus.^34,35 

The observed mild but statistically significant global association of height with CH and CRF suggests that corneal structural development relies on genes involved in connective tissue maturation throughout the body. This result provides a quantitative genomic context for the phenotypic associations observed with several corneal conditions, including keratoconus.^1,36–38 In addition to keratoconus itself, several keratoconus-related endophenotypes are associated with overall body habitus. For instance, one recent study identified a correlation between reduced corneal refractive power and greater body height in a European population.¹¹ Other studies have demonstrated associations between keratoconus incidence and both higher BMI and lower body height.^39,40 In addition to providing evidence of a common genomic architecture for structural corneal traits and body habitus, local genetic covariance analysis also allows for the identification of genomic regions driving this association. While previous GWAS have identified dozens of regions associated with keratoconus endophenotypes, including genes involved in collagen synthesis,¹⁷ cell differentiation,^15,41 and extracellular matrix formation,⁴² such studies collectively explain less than 20% of known heritability for these traits. This discrepancy is thought to derive from the abundant role of common, low-effect variants in corneal biomechanical development, an ideal setting for genetic correlation analysis.¹⁵ 

The genes contained within the 68 regions of high covariance facilitate a variety of processes with experimental or associative links to corneal ectasia. Cadherin signaling, which is broadly involved in connective tissue maturation, represented the single most enriched functional pathway among genes in correlated regions, providing additional evidence for its importance in corneal function specifically. For instance, several cadherin-family proteins, such as cadherin 11 and desmoglein 1, are known biomarkers for keratoconus, and N-cadherin is critical for corneal epithelial cell maturation.^43–45 We also identified significant enrichment of Wnt signaling genes, which likewise mediate collagen dysfunction in the corneal epithelium and axial skeletal connective tissue.^46,47 Beta-catenin, a central Wnt signaling hub, has been implicated as a potentially pathogenic mechanotransducer in keratoconus, and Wnt-protein sequence variants are associated with keratoconus risk.^48,49 Finally, we identified three genes, NOX4, FNDC3B, and ADAMTS17, present in regions of significant local covariance that also have known associations with both keratoconus and body height. In particular, NOX4 mediates oxidative stress response, which is dysregulated in corneal ectasias, whereas FNDC3B is broadly involved in cellular differentiation.^50–52 

We also observed enrichment of gene sets interacting with metal cations such as zinc, copper, and calcium and the selective estrogen receptor modulator tamoxifen. Although the therapeutic potential of such compounds in corneal ectasia remains to be established, these results provide additional support for the involvement of metal ion homeostasis and sex hormone signaling in keratoconus.^53–57 In particular, metal ion homeostasis, including of zinc, copper, and calcium specifically, is critical for extracellular matrix development and other supportive structures, providing a potential explanation for its shared role in both height and keratoconus endophenotypes.^58–60 Such associations may be therapeutically relevant, because a recent Phase 1/2a clinical trial of copper eyedrops showed a significant reduction in corneal steepness in keratoconus patients (iVeena Delivery Systems. Safety and preliminary efficacy of IVMED-80 eye drops in keratoconus patients. ClinicalTrials.gov identifier NCT05241145.) 

Our analysis has several limitations. The study population is predominantly of European ancestry, which may limit the generalizability of the results. We also note that although LD regression analysis identifies correlative relationships between genomic landscapes, it does not assess causal relationships, is limited in its ability to infer significance at the individual gene level, and may be susceptible to bias because of phenotype-specific mating patterns.⁶¹ Although our downstream pathway and expression analyses, as well as a lack of evidence for assortative mating in corneal ectasia, support the biological relevance of our genetic correlation findings, our analysis relies on the standard, imperfect assumptions of genetic correlation analysis. Furthermore, several groups have reported paradoxically inverted directionality findings from global genetic correlation analysis in comparison to corresponding phenotypic associations.^62,63 Although the theoretical basis for this phenomenon remains poorly understood, it nonetheless limits our ability to make assumptions about the direction of a phenotypic association based on genetic correlations. Finally, we acknowledge that ρ-HESS is inherently an estimation tool. Other approaches to assess local genetic correlation, such as LOGODetect (which implements a scan statistic search⁶⁴), SUPERGNOVA (which uses a random effects model⁶²), and LAVA (which uses partial correlation and multiple regression⁶⁵), could provide additional insights into local correlation patterns. Ultimately, given the associative nature of genetic correlation analysis, knockout and overexpression studies will be necessary to make conclusions about causality. 

Corneal structural abnormalities underly ectatic disease and yet remain poorly understood at a genetic level. Our results suggest that corneal development shares a common genomic architecture with gross body habitus, informing the design of both clinical biomarkers and therapeutic interventions. Further investigations, including in vivo interrogation of candidate genes, will provide additional insight into this relationship. 

Acknowledgments

Supported by an unrestricted Challenge Grant from Research to Prevent Blindness; NEI R01EY015473, EY032599, NEI R01EY023242, R01EY032960, R21EY033961, R21EY028671, and P30EY031631. 

Disclosure: H.C. Cousins, None; C.C. Cousins, None; G. Valluru, None; R.B. Altman, None; Y. Liu, None; L.R. Pasquale, Twenty Twenty (C), Character Bio (C); S. Ahmad, None