The GERA cohort comprises 110,266 adult men and women who are consented participants in the Research Program on Genes, Environment, and Health, established for members of the Kaiser Permanente Medical Care Plan, Northern California Region (KPNC) (Banda et al., 2015; Kvale et al., 2015). The current study population included 19,420 GERA participants, 18 years and older, who were of non-Hispanic white, Hispanic/Latino, Asian, or African American race/ethnicity, and who had at least one recorded AL measurement on both eyes during the same visit (Table 1). In KPNC, biometry measurements, including AL, are collected using the Haag-Streit Lenstar 900 device prior to cataract surgery. The mean (in millimeters) AL of an individuals’ two eyes was used for the analyses. All study procedures were approved by the Institutional Review Board of the Kaiser Foundation Research Institute.

GERA DNA samples were genotyped at the Genomics Core Facility of the University of California, San Francisco (UCSF) on four ancestry-specific Affymetrix Axiom arrays (Affymetrix, Santa Clara, CA, USA) optimized for individuals of European, Latino, East Asian, and African American ancestry (Hoffmann et al., 2011a; Hoffmann et al., 2011b). Genotype QC (quality control) procedures were performed on an array-wise basis (Kvale et al., 2015), as follows: SNPs with initial genotyping call rate ≥97%, allele frequency difference (≤0.15) between males and females for autosomal markers, and genotype concordance rate (>0.75) across duplicate samples were included.

Imputation was also conducted on an array-wise basis. For imputation, we additionally removed variants with call rates <90%. Genotypes were pre-phased with Eagle (Loh et al., 2016) v2.3.2, and then imputed with Minimac3 (Das et al., 2016) v2.0.1, using two reference panels: variants were preferred if present in the EGA release of the Haplotype Reference Consortium (HRC; n = 27,165; no indels) reference panel (McCarthy et al., 2016), and from the 1000 Genomes Project Phase III release if not (n = 2,504; including indels) (Birney and Soranzo, 2015). As a QC metric, we used the info r ² from IMPUTE2, which is an estimate of the correlation of the imputed genotype to the true genotype. We considered variants with a good imputation score if r ² ≥ 0.7.

Each of the four GERA ethnic groups (non-Hispanic white, Hispanic/Latino, East Asian, and African-American) were first analyzed individually. We performed a linear regression of AL and each SNP using PLINK (Chang et al., 2015) v1.90 (www.cog-genomics.org/plink/1.9/) with the following covariates: age, sex, and ancestry principal components (PCs). The top ten ancestry PCs and the percentage of Ashkenazi (ASHK) ancestry were included as covariates for the non-Hispanic white ethnic group, while the top six ancestry PCs were included for the three other ethnic groups (Banda et al., 2015). The genomic inflation factor λ was calculated for each GWAS analysis to assess inflation due to population stratification (Yang et al., 2011). Data from each genetic variant were modeled using additive dosages to account for the uncertainty of imputation (Huang et al., 2009). The GWAS analyses were also conducted using a recent approach accounting for relatedness that-fits a whole-genome regression model, implemented in REGENIEv2.0.2 (Mbatchou et al., 2021). This new machine-learning method requires only local segments of the genotype matrix to be loaded in memory, making it faster and more efficient in memory use compared to other methods.

