Association of Novel Loci With Keratoconus Susceptibility in a Chinese Genome-Wide Association Study

Abstract

Purpose: Keratoconus (KC) is a progressive corneal disease that can lead to corneal blindness if not properly managed. The purpose of this study was to identify genetic associations with KC in China and to investigate whether these genetic variants are associated with corneal thickness and corneal curvature in KC cases.

Methods: A genome-wide association study was conducted on 853 patients with KC and 6248 controls. The KC cases were genotyped with the Illumina Infinium Human Asian Screening Array BeadChip, and the controls were genotyped with the Illumina Infinium Human Global Screening Array BeadChip. Genetic associations with KC, as well as correlations between the positive variants and corneal parameters including central corneal thickness (CCT) and mean keratometry (Km), were compared using PLINK version 1.90.

Results: Our present study identified four single-nucleotide polymorphisms (SNPs) within four risk loci (PTGER3: rs2300163, EYA1: rs1077435, ASS1: rs141365191, and CHTF8: rs3743680) associated with KC in Chinese patients that reached genome-wide significance. Among the identified SNPs with P < 1.00 × 10⁻⁴, seven SNPs (FOSL2-PLB1: rs12622211, RXRA-COL5A1: rs3118515, rs3132306, rs1536482, rs3118520, KAT6B: rs192187772, RAP2A-IPO5: rs41361245) were observed to be associated with CCT, and one SNP (USP13: rs6767552) was found to be associated with Km.

Conclusions: In the first genome-wide association study of KC with a relatively large study population in China, we identified four SNPs in four risk loci associated with the disease. The findings enriched the understanding of genetic susceptibility to KC and provided new insights into the genetic etiology of the disease.

Keratoconus (KC) is a complex corneal disorder characterized by progressive corneal thinning and protrusion.¹ Generally, the disease occurs in adolescence, with progression into the third or fourth decade of life.² The estimated prevalence of KC varies among different ethnicities and populations, with a higher prevalence in the Asian population.³ Both sexes are affected,⁴ but the percentage of males is higher than that in females in central China.⁵ Because the onset of KC is insidious and its progression is irreversible, early diagnosis and proper treatment are essential. Notably, understanding the etiology and pathogenesis of a disease is highly important for its diagnosis and treatment. Therefore the need to explore the etiology and pathogenesis of KC is urgent. 

Currently, the etiology and pathogenesis of KC remain unclear. Both twin studies and familial aggregation studies have indicated a genetic component of the disease.^6,7 A twin study by Tuft et al.⁸ identified that the severity of KC was more concordant in monozygotic twins than in dizygotic twins, suggesting a genetic contribution to the pathogenesis of the disease. Familial aggregation studies revealed a higher prevalence of KC in first-degree relatives of KC patients,^9,10 indicating a strong genetic predisposition of the disease. To date, multiple studies have been conducted on exploring the genetic mechanism of KC. Several candidate genes or variants associated with KC have subsequently been identified.^11–13 However, some genes or variants exhibited different associations in different populations, indicating genetic heterogeneity in KC. 

Genome-wide association study (GWAS), an effective method for exploring genotype–phenotype associations, has been widely applied to search for susceptibility loci in multiple disorders, including KC. A multicenter GWAS identified a novel locus in the patatin-like phospholipase domain–containing 2 (PNPLA2) gene in white patients with KC.¹⁴ The variants in RAB3 GTPase activating protein catalytic subunit 1 (RAB3GAP1) gene¹⁵ and hepatocyte growth factor (HGF) gene¹⁶ have been reported to be associated with KC in GWAS of the Caucasian population. In addition, several susceptibility loci associated with KC, such as variants in the stonin 2 (STON2) gene¹⁷ and collagen type V alpha 1 chain (COL5A1) gene,¹⁸ were identified through large GWASs of corneal thickness, indicating that these loci might result in KC by influencing corneal thickness. Considering that KC is a complex disease with strong genetic heterogeneity and that no GWAS has been performed in Chinese patients with KC, it is of great significance to conduct GWAS of KC in China to explore the genetic etiology of KC in Chinese patients. 

In the present study, we conducted the first GWAS of KC with a relatively large study population (Supplementary Table S1) in China, aiming to identify genetic susceptibility loci associated with KC in Chinese patients. Besides, the correlations between positive variants and corneal parameters including central corneal thickness (CCT) and mean keratometry (Km) in KC patients were also explored. The study might provide novel insights into the genetic etiology of KC in Chinese patients. 

