Corneal curvature (CC) measures the steepness of the cornea and is an important parameter for clinically diseases such as astigmatism and myopia. Despite the high heritability of CC, only two associated genes have been discovered to date. We performed a three-stage genome-wide association study meta-analysis in 12 660 Asian individuals. Our Stage 1 was done in multiethnic cohorts comprising 7440 individuals, followed by a Stage 2 replication in 2473 Chinese and Stage 3 in 2747 Japanese. The SNP array genotype data were imputed up to the 1000 Genomes Project Phase 1 cosmopolitan panel. The SNP association with the radii of CC was investigated in the linear regression model with the adjustment of age, gender and principal components. In addition to the known genes, MTOR (also known as FRAP1) and PDGFRA, we discovered two novel genes associated with CC: CMPK1 (rs17103186, P = 3.3 × 10−12) and RBP3 (rs11204213 [Val884Met], P = 1.1 × 10−13). The missense RBP3 SNP, rs11204213, was also associated with axial length (AL) (P = 4.2 × 10−6) and had larger effects on both CC and AL compared with other SNPs. The index SNPs at the four indicated loci explained 1.9% of CC variance across the Stages 1 and 2 cohorts, while 33.8% of CC variance was explained by the genome-wide imputation data. We identified two novel genes influencing CC, which are related to either corneal shape or eye size. This study provides additional insights into genetic architecture of corneal shape.

INTRODUCTION

Corneal curvature (CC) is an ocular biometric parameter that measures the steepness of the cornea. Abnormalities of the CC is related to clinically important diseases such as corneal astigmatism and keratoconus (1). Moreover, CC is also a key determinant of refractive error and myopia (2).

Both family and twins studies have shown that CC has highly heritability; 48–95% of interpersonal CC variation is due to genetic factors (3,4). However, limited work has been done in seeking for genetic variants that influence CC. A few genome-wide association studies (GWASs) have been conducted. In our previous meta-analysis of GWAS in Singaporean Asians, common variants in MTOR (also known as FRAP1) and PDGFRA were identified as possible markers for CC. The association of PDGFRA was later confirmed in two independent cohorts of European ancestry in UK (5) and Australia (6). In addition to its influences on CC, PDGFRA was associated with corneal astigmatism and ocular axial length (AL) (5,7), and was thus regarded as a regulatory gene of cornea shape and eye size. However, neither of the two European studies replicated the association of MTOR.

In this study, we conducted a GWAS meta-analysis of CC in 7440 multiethnic Asians and recruited three more cohorts (including 2473 Chinese and 2747 Japanese) in the current analysis to discover new genetic loci for CC in Asians.

RESULTS

Genome-wide association discovery

The amounts of genomic variants across the cohorts in the post-quality controls (QCs) SNP array data ranged from 496 611 to 940 294. The 1000 Genomes imputation data contain ∼38 million genomic variants. However, the amounts of variants available for meta-analysis remained some 7 million after the postassociation QC (Table 1).

Table 1.

Characteristics of the study participants

Study N Age Gender CC AL Corneal astigmatism No. of genetic variants
 
λGCc 
Mean (SD) Males/females Mean (SD) N/mean/SD Case/control Array SNPa Imputed SNPb 
Stage 1 
SP2-1M 851 46.7 (10.2) 543/308 7.75 (0.28) – 503/312 940 294 7 325 522 1.006 
SP2-610 1090 48.5 (11.4) 250/840 7.70 (0.28) – 542/333 538 961 7 044 452 1.017 
SP2-550 325 49.6 (12.7) 246/79 7.78 (0.25) – 186/140 496 611 6 960 927 0.982 
SiMES 2138 57.6 (10.7) 1054/1084 7.66 (0.25) 2141/23.57/1.04 1018/1220 555 718 7 258 467 1.015 
SINDI 2124 55.9 (8.8) 1091/1033 7.62 (0.27) 2120/23.41/1.08 825/1314 556 499 7 874 872 0.990 
SCORM 912 10.8 (0.8) 473/439 7.76 (0.24) 926/24.13/1.12 746/167 507 125 6 993 807 1.008 
Subtotal 7440     5593/4427    
Stage 2 
STARS 742 38.9 (5.3) 379/363 7.68 (0.28) 745/24.64/1.51 596/244 530 923 7 040 338 1.012 
SCES 1731 57.6 (9.0) 893/838 7.66 (0.25) 1720/23.95/1.31 1177/697 535 781 7 038 767 1.011 
Subtotal 2473     1173/941    
Stage 3 
Nagahama 2747 52.0 (13.8) 933/1814 7.68 (0.26) –     
Study N Age Gender CC AL Corneal astigmatism No. of genetic variants
 
λGCc 
Mean (SD) Males/females Mean (SD) N/mean/SD Case/control Array SNPa Imputed SNPb 
Stage 1 
SP2-1M 851 46.7 (10.2) 543/308 7.75 (0.28) – 503/312 940 294 7 325 522 1.006 
SP2-610 1090 48.5 (11.4) 250/840 7.70 (0.28) – 542/333 538 961 7 044 452 1.017 
SP2-550 325 49.6 (12.7) 246/79 7.78 (0.25) – 186/140 496 611 6 960 927 0.982 
SiMES 2138 57.6 (10.7) 1054/1084 7.66 (0.25) 2141/23.57/1.04 1018/1220 555 718 7 258 467 1.015 
SINDI 2124 55.9 (8.8) 1091/1033 7.62 (0.27) 2120/23.41/1.08 825/1314 556 499 7 874 872 0.990 
SCORM 912 10.8 (0.8) 473/439 7.76 (0.24) 926/24.13/1.12 746/167 507 125 6 993 807 1.008 
Subtotal 7440     5593/4427    
Stage 2 
STARS 742 38.9 (5.3) 379/363 7.68 (0.28) 745/24.64/1.51 596/244 530 923 7 040 338 1.012 
SCES 1731 57.6 (9.0) 893/838 7.66 (0.25) 1720/23.95/1.31 1177/697 535 781 7 038 767 1.011 
Subtotal 2473     1173/941    
Stage 3 
Nagahama 2747 52.0 (13.8) 933/1814 7.68 (0.26) –     