We then performed a meta-analysis of AL in GERA by combining the results of the four ethnic groups using the package “meta” of R (https://www.R-project.org) (Balduzzi et al., 2019). Fixed effects and random effects summary estimates were calculated for an additive model. Heterogeneity index, I ² (0%–100%) and p-value for Cochrane’s Q statistic among ethnic groups were assessed. For each locus, the lead SNP was defined as the most significant variant within a 2 Mb window, and novel loci were defined as those that were located over 1 Mb apart from any previously reported locus.

To test the 5 novel AL-associated SNPs identified in the current study for replication, we evaluated associations in a subset of CREAM participants with AL measurement (Fan et al., 2020). GWAS summary statistics for AL from the study of Fan et al. (Fan et al., 2020), consisting of 10,851 individuals of European and Asian descent from 9 studies, were used for those replication analyses. Among the 5 novel AL-associated SNPs identified in the current study, 2 were available in the CREAM GWAS summary statistics; thus, the Bonferroni threshold was set up at a p ≤ 0.025 (to account for a total of 2 SNPs tested) for associations that replicated. However, the 5 novel loci were visually inspected by generating regional plots using the European descent CREAM dataset. To generate those regional plots, we used LocusZoom web-based plotting tool (Pruim et al., 2010; Boughton et al., 2021).

To potentially identify independent signals within the 15 identified genomic regions, we performed a multi-SNP-based conditional and joint association analysis (COJO) (Yang et al., 2012), which is implemented in the Genome-wide Complex Trait Analysis (GCTA) integrative tool (Yang et al., 2013). This COJO analysis was conducted on the GERA non-Hispanic white ethnic group GWAS analysis results. To calculate linkage disequilibrium (LD) patterns we used 10,000 random samples from GERA non-Hispanic white ethnic group as a reference panel. For this COJO analysis, a p-value less than 5 × 10⁻⁸ was considered significant.

Array heritability estimate was obtained for AL in GERA non-Hispanic white people (the largest ethnic group of GERA) using LD score regression (Zheng et al., 2017).

To prioritize variants within the 15 AL-associated loci identified in the GERA non-Hispanic white ethnic group GWAS analysis, we used a Bayesian approach (CAVIARBF) (Chen et al., 2015). For each of the 15 associated signals, we computed each variant’s capacity to explain the identified signal within a 2 Mb window (±1.0 Mb with respect to the lead SNP). Then, the smallest set of SNPs that included the causal variant with 95% probability (95% credible set) was derived. For this CAVIARBF analysis, we used 10,000 random samples from the GERA non-Hispanic white ethnic group as a reference panel to calculate LD patterns.

To prioritize genes and biological pathways, we conducted gene-based and pathways association analysis using the Versatile Gene-based Association Study - 2 version 2 (VEGAS2v02) web platform (Mishra and Macgregor, 2015). We first performed a gene-based association analysis on the GERA non-Hispanic white ethnic group GWAS analysis results using the default ‘-top 100’ test that uses all (100%) variants assigned to a gene to compute gene-based p-value. We used 1000 Genomes phase 3 European population as the reference panel. As 20,950 genes were tested, the p-value adjusted for Bonferroni correction was set as p < 2.28 × 10⁻⁶ (0.05/20,950).

Second, we performed a pathway association analysis based on VEGAS2 gene-based p-values. We tested enrichment of the genes defined by VEGAS2 in 9,734 pathways or gene-sets derived from the Biosystem’s database (https://vegas2.qimrberghofer.edu.au/biosystems20160324.vegas2pathSYM). We adopted the resampling approach to perform this pathway analysis using VEGAS2 derived gene-based p-values considering the default ‘−10 kbloc’ parameter as previously described (Iglesias et al., 2018). As 9,734 pathways or gene-sets were tested, the p-value adjusted for Bonferroni correction was set as p < 5.14 × 10⁻⁶ (0.05/9,734).

We evaluated the associations of the 17 lead AL-associated SNPs (sixteen from the multiethnic meta-analysis and an additional one from the non-Hispanic white ethnic group GWA analysis) with MSE by linear regression under an additive model, and adjusting for age, sex, and ancestry PCs. As previously described (Hysi et al., 2020; Choquet et al., 2022), GERA individuals had at least one spherical equivalent value measured during routine eye examinations. Most subjects had multiple measures for both eyes. Spherical equivalent was calculated as the sphere + (cylinder/2). The spherical equivalent was selected from the first documented refraction assessment, and the mean of both eyes was used. After excluding participants with histories of cataract surgery (in either eye), refractive surgery, keratitis or corneal diseases, our GERA sample for this analysis included 72,388 individuals with spherical equivalent measurement from four ethnic groups (Supplementary Table S1). As a note, among those with spherical equivalent measurement, 15,503 (79.8%) were included in the GERA AL sample.

LD score regressions (Bulik-Sullivan et al., 2015) were conducted, using the LDSC v1.0.1 command line tool (https://github.com/bulik/ldsc), to estimate the genome-wide genetic correlations (r _g) between AL and MSE and between AL and myopia. We used as input data GWAS summary statistics from our previous GWASs of MSE and myopia conducted in the GERA cohort (Choquet et al., 2022) (including 59,094 individuals with RE measurement, and 19,540 myopia cases and 36,487 controls of non-Hispanic white ethnic group). As previously described (Choquet et al., 2022), in GERA, myopia cases were defined as having a MSE ≤ -0.75 diopters (D) and controls as having a MSE > -0.75 D, corresponding to the definition of ‘any myopia’, as previously used (Hysi et al., 2020). This more myopic threshold (compared to the well-accepted threshold for myopia of MSE ≤ −0.50 D), increases the probability that only ‘‘true’’ myopes are included (Flitcroft et al., 2019).

The GERA cohort is an unselected cohort of adult members of the KPNC integrated healthcare delivery system, with ongoing longitudinal records from vision examinations. For this study, our GERA sample consisted of 19,420 individuals from 4 ethnic groups (85.1% non-Hispanic white ethnic group, 6.2% Hispanic/Latino, 6.2% East Asian, and 2.5% African American) with a measured AL (Table 1). In our GERA sample, East Asian individuals had higher ALs on average than other groups, consistent with previous studies (Ip et al., 2007; Knight et al., 2012).

We identified 16 loci that reached genome-wide significance in the GERA meta-analysis, of which 5 were novel: SLC25A12, near BMP3, RGR, RBFOX1, and MYO5B (Figure 1 ; Supplementary Table S2 ; Table 2). The genomic inflation factor, λ, of 1.105, is reasonable for a sample of this size (Yang et al., 2011) (Supplementary Figure S1). Regional plots of the association signals at the 5 novel loci are presented in Supplementary Figure S2. These 16 associations with AL were also examined in each individual ethnic group (Supplementary Table S2) and no other SNPs at the 16 loci identified in the multiethnic GWA meta-analysis reached genome-wide level of significance in East Asians, Hispanic/Latinos, or African Americans from the GERA cohort. While the GWAS results generated using REGENIE (Mbatchou et al., 2021) were similar compared to the results generated using PLINK, we found that MYO5B was no longer genome-wide significantly associated with AL using REGENIE (lead SNP rs55754534 p-value = 3.0 × 10⁻⁷) (Supplementary Figures S3, S4 ; Supplementary Table S3). Conducting a GWAS of AL in GERA non-Hispanic white people resulted in the identification of one additional genome-wide significant locus, PRSS56, which has been recently reported to be associated with AL in Japanese populations (Fuse et al., 2022) (Supplementary Figures S5, S6 ; Supplementary Table S4).

Of the 5 novel lead SNPs, two were available in the CREAM dataset. While lead SNP rs1353386 at C4orf22-BMP3 replicated at a Bonferroni level of significance (p ≤ 0.025, 0.05/2 available), SNP rs55754534 at MYO5B was nominally associated with AL (p = 0.030) in the combined (European and Asian ancestry) CREAM subset (Table 2). In addition, we generated the regional plots at the 5 novel loci using the CREAM European ancestry subset (Supplementary Figure S7), and reported association results for the lead SNPs in CREAM (even if those are different than lead SNPs in GERA) in Supplementary Table S5. Those CREAM results provide evidence of replication at the novel AL identified loci.

Conditional (COJO) analysis in the GERA non-Hispanic white ethnic group revealed eight additional independent signals within the identified genomic regions, including at the CD55, LYPLAL1-SLC30A10, HAT1, PRSS56, C4orf22-BMP3, RGR-LOC170425, RDH5, and MYO5B-MIR4320 loci (Supplementary Table S6).

We then estimated the array heritability for AL in the GERA non-Hispanic white ethnic group (the largest group of individuals from GERA) using LD score regression (Zheng et al., 2017), and found an overall heritability estimate of 0.37 (s.e. = 0.04).

To prioritize variants within the 15 AL-associated genomic regions identified in the GERA non-Hispanic white ethnic group GWAS analysis, we computed each variant’s ability to explain the observed signal and derived the smallest set of variants that included the causal variant with 95% probability (Chen et al., 2015). The 15 credible sets, corresponding to each of the 15 AL-associated loci, contained from one to 74 variants (355 total variants, Supplementary Table S7). Interestingly, two sets included a unique variant (i.e., KCNQ5-MIR4282 rs7744813, and LAMA2 rs12193446 with 99.6%, and 99.7% posterior probability of being the causal variants, respectively), suggesting that those variants may be the true causal variants.

To prioritize genes within the 15 AL-associated genomic regions identified in the GERA non-Hispanic white ethnic group GWAS analysis, we performed a gene-based association analysis using the Versatile Gene-based Association Study (VEGAS2v02) integrative tool (Mishra and Macgregor, 2015). We identified 11 significant genes for AL after correcting for multiple testing (Bonferroni correction was set as p < 2.28 × 10⁻⁶ (0.05/21,950 genes tested)), including BMP4, PRSS56, and RGR (Supplementary Table S8).

We then conducted a pathway analysis using VEGAS2v02 and based on the gene-based association results to assess enrichment in 9,734 pathways (or gene-sets) derived from the Biosystem’s database. We found that pathways involving the cranial skeleton system development as well as the elastic fiber formation were significantly enriched after correcting for multiple testing (Bonferroni correction was set as p < 5.14 × 10⁻⁶ (0.05/9,734 pathways tested)) (Supplementary Table S9).

As previous studies have reported that AL is significantly higher in eyes with myopia compared to normal eyes (Adams, 1987; Grosvenor, 1988; Tideman et al., 2018), we evaluated the effect estimates of the 17 AL-SNPs identified in the current study between AL and MSE in GERA. Our GERA MSE sample consisted of 72,388 participants (Supplementary Table S1). All of the 17 AL-associated SNPs were significantly associated with MSE after Bonferroni correction (p < 0.05/17 = 2.94 × 10⁻³), including 15 at genome-wide level of significance (Supplementary Table S10). The Pearson correlation coefficient between AL effect size (beta) and the MSE effect size (beta) was −0.86 (p = 1.23 × 10⁻⁵; Figure 2).

We also performed genome-wide genetic correlation analyses to quantify genetic overlap between AL and MSE or myopia risk in GERA non-Hispanic white people. We found that AL was genetically correlated with MSE (r_g = −0.83; SE, 0.04; p = 1.95 × 10⁻⁸⁹) and myopia (r_g = 0.80; SE, 0.05; p = 2.84 × 10⁻⁵⁵). These results suggest considerable shared genetic influences between AL and myopic refractive error.