Scleral HIF-1α is a prominent regulatory candidate for genetic and environmental interactions in human myopia pathogenesis.

BACKGROUND: Myopia is a good model for understanding the interaction between genetics and environmental stimuli. Here we dissect the biological processes affecting myopia progression.
METHODS: Human Genetic Analyses: (1) gene set analysis (GSA) of new genome wide association study (GWAS) data for 593 individuals with high myopia (refraction ≤ -6 diopters [D]); (2) over-representation analysis (ORA) of 196 genes with de novo mutations, identified by whole genome sequencing of 45 high-myopia trio families, and (3) ORA of 284 previously reported myopia risk genes. Contributions of the enriched signaling pathways in mediating the genetic and environmental interactions during myopia development were investigated in vivo and in vitro.
RESULTS: All three genetic analyses showed significant enrichment of four KEGG signaling pathways, including amphetamine addiction, extracellular matrix (ECM) receptor interaction, neuroactive ligand-receptor interaction, and regulation of actin cytoskeleton pathways. In individuals with extremely high myopia (refraction ≤ -10 D), the GSA of GWAS data revealed significant enrichment of the HIF-1α signaling pathway. Using human scleral fibroblasts, silencing the key nodal genes within protein-protein interaction networks for the enriched pathways antagonized the hypoxia-induced increase in myofibroblast transdifferentiation. In mice, scleral HIF-1α downregulation led to hyperopia, whereas upregulation resulted in myopia. In human subjects, near work, a risk factor for myopia, significantly decreased choroidal blood perfusion, which might cause scleral hypoxia.
INTERPRETATION: Our study implicated the HIF-1α signaling pathway in promoting human myopia through mediating interactions between genetic and environmental factors. FUNDING: National Natural Science Foundation of China grants; Natural Science Foundation of Zhejiang Province.

Evidence before this study

Myopia is caused by a combination of genetic and environmental factors; however, the mechanisms underlying their interactions remain unclear. Even though it is widely recognized that excessive near work is a risk factor for human myopia, the underlying mechanisms controlling this process are poorly understood. We previously proposed that scleral hypoxia leads to scleral extracellular matrix (ECM) remodeling, which in turn, is followed by myopia development in experimenal myopia models. However, the possibility that this stress contributes to human myopia needs to be further investigated.

Added value of this study

The results show that the model previously generated to account for myopia development in mice and guinea pigs is indeed translatable to humans. The derived human model postulates that genetic variations in the hypoxia-inducible factor-1α (HIF-1α) signaling pathway likely promote the development of extremely high myopia (refraction ≤-10 diopter [D]). High myopia with -10 D ≤ refraction ≤ -6 D is likely to be caused by the combined effects of environmental factors and genetic alterations in the regulation of the actin cytoskeleton and/or ECM receptor interaction pathways. Furthermore, our study highlights the involvement of the scleral HIF-1α signaling pathway in promoting human myopia through modulating interactions between genetic and environmental factors. Moreover, our results suggest that near work significantly reduces the choroidal blood perfusion, which might lead to scleral hypoxia and ECM remodeling. Taken together, this well-grounded human model provides new insight that may help identify novel targets for suppressing myopia development in a clinical setting.

Implications of all the available evidence

Our results, taken together with previous findings, show that scleral hypoxia promotes myopia development, not only in the experimental animal models but also in humans, by modulating the interaction between genetic and environmental factors.

Alt-text: Unlabelled box

Introduction

Myopia is the most prevalent refractive error and a leading cause of visual impairment worldwide [1,2]. In the vast majority of cases, this condition is characterized by excessive increases in ocular axial length (AL) and vitreous chamber depth (VCD). Myopia of −6.00 diopters (D) or worse, called high myopia, often causes visual impairment due to complications such as posterior staphyloma, choroidal neovascularization, and retinal detachment [3]. However, the underlying etiology and pathogenesis of myopia require clarification. Both genetic anomalies [4,5], and environmental risk factors, e.g., intensive near work such as reading/writing [6], insufficient time spent outdoors [7], and extended exposure to low levels of illumination [8], are implicated in myopia development. Thus, identifying interactions between gentic and environmental factors is a viable approach for improving our understanding of this complex condition.

In recent decades, many genetic factors that contribute to myopia development have been identified. Family-based exome sequencing, meta-analyses, and genome-wide association studies (GWAS) of unrelated individuals have led to the identification of 284 candidate myopia risk genes [9,10]. Tedja et al. conducted a meta-analysis of GWAS data of 160,420 individuals from multiple ancestries with quantitative information on refraction-affecting related functions [9]. Their findings are consistent with the hypothesis that myopia is mediated by a vision-dependent retina-to-sclera signaling cascade. These findings implicate environmental factors in triggering or aggravating myopia development, but provide only limited insight into how genetic and environmental factors interact with one another to modulate this response.

The universality of putative genetic and environmental interactions and causal pathways shared across species can be tested in animal models. Tkatchenko et al. demonstrated interactions between amyloid precursor-like protein 2 (APLP2) and environmental factors during myopia development in humans and mice [11]. However, the biological mechanisms involved in these gene-environment interactions need further investigation. One possibility is that myopia results from abnormalities in complex gene networks. System network analysis has been successful at identifying the causes of various complex diseases [12], [13], [14], [15]. This approach does not merely characterize the roles of individual genes in contributing to a disease phenotype, but it also can provide insights into the combined phenotypic effects of sets of genes, even if the contribution of any one set is relatively small. Each gene set can form transcriptional networks with other gene sets and modulate disease risk in ways that might not be obvious from considering single genes or sets of genes.

We provide here an analysis of the combined effects of genes and the environment on human myopia development. We first performed gene set analysis (GSA) for common genetic variants within each of 321 human regulatory pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) across independent cohorts, which included a high myopia (refraction <-6 D) population in Wenzhou. We then applied over-representation analyses (ORA) in protein-protein interaction (PPI) networks constructed for all the KEGG pathways (KEGG-PPI networks, Fig. 1a) for both the 284 previously reported myopia risk genes and the newly identified genes with de novo mutations in 45 trio families (two non-myopic parents and one high-myopia offspring). The results revealed the critical roles of two KEGG-PPI networks in myopia development. One of the KEGG-PPI networks is centered on the extracellular matrix (ECM)-receptor interaction pathway, and the other is centered on the regulation of the actin cytoskeleton pathway.

Fig. 1

Study design for discovering genetically directed biological processes underlying myopia development. (a) Three different genetic analyses were performed, including GSA of GWAS data obtained from the Wenzhou study, whole genome sequencing for 45 trio families of high myopia to identify functional de novo mutations, and an assemblage of 284 previously reported myopia risk genes. (b) To determine the enrichment of these genetic factors in biological processes, we applied GSA to detect the enrichment of minor effects of genetic polymorphisms in all 321 human regulatory KEGG pathways. Furthermore, we expanded the KEGG pathways by combining the protein-protein interaction (PPI) evidence and constructed the KEGG-PPI networks with the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database. Then, we also applied the ORA strategy to detect the over-representation of 196 genes with functional DNMs and 284 myopia risk genes in KEGG-PPI networks. (c) Numbers in overlapping regions indicate the number of KEGG pathways in common, between and among the analyzed groups. Among all significant gene sets (P<0.05), the same four KEGG pathways were identified as being significant, by analyses of three different data sets: Wenzhou, GWAS data in Wenzhou; DNM196, 196 genes with functional DNMs; and GWAS284, 284 previously reported myopia risk genes. GSA: gene set analysis; GWAS: genome wide association; PPI: protein-protein interaction; ORA: over-represatation analysis.