Material and Methods

Study Population

A total of 958 patients with KC from the Chinese keratoconus (CKC) cohort study and 7186 controls were included in the current study. All subjects were from China and self-reported Chinese. The KC cases were recruited from Henan Eye Hospital and diagnosed according to the following criteria: corneal topography revealing an asymmetric bowtie pattern with or without skewed axes or keratoconus signs detected by slit lamp examination, such as localized stromal thinning, conical protrusion, Vogt's striae, Fleischer's ring or anterior stromal scar.² Patients with KC who had syndromic diseases (e.g., Down syndrome, Ehlers-Danlos syndrome, Leber congenital amaurosis), as well as patients concomitant with corneal dystrophy were excluded. The healthy controls (n = 7186) were recruited from the Chinese Han population through collaboration with multiple hospitals in China.¹⁹ The controls were clinically assessed as being without systemic disorders, autoimmune disorders, or family history of autoimmune disorders (including first-, second-, and third-degree relatives). The study was approved by the Ethics Committee of Henan Eye Hospital and the First Affiliated Hospital of Anhui Medical University in accordance with the guidelines of the Declaration of Helsinki. In addition, written informed consent was obtained from each participant or his/her legal guardians before participation. 

Genotyping and Quality Control

Genomic DNA was extracted from peripheral blood samples according to the manufacturer's recommendations. As shown in Figure 1, the KC cases were genotyped with the Illumina Infinium Human Asian Screening Array (ASA) BeadChip, which is introduced by Illumina in 2017 for East and Southeast Asian populations, and contains a total of 659184 single-nucleotide polymorphisms (SNPs). The controls were genotyped with the Illumina Infinium Human Global Screening Array (GSA) BeadChip, which is launched by Illumina in 2016. The GSA BeadChip is a global population genotyping chip designed according to the Phase III of the 1000 Genomes Project, featuring 700078 SNPs. In our prior study, the controls were genotyped using the GSA BeadChip.¹⁹ For the current study, in an effort to reduce expenses, only the patients underwent genotyping with the ASA BeadChip. All genotype analyses were performed using Genome Studio v2.0 and PLINK version 1.90 software. To integrate these distinct chip datasets, we initially aligned both datasets to the forward strand as per their corresponding manifest files (GSAMD-24v1-0_20011747_A1.csv and ASA-24v1-0_E1.csv, GRCh37). Subsequently, after excluding 7387 SNPs lacking chromosome identifiers in the GSA chip, both datasets were cross-referenced with the 1000 Genomes Phase III dataset to eliminate mismatched SNPs (33 variants from ASA and 790 from GSA were excluded). After this, the patient and control datasets were merged, and two variants with three or more alleles were eliminated before proceeding to the quality control assessment. SNPs on sex chromosomes, major histocompatibility complex regions and mitochondria were excluded. Then, samples with SNP call rates less than 95% and individuals closely related based on estimated identity-by-descent (PI_HAT > 0.25) were excluded from further analysis. Samples with abnormal heterozygosity (> mean ± 3SD) or genotyping call rate <95% were also excluded. 

Figure 1.

 

The study design of identifying associated loci with keratoconus in a Chinese population.

Principal Component Analysis and GWAS

After the quality control procedures, principal component analysis was performed to investigate the population structure. Twenty principal components were extracted from the variance-standardized relationship matrix. In the PC1 versus PC2 coordinate, the average distance from the centroid was 0.014, with a standard deviation of 0.007. We used the average distance plus three times SD as the threshold value (0.035) to identify outliers, and 1043 samples were removed. Finally, 853 KC cases and 6248 controls were left for the GWAS, and the first to three principal components were used as covariates in the subsequent association analysis. SNPs were excluded if they had a call rate <95% in cases or controls, minor allele frequency <1% in the population, or deviation from Hardy-Weinberg equilibrium in the controls (P_hwe < 1.00 × 10⁻⁴). Conditional analysis was carried out to determine independent association signals. For the identified variants reaching genome-wide significance, regional association plots were generated using LocusZoom (http://locuszoom.sph.umich.edu//). 

Correlation Analysis

Linear regression analyses were conducted in KC cases to explore the correlations between identified variants and corneal parameters including CCT and Km. Both parameters were recorded from Pentacam HR (Oculus, Lynnwood, WA, USA), which is a widely used Scheimpflug-based tomographic device. The parameters in the worse eye were used to conduct the analysis. 

Results

Clinical Characteristics of the Study Population

As shown in Table 1, a total of 853 patients with KC and 6248 controls were included in the GWAS. The average age of KC patients was 21.91 ± 6.74 years. There were 647 (75.8%) male patients and 206 (24.2%) female patients included in this study. The mean CCT of the patients was 436.46 ± 58.47 µm, and the mean Km was 53.29 ± 8.48D. In the control group, there were 3042 (48.7%) males and 3206 (51.3%) females. The average age of the individuals in the control group was 33.52 ± 15.42 years. 

Table 1.

 

Clinical Characteristics of the Study Population

 

Genome-Wide Association Analysis

As shown in Supplementary Figure S1, no abnormal deviation was found among the samples, and the cases and controls were basically matched (PC1-2, PC1-3, PC1-4, PC1-5). The quantile–quantile plot is shown in Supplementary Figure S2. The genomic inflation factor (λ_GC) was 1.087, indicating no impact of population stratification. A Manhattan plot for the association analysis with sex and the first three principal components as covariates is presented in Figure 2. As shown in Supplementary Table S2, 53 SNPs at 44 loci had P values < 1.00 × 10⁻⁴. Among these loci, five were previously reported to be associated with KC, namely, 3q26.31, 9q34.3, 13q32.2, 15q23, and 17p11.2. Among the identified SNPs with P < 1.00 × 10⁻⁴, a total of four SNPs were previously reported to be associated with KC including rs4894535 (fibronectin type III domain containing 3B, FNDC3B, P = 1.26 × 10⁻⁵), rs1536482 (retinoid X receptor alpha - collagen type V alpha 1 chain, RXRA-COL5A1, P = 1.60 × 10⁻⁷), rs12913547 (SMAD family member 3, SMAD3, P = 4.67 × 10⁻⁶), and rs2228100 (aldehyde dehydrogenase 3 family member A1, ALDH3A1, P = 1.27 × 10⁻⁶). In addition, the results also revealed other SNPs located on previously reported loci with P value < 1.00 × 10⁻⁴. These findings suggested that the genotype data in the current study were highly reliable. 

Figure 2.

 

The Manhattan plot for genome-wide associations with keratoconus. (Genome-wide associations results with sex and the first three principal components as covariates were plotted for each chromosome. The x axis represents the genome in physical order. The y axis shows −log10 P values for all SNPs. The red line represents the genome-wide significance threshold of a P value of 5.00 × 10⁻⁸).

Totally, there were four SNPs in four loci that reached genome-wide significance (Table 2). The significant SNPs at each locus were rs2300163 (1p31.1, prostaglandin E receptor 3, PTGER3, P = 7.70 × 10⁻¹⁰), rs1077435 (8q21.11, EYA transcriptional coactivator and phosphatase 1, EYA1, P = 3.03 × 10⁻⁸), rs141365191 (9q34.2, argininosuccinate synthase 1, ASS1, P = 3.05 × 10⁻⁸), and rs3743680 (16q22.1, chromosome transmission fidelity factor 8, CHTF8, P = 7.81 × 10⁻¹⁰). The regional plots corresponding to the SNPs that reached genome-wide significance were shown in Figure 3. 

Table 2.

 

Association Results for SNPs With P Value < 5 × 10⁻⁸

 

Figure 3.

 

Regional plots for keratoconus-related variants with genome-wide significance. (The P values of SNPs (shown as –log10 P values in y-axis) were plotted against their map positions (x-axis) with Locus Zoom (http://locuszoom.sph.umich.edu//). The selected SNP was in purple, and its linkage disequilibrium values (r²) with nearby SNPs are indicated by different colors. (A) rs2300163 ± 1000 kb; (B) rs1077435 ± 1000 kb; (C) rs141365191 ± 1000kb; (D) rs3743680 ± 1000 kb).

Correlation of the Identified Variants With Corneal Parameters

To explore whether the identified variants were correlated to corneal clinical characteristics such as corneal thickness or corneal curvature in KC cases, linear regression analyses were conducted on the identified variants and corresponding clinical parameters. As shown in Table 3, no SNPs with genome-wide significance were associated with CCT or Km (all P > 0.05). For identified SNPs with P < 1.00 × 10⁻⁴, seven SNPs (FOS like 2, AP-1 transcription factor subunit–phospholipase B1, FOSL2-PLB1: rs12622211, RXRA-COL5A1: rs3118515, rs3132306, rs1536482, rs3118520, lysine acetyltransferase 6B, KAT6B: rs192187772, RAP2A, member of RAS oncogene family–importin 5, RAP2A-IPO5: rs41361245) exhibited a nominal association with CCT (all P < 0.05), indicating that these SNPs might be associated with KC by influencing corneal thickness. Additionally, rs6767552 located in ubiquitin-specific peptidase 13 (USP13) gene was nominally associated with Km (P < 0.05), indicating that the SNP might be associated with KC by influencing corneal curvature. For other SNPs, there were no significant associations with CCT or Km in our study population. 

Table 3.

 

Correlation Between Identified SNPs and Corneal Parameters in Keratoconus

 

Discussion

KC is a complex heterogeneous disease with genetic factors involved in its pathogenesis.²⁰ In this study, we performed the first GWAS of KC in China, and identified four SNPs within four risk loci (PTGER3: rs2300163, EYA1: rs1077435, ASS1: rs141365191, CHTF8: rs3743680) associated with KC in Chinese patients that reached genome-wide significance. Among the identified loci with P < 1.00 × 10⁻⁴, there were 39 new risk loci and five previously reported loci. In addition, among the SNPs with P < 1.00 × 10⁻⁴, seven SNPs (FOSL2-PLB1: rs12622211, RXRA-COL5A1: rs3118515, rs3132306, rs1536482, rs3118520, KAT6B: rs192187772, RAP2A-IPO5: rs41361245) were observed to be associated with CCT in KC cases, and the rs6767552 located in USP13 gene was found to be associated with Km, indicating that these SNPs might result in KC by influencing corneal thickness and corneal curvature respectively. 

As reported in previous studies, KC is considered as a complex disorder with a strong genetic predisposition.^3,20 Currently, many researchers are committed to exploring the genetic etiology of KC and have identified a large number of variants associated with the disease via various study strategies.²¹ As an effective method to identify genetic risk loci, GWAS has also been applied in KC. In the first GWAS of KC, Li et al.¹⁵ included 222 Caucasian KC cases and 3324 Caucasian controls in the discovery stage and finally identified a novel potential KC locus at 2q21.3, containing a candidate gene RAB3GAP1. McComish et al.¹⁴ conducted a relatively large GWAS with 522 KC patients and 665 control participants in the discovery stage and identified a genome-wide significant locus for KC in PNPLA2. The largest GWAS of KC conducted by Hardcastle et al.²² identified 36 genomic loci associated with KC, and implicated both dysregulation of corneal collagen matrix integrity and cell differentiation pathways as primary disease-causing mechanisms. However, because of the strong genetic heterogeneity of KC and the limited studies in China, it is also of great importance to conduct GWAS in Chinese patients with KC. In the present study, we conducted the first GWAS of KC in a relatively large population in China with Illumina BeadChip screening arrays, containing 853 Chinese KC patients and 6248 controls. Finally, four SNPs located in four risk loci were found to be associated with KC at genome-wide significance. 

Among the variants with genome-wide significance in the present study, the most significant variant at 1p31.1 was rs2300163 in PTGER3 gene. The T allele at rs2300163 was associated with an increased risk of KC in our study. The prostaglandin E receptor 3, encoded by PTGER3 gene, is one of four receptors for prostaglandin E2 and has many biological functions including inflammatory response.^23,24 Although KC has been classically defined as a non-inflammatory condition, several recent studies, as well as our previous study, have indicated an association between KC and inflammation.^3,25–27 Therefore we speculated that PTGER3 might be related to KC via its involvement in inflammation. 

At 8q21.11, rs1077435 in the EYA1 gene was shown to be associated with KC at genome-wide significance, and the A allele of the variant was revealed to be associated with an increased risk of KC. The EYA1 gene encodes a member of the eyes absent family of proteins that might play a role in eye anomalies. A previous study conducted by Azuma et al.²⁸ identified novel mutations of EYA1 gene in patients who had congenital cataracts and ocular anterior segment anomalies, indicating a potential role of EYA1 in eye morphogenesis. As is known to us all, the VSX1 gene is a genetic susceptibility gene for KC, and it might be related to KC via its involvement in corneal development.²⁹ Similarly, we inferred that the EYA1 gene might influence the occurrence of KC by affecting eye morphogenesis. 

The genome-wide significant variant identified at 9q34.2 was rs141365191 in ASS1 gene. The A allele at rs141365191 was observed to be associated with an increased risk of KC. Argininosuccinate synthetase 1, encoded by the ASS1 gene, is a rate-limiting enzyme that catalyzes the biosynthesis of arginine.³⁰ Interestingly, the suppressed arginase activity was detected in KC-keratocytes.³¹ Moreover, another study found that arginine supplementation could promote extracellular matrix and metabolic changes in KC, indicating a potential role and therapeutic applications of arginine in the disease.³² As a rate-limiting enzyme in the biosynthesis of arginine, the ASS1 might be involved in the pathogenesis of KC via mediation the biosynthesis of arginine. 

Another variant with genome-wide significance was rs3743680 in CHTF8 gene. The T allele at rs3743680 exhibited an increased risk of KC in our study population. The CHTF8 gene encodes a short protein that forms part of the Ctf18 replication factor C complex and has been demonstrated to be involved in DNA replication and repair.³³ Although the pathogenesis of KC has not been precisely defined, oxidative stress is considered to play an important role in the development of the disease.³⁴ The exposure of the cornea to endogenous and exogenous reactive oxygen species can result in various types of molecular damage, including DNA damage.³⁵ Therefore the accurate DNA repair mechanism is essential for maintaining genome integrity, as well as the normal structure and function of the cornea. In addition, some genes related to DNA repair have also been identified to be associated with KC.^36,37 Consequently, we speculated that the CHTF8 might be involved in KC by influencing the process of DNA repair. 

As is known to us all, the progressive corneal thinning and protrusion are typical corneal clinical characteristics of KC. To explore whether the identified SNPs were related to corneal clinical characteristics in KC cases, correlations between the identified SNPs and corneal clinical characteristics such as corneal thickness or corneal curvature were investigated. Consequently, the variants with genome-wide significance showed no correlation to CCT or Km. However, seven SNPs with P < 1.00 × 10⁻⁴ (FOSL2-PLB1: rs12622211, RXRA-COL5A1: rs3118515, rs3132306, rs1536482, rs3118520, KAT6B: rs192187772, RAP2A-IPO5: rs41361245) were identified to be correlated to CCT in our study, indicating a potential connection between corneal thinning and KC. The findings were consistent with several previous studies that demonstrated associations of variants in RXRA-COL5A1 (rs3118515, rs3132306, rs1536482) with CCT.^18,38,39 In addition, rs6767552 in USP13 gene was found to be associated with Km, indicating that the variants might result in KC by influencing corneal curvature. 

Several limitations should be noted in the current study. First, the study included only 853 patients with KC and 6248 controls to identify KC-related variants in the analysis. Considering the limited prevalence of KC and difficulty in collecting blood samples, the results were not verified in a replication population. Future studies should be conducted to verify the findings. Second, the molecular mechanism of identified genes or variants was not explored in the present study. Further studies are warranted to clarify the pathogenesis of KC caused by those genes or variants. 

Conclusions

In conclusion, the first GWAS of KC in China identified four SNPs in four risk loci associated with the disease. Moreover, seven SNPs with P < 1.00 × 10⁻⁴ (FOSL2-PLB1: rs12622211, RXRA-COL5A1: rs3118515, rs3132306, rs1536482, rs3118520, KAT6B: rs192187772, RAP2A-IPO5: rs41361245) were observed to be associated with CCT in KC cases, and one SNP (USP13: rs6767552) was found to be associated with Km, suggesting a potential role of these variants in KC by influencing corneal thickness and corneal curvature. The findings provided convincing evidences for the genetic etiology of KC, and further studies would be warranted to investigate the molecular mechanisms underlying the identified variants in the pathogenesis of the disease. 

Acknowledgments

Supported by the Henan Young Health Science and Technology Innovation Outstanding Program (No. YXKC2020023), Special Program for Basic Research of Henan Eye Hospital (No. 20JCZD003), Youth Special Program for Basic Research of Henan Eye Hospital (No. 21JCQN006, 21JCQN008). 

Disclosure: L. Xu, None; X. Zheng, None; S. Yin, None; K. Yang, None; Q. Fan, None; Y. Gu, None; Y. Yuan, None; C. Yin, None; Y. Zang, None; C. Pang, None; L. Sun, None; S. Ren, None