CC, corneal curvature; AL, axial length; N, the number of samples included in this study; SD, standard deviation; SP2, Singapore Prospective Study Program; SiMES, Singapore Malay Eye Study; SINDI, Singapore Indian Eye Study; SCORM, Singapore Cohort Study of the Risk Factors for Myopia; STARS, Strabismus, Amblyopia and Refractive Error Study in Singapore; SCES, Singapore Chinese Eye Study. Nagahama, Nagahama Prospective Genome Cohort for Comprehensive Human Bioscience.

aThe number of array SNPs post-QC is listed.

bThe imputed variants was filtered by MAF ≥ 0.01, call rate ≥0.95, Hardy–Weinberg equilibrium P-value > 1 × 10−6 and imputation quality ≥0.5.

cThe inflation factor of each study in Stage 1 was calculated by genomic control.

Of the eight cohorts across Stages 1 and 2, the maximal inflation factor was 1.017, and Stage 1 meta-analyzed results had an inflation factor of 1.000, indicating there was no obvious confounding in our analysis (Table 1). The quantile–quantile (Q–Q) plot showed strong departure from the null hypothesis of no association between the genomic variants and CC measurements. This trend persisted to some extent, even after we removed the variants near the known associated loci, MTOR and PDGFRA (Supplementary Material, Fig. S1). With the Stages 1 and 2 cohorts alone (n = 9913), we have 74 or 99% power to detect the variants that explained 0.5 or 1% of phenotype variance.

The GWAS meta-analysis in Stage 1 revealed three loci with genome-wide significance (genomic controlled P≤ 5 × 10−8, Fig. 1). Of these, MTOR and PDFGRA have been known to be associated with CC (5,6,8). We calculated the linkage disequilibrium (LD) using the ASN panel (CHB, JPT and CHS) of 1000 Genomes Phase 1 haplotypes. The index SNP near MTOR gene, rs74225573, was in LD (r2 = 0.87; D′ = 0.93) with the known index SNP rs17036350, so was the index SNP near PDFGRA gene rs1800813 with the known index rs2114039 (r2 = 0.68; D′ = 1). In addition, we identified one novel genome-wide significant loci near CMPK1 (rs60078183, P = 7.1 × 10−9, Fig. 1 and Table 2).

Table 2.

Association of the index SNPs with corneal curvature

SNP CHR BP Gene Allelesa Stageb N EAF Effect StdErr P-value 
rs74225573 11 280 994 MTOR C/T 7440 0.19 0.16 0.02 7.4E−14 
     2473 0.13 0.19 0.04 8.8E−06 
     1+2 9913 0.18 0.16 0.02 5.7E−18 
rs12139042 11 167 146 MTOR A/G 7440 0.12 0.15 0.02 3.4E−09 
     2473 0.08 0.18 0.05 5.2E−04 
     1+2 9913 0.11 0.15 0.02 1.0E−11 
rs17103186 47 854 514 CMPK1 C/T 7436 0.82 0.11 0.02 4.1E−07 
     2473 0.80 0.12 0.03 6.0E−04 
     1+2 9909 0.81 0.11 0.02 1.2E−09 
     2747 0.86 0.13 0.04 8.2E−04 
     1+2+3 12331 0.82 0.11 0.02 3.3E−12 
rs60078183 47 857 307 CMPK1 G/A 7440 0.82 0.13 0.02 7.1E−09 
     2473 0.81 0.13 0.04 2.2E−04 
     1+2 9913 0.82 0.13 0.02 8.2E−12 
rs1800813 55 094 467 PDGFRA G/A 7440 0.76 0.12 0.02 3.3E−10 
     2473 0.81 0.02 0.04 5.9E−01 
     1+2 9913 0.77 0.10 0.02 7.5E−09 
rs2280823 1 897 657 ARHGEF10 G/A 7437 0.22 0.09 0.02 1.4E−06 
     2472 0.24 0.05 0.03 1.1E−01 
     1+2 9909 0.22 0.08 0.02 7.3E−07 
rs117523291 105 907 982 8q22.3c G/C 5316 0.97 0.26 0.05 1.6E−06 
     2473 0.97 0.04 0.08 5.7E−01 
     1+2 7789 0.97 0.19 0.04 2.2E−05 
rs11204213 10 48 388 228 RBP3 T/C 5316 0.04 0.27 0.05 3.6E−07 
     2473 0.04 0.33 0.08 1.7E−05 
     1+2 7789 0.04 0.29 0.04 3.8E−11 
     2747 0.02 0.29 0.09 9.9E−04 
     1+2+3 10211 0.04 0.29 0.04 1.1E−13 
rs7183510 15 25 129 974 SNRPN G/T 7440 0.16 0.15 0.03 2.5E−07 
     2473 0.04 −0.15 0.08 4.7E−02 
     1+2 9913 0.14 0.11 0.03 3.3E−05 
SNP CHR BP Gene Allelesa Stageb N EAF Effect StdErr P-value 
rs74225573 11 280 994 MTOR C/T 7440 0.19 0.16 0.02 7.4E−14 
     2473 0.13 0.19 0.04 8.8E−06 
     1+2 9913 0.18 0.16 0.02 5.7E−18 
rs12139042 11 167 146 MTOR A/G 7440 0.12 0.15 0.02 3.4E−09 
     2473 0.08 0.18 0.05 5.2E−04 
     1+2 9913 0.11 0.15 0.02 1.0E−11 
rs17103186 47 854 514 CMPK1 C/T 7436 0.82 0.11 0.02 4.1E−07 
     2473 0.80 0.12 0.03 6.0E−04 
     1+2 9909 0.81 0.11 0.02 1.2E−09 
     2747 0.86 0.13 0.04 8.2E−04 
     1+2+3 12331 0.82 0.11 0.02 3.3E−12 
rs60078183 47 857 307 CMPK1 G/A 7440 0.82 0.13 0.02 7.1E−09 
     2473 0.81 0.13 0.04 2.2E−04 
     1+2 9913 0.82 0.13 0.02 8.2E−12 
rs1800813 55 094 467 PDGFRA G/A 7440 0.76 0.12 0.02 3.3E−10 
     2473 0.81 0.02 0.04 5.9E−01 
     1+2 9913 0.77 0.10 0.02 7.5E−09 
rs2280823 1 897 657 ARHGEF10 G/A 7437 0.22 0.09 0.02 1.4E−06 
     2472 0.24 0.05 0.03 1.1E−01 
     1+2 9909 0.22 0.08 0.02 7.3E−07 
rs117523291 105 907 982 8q22.3c G/C 5316 0.97 0.26 0.05 1.6E−06 
     2473 0.97 0.04 0.08 5.7E−01 
     1+2 7789 0.97 0.19 0.04 2.2E−05 
rs11204213 10 48 388 228 RBP3 T/C 5316 0.04 0.27 0.05 3.6E−07 
     2473 0.04 0.33 0.08 1.7E−05 
     1+2 7789 0.04 0.29 0.04 3.8E−11 
     2747 0.02 0.29 0.09 9.9E−04 
     1+2+3 10211 0.04 0.29 0.04 1.1E−13 
rs7183510 15 25 129 974 SNRPN G/T 7440 0.16 0.15 0.03 2.5E−07 
     2473 0.04 −0.15 0.08 4.7E−02 
     1+2 9913 0.14 0.11 0.03 3.3E−05 

CHR, chromosome; BP, position based on GRCh human genome build 37; N, the number of individuals with valid information at each locus; EAF, effect allele frequency; StdErr, the standard error of the effect.

The effect sizes were based on the normalized CC values.

aAlleles were given by the effect allele/other allele.

bStage refers to the meta-analysis of Stages 1, 2 and/or 3.

cThe cytoband location is presented since there is no gene in the 1 Mb region around the index SNP.

Figure 1.

Manhattan plots of the genome-wide association meta-analysis in Stage 1 (n = 7440). The −log 10 of the genomic controlled P-vales are plotted against the genomic positions. The genome-wide significant loci in this stage are labeled in brown. In addition to the index SNPs at the brown loci, the index SNPs in suggestively significant loci (P < 5 × 10−6, labeled in aquamarine) were carried forward to the Stage 2. The dashed line represents P < 5 × 10−8; the dotted line represents P < 5 × 10−6.

Figure 1.

Manhattan plots of the genome-wide association meta-analysis in Stage 1 (n = 7440). The −log 10 of the genomic controlled P-vales are plotted against the genomic positions. The genome-wide significant loci in this stage are labeled in brown. In addition to the index SNPs at the brown loci, the index SNPs in suggestively significant loci (P < 5 × 10−6, labeled in aquamarine) were carried forward to the Stage 2. The dashed line represents P < 5 × 10−8; the dotted line represents P < 5 × 10−6.

Replication

In Stage 2, in addition to the three index SNPs at the two known loci and CMPK1, we further selected five index SNPs with P < 5 × 10−6 for replication, including rs17103186 at CMPK1, rs2280823 at ARHGEF10, rs117523291 at 8q22.3, rs11204213 at RBP3 and rs7183510 at SNRPN gene. In addition to rs60078183, rs17103186 at the CMPK1 loci was followed up because it was directly genotyped across all the cohorts in Stage 1.

The associations of SNPs in MTOR, CMPK1 and RBP3 with CC were further validated in Stage 2 (P < 0.001 in Stage 2 and P < 5 × 10−8 in meta-analysis of Stages 1 and 2; Table 2), the direction of effects being consistent with those in Stage 1. However, the association of the index SNP rs1800813 in PDGFRA was not significant in Stage 2, despite in the same direction with a smaller effect. Nevertheless, this SNP remained genome-wide significant in the meta-analysis of Stage 1 and Stage 2 (P = 7.5 × 10−9). The regional association plots of the four genome-wide significant loci are presented in Figure 2. Unfortunately, the associations of index SNPs in the other three loci (ARHGEF10, 8q22.3 and SNRPN) were either not significant (P > 0.05) or in the opposite direction (rs7183510 at SNRPN) in Stage 2.

Figure 2.

Regional association plots of the associated loci in the meta-analysis of Stages 1 and 2 (n = 9913). The −log 10 of the genomic controlled P-values are plotted against the genomic positions. The variants within the 500 kb regions around (A) MTOR (rs74225573), (B) CMPK1 (rs17103186), (C) PDGFRA (rs1800813) and (D) RBP3 (rs11204213) gene are plotted. The round dots represent the Stage 1 signals, while the square dots represent the meta-analysis of Stage 1 and 2.

Figure 2.

Regional association plots of the associated loci in the meta-analysis of Stages 1 and 2 (n = 9913). The −log 10 of the genomic controlled P-values are plotted against the genomic positions. The variants within the 500 kb regions around (A) MTOR (rs74225573), (B) CMPK1 (rs17103186), (C) PDGFRA (rs1800813) and (D) RBP3 (rs11204213) gene are plotted. The round dots represent the Stage 1 signals, while the square dots represent the meta-analysis of Stage 1 and 2.

The novel associations of CMPK1 (rs17103186) and RBP3 (rs11204213) with CC were further strongly supported in Stage 3 (both P < 0.001). The meta-analysis P-values of the Stages 1, 2 and 3 together were 3.3 × 10−12 and 1.1 × 10−13, respectively (Table 2). The major allele of the index SNP at CMPK1 and the minor allele of the index SNP at RBP3 were associated with increased radius of CC.

Conditional association analysis at the indicated loci

With the aim of identifying any secondary association signal at each indicated locus, we conducted conditional association tests adjusted for the allele dosage of the most significant index SNP at each locus. Meta-analysis of Stages 1 and 2 conditional association results did not reveal any secondary signal at the four associated loci (Supplementary Material, Fig. S2).

CC variance explained by SNPs

The explained phenotype variance was first estimated from genome-wide imputation data. Since no secondary signal was found in the four indicated loci, we also used the most significant index SNPs of the identified loci to estimate the phenotypic variance explained by the genetic variants. Overall, the index SNPs at the four indicated loci explained 1.9% of CC variance across the Stages 1 and 2 cohorts, while 33.8% of CC variance was explained by the genome-wide imputation data (Supplementary Material, Table S1). Hence, the identified loci explained 5.5% of the genome-wide SNP heritability. We further estimated the variance explained by the two novel index SNPs in Stage 3. CMPK1 rs17103186 and RBP3 rs11204213 explained 0.25 and 0.26% of CC variance in this stage.

Associations with AL, corneal astigmatism and spherical equivalent

We further evaluated the effects of CC-associated variants on other clinically related ocular traits, including AL, corneal astigmatism and spherical equivalent (SE) in the study cohorts involved in Stages 1 and 2. The PDGFRA index SNP (rs1800813) was significantly associated with corneal astigmatism (P = 2.6 × 10−6, Supplementary Material, Table S2). This was consistent with the previous report by Fan et al. (7), although the reported index SNPs were not the same. The identified CC SNPs in the CMPK1, PDGFRA and RBP3 loci showed some evidence of association with AL (P < 0.05), but none of them were significantly associated with SE, inferring that these SNP are likely ‘eye size’-determining variants, rather than ‘refractive error’-determining variants.

Gene- and pathway-based analysis

In the gene-based test, we tested 23 436 genes from the HUGO Gene Nomenclature Committee (HGNC; http://www.genenames.org). A total of five loci were identified genome-wide significantly (defined as Pgene based < 2.2 × 10−5, see Materials and Methods for details) associated with CC (Supplementary Material, Table S3). Of them, four were consistently with the results of SNP-based analysis. Moreover, we identified one more locus, WNT7B to be marginally associated with CC (Pgene based = 2.0 × 10−5). We further scanned 4850 pathways/gene sets from the curated gene sets in the Molecular Signatures Database (MSigDB v4.0, http://www.broadinstitute.org/gsea/msigdb/index.jsp). None of these pathways was significantly associated with CC (all Ppathway > 2.1 × 10−5).

DISCUSSION

In this GWAS meta-analysis of CC in 12 660 Asians, we identified two novel loci, CMPK1 and RBP3 for CC, and provided further supportive evidence for the two previously known loci, MTOR and PDGFRA. Although the association of MTOR has not been replicated in non-Asian populations, our Stage 2 cohorts, which were not involved in our previous meta-analysis (8), strongly supported this association (P = 8.8 × 10−6).

Our power to discover associations with weak effects was enhanced by imputation up to the 1000 Genomes cosmopolitan reference panel. The most significant hit near CMPK1 was rs60078183, which was not present in HapMap panels. In fact, it was the only genome-wide significant hit near CMPK1 in Stage 1. The association of CMPK1 and CC was further supported by the evidence from rs17103186, which was directly genotyped across the three stages in this study.

A common missense SNP rs11204213 in exon 2 of RBP3 gene was associated with CC for the first time. It possesses a larger effect as compared with other SNPs in this study (Table 2), and explained 0.5% of the trait variance. rs11204213 is an East Asian specific mutation with the minor allele frequency (MAF) found to be 0.045 in 1000 Genomes ASN panel but 0.001 in EUR panel (only one minor allele was found in 758 individuals of European ancestry). It was genotyped by SNP array in the SP2-1M cohort in Stage 1 and the Nagahama study in Stage 3, and was imputed with high confidence in other cohorts included in the meta-analysis of this SNP. The imputation quality ranging from 0.73 to 0.92 across the meta-analyzed cohorts. In SINDI (Indians), rs11204213 was imputed with low confidence and frequency (score_info = 0.46, MAF = 0.004), and thus SINDI was excluded from the meta-analysis of this SNP. The association of RBP3 gene with CC was strongly supported by our Stage 3 result. As a result, we validated the associations of the common variants at the two novel loci with CC in over 10 000 individuals.

Furthermore, RBP3 gene was also significantly associated with AL in our analysis. Each copy of the rare allele rs11204213-T significantly increased AL by 0.24 mm. This effect is, again, apparently larger than those of other SNPs in previous studies (9,10). Taking together, RBP3 gene correlates strongly with CC and AL and thus may be involved in the development of eye size and shape.

In contrast, our results showed that the PDGFRA gene was also associated with corneal astigmatism, but not with AL. This suggests that the effect of PDGFRA on ocular biometrics lies in the irregularity of the cornea, and not the length of the eyeball.

We investigated the tissue-specific expression spectrum of RBP3 in human and mouse by querying tissue-specific gene expression databases, TiGER (Supplementary Material, Fig. S3) (11) and TISGED (12). Data showed that RBP3 gene was specifically expressed in human eyes and mouse retina. The specificity measurement of RBP3 expression in TISGED was 0.99. Interestingly, the expression level of RBP3 gene has been found to be tuned down in the early stage of diabetic retinopathy (13), indicating the potential role of RBP3 in protecting against the neuroretinal degeneration.

RBP3 (retinoid-binding protein 3) was believed to be expressed specifically in retina, and is responsible for the transportation of retinoid between retinal pigment epithelium and photoreceptors (14,15). In human, rs11204213 induces the amino acid change of Val884Met, a residual in the third module among the four tandem homologous modules. Protein sequence blast (Supplementary Material, Fig. S4) revealed that the amino acid 884 is homologous to the 863 amino acid in Xenopus lavis RBP3 protein. The prediction of the three-dimensional structure (http://raptorx.uchicago.edu/) of the third module illustrated clearly the position of this amino acid in the 3D structure. Unfortunately, it is not located in the putative retinoid-binding site (Supplementary Material, Fig. S5) (16). Moreover, the mutated protein amino acid sequence did not disrupt the predicted protein folding (data not shown).

Using Encyclopedia of DNA Elements data (17) and HaploReg (18), we identified variants in LD with the index SNPs (r2 > 0.8 in the 1000 Genomes ASN individuals) within each of the identified CC loci and applied functional annotations relevant to the regulation of transcription (Supplementary Material, Table S4). The index SNP (rs17103186) at the CMPK1 locus lies in a region with potential regulatory activities, such as DNase I hypersensitivity, modification of histone marks and binding of transcription factors, and the missense index SNP of RBP3 (rs11204213) is located within a DNase I hypersensitive site, thus suggesting their roles in regulating transcription.

All the genes we identified in SNP-based association analysis were further supported by the gene-based association analysis (Supplementary Material, Table S3). WNT7B was also marginally associated with CC in our gene-based analysis; however, there has been little evidence of this gene's connection with eye development or ocular diseases in the literature.

In conclusion, our study discovered two novel genes associated with CC in Asian populations and expanded our knowledge of genetic factors underlying this ocular trait associated with important clinical diseases such as corneal astigmatism and myopia. In addition to the statistical evidences from >12 000 samples, the top SNP of RBP3 leads to a missense mutation and the gene has eye tissue-specific expression and may be the potential targets of the functional studies in the future. Moreover, the missense RBP3 SNP has larger effects on CC and AL, as compared with other SNPs. This may be of particular interest in the growth of eyeball or the development of related phenotypes such as myopia.

MATERIALS AND METHODS

Study design

We performed a three-stage meta-analysis (Table 1). In Stage 1, we conducted a trans-ethnic GWAS meta-analysis in four Asian cohorts (SP2, SiMES, SINDI and SCORM), which involved 7440 individuals, comprising 3178 Chinese, 2138 Malays and 2124 Asian Indians. Although individuals in SCORM were school children while the other cohorts were adults, previous studies supported the stability of CC from an early age (19), and the distribution of CC in SCROM was similar to the adults in the other cohorts (Table 1). The positive associated signals of P <5 × 10−6 from Stage 1 were brought forward for replication in two independent cohorts (STARS and SCES) in Stage 2, comprising 2473 Chinese individuals. In Stage 3, we further tested the index SNPs at the novel loci which were replicated in Stage 2, in 2747 Japanese from the Nagahama study, using the SNP array genotypes. Overall, this GWAS meta-analysis included 12 660 Asian individuals, comprising 5651 Chinese, 2747 Japanese, 2138 Malays and 2124 Asian Indians. The detailed study descriptions are presented in Supplementary Material. All participant studies were population based.

Phenotype measurement and definition

All participating cohorts used similar protocols for keratometry and other ocular biometric measurements. The protocols have been described in detail elsewhere (7,9). In brief, in STARS, SiMES, SCES and SINDI, the horizontal and vertical radii of CC, as well as AL were determined using a non-contact partial coherence laser interferometry (IOLMaster, Carl Zeiss Meditec AG, Jena, Germany). Refraction in all cohorts and CC in SP2 and SCORM subjects were measured with an Auto Ref-Keratometer (RK-5 Autorefractor Keratometer; Canon Inc., Tokyo, Japan). For SCORM subjects, A-scan ultrasound biometry (Echoscan US-800, Nidek Co, Tokyo, Japan) was employed for the measurement of AL. AL measurement was not available for SP2. For SiMES, SCES and SINDI, samples were excluded if they have undergone a cataract surgery before.

Both AL and CC were reported in millimeter. Corneal power in diopters was calculated as 337.5 divided by CC in mm. Corneal astigmatism was calculated as the difference in corneal cylinder power between the flattest and steepest meridians. SE was calculated according to the standard formula (SE = sphere + 1/2 cylinder, in diopter). Corneal astigmatism cases were defined as those who have corneal cylinder power difference ≤ −0.75 D, while controls were defined as corneal cylinder power difference between −0.75 and 0 D. The means of CC, SE, AL and corneal cylinder power from individuals' two eyes were used for analysis, while the means of the readings from one eye were used when the readings from another eye is unavailable.

SNP array genotyping and QC

We genotyped total 3072 SiMES, 2953 SINDI, 1952 SCES and 1451 STARS individuals using Illumina HumanHap 610Quad arrays (San Diego, CA, USA). For SCORM, 1116 samples were genotyped on Illumina HumanHap 550Duo array. For SP2, 1467 samples were genotyped on the HumanHap 610Quad (SP2-610), 1016 samples on the HumanHap 1Mduov3 (SP2-1M) and 535 samples on the HumanHap 550Duo (SP2-550). For this study, genotyping in the Nagahama study was done using Illumina HumanOmni2.5M Arrays.

We performed strict QC (7,9). SNPs that cannot be mapped to the forward strand of Genome Reference Consortium human genome build 37 (GRCh37) or have inconsistent allele/position with the variant annotation of the 1000 Genomes Project Phase 1 (March 2012 release) were removed. Within each cohort, a preliminary SNP QC was conducted to keep the high-quality SNPs for the successive sample QC. We removed the monomorphic, non-autosomal and poorly called (missing rete > 5%) SNPs, as well as those showed significant Hardy–Weinberg Equilibrium (HWE; P ≤ 10−6). The sample QC was based on the remaining SNPs. We removed samples of high missing rate (>5%), excessive heterozygosity, cryptic relatedness, discordant ethnic membership or discrepancy between the genetically inferred and clinically recorded gender information. Furthermore, we also excluded the STARS samples with >1% SNPs showing Mendelian error. We calculated the SNP quality metrics again on the basis of the post-QC samples. Finally, the SNPs were qualified with the same criteria as in the preliminary SNP QC.

Imputation

The 1000 Genomes Phase 1 haplotypes were inferred from the whole-genome/exome sequencing data of 1092 individuals from 14 populations around the world (20). To increase our chance to detect the association signals, we imputed the missing variants for all the participant studies separately using the cosmopolitan reference panel. Specifically, the imputation was done in a pre-phase strategy. The SNP array genotype data of each study were split into 22 autosomes and phased into haplotypes by SHAPEIT version 2.644 (21). To facilitate parallel computation, the haplotypes were divided into chunks of ∼2500 SNPs with 250 buffer SNPs at each flank. The indels on Illumina HumanHap 1M Duo array were also included in these chunks. Finally, the chunks were imputed up to the reference panel by minimac (2012-11-16 release) (22). The imputation resulted in over 38 million variants.

SNP-based association testing and meta-analysis

Within each cohort, CC was normalized by the rank-based inverse-normal transformation in the aim to minimize the effect of different measure methods among these studies. The quantitative normalized CC values and raw AL and SE values, as well as the binary case–control status of corneal astigmatism were tested for the association with SNP genotypes. The associations with SNP genotypes were evaluated using linear (for CC and AL) and logistic (for corneal astigmatism) regression analysis assuming additive effect models using SNPTEST (version 2.4.1) (23). In the regression model, age and sex were included as covariates for analysis of CC, SE and corneal astigmatism, while age, sex, height and education level were included as covariates for AL. To minimize the confounding by population structure in each cohort, we calculated the principal components (PCs) using the post-QC common autosomal SNPs of MAF ≥ 0.01 (24). Obvious population structures were found in SiMES (Malays) and SINDI (Indians), but not in the cohorts of Chinese ancestry (Supplementary Material, Fig. S6 and Table S5). Hence, we additionally adjusted the first five PCs for SiMES and SINDI in the regression models for CC and corneal astigmatism. The PCs included for AL were described previously (9). Both post-QC SNP array genotypes and the imputed genotypes were tested, and we reported the association testing results using the SNP array data wherever available. The conditional analysis was done in a similar way, except adjusting for the allele dosage of the index SNP at each locus in addition.

In the meta-analysis of CC in Stages 1 and 2, we conducted a postassociation QC of the variants in each study. The variants of HWE P ≤ 1 × 10−6, MAF < 0.01, call rate < 0.95 were excluded. The imputed variants of SNPTEST score_info <0.5 were also excluded. We combined the association evidence across our cohorts by the fixed effect inverse-variance meta-analysis implemented in METAL (25). The genomic control correction was applied by METAL to the test statistic before combining the results together (26).

Gene- and pathway-based analysis

The gene- and pathway-based genome-wide association test was carried out using KGG version 2.5 (27,28). The SNP-based association P-values were genomic controlled and mapped to the genes in a putative list. SNPs within the 5 kb flanking regions of each gene were also considered as within the gene region. Gene-based P-values were calculated by the GATES algorithm as implemented in KGG software. A sequential Bonferroni-type procedure (29) was employed for multiple testing correction for 23 436 of the 32 630 genes from the HGNC, and genome-wide significance was defined as P < 2.24 × 10−5 by the software. The gene-based P-values were then taken forward to the pathway-based association analysis, in which the genes were grouped into predefined gene sets or pathways in the MSigDB database. The HYST algorithm in KGG was employed for pathway-based association analysis, and significance was defined as P < 2.08 × 10−5 (29).

Explained phenotype variance

We calculated the proportion of the phenotype variance explained by the SNPs using GCTA version 1.24 (30,31). Within each cohort in Stage 1 and Stage 2, the explained variance was estimated using the restricted maximal likelihood method. The overall phenotype variance explained was calculated from a sample-size weighted average across all the cohorts. For the estimates in Stage 3, a linear regression model in R software (http://www.r-project.org) was employed.

SUPPLEMENTARY MATERIAL

Supplementary Material is available at HMG online.

FUNDING

SP2 acknowledges the support of the Yong Loo Lin School of Medicine, the National University Health System and the Life Sciences Institute from the National University of Singapore. We also acknowledge the support from the National Research Foundation (NRF-RF-2010-05), the Biomedical Research Council (BMRC), Singapore under the individual research grants scheme and the National Medical Research Council (NMRC), Singapore under the individual research grant and the clinician scientist award schemes.

The Singapore Chinese Eye Study (SCES), Singapore Malay Eye Study (SiMES) and Singapore Indian Eye Study (SINDI) were supported by the NMRC, Singapore (grants 0796/2003, IRG07nov013, IRG09nov014, NMRC 1176/2008, STaR/0003/2008, CG/SERI/2010) and BMRC, Singapore (08/1/35/19/550 and 09/1/35/19/616).

The Singapore Cohort Study of the Risk Factors for Myopia (SCORM) was supported by a BMRC grant (06/1/21/19/466).

The Strabismus, Amblyopia and Refractive Error Study of Preschool Children (STARS) was supported by an NMRC grant (1176/2008). C.-Y. C. is supported by an award from NMRC (CSA/033/2012). The Singapore Tissue Network and the Genome Institute of Singapore, Agency for Science, Technology and Research, Singapore provided services for tissue archival and genotyping, respectively.

ACKNOWLEDGEMENTS

Nagahama study also acknowledges the efforts of Takeo Nakayama (Department of Health Informatics, Kyoto University School of Public Health, Kyoto, Japan), Akihiro Sekine (Department of Genome Informatics, Kyoto University School of Public Health, Kyoto, Japan), Shinji Kosugi (Department of Medical Ethics, Kyoto University School of Public Health, Kyoto, Japan), Yasuharu Tabara (Center for Genomic Medicine, Graduate School of Medicine, Kyoto University) and Fumihiko Matsuda (Center for Genomic Medicine, Graduate School of Medicine, Kyoto University).

Conflict of Intereststatement. None declared.

REFERENCES

1
Rabinowitz
Y.S.
Keratoconus
Surv. Ophthalmol.
 
1998
42
297
319
2
Carney
L.G.
Mainstone
J.C.
Henderson
B.A.
Corneal topography and myopia. A cross-sectional study
Invest. Ophthalmol. Vis. Sci.
 
1997
38
311
320
3
Klein
A.P.
Suktitipat
B.
Duggal
P.
Lee
K.E.
Klein
R.
Bailey-Wilson
J.E.
Klein
B.E.
Heritability analysis of spherical equivalent, axial length, corneal curvature, and anterior chamber depth in the Beaver Dam Eye Study
Arch. Ophthalmol.
 
2009
127
649
655
4
Cagigrigoriu
A.
Gregori
D.
Cortassa
F.
Catena
F.
Marra
A.
Heritability of corneal curvature and astigmatism: a videokeratographic child-parent comparison study
Cornea
 
2007
26
907
912
5
Guggenheim
J.A.
McMahon
G.
Kemp
J.P.
Akhtar
S.
St Pourcain
B.
Northstone
K.
Ring
S.M.
Evans
D.M.
Smith
G.D.
Timpson
N.J.
et al.  
A genome-wide association study for corneal curvature identifies the platelet-derived growth factor receptor alpha gene as a quantitative trait locus for eye size in white Europeans
Mol. Vis.
 
2013
19
243
253
6
Mishra
A.
Yazar
S.
Hewitt
A.W.
Mountain
J.A.
Ang
W.
Pennell
C.E.
Martin
N.G.
Montgomery
G.W.
Hammond
C.J.
Young
T.L.
et al.  
Genetic variants near PDGFRA are associated with corneal curvature in Australians
Invest. Ophthalmol. Vis. Sci.
 
2012
53
7131
7136
7
Fan
Q.
Zhou
X.
Khor
C.C.
Cheng
C.Y.
Goh
L.K.
Sim
X.
Tay
W.T.
Li
Y.J.
Ong
R.T.
Suo
C.
et al.  
Genome-wide meta-analysis of five Asian cohorts identifies PDGFRA as a susceptibility locus for corneal astigmatism
PLoS Genet.
 
2011
7
e1002402
8
Han
S.
Chen
P.
Fan
Q.
Khor
C.C.
Sim
X.
Tay
W.T.
Ong
R.T.
Suo
C.
Goh
L.K.
Lavanya
R.
et al.  
Association of variants in FRAP1 and PDGFRA with corneal curvature in Asian populations from Singapore
Hum. Mol. Genet.
 
2011
20
3693
3698
9
Fan
Q.
Barathi
V.A.
Cheng
C.Y.
Zhou
X.
Meguro
A.
Nakata
I.
Khor
C.C.
Goh
L.K.
Li
Y.J.
Lim
W.
et al.  
Genetic variants on chromosome 1q41 influence ocular axial length and high myopia
PLoS Genet.
 
2012
8
e1002753
10
Cheng
C.Y.
Schache
M.
Ikram
M.K.
Young
T.L.
Guggenheim
J.A.
Vitart
V.
Macgregor
S.
Verhoeven
V.J.
Barathi
V.A.
Liao
J.
et al.  
Nine loci for ocular axial length identified through genome-wide association studies, including shared loci with refractive error
Am. J. Hum. Genet.
 
2013
93
264
277
11
Liu
X.
Yu
X.
Zack
D.J.
Zhu
H.
Qian
J.
TiGER: a database for tissue-specific gene expression and regulation
BMC Bioinformatics
 
2008
9
271
12
Xiao
S.J.
Zhang
C.
Zou
Q.
Ji
Z.L.
TiSGeD: a database for tissue-specific genes
Bioinformatics
 
2010
26
1273
1275
13
Garcia-Ramirez
M.
Hernandez
C.
Villarroel
M.
Canals
F.
Alonso
M.A.
Fortuny
R.
Masmiquel
L.
Navarro
A.
Garcia-Arumi
J.
Simo
R.
Interphotoreceptor retinoid-binding protein (IRBP) is downregulated at early stages of diabetic retinopathy
Diabetologia
 
2009
52
2633
2641
14
Gonzalez-Fernandez
F.
Interphotoreceptor retinoid-binding protein – an old gene for new eyes
Vis. Res.
 
2003
43
3021
3036
15
Wu
Q.
Blakeley
L.R.
Cornwall
M.C.
Crouch
R.K.
Wiggert
B.N.
Koutalos
Y.
Interphotoreceptor retinoid-binding protein is the physiologically relevant carrier that removes retinol from rod photoreceptor outer segments
Biochemistry
 
2007
46
8669
8679
16
Gonzalez-Fernandez
F.
Baer
C.A.
Ghosh
D.
Module structure of interphotoreceptor retinoid-binding protein (IRBP) may provide bases for its complex role in the visual cycle - structure/function study of Xenopus IRBP
BMC Biochem.
 
2007
8
15
17
Encode Project Consortium
The ENCODE (ENCyclopedia Of DNA Elements) Project
Science
 
2004
306
636
640
18
Ward
L.D.
Kellis
M.
HaploReg: a resource for exploring chromatin states, conservation, and regulatory motif alterations within sets of genetically linked variants
Nucleic Acids Res.
 
2012
40
D930
D934
19
Wong
H.B.
Machin
D.
Tan
S.B.
Wong
T.Y.
Saw
S.M.
Ocular component growth curves among Singaporean children with different refractive error status
Invest. Ophthalmol. Vis. Sci.
 
2010
51
1341
1347
20
Abecasis
G.R.
Auton
A.
Brooks
L.D.
DePristo
M.A.
Durbin
R.M.
Handsaker
R.E.
Kang
H.M.
Marth
G.T.
McVean
G.A.
An integrated map of genetic variation from 1,092 human genomes
Nature
 
2012
491
56
65
21
Delaneau
O.
Zagury
J.F.
Marchini
J.
Improved whole-chromosome phasing for disease and population genetic studies
Nat. Methods
 
2013
10
5
6
22
Howie
B.
Fuchsberger
C.
Stephens
M.
Marchini
J.
Abecasis
G.R.
Fast and accurate genotype imputation in genome-wide association studies through pre-phasing
Nat. Genet.
 
2012
44
955
959
23
Marchini
J.
Howie
B.
Myers
S.
McVean
G.
Donnelly
P.
A new multipoint method for genome-wide association studies by imputation of genotypes
Nat. Genet.
 
2007
39
906
913
24
Price
A.L.
Patterson
N.J.
Plenge
R.M.
Weinblatt
M.E.
Shadick
N.A.
Reich
D.
Principal components analysis corrects for stratification in genome-wide association studies
Nat. Genet.
 
2006
38
904
909
25
Willer
C.J.
Li
Y.
Abecasis
G.R.
METAL: fast and efficient meta-analysis of genomewide association scans
Bioinformatics
 
2010
26
2190
2191
26
Devlin
B.
Roeder
K.
Genomic control for association studies
Biometrics
 
1999
55
997
1004
27
Li
M.X.
Sham
P.C.
Cherny
S.S.
Song
Y.Q.
A knowledge-based weighting framework to boost the power of genome-wide association studies
PLoS ONE
 
2010
5
e14480
28
Li
M.X.
Gui
H.S.
Kwan
J.S.
Sham
P.C.
GATES: a rapid and powerful gene-based association test using extended Simes procedure
Am. J. Hum. Genet.
 
2011
88
283
293
29
Benjamini
Y.
Hochberg
Y.
Controlling the false discovery rate: a practical and powerful approach to multiple testing
J. R. Stat. Soc. Ser. B – Methodol.
 
1995
57
289
300
30
Yang
J.
Benyamin
B.
McEvoy
B.P.
Gordon
S.
Henders
A.K.
Nyholt
D.R.
Madden
P.A.
Heath
A.C.
Martin
N.G.
Montgomery
G.W.
et al.  
Common SNPs explain a large proportion of the heritability for human height
Nat. Genet.
 
2010
42
565
569
31
Yang
J.
Lee
S.H.
Goddard
M.E.
Visscher
P.M.
GCTA: a tool for genome-wide complex trait analysis
Am. J. Hum. Genet.
 
2011
88
76
82

Author notes

These authors contributed equally to this work.