We then performed similar studies based on individuals within a subset of the Wenzhou cohort who had extremely high myopia (refraction ≤ -10 D). Within that subset, there was a significant enrichment of hypoxia-inducible factor-1α (HIF-1α) signaling pathway. To pursue the possible interaction between scleral hypoxia and the ECM receptor interaction pathway and the regulation of the actin cytoskeleton pathway, we used both in vivo and in vitro approaches to assess the possible roles of the PPIs in response to visual environmental stress during myopia development. In vivo, we assessed the effects of scleral HIF-1α knock-down on form deprivation (FD)-induced myopia development in mice. In vitro, we determined if silencing the key nodes of the interaction networks could antagonized the biological effects of hypoxia on human scleral fibroblasts (HSFs). Overall, the present study reveals, in both a mouse model and humans, the presence of KEGG-PPI networks that could mediate the effects of scleral hypoxia on myopia development. Among the three KEGG-PPI networks we investigated in detail, scleral HIF-1α is a prominent regulatory candidate for genetic and environmental interactions.

Materials and methods

Ethics approval and consent to participate

Studies involving human research participants were approved by the Institutional Review Board at all institutions and organizations, i.e., Wenzhou Medical University (Wenzhou, Zhejiang, China Approval number: KYK-2015-02 and KYK-2018-13), WeGene, Inc. (Shenzhen, China), Zhongshan Ophthalmic Center, Sun Yat-sen University (Guangzhou, China), and Sichuan Provincial Key Laboratory for Human Disease Gene Study (Sichuan Provincial People's Hospital, University of Electronic Science and Technology of China, Chengdu, China). All participants of the three studies were treated in accordance with the tenets of the Declaration of Helsinki, and each signed informed consent.

All animal experiments were approved by the Animal Care and Ethics Committee at Wenzhou Medical University (Approval number: KYK-2018-10) and followed the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Visual Research (https://www.arvo.org/About/policies/statement-for-the-use-of-animals-in-ophthalmic-and-vision-research/). Wild-type C57BL/6 mice and Hif-1α mice were bred in the animal breeding unit at Wenzhou Medical University. They were raised in standard mouse cages (24 × 18 × 13 cm) at 22 ± 2°C, with a 12-h light/12-h dark cycle (light from 8 a.m. to 8 p.m., with luminance of approximately 100–200 lux) and ad libitum access to food and water.

Samples for genome-wide genotyping and whole genome sequencing

For the case-control GWAS, 1,204 non-myopic control subjects and 593 high-myopia subjects were enrolled. Of the non-myopic controls, 271 were from Wenzhou Medical University, and 933 were customers of WeGene, Inc. (Shenzhen, China). For the non-myopic and high myopic participants from Wenzhou, detailed clinical information was available regarding ocular, genetic, or systemic connective tissue disorders, and none of the control participants had a family history of high myopia. Each of the Wenzhou participants underwent complete ophthalmic examinations, including measurement of visual acuity (Topcon RM-8800; Topcon Corp., Tokyo, Japan), AL (Zeiss IOL Master; Carl Zeiss Meditec, Jena, Germany), and spherical equivalent refractive error (KOH3 Keratometer; Nikon, Tokyo, Japan), as well as fundus photography (Canon CR6–45NM Fundus Camera; Canon Inc., Tokyo, Japan). Individuals with systemic diseases or ocular conditions such as cataract or retinal detachment were excluded from the current study because these conditions could alter refraction and skew the data. For the 271 Wenzhou control subjects, the refraction in both eyes was -0.50 D to + 2.00 D. Similar criteria, along with ALs < 24.00 mm, were used to select the participants from WeGene, Inc. For the high myopes, at least one eye had a refraction ≤ -6.00 D and/or AL > 26.00 mm (Supplementary Table 1). Among the 593 high myopes, 369 had extremely high myopia (refraction ≤ -10.00 D).

Zhongshan Ophthalmic Center, Sun Yat-sen University provided a cohort of 139 high myopes (Guangzhou cohort). Sichuan Provincial Key Laboratory for Human Disease Gene Study provided a cohort of 949 high myopes and 2,606 non-myopic controls (Sichuan cohort) .

We also enrolled 45 trio families and the offspring (probands) included 21 females and 24 males, with a mean ± standard deviation age of 17.2 ± 9.3 years. The left and right eye spherical refractions were -10.92 ± -4.54 D and -10.75 ± -4.54 D respectively, and the ALs were 27.99 ± 2.34 mm and 27.98 ± 2.39 mm, respectively.

Genotyping and GSA for GWAS data

DNA samples obtained from individuals enrolled from Wenzhou Medical University were genotyped using the Illumina CoreExome Array (Illumina, San Diego, CA, USA). To combine this dataset with genotype data derived using the Affymetrix Customized Array (Affymetrix, Santa Clara, CA, USA) for the 933 controls provided by WeGene, Inc., additional genotypes were imputed with MACH software (version 1.0.18) using whole genome sequencing (WGS) data of Beijing Chinese Han and Southern Chinese Han individuals from the 1000 Genomes Project (Phase III) as a reference panel [16]. The single nucleotide polymorphism (SNP) locations in this analysis referred to the human reference genome (hg19/GRCH37). All included participants had no more than 5% missing genotypes (call rate >0.95), and the estimated pair-wise kinship coefficients were <0.0221. The coefficients represent the probability that alleles sampled at random from each individual are identical by descent (IBD). Pair-wise kinship coefficients <0.0221 mean that the relationship between the two tested individuals is further than 4th degree relatives, ensuring that all participants were probably unrelated. In further analyses, the SNP genotyping rate was >99% in all individuals, and those with a minor allele frequency (<5%) or a significant departure from the Hardy–Weinberg equilibrium in controls (P≤0.001) were excluded. Finally, 3,825,323 SNPs were included in the analysis of the Wenzhou cohort.

For the Wenzhou cohorts, we estimated the power of single variant GWAS using Genetic Association Study Power Calculator (http://csg.sph.umich.edu/abecasis/gas_power_calculator/). The current sample size in this cohort was 593 cases vs 1,204 controls. The prevalence of high myopia in Chinese University Students is about 20% [17]. Previously, disease risk allele frequencies documented in the GWAS catalog showed a broad range, from 0.1 to 0.9, and we used 0.5 in the power calculation model. Moreover, the highest odds ratio for those risk alleles was no more than 1.3, which was the cutoff for the genotype relative risk. Finally, the expected power for this study based on the multiplicative model was 0.766 at a significance level of P-values <10−4. No sites were identified according to the standard P-value threshold of 5 × 10−8 for single variant GWAS. We therefore applied the GSA strategy to examine the combined effects of SNPs within each pathway. This strategy achieved by MAGMA (v1.07b) can increase the power in comparison with PLINK by 20% in studies on based on 500 cases and 500 controls [18].

Single variant GWAS analysis was performed using PLINK software (version v1.90b3.44, parameters: –adjust –allow-no-sex –assoc —ci 0.95 –geno 0.01 –hwe 0.001 –maf 0.05 –mind 0.01) in the Wenzhou and Sichuan cohorts, with genomic control for population stratification (Genomic inflation est. lambda = 1.00336, median chisq) and Bonferroni correction for multiple testing corrections of associated P-values. Prior to GSA, the curated gene sets for all KEGG pathways (February 2019) (http://www.kegg.jp) were extracted by R package KEGGREST (Tenenbaum D 2019. KEGGREST: Client-side REST access to KEGG. R package version 1.26.1) via KEGG API (application programming interface). In all, curated gene lists for 321 out of the 330 pathways were provided in the database. In GSA conducted by MAGMA (v1.07b), genotyping data from the Wenzhou cohort was the reference specified by –bfile flag in the command for linkage disequilibrium estimation between SNPs. The files containing all SNP P-values in the single variant GWAS served as the input specified by –pval flag with the N parameter indicating the cohort size. Gene-level analysis for contribution of members within each pathway used the snp-wise=mean gene analysis model. The GSA identified pathways of significance (gene-set P-value <0.05, Z-test, MAGMA), as well as genes in these pathways that may contribute (gene P-value <0.05, F-test, MAGMA) to the significance of pathways.

In the Guangzhou cohort of 139 high myopia subjects without a control group, we used linear regression analysis by PLINK software (version v1.90b3.44) to analyze AL as the quantitative trait based on 7,705,601 SNPs for all of the KEGG pathways.

WGS and de novo mutation (DNM) calling

Paired-end (2 × 150 bp) WGS was performed for the 45 trio families (135 subjects total) using the Illumina HiSeq Xten sequencing system (Illumina). The average sequencing coverage was 30 ×, yielding an output of 90 GB raw data. Sequencing reads with a quality lower than Q20 were removed using the FastQC program (v0.10.1) (www.bioinformatics.babraham.ac.uk/projects/fastqc/). For each read, the 3’/5’ adapter sequences were trimmed using the Cutadapt program (v1.9) (https://cutadapt.readthedocs.io/en/stable/). Burrows-Wheeler Aligner software (v0.5.9) mapped the trimmed sequencing reads to the human reference genome (hg19/GRCH37), while SAMtools (v1.0, http://www.htslib.org/download/) was used to remove polymerase chain reaction (PCR) duplicates. Local realignments of reads, quality recalibration, and germline variant calling (single nucleotide variations, SNVs;structure variaions, SVs) were performed using the Genome Analysis Toolkit following their best practice guidance (GATK, v2.7.4, https://gatk.broadinstitute.org/). For non-synonymous variants, we applied filter-based annotation in ANNOVAR software (2016Feb01, https://doc-openbio.readthedocs.io/projects/annovar/) against LJB* (dbNSFP) (dbnsfp35c) and pre-computed functional importance scores (SIFT, Polyphen2, and MutationTester) for all possible non-synonymous mutations in the human genome. After removal of known variants, candidates meeting any of the following criteria were defined as deleterious: 1) SIFT scores <=0.05; 2) Polyphen2 score >=0.957; 3) predictions A ("disease_causing_automatic”) and D ("disease_causing") in MutationTester. Annotation of SVs was performed using Meerkat (0.189, http://compbio.med.harvard.edu/Meerkat/). Candidate functional DNMs in offspring (SNVs, such as missense, 3’ and 5’ UTRs, and intronic; SVs, within exonic regions) from each trio family were detected by the Denovogear program (version 1.1.1, https://github.com/ultimatesource/denovogear), based on realigned bam files. We limited the mutations to the strict mask regions of the 1000 Genomes Project (Phase III) (22). Sequencing data have been deposited in the Beijing Institute of Genomics (BIG) Data Center database http://bigd.big.ac.cn/bioproject/.

ORA of previously reported myopia risk genes and genes with de novo variants in candidate KEGG-PPI networks

In addition to genes with DNMs, we queried the GWAS catalog (October 2017) [19] for candidates using the keyword “Myopia” and combined them with other risk genes from the International Consortium for Refractive Error and Myopia (CREAM) as myopia risk genes [9]. Before performing ORA for these two gene sets, we first constructed the PPI networks for all aforementioned 321 KEGG pathways. They were expanded to incorporate genes that directly interacted with the original pathway members with an interaction score > 0.7, according to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) protein-protein interaction (PPI) database (https://string-db.org). Interactions with high confidence (interaction score > 0.9) for PPI networks were used in network visualization. For the node genes, we performed Gene Ontology (GO) enrichment analysis to examine their over-representation in cellular components (http://geneontology.org/).

To identify the over-represented PPI networks in the myopia risk gene set (GeneSetGWAS) and the DNM gene set identified in WGS for trio families (GeneSetDNM), we first obtained the a priori distribution for the number of members (N) belonging to each network for these pathways in a random gene set (GeneSetR) with the same size (M) as GeneSetGWAS or GeneSetDNM. This was achieved by sampling M genes with replacements from all protein-coding genes (approximately 20,000, repeated 10,000 times) and determining the cumulative probability of observing N within each KEGG PPI network in a GeneSetR. The significance of the over-represented KEGG PPI network was the probability value (P) to find N of this network in GeneSetGWAS or GeneSetDNM.

Induction of monocular form-deprivation myopia (FDM) in mice and ocular biometric measurements

FDM was induced by occluding the right eye with a handmade white translucent occluder as previously described [20]. It was attached to the fur around the eye with polystyrene glue to prevent exposure of the eye to any direct unfiltered light. A collar made of thin plastic was fitted around the neck to prevent the mice from removing the occluder. Ocular refraction was measured in a dark room using an eccentric infrared photorefractor, with each measurement repeated three times to obtain the final mean refraction. Ocular biometric parameters, including VCD, corneal radius of curvature, and AL were measured using a custom-made optical coherence tomography (OCT) apparatus [20]. Each eye was scanned at least three times, and the mean measurements were used as the final values.

Design of animal experiments

Alterations in gene expression during myopia progression were assessed in 4-week-old wild-type male mice (interocular refraction difference, <3.00 D). Because the interocular difference in refraction for each mouse is different, the animals were randomly divided into two groups. The right eyes of mice in one group were form deprived for two days (labeled as FD-T), and the untreated fellow eyes (FD-F) served as the control. While the age-matched mice in the other group were not form deprived, and they served as normal controls (NC). Scleral tissues from these mice were subsequently used to evaluate changes in scleral mRNA expression levels of the enriched candidate genes.

Effects of scleral Hif-1α knock-down on normal refractive development were evaluated in 39 four-week-old Hif-1α mice (B6.129-Hif-1α/J, stock No: 007561, The Jackson Laboratory, Bar Harbor, ME, USA) [21]. Both male and female Hif-1α mice were assigned to three groups: (1) sub-Tenon's capsule injection of an AAV8-packaged Cre overexpression vector (AAV8-Cre, n=16); (2) injection of an AAV8-packaged empty vector (AAV8-Vector Control, n=13); or (3) non-injected controls (Uninjected Control, n=10). Ocular measurements were recorded before, and at two and four weeks after injection of AAV8-Cre, AAV8-Vector, or uninjected control treatment (weeks four, six and eight).

Effects of scleral Hif-1α knock-down on myopia development were evaluated in 44 four-week-old Hif-1α mice. The mice were assigned to three groups: (1) AAV8-Cre overexpression vector plus FD (AAV8-Cre+FD; n=19); (2) AAV8-empty package plus FD (AAV8-Vector+FD, n=14); or (3) non-injected controls subjected to FD (Uninjected Control+FD, n=11). FD was initiated one week after the injections and continued for the subsequent two weeks. Refraction and ocular biometric parameters were measured before and at two weeks after beginning FD (weeks five and seven).

Tissue RNA extraction, cDNA preparation, and quantitative real-time PCR (qRT-PCR)

At the end of the different treatment periods, the mice were euthanized and their eyes were removed. The scleras and corneas were collected by dissection and each was homogenized separately in a ball mill. Total RNA was extracted using the RNeasy Fibrous Tissue Mini Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. Total retinal RNA was extracted using TRIZOLTM reagent (Invitrogen, Waltham, MA, USA) according to the manufacturer’s protocol. After treatment with RQ1 RNase-Free DNase (Promega, Madison, WI, USA), the RNAs were reverse-transcribed, using random primers and M-MLV reverse transcriptase (Promega) to synthesize cDNA. qRT-PCR was performed on an ABI ViiA 7 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) using specific primers and the Power SYBR Green PCR Master Mix (Applied Biosystems). The expression level of each target mRNA, relative to that of 18S rRNA, was determined using threshold cycle values and the 2−ΔΔCt method [22]. The primer sequences were proivded in Supplementary Table 2.

Experimental hypoxia in cultured HSFs

To determine if myopia risk genes could mediate hypoxia-induced myofibroblast transdifferentiation and inhibition of collagen production, HSF cultures were established as previously described [23]. HSFs were cultured in a 5% CO2 humidified incubator at 37°C in Dulbecco's Modified Eagle's Medium (DMEM, Gibco, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (Gibco) and 2 mM GlutaMAXTM (Gibco). The cells were seeded in 6-well culture plates, at a density of 2.5 × 105 cells/mL for 24 hours, and then transfected with small interfering RNAs (siRNAs) using Lipofectamine® RNAiMAX reagent (LipoRNAiMax, Invitrogen), according to the manufacturer's instructions. Sequences of siRNAs targeting the human genes tested (APLP2, CLUSTERIN, GAS6, SPTBN1, and STK11) are shown in Supplementary Table 3. Cells transfected with scrambled-sequence normal control (NC) siRNA served as controls. After incubation for 36 hours, siRNA-transfected cells were exposed to 1% O2 for 12 hours in a three-gas incubator (Thermo Fisher Scientific, Waltham, MA, USA). siRNA-transfected cells incubated under normoxia (21% O2) served as controls. At the end of the designated periods, the cells were harvested, and quantities of HIF-1A, α-SMA, paxillin, and COL1A1 proteins were determined by western blotting.

Western blot analysis

Mouse scleral tissues (pools of two scleras, dissected free as described above) and HSFs were homogenized in RadioImmunoPrecipitation Assay buffer (RIPA Lysis Buffer, Beyotime, Shanghai, China) supplemented with a protease inhibitor cocktail (Roche, Grenzach-Wyhlen, Germany) and 1 mM phenylmethanesulfonyl fluoride (PMSF, Beyotime). Total protein was extracted and the concentration was determined using an Enhanced BCA Protein Assay Kit (Beyotime).

Equal amounts (50 µg) of total protein from mouse scleral samples, or 30 µg of total protein from HSFs, were loaded onto a 10% sodium dodecyl sulfate-polyacrylamide gel, separated by electrophoresis, and then electro-transferred onto a nitrocellulose membrane (Millipore, Billerica, MA, USA). After blocking with 5% non-fat milk for 2 hours at room temperature, the membranes were incubated with primary antibodies overnight at 4°C. Primary antibodies against α-SMA (Cat# ab5694, 1:500, RRID: AB_2223021, Abcam, Cambridge, MA, USA), α-tubulin (Cat# ab52866 1:2000, RRID: AB_869989, Abcam), COL1A1 (Cat# ab88147, 1:2000, RRID: AB_2081873, Abcam), CLUSTERIN (Cat# 34642, 1:1000, RRID:AB_2799057, Cell Signaling Technology, Inc, Danvers, MA, USA), GAS6 (Cat# 67202, 1:2000, RRID:AB_2799720, Cell Signaling Technology), HIF-1A (Cat# 610958, 1:800, RRID:AB_398271, BD Bioscience, San Diego, CA, USA), paxillin (Cat# ab32084 1:500, RRID: AB_779033, Abcam), SPTBN1 Cat# ab124888, 1:1000, RRID: AB_10973548, Abcam), and STK11 (Cat# ab15095, 1:1000, RRID: AB_301641, Abcam) were used. After washing three times for 5 minutes each with tris-buffered saline containing Tween detergent (TBST, 10 mM Tris-HCl, pH 7.2-7.4, 150 mM NaCl, and 0.1% Tween-20), the membranes were incubated for 2 hours at room temperature with 1:2000 IRDye 800CW goat anti-rabbit IgG (Cat# 926-32211, RRID: AB_621843, LI-COR Biosciences, Lincoln, NE, USA,) or 1:2000 IRDye 800CW goat anti-mouse IgG (Cat# 926-32211, RRID: AB_621842, LI-COR Biosciences) antibodies. The membranes were washed again three times in TBST, followed by visualization of protein bands using the Odyssey Infrared Imaging System (LI-COR Biosciences). Densitometric analysis of protein bands was performed using ImageJ software (Version 1.48v, National Institutes of Health [NIH], Bethesda, USA; https://imagej.nih.gov/ij). Values were normalized to those of the α-tubulin loading control. All western blots shown are representative of at least three independent experiments.

Measurement of human choroidal thickness (ChT) and choroidal blood perfusion (ChBP) during accommodation

ChBP in the left eye was measured in 16 non-smoking subjects, age 22 to 28 years (24.9±1.7, mean ± standard deviation), with a spherical equivalent refraction of -2.74 ± 1.45 D and a best corrected visual acuity of 0.0 logMAR or better in both eyes (Supplementary Table 4). There is a close correlation between theChT and refractive diopters in humans [24], and ChT decreases with the severity of myopia [25,26]. To exclude the effect of refractive diopter on the accommodative response and the ChT at baseline, we excluded subjects with astigmatism or anisometropia of >-1.00 D.

The subjects first watched television at a distance of 5 m for 20 min, in a room with mean ambient illuminance of about 100 lux. Then ChT and ChBP were determined using optical coherence tomography angiography (OCTA, Spectralis HRA+OCT System, Heidelberg Engineering, Baden-Württemberg, Germany) with a Badal system [27,28]. This system was composed of a movable target and a fixed positive power lens that was placed at its focal distance away from the eye. This system was used to keep the fixed target at the same size, independent of the target position. Because accommodation occurs in binocular synchronization, when an accommodative stimulus is induced in the right eye using the Badal system, the response of the left eye to the stimulus is the same as the response of the right eye [29]. ChT and ChBP of the left eye were measured, respectively. The control baseline was an accommodative stimulus of 0 D. First, the right eye of each subject was exposed to a 0 D stimulus (the fixation target was placed at the farthest distance away from the right eye), while ChT and ChBP were measured in the left eye. Then, a 6 D hyperopic defocus was imposed upon the right eye by moving the fixation target 3.375 cm closer to the eye, and the ChT and ChBP were again measured in the left eye. All subjects were instructed to focus on the accommodative target and maintain focus as sharp as possible during measurements. All procedures were performed between 3 and 6 p.m.

The percentage of blood perfusion in the choriocapillaris, which provides an estimate of ChBP, was measured by OCTA (scan area of 10° × 10° [3 × 3 mm] centered on the fovea) using the enhanced depth imaging mode and the follow-up modes. The choriocapillaris was defined as all layers from Bruch's membrane to 20 μm behind it, according to the OCTA instrument instructions. Projection artifacts from superficial retinal vessels can confound the estimation of true choriocapillaris flow. Therefore we excluded these projection artifacts by utilizing our own previously described custom MATLAB program [30]. Briefly, we used a linear model to express the relationship between the choriocapillaris and superficial retinal vessels as follows:where I is the flow signal of choriocapillaris obtained from OCTA, Iis the flow signal of the superficial retinal vessels, I is the true flow signal of choriocapillaris, a is the formula coefficient, and b is background noise.

The mean signal values of the foveal avascular zone of the superficial retina were measured and taken as the background noise. It served as the threshold for the retinal layer and choriocapillaris layer. The signal values of the 2.5 mm annular area of the superficial retinal layer and choriocapillaris layer were calculated, respectively. The formula coefficient a is the ratio of the signal value of the choriocapillaris layer to the signal value of the retinal superficial annular area. Thus, when I is greater than the sum of aI and the background noise, it is considered to be the true flow signal of choriocapillaris. The vascular percentage of the choriocapillaris was analyzed in 2.5-mm diameter of en face OCTA images. The choroid layer was defined as extending from the exterior surface of the retinal pigment epithelium to the interior surface of the sclera. ChT was measured manually at the fovea and at 500, 1000, and 1500 µm away from the fovea. The average ChT at all these positions was taken as the sample value.

Statistics for studies based on the experimental models

Data were analyzed using GraphPad Prism software (Version 8.3.0; GraphPad Inc., San Diego, CA, USA). The Shapiro-Wilk normality test was used to analyze the distribution of all datasets. Differences between two groups were assessed using two-tailed Student's t-tests. Multiple comparisons were performed using one-way ANOVA with Bonferroni's post hoc tests (for normally distributed data), Kruskal-Wallis non-parametric test with Dunn's post hoc test (for non-normally distributed data), or two-way ANOVA with Bonferroni's post hoc tests. A two-way repeated-measure ANOVA (RM-ANOVA) was used to assess the effect of scleral Hif-1α knock-down on normal refractive development in mice at different times. All investigators were routinely blinded to the group allocation during data collection and analysis procedures. P-values <0.05 were considered statistically significant. Other descriptive statistics and tests are provided in the text or figure legends or tables they support.

Availability of data and materials

The raw genotype and phenotype data of the subjects enrolled in Wenzhou Medical University are available through application at https://www.wmubiobank.org. Sequencing data have been deposited in the BIG Data Center database (http://bigd.big.ac.cn/bioproject/, Accession number: subCRA001412). All other data that support our findings are available from the corresponding authors upon request.

Results

Biological processes mediating high myopia: Analysis of human genetic data

In the Wenzhou cohort, the single variant GWAS results did not identify any single sites reaching the 10−8 significance level (P < 10−8). Nevertheless, numerous associated SNP sites were present with moderate significance (Supplementary Fig. 1a-c). To evaluate the roles of diverse biological processes in myopia development, especially those that were affected by most of the common genetic polymorphisms within gene members of these pathways, we performed a GSA for all human regulatory KEGG pathways available at that time. Genetic polymorphisms with moderate or minor associations were over-represented in 35 KEGG pathways (P<0.05, z-test, Supplementary Table 5, Supplementary Fig. 1a), suggesting that they might be functionally important in myopia development.

Rare genetic variants that predispose myopia development can also influence biological functions in ways similar to those caused by common genetic polymorphisms. To assess this possibility, we performed WGS of 45 high myopia trio families and identified 2,024 DNMs with a mutation rate of 0.75 × 10−8. Each proband had an average of 45 DNMs. In all, 196 genes with canditate functional DNMs were identified (Materials and Methods, Supplementary Table 6). Among them, only KCNMA1 was reported previously as a high myopia risk gene [31]. In addition, we obtained 284 previously reported myopia risk genes (Supplementary Table 7) from the International Consortium for Refractive Error and Myopia (CREAM) [9] and the GWAS catalogue (October 2017) [19].

ORA of the 196 genes with DNMs and the 284 previously reported risk genes in the KEGG-PPI networks (Fig. 1b) revealed that 41 over-represented KEGG-PPI networks were identified for genes with DNMs (Supplementary Table 8), and 138 for the myopia risk genes (P<0.05, ORA, Supplementary Table 9). Notably, there were four significant KEGG pathways shared across these datasets, i.e., the amphetamine addiction, ECM-receptor interaction, neuroactive ligand-receptor interaction, and regulation of actin cytoskeleton pathways (Fig. 1c), as candidates for further investigation.

In a recent study, we provided evidence that scleral hypoxia activates the HIF-1α signaling pathway, which in turn, induces myopia development [32]. This response to hypoxia included promotion of myofibroblast transdifferentiation through stimulation of the actin cytoskeleton pathway and collagen remodeling via regulation of the ECM receptor interaction pathway. In the current study, both of these pathways had significant enrichment of genetic polymorphisms that contributed to myopia development in the Wenzhou cohort. The ECM-receptor interaction pathway was validated in the two independent cohorts. It was significantly enriched both the Guangzhou cohort (P<0.05, z-test, Supplementary Table 10) and the the Sichuan cohort (Table 1). In the Wenzhou cohort, over-representation of genetic variations with moderate associations in the HIF-1α signaling pathway were solely associated with extremely high myopia (refraction of ≤ -10 D; P=0.031, z-test, Table 1; Supplementary Table 11), rather than with moderately high myopia (refraction of -6 D to -10 D; P=0.12, z-test, Table 1; Supplementary Table 12). Myopia risk genes were also slightly over-represented in the KEGG-PPI network for the HIF-1α signaling pathway (P=0.049, z-test, ORA, Table 1, Supplementary Table 9). Additionally, over half of the probands in trio families (27/45) had de novo mutations with functional changes detected in KEGG-PPI networks for HIF-1α signaling pathway, regulation of the actin cytoskeleton pathway, and the ECM receptor interaction pathway (Supplementary Table 13).

Table 1

Enrichment analyses of myopia-associated genetic polymorphisms or genes in gene sets for three candidate pathways.

aKEGG Pathways
		Significant pathways in the highly myopic Han populations
		GWASWenzhou		GuangzhouCohort		SichuanCohort		GWAS Wenzhou (refraction ≤-10 D)		GWAS Wenzhou (-10 D ≤refraction≤-6 D)		cGenes with DNMs		dPreviously reported myopia risk genes
		P value	bBeta	P value	bBeta	P value	bBeta	P value	bBeta	P value	bBeta	eGenes	P value	eGenes	P value
Candidate pathways	Hif-1α signaling pathway	0.41	0.17	0.73	-0.054	0.4	0.024	0.031	0.207	0.12	0.13	34	0.62	63	0.049
	Regulation of actin cytoskeleton	0.037	0.4	0.9	-0.077	0.031	0.124	0.149	0.081	0.26	0.05	42	0.047	65	0.003
	ECM receptor interaction	0.007	0.2	0.008	0.23	0.46	0.011	0.07	0.181	0.03	0.22	17	0.017	27	0.001
fNegative controls	Purine metabolism	0.841	-0.064	0.96	-0.139	0.36	-0.028	0.59	0.0226	0.59	-0.02	17	0.05	15	0.63
	Cell cycle	0.458	0.007	0.14	0.0814	0.68	-0.039	0.85	-0.107	0.86	-0.02	28	0.77	49	0.32

GWAS Wenzhou, genome-wide association studies of 1,204 controls, 594 high myopia subjects in Wenzhou; Guangzhou Cohort, 139 pathological myopia subjects in Guangzhou; Sichuan Cohort, 2,606 controls, 949 high myopia subjects in Sichuan; GWAS Wenzhou (refraction≤-10 D), genome-wide association studies of 1,204 controls and 369 extremely high myopia subjects; GWAS Wenzhou (-10 D ≤refraction≤-6 D), genome-wide association studies of 1,204 controls and 225 high myopia subjects (-10 D ≤refraction≤ -6 D); DNMs, de novo mutations identified in 45 trio families;

KEGG pathways: Kyoto Encyclopedia of Genes and Genomes pathway.

Beta, the mean association beta-coefficient of genes in the gene set.

Genes with DNMs: genes with de novo mutations.

284 candidate risk genes revealed by GWAS (genome-wide association studies).

Genes, number of genes that are members of KEGG-PPI networks.

Negative control: Purine metabolism and cell cycle pathways, each with a gene set size similar to the three candidate pathways but are not associated with myopia development. P-values <0.05 shown in italic fonts.

The protein products of genes in these three KEGG pathways are predicted to substantially interact with one another (Fig. 2). In the network, 61 out of the 284 myopia risk genes had more than one protein-protein interaction (Fig. 2a). Additionally, 32 out of the 196 DNM genes in the 45 high myopia trio families (Fig. 2c) and 24 genes contributing most to the significance of three KEGG pathways in the Wenzhou and Guangzhou cohorts (as revealed by GSA) also had more than one protein-protein interaction (Fig. 2c). Furthermore, the node genes with DNMs and the prevously reported myopia risk genes were overrepresented in cellular components such as cell periphery and extracellular region parts (Supplementary Tables 14 and 15).

Fig. 2

Protein-protein interaction (PPI) network constructed for myopia-associated genes in candidate pathways. (a) Myopia risk genes were identified by previous genome-wide association study (GWAS) results. (b) Genes with de novo mutations in whole-genome sequencing of 45 trio families of high myopia. (c) Genes with P-values <0.05 in gene set analysis (GSA) results for three candidate pathways in the Wenzhou and Guangzhou studies. The networks were constructed based on interactions extracted from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) and visualized using Cytoscape (https://cytoscape.org/).

Studies with experimental models

Scleral differential expression of key nodal genes within three KEGG-PPI networks in myopic mice

To evaluate the potential contributions by each of the three KEGG-PPI networks to myopia development, we measured the scleral mRNA expression levels of key network nodes in myopia-induced eyes. We first extracted 171 genes that were involved in KEGG-PPI networks (Fig. 2a–c and Supplementary Table 16). Fifty-five key nodes, i.e., genes with a connection degree >3 (shown in bold fonts in Supplementary Table 16), were chosen for further study. After two days of FD, we measured the scleral mRNA expression levels of these 55 genes, as well as that of APLP2, which was earlier reported to be involved in mediating genetic and environmental interactions during myopia development [11].

After two days of FD, qRT-PCR analysis showed that scleral mRNA expression levels in FD-T eyes were significantly lower for a few candidate genes, i.e., Gpc6, Sptbn1, and Stk11, compared with those in the FD-F or NC eyes (Fig. 3a–d). For other genes, i.e., Bmp4, clusterin, and Gas6, their scleral expression levels were significantly higher in FD-T eyes than in FD-F or NC eyes.

Fig. 3

Differential expression of the key nodal genes in the three KEGG-PPI networks in the scleras of FD-T and FD-F eyes after two days of monocular form deprivation, and in non-treated control (NC) mice. Panels a-d represent the relative mRNA expression levels for different sets of genes. Expression levels were normalized to levels in the NC group. Each data point represents an independent mouse. * P<0.05, ** P<0.01, *** P<0.001; one-way ANOVA with Bonferroni's post hoc tests (for normally distributed data), Kruskal-Wallis non-parametric test with Dunn's post hoc test (for non-normally distributed data). The detailed statistical data, including the estimated differences, the appropriate 95% confidence intervals for the differences, and P-values were given in Supplementary Table 17.

Silencing the key nodes antagonized the effects of hypoxia on HSFs

To determine if silencing of the differentially expressed key nodes in these networks can reverse the biological effects of hypoxia on HSFs, we established an in vitro model system. HSFs were transfected with siRNAs that selectively targeted each of the key nodes: genes APLP2, CLUSTERIN, GAS6, SPTBN1, and STK11. BMP4, and GPC6 were not included because highly effective siRNA sequences were unavailable.

After 36 hours of siRNA transfection, HSFs were incubated under a hypoxic condition (i.e., 1% oxygen) for another 12 hours. Consistent with our previous results [32], western blot analysis showed that the level of type I collagen α1 (COL1A1) protein decreased, whereas the levels of HIF-1A, the focal adhesion protein paxillin, and the myofibroblast marker α-smooth muscle actin (α-SMA) protein expression were significantly increased (Fig. 4a–f; Supplementary Fig. 2a-b, and Supplementary Fig. 3a-b). Knock-down of gene GAS6 expression significantly suppressed the hypoxia-induced upregulation of paxillin and α-SMA proteins, but had no significant effect on hypoxia-induced declines in collagen expression levels (Fig. 4a, c–f). Thus, reducing the expression of the key nodal gene GAS6 suppressed the hypoxia-induced transdifferentiation of HSFs. This inhibitory effect was not achieved by silencing the expression of other nodal genes, i.e., APLP2, CLUSTERIN, SPTBN1, and STK11 (Supplementary Fig. 2a-b, and Supplementary Fig. 3a-b). CLUSTERIN expression was also significantly inhibited by hypoxic treatment (Supplementary Fig. 2b).

Fig. 4

Loss of GAS6 function blunts on hypoxia-dependent events in HSFs. (a) HSFs transfected with siRNAs targeting GAS6 were treated with 1% O2, and protein levels of HIF-1A, GAS6, COL1Α1, paxillin, and α-SMA proteins were analyzed by western blotting. (b-f) The bar graphs represent relative protein levels of HIF-1A (b), GAS6 (c), COL1Α1 (d), paxillin (e), and α-SMA (f), respectively (n=3). Data are expressed as means ± standard error of the means. Mock: blank control (transfection reagent only); NC, normal control (cells transfected with scrambled normal control siRNA). α-tubulin was the loading control. **P<0.05; ***P<0.001; two-way ANOVA with Bonferroni's post hoc tests. The detailed statistical data, including the estimated differences, the appropriate 95% confidence intervals for the differences, and P-values were given in Supplementary Table 17.

Overall, these results demonstrate that key nodes in these KEGG-PPI networks are likely to play critical roles in mediating interactions between genetic (gene expression) and environmental (microenvironmental oxygen levels) during myopia development.

HIF-1α signaling pathway: Master regulator of the KEGG-PPI networks

Hypoxic conditions commonly activate the HIF-1α signaling pathway, which in turn plays a critical role in mediating cellular responses to this imposed stress (24). Therefore, we hypothesized this signaling pathway activates the other two scleral KEGG-PPI networks that underlie myopia development, specifically, the regulation of actin cytoskeleton and ECM receptor interaction pathways. To test this hypothesis, we established a scleral Hif-1α knock-down mouse model.

We assessed the role of the HIF-1α signaling pathway in regulating the KEGG-PPI networks and ocular development in two stages. In the first stage, we determined the effect of scleral HIF-1α knock-down on refractive development under an unobstructed normal visual condition (non-deprived). At four weeks after AAV8-Cre injection, scleral Hif-1α mRNA and protein levels were significantly lower than in the AAV8-Vector Control eyes (Supplementary Fig. 4a and Fig. 5a, b). This procedure selectively reduced scleral Hif-1α gene expression, whereas there was no significant difference between Hif-1α mRNA levels in the cornea and retina of AAV8-Cre eyes, and those in AAV8-Vector Control eyes (Supplementary Fig. 4b, c). In addtion, AAV8-Cre had no statistically significant effects on content of the scleral COL1α1 and α-SMA proteins in the sclera (Fig. 5a, b). Hyperopia developed in all three treatment groups (Uninjected Control, AAV8-Vector Control, and AAV8-Cre) with ageing, but AAV8-Cre eyes were more hyperopic than control eyes (Fig. 5c, d). These differences reached significance by four weeks after AAV8-Cre injection (Fig. 5c, d). AAV8-Cre eyes were slightly more hyperopic than AAV8-Vector Control eyes, but the difference was not significant (Fig. 5c, d). In parallel with changes in refraction, AL increased less in AAV8-Cre–injected eyes than in AAV8-Vector Control eyes, with the difference becoming significant by four weeks after AAV8 injection (Fig. 5e). On the other hand, the increase in VCD was not significantly different between the AAV-Cre eyes and the AAV-Vector Control eyes (Fig. 5f).

Fig. 5

Scleral HIF-1α knock-down shifts refraction towards hyperopia in normal eyes under a normal viewing environment. The right eye of each Hif-1α mouse (Right) was treated with AAV8-Cre or AAV8-Vector Control; the left eye (Left) served as untreated fellow control. Data from age-matched Hif-1α mice that received no injections (Control) are shown for comparison. (a) Levels of scleral HIF-1α, COL1α1, and α-SMA proteins, were determined by western blot analysis four weeks after AAV8 injection; α-tubulin was the loading control. (b) Densitometric quantification of western blot results. Each data point represents one scleral sample pooled from two mice. (c) Refraction of Left and Right eyes in the three groups before and two weeks and four weeks after AAV8 injection. Box-plot diagrams showing the distribution of the data. Lines within the boxes indicate medians; bars, range; black dots, outliers. (d-f) Inter-ocular differences in refraction (d), AL (e), and VCD (f). Data are expressed as means ± standard errors of the mean, except as indicated otherwise (c). Data in b were assessed by one-way ANOVA with Bonferroni's post hoc tests; n=3 and 4 for AAV8-Vector Control and AAV8-Cre respectively. Data in c-f were assessed by two-way RM-ANOVA with Bonferroni's post hoc tests; n=10, 13, and 16 mice for Control, AAV8-Vector Control, and AAV8-Cre respectively. *P<0.05, **P<0.01. The detailed statistical data, including the estimated differences, the appropriate 95% confidence intervals for the differences, and P-values were given in Supplementary Table 17. D, diopter.

In the second stage of this analysis, we assessed the effect of scleral Hif-1α knock-down on myopia development in a visually obstructed environment (i.e., form deprived, FD). After two weeks of FD, the HIF-1α protein level was significantly higher in FD-T eyes than in FD-F eyes in the AAV8-Vector group (Fig. 6a, b), as previously described [32]. However, this increase was suppressed in FD-T eyes of AAV8-Cre mice, compared with FD-F eyes of AAV8-Cre mice or FD-T eyes of AAV8-Vector mice (Fig. 6a, b). To determine if loss of scleral HIF-1α function affected myofibroblast transdifferentiation and ECM remodeling, we measured scleral COL1α1 and α-SMA protein levels. After two weeks of FD, the levels of COL1α1 and α-SMA protein in the AAV8-Vector Control and AAV8-Cre groups were not significantly different between the two groups (Fig. 6a, b). After two weeks of FD, significant myopia was induced in the FD-T eyes compared to the respective FD-F eyes in the uninjected control mice and in the AAV8-Vector groups, but not in the AAV8-Cre group (Fig. 6c). Myopia development was significantly suppressed in parallel with the reduction in scleral HIF-1α protein level (Fig. 6d). However, increases in AL and VCD in AAV8-Cre mice and in AAV8-Vector mice were not statistically significant (Fig. 6e, f).

Fig. 6

Scleral HIF-1α knock-down inhibits myopia development in form deprived eyes. The sub-Tenon's space in the right eye of each Hif-1α mouse was injected with AAV8-Cre or AAV8-Vector Control. After one week, monocular form deprivation (FD) was imposed for two additional weeks in the right eye (FD-T), while the left eye served as non-FD fellow-eye control (FD-F). A separate group of form deprived, non-injected eyes served as uninjected FD-only controls (Control+FD). (a) Western blot analysis of scleral levels of proteins HIF-1α, α-SMA, and COL1α1. α-tubulin was the loading control. (b) Densitometric quantification of western blot results Each data point represents one scleral sample pooled from two mice. (c) Refraction of FD-T and FD-F eyes in the three groups. (d-f) Interocular differences (FD-T eye minus FD-F eye) in refraction (d), AL (e), and VCD (f), before and after two weeks of FD. Data are expressed as means ± standard error of the means (b). Box-plot diagrams showing the distribution of the data; lines within the boxes indicate medians; bars, range; black dots, outliers (c-f). Data in b were assessed by one-way ANOVA with Bonferroni's post hoc test (for normally distributed data), or Kruskal-Wallis with Dunn's post hoc test (for non-normally distributed data); n=4 for each group. Data in c-f were assessed by two-way RM-ANOVA with Bonferroni's post hoc tests; n=11, 14, 19 mice, for Control, AAV8-Vector, and AAV8-Cre, respectively; *P<0.05, **P<0.01. The detailed statistical data, including the estimated differences, the appropriate 95% confidence intervals for the differences, and P-values were given in Supplementary Table 17. D, diopter.

Taken together, these results indicate that the knock-down of scleral HIF-1α shifted the refraction towards hyperopia during normal refraction development, and it suppressed myopia development.

Accommodation thins human ChT and reduces ChBP

Our recent study revealed that environmental manipulation, specifically the induction of myopia by FD or negative lens wear, caused choroidal thinning and blood perfusion reduction. These changes likely cause scleral hypoxia and thus contribute to myopia development [33]. Having identified the essential role of scleral hypoxia in the experimental animal myopia model, our next step was to clarify its potential role in human myopia. Near-work–induced accommodation (focusing the lens of the eye) is widely accepted as a risk factor for myopia development in humans. Therefore, we enrolled 16 volunteers and determined if this environmental stressor affected ChT and ChBP. The subjects first watched television at a distance of 5 meters for 20 minutes; then a 6 D accommodative stimulus (lens) was imposed, and ChT and ChBP were determined by OCTA with a Badal system (Fig. 7a, b). This powerful accommodative stimulus significantly thinned the ChT (Fig. 7c) and reduced the ChBP (Fig. 7d), compared with the initial 0 D stimulus. It would be reasonable to assume that the decline in ChBP during accommodation reduces the oxygen supply to the sclera, rendering it relatively hypoxic.

Fig. 7

Effects of accommodation on choroidal thickness (ChT) and choroidal blood perfusion (ChBP). (a) Schematic diagram of OCTA measurement, under 0 D accommodative stimulus. The white symbol “+” on black background served as the fixation target, and depending upon viewing distance, provided either 0 D defocus, or 6 D hyperopic defocus, causing linked accommodative responses in both eyes. When the stimulus was presented to the right eye, ChT and ChBP of the left eye were examined by OCTA. (b) Schematic diagram of OCTA detection under a 6 D accommodative stimulus. When moving the fixation target, L1 did not move, but L2 moved together with the target. Meanwhile, the left eye turned inward, and the OCTA position was adjusted to maintain a clear, well-focused fundus image. L1: Badal lens, f=75 mm; L2: Badal lens, f=50 mm. (c) ChT and (d) ChBP after 0 D and 6 D accommodation stimuli. Data are expressed as mean ± standard error of the mean. n=16; ** P<0.01; paired two-tailed t-tests. The detailed statistical data, including the estimated differences, the appropriate 95% confidence intervals for the differences, and P-values were given in Supplementary Table 17. D, diopter.

Discussion

A recent genetic analysis of a large cohort supported the notion that a light-dependent signaling cascade within the neural retina and non-neuronal components at the back of the eye contributes to mediating genetic control of human eye growth and myopia development [9]. However, interpreting the GWAS data, to gain insight into the underlying pathophysiological processes, remains a daunting task. In the present study, we first identified significantly enriched KEGG pathways by using GSA of the GWAS data in high myopia populations in Wenzhou. Our newly identified genes with DNMs and the previously reported myopia risk genes may not belong to the KEGG pathways, but instead might modulate functions in these pathways via PPIs. Accordingly, we expanded the KEGG-PPI networks to identify the enrichment of both the genes with DNMs and the myopia risk genes. This approach substaniated that the KEGG-PPI network likely has a potential role in myopia development. Background gene set selection in ORA analysis for GWAS and DNM genes may introduce potential bias in functional interpretation of the results [34]. Here, we chose all protein coding genes as the background with the assumption that they should have interactions with at least one of the other proteins, while the network construction was in accordance with the PPI evidence in the STRING database. Current experimental procedures may fail to detect some moderate or subtle interactions. Therefore, our constructed network may represent a subset within the true entire PPI network for the pathway. A further consideration is that false positives may exist in both GWAS [35] and DNM gene sets [36], thus potentially biasing the ORA results. Hence we chose commonly enriched pathways in all of these analyses, for further investigation.

Within the human KEGG-PPI networks, four pathways, i.e., amphetamine addiction, neuroactive ligand-receptor interaction, regulation of actin cytoskeleton, and ECM receptor interaction pathways were identified as having significant potential to regulate the development of myopia. Specially, in the Wenzhou cohort, over-representation of genetic variations with moderate associations in the HIF-1α signaling pathway were especially important because they were linked solely with extremely high myopia. Our previous study in experimental myopia models showed that scleral hypoxia activates the HIF-1α signaling pathway, which, in turn, contributes to myopia development by promoting scleral myofibroblast transdifferentiation via the regulation of actin cytoskeleton pathway and ECM remodeling via the ECM receptor interaction pathway [32]. Given the known cellular events of scleral ECM remodeling [32,37], we propose that myopia development is initiated by activation of the HIF-1α signaling pathway. In experimental myopia models, that pathway can be triggered by an environmental stress, and in humans, genetic alterations can initiate it, leading to myopia.

In addition to the ECM-receptor interaction and the regulation of actin cytoskeleton pathways, KEGG-PPI networks in the amphetamine addiction and neuroactive ligand-receptor interaction pathways also displayed significant enrichment in all of the genetic analyses. Both of them may contribute to myopia development by modifying retinal dopamine signaling, which plays an essential role in ocular development [38]. This supposition is supported by various studies. Specifically, the dopamine D2-like receptor is a potential target of the neuroactive ligand-receptor interaction pathway [39,40]. Additionally, the amphetamine addiction pathway could in turn regulate dopamine release [41], and thereby regulate ocular growth. Clearly, further studies are required to evaluate the roles played by these and other signaling pathways in controlling myopia onset and progression.

Although our in vivo studies were conducted mainly in experimental mouse models, the candidate risk genes used for genetic enrichments were derived from the human GWAS data for myopia and refractive error. Alterations in the expression of risk genes during FDM in mice further suggest a critical role for the common genes and pathways in myopia development. Additionally, the relevance of this study is supported by our use of HSFs to clarify the effects of hypoxia on changes in cellular phenotype and collagen expression. We showed that hypoxia-induced myofibroblast transdifferentiation was alleviated by silencing GAS6, a key node of the PPI interaction network in HSFs. These findings suggest that the PPI network might be involved in mediating genetic and environmental interactions underlying myopia development in both the mouse myopia model and in humans.

To investigate the potential of GAS6 knock-down to reverse the hypoxia-induced mouse myopia, additional studies are needed to assess the effect of scleral GAS6 knock-down on normal refractive development and myopia development. On the other hand, loss of APLP2 function did not reverse hypoxia-induced increase in myofibroblast transdifferentiation was not observed in APLP2-silenced cells. This negative effect supports the hypothesis put forth in the previous study that APLP2 mediates genetic and environmental interactions by modulating visual image processing and neurotransmitter release in the retina [11]. Furthermore, this negative effect also applies to in CLUSTERIN, SPTBN1, and STK11-silenced cells, suggesting that these genes have only minor effects on myopia development. Another outcome of our findings that are supportive of our hypothesis is that myopia is caused by abnormalities in complex gene networks, rather than in a single gene. Thus, an alteration in any one of these genes alone is not sufficient to cause the phenotypic change. Therefore, the potential for combined interactions in the KEGG-PPI network underlying myopia development needs further investigation.

The onset of scleral hypoxia could be due to a decline in ChBP, which is the highly vascularized layer sandwiched between the retina and the sclera. This condition could occur as a consequence of the choroidal thinning that is observed consistently during myopia development in humans [42], [43], [44], as well as in experimental animal models [45,46]. However, the causal mechanisms, biological effects, and functional significance of choroidal thinning in myopia development remain largely unknown [47]. Our previous results showed that FD and negative-lens induction of myopia caused thinning and reduced blood perfusion in the choroid. It is conceivable that such effects can result in scleral hypoxia [33]. Epidemiological surveys suggest that intensive near work, such as writing, reading, and computer use, is a risk factor for humans myopia [2]; however, some studies disagree because the authors contended that extensive near work is not associated with myopia progression [48,49]. The role of near work in inducing myopia progression remains unclear because there are no experimental studies assessing the association between the near work and myopia progression. In the current study, we demonstrated that near work decreased both the accommodative response and ChBF, We propose that this effect is large enough to reduce the available oxygen supply, and as a result, create a hypoxic environment in the adjacent scleral tissue. Our results are supportive of the notion that reducing time spent performing near work during childhood, when most eye growth occurs, may reduce the risk of myopia development. Further studies are needed to investigate if other environmental factors, e.g., bright light and increasing time spent outdoors, inhibit myopia development by suppression of the signaling pathway controlling biochemical pathway networks linked by PPIs that in turn reuduce ChBF, increase scleral hypoxia, myofibroblast transdifferentiation, and ultimately, ECM remodeling.

Eye growth is regulated by the interplay between different homeostatic control mechanisms that are modulated by changes in retinal signal transduction pathway activity. Visual stimuli trigger cascades of retinal signaling processes. e.g., responses induced by neurotransmitter release that ultimately induce biochemical changes that promote or inhibit scleral ECM remodeling [50]. Excessive or reduced ocular growth causes images of distant objects to be focused either in front of or behind the photoreceptors, resulting in myopia or hyperopia, respectively. Previous studies showed that changes in choroidal thickness move the retina either closer or further away from the optical focal point determined by the combined refractive power of the cornea and lens [50,51]. Specifically, increased choroidal thickness moves the focal plane (hyperopia), in front of the retina whereas choroidal thinning moves focal plane behind the retina (myopia) [51]. In a previous study, we have shown that FD upregulated scleral HIF-1α expression and promoted myopia development, wheras inhibition of this process by anti-hypoxia drugs slowed myopia progression [32]. In the present study, we confirm HIF-1α involveent by showing that loss of its function suppressed myopia progression, while hyperopia developed in scleral HIF-1α knock-down mice under a normal visual environment. These findings document the importance of HIF-1α in regulating scleral homeostasis in the mouse model, and provide further support for the HIF-1α signaling pathway as the master regulator of the PPI network involved in myopia development in humans and experimental animal models. Furthermore, we propose that structural changes in the choroid cause changes in ChBP that are sufficient changes to modulate oxygenation of the scleral microenvironment.

In Fig. 8, we summarize the the work flow along with the key findings of the current study. Based on these data, a novel model is proposed that accounts for how human myopia develops. In human myopia, genetic variations underlying aberrant HIF-1α signaling pathway events are likely to be associated with the development of extremely high myopia (refraction ≤-10 D). Lack of enrichment in the HIF-1α signaling pathway in cases of non-extremely high myopia (-10 D ≤ refraction ≤ -6 D) may be due to the hypoxia that is caused by enoviromental stress, such as near work, reading, and insufficient light exposure in humans [52,53]. This hypothesis is supported by our findings that near-work induced accomodation significanly reduced the ChBP, which likely cause scleral hypoxia in humans. Such an association suggests that environmental conditions can augment changes induced by losses in risk genes function. Such possible interactions posit that high myopia is caused by contributions from both environmental and genetic factors. Our hypothesis is also supported by heritability and familial aggregation studies demonstrating that the prevalence of extremly high refractive errors (i.e. myopia or hyperopia) is more likely to be influenced by genetic factors, e.g., the genetic alteration in HIF-1α signaling pathway in our current study, than the milder forms of these refractive errors, for which the causes may be mainly environmental [54,55]. Therefore, a genetic variation in the HIF-1α signaling pathway might serve as a marker for identifying human extremely high myopia.

Fig. 8

Recapitulation of findings and derived model describing how scleral hypoxia induces myopia. Previous known data are shown in blue font. ChT: Choroidal thickness; ChBP: Choroidal blood perfusion; DNMs: De novo mutations; ECM: Extracellular matrix; GSA: Gene set analysis; GWAS: Genome wide association; KEGG: Kyoto Encyclopedia of Genes and Genomes; ORA: Over-representation analysis; PPI: Protein-protein interaction; WGS: Whole genome sequencing; D: diopter (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

We acknowledge several limitations of our study. First, although we have attempted to combine evidences from multiple study designs, e.g. diverse cohorts, trio families and previous GWAS data, larger population studies are required to further validate the commonly enriched KEGG pathway. Second, a detailed analysis of the contributions by KEGG-PPIs to myopic individuals with the quantitative trait of refraction is needed. Myopia is a complex phenotype in which diverse SNPs only have moderate or minor effects on myopia development [9]. On the other hand, high myopia is a more extreme condition in which genetic factors may contribute more to disease pathogenesis [56]. Third, further in vivo validation for key nodes of KEGG-PPIs in mediating genetic and environmental interaction is needed. Fourth, new experimental evidences of protein-protein interactions keep accumulating and PPI networks will be require continuous updating to in future studies.

Conclusions

Collectively, we hypothesize that scleral hypoxia triggers what we will now call the “Myopia - Scleral Hypoxia-dependent Pathway for modulating Myofibroblast Transdifferentiation and ECM Remodeling”. Specifically, this cascade of signal transduction pathways is a characteristic of myopia development in both mice and humans. Furthermore, our results highlight that hypoxia is a key regulator of the genetic and environmental interactions underlying the development and progression of myopia. The results of this study thus provide a novel working model of myopia development, and a foundation for exploring additional gene networks affecting ocular development. This approach to analyzing genetic and environmental interactions is especially useful for a complex condition such as myopia, in which the etiology cannot be explained by alterations of a single gene or environmental factor.