All animal experiments were approved by the Atlanta Veterans Affairs Institutional Animal Care and Use Committee (protocol V015-14). Male C57BL/6J mice were housed in a standard 12-12 hour light-dark cycle. Food and water were provided ad libitum. From postnatal day (P) 28, mice were exposed to FDM by surgically attaching a head-mounted diffuser goggle over the right eye (17). The left eye remained uncovered and served as an internal control. Refractive error, corneal curvature and axial length were measured at P28, P35, and P42. After dilating the eyes with 1% tropicamide, mice were anesthetized (ketamine 80 mg/kg; xylazine 16 mg/kg). Refractive error was measured with an automated infrared eccentric photorefractor (5), corneal radius of curvature was measured using an automated keratometer (18) and axial length with a spectral domain optical coherence tomography system (Bioptigen, Durham, NC, USA). Myopic shift was calculated as the refractive error difference between the goggled and uncovered eye.

At P42, mice were sacrificed by cervical dislocation during the light phase; their retinal tissue was immediately collected, flash-frozen on dry ice, and stored at −80°C until further processing. Total RNA was extracted using TRIzol and RNeasy Micro kit RNA extraction methods according to the manufacturers' protocols. RNA quantity was measured with a NanoDrop 1000 (Thermo Fisher Scientific), RNA quality was analyzed with the Agilent 2100 Bioanalyzer using the Pico chip (Agilent Technologies), and samples with an RNA integrity number (RIN) >7.5 were used for RNA sequencing (RNA-Seq). RNA was submitted to Emory Integrated Genomics Core. Following poly-A enrichment, 50-base paired-end libraries were prepared and sequenced on the Illumina HiSeq Sequencing System at 50 M reads per sample. The RNA-Seq data of this study are available in the Gene Expression Omnibus repository, accession number GSE284642, and BioProject, accession number PRJNA1200000.

FASTQ files were uploaded to the Galaxy web platform (19). Read quality was analyzed with FastQC (20). Reads filtered for low-quality reads and trimmed using Trimmomatic (21). Transcript abundance was quantified with Salmon (22). Count normalization and differential analysis were conducted with DESeq2 (23) in R (24) and an unadjusted p value of < 0.05 was considered statistically significant, as this analysis is exploratory in nature with an aim to suggest further hypotheses. Differentially regulated pathways and cellular functions were further analyzed using the KEGG database (25) and GO terms (26) and visualized with pathfindR (27) and ggplot2 (28).

A search was conducted on the PubMed, Gene Expression Omnibus (GEO) and BioProject databases using the following search statement: retina^* AND (transcriptom^* OR “RNA seq^*” OR “RNA-seq^*” OR microarray) AND (myopia OR “short^*sighted^*” OR “refractive error” OR “refractive development” OR “ocular growth” OR “eye growth” OR “experimental myopia” OR “lens-induced myopia” OR “image defocus” OR “form^*deprivation”). All searches were conducted on June 5, 2024. The results were screened for studies that met the following inclusion criteria:

With this search, six RNA-Seq datasets were identified (Table 1). The unprocessed data files of identified studies in the FASTQ format were located. In two instances, the studies filled the inclusion criteria, yet the raw data was not found nor accessed after contacting the corresponding author. Identified studies used either the chick or mouse model of myopia, and the samples from different species were analyzed separately.

The FASTQ files were imported to the Galaxy platform (19), filtered for low-quality reads, and trimmed using Trimmomatic (21), specifying the adapter parameters to what was used in each experiment. The quality of individual reads was analyzed with FastQC (20) before and after the trimming process. Samples were excluded if the percentage of duplicates, read count, or GC content was >2 standard deviations from the mean of the respective study. Transcript abundance was quantified with Salmon (22), estimates were aggregated to the gene level.

Non-normalized untransformed count matrices of each study were obtained using tximport package (v 1.28.0) (29) in R. Next, ComBat-seq (30) from the sva package (v 3.48.0) (31) was used to remove batch effects between studies, specifying the study identifier as the batch variable and the treatment group (experimental myopia vs. control) as the grouping variable. Manipulations used for myopia induction were not separated in this meta-analysis. Low-abundance transcripts were eliminated by retaining only genes with >10 counts in the number of samples of the smallest experimental group. Differential gene expression between control and experimental myopia samples was conducted using the DESeq2 package (23). An unadjusted p value of < 0.05 was considered statistically significant, which does increase type I error risk, but as the analysis has inherently high biological variability, a more stringent criteria could eliminate potentially interesting and relevant findings. Transcripts were annotated using the Mus musculus GRCm38 Ensembl v91 release or the Gallus gallus GRCg7b Ensembl v112 release genome assemblies. Differentially expressed pathways and cellular functions were further analyzed using the KEGG database (25), and visualized with pathfindR (27) and ggplot2 (28). Principal component analyses (PCA) of experiments were performed on the variance-stabilized counts of the RNA-Seq data using DESeq2 (23). PCA coordinates were extracted, considering treatment and batch as grouping factors.

Genes associated with refractive errors from genome-wide association studies (GWAS) were retrieved from the GWAS Catalog (32). Single nucleotide variants (SNVs) and their corresponding mapped genes were obtained for traits “myopia” (Experimental Factor Ontology [EFO] trait HP_0000545) and “abnormality of refraction” (EFO trait HP_0000539). Pathways enriched for genes associated with a genetic predisposition to myopia were determined using the KEGG database (25) and visualized with pathfindR (27) and ggplot2 (28).

Form deprivation was induced in three male mice from postnatal day (P) 28 by placing a translucent goggle in front of the right eye with a head-mounted pedestal, while the contralateral left eye remained uncovered and served as a control eye (17). By P42, the form-deprived eyes developed myopia, as their refractive error was on average 1.37 ± 0.87 D [mean ± standard deviation (SD)] more myopic than the contralateral uncovered eyes (P28 vs. P42 p = 0.008, one-way ANOVA with Šidák's correction for multiple comparisons; Figure 1A). We did not observe statistically significant changes in axial length and corneal radius of curvature between the control and form-deprived eyes (two-way repeated measures ANOVA, age x eye interaction effect p > 0.05; Figures 1B, C), which is not uncommon as the magnitude of changes is very small (33, 34).

The retinas of these animals were submitted for RNA-Seq analysis. Upon quality control, one sample was removed as its read duplication level was >2 SD from the mean of the study. Differential gene expression analysis showed that 235 transcripts were differentially expressed between the control and experimental retinas (Figure 2A and Supplementary Table 1). Among the most highly upregulated transcripts were several crystallins (Cry), such as Cryaa, Cryba1, Cryba2, Cryba4, Crybb1, Crybb2 and Crygs. Upregulation of different crystallin transcripts has been also demonstrated previously in chick compound retina/RPE/choroid tissue in response to FDM (35) and in chick retinas after FDM and LIM induction (36). Among the downregulated genes were paraoxonase 1 (Pon1), an antioxidative protein associated with lower activity in AMD (37) and potentially linked to oxidative stress, which has been reported in myopia (38); and retinoic acid early transcript 1E (Raet1e), which is intriguing as retinoic acid signaling has been implicated in myopia (39).

The genes were further analyzed for enriched pathways using the KEGG database and GO biological process terms. The enriched pathways included various cell signaling pathways, e.g. thyroid hormone and Ras signaling pathways. Thyroid hormone signaling has also been implicated in prior studies investigating the retinal response to experimental myopia (10, 40). A key gene in the Ras signaling pathway, RASGRF1, has been associated with myopia (41), and there is evidence that its downstream pathway, the MAPK pathway, is differentially regulated in myopic retinas (11). They also included the regulation of the BMP signaling pathway, which has been associated with experimental myopia in previous studies in the retina and choroid (42), as well as the sclera (43, 44) (Figures 2B, C, Supplementary Tables 2, 3).

To increase the number of experimental samples and thereby increase the statistical power of our analysis, while simultaneously reducing the variation induced by conducting experiments in different laboratories in different species and using slightly different RNA-Seq protocols, we sought to combine and analyze all similar experiments conducted to date. We performed systematic searches in PubMed, GEO, and BioProject to identify studies employing transcriptome profiling to study retinal gene expression in response to experimental myopia in animal models. These searches yielded 71, 42, and 8 results, respectively. After screening the articles for inclusion criteria (see Materials and methods), we determined six suitable studies: two using mice (8, 9), three using chicks (7, 10, 45), and one using marmosets (6) as experimental animals. Four of the studies were included in the meta-analysis, as well as the RNA-Seq experiment presented earlier (Table 1).

There were several global differences between the studies conducted in mouse and chick models regarding the developmental stage of animals and the duration of myopia induction. The experiments conducted in mice used juvenile animals aged between 3 and 4 weeks, while chicks were as young as newly hatched to 7 days old. In mice, myopia was induced for a total of 2 to 4 weeks, while in chicks, the induction was limited to 4–48 h. Four of the five studies induced myopia with form deprivation, and one study in chicks used a −10D lens. Four studies used the contralateral eye as a control eye, while one used age-matched untreated animals. Three of the studies used only male animals, one used male and female, and the sex of animals in one study was not collected (see an overview in Table 1).

Since the chick and mouse samples differed not only in the experimental animal used but also in the developmental stage and myopia induction length, we conducted two separate meta-analyses, one including samples from mouse studies and one including all samples from chick studies. In these two sets of samples, we analyzed the effect of experimental myopia on the retinal transcriptome. A total of 9 control and 8 experimental myopia samples were identified using the mouse model, and 49 control and 50 experimental myopia samples in the chick model. After removing samples with outliers regarding sequencing read quality metrics (see Materials and methods; Supplementary Tables 4, 5), a total of nine control and seven myopia mouse samples, and 46 control and 45 myopia chick samples were included in the downstream analysis.

Samples from all included studies underwent the same analysis pipeline to ensure the results were directly comparable. Principal component analyses (PCA) of the individual studies did not indicate a clear clustering of the experimental samples by myopia treatment (Supplementary Figure 1). In studies that used the contralateral eye as a control, the samples from the same experimental animal were clustered together (Supplementary Figures 1B, C, E). In the only study using experimental animals of both sexes, there was clear clustering based on sex (Supplementary Figure 1D).

Next, we combined the two studies conducted in mice and the three studies conducted in chicks separately by applying batch correction and normalization using ComBat-seq (30). Before batch correction and normalization, we observed the expected batch effects, revealed by the sample distribution in PCA (Figures 3A, C). After data processing, the samples were clustered significantly closer, while a certain level of segregation between samples remained (Figures 3B, D).

The meta-analysis of 16 mouse retinal samples identified 427 differentially regulated genes (unadjusted p value < 0.05, Figure 4A, Supplementary Table 6). Among the enriched KEGG pathways, significant changes were observed in signaling pathways such as TGF-beta (Figures 4B, C, Supplementary Table 7), which is critical in ocular growth and myopia development (42, 46). In addition, among the enriched pathways were dopaminergic, GABAergic and serotonergic synapse and circadian entrainment. The genes contributing to these pathways largely overlapped and included G-protein gamma subunit genes Gng3 and Gng4 (Figure 4C), which contribute to various receptor signaling pathways (47).

A significantly larger number of samples, a total of 91, were included in the meta-analysis of chick retinal responses to experimental myopia. We found that 1,110 genes were differentially regulated (Figure 5A, Supplementary Table 8), a larger proportion of them (54%), with highest changes in significance, were downregulated. Among those genes was vasoactive intestinal polypeptide (VIP), which is expressed in a subset of amacrine cells (48), and its downregulation in response to myopia has been shown in macaque retinas (48), chick retinas (45) and chick retina/RPE complex (49). Furthermore, intravitreal administration of VIP reduced the magnitude of FDM, while VIP antagonists abolished FDM development in chicks (50). The most significantly downregulated gene was urotensin-2B (UTS2B; Benjamini-Hochberg adjusted p value 1.05E-35). UTS2B encodes for urotensin II-related peptide (URP), and the downregulation of prepro-URP has been reported previously in chick retina (45) and retina/RPE complex (49). In addition to its vasoactive properties, URP is also known to induce cell proliferation (51). To our knowledge, the effect of URP in the retina is not well understood. Also among the top downregulated genes was brain-derived neurotrophic factor (BDNF), which has been established as neuroprotective agent in the retina (52–54). While the exact role of BDNF in myopia has not been clearly identified to date, lower levels of BDNF have been found in the aqueous humor of myopic individuals (55), and polymorphisms in the noncoding RNA gene BDNF-AS, the antisense RNA of BDNF, have been associated with myopia (56).

In line with a higher number of differentially expressed genes, the gene set was enriched for more pathways (Figures 5B, C, Supplementary Table 9). Similarly to the mouse model, the TGF-beta pathway and circadian entrainment pathways were prominent. Among the pathways that were overrepresented in the genes changing in chick retinas, but not in mouse retinas, were the HIF-1 signaling pathway as well as glucagon signaling.

Next, we sought to understand to what extent the genes and pathways differentially regulated in the mouse and chick retinas overlap. Furthermore, we wanted to determine how the findings from experimental myopia models compare to the knowledge about refractive error genetics in humans. To this end, we analyzed the results of the meta-analyses in the context of genes implicated in refractive errors in GWAS studies (Supplementary Table 10), defined as being mapped in the proximity of genetic loci associated with refractive error and myopia. Regarding gene-level changes, we found that one gene, follistatin (Fst), was downregulated in mouse and chick retinas, and has also been implicated in refractive error in humans (Figure 6A, Supplementary Table 11). The impact of light perception and photoreceptors also seems to be conserved in the mechanisms. In particular, we observed that the gene encoding for melanopsin (Opn4) was downregulated in both mouse and chick retinas (Figure 6A). In addition, genes associated with photoreceptor function and phototransduction, such as Vsx1 and Rdh5, were changing in chick retinas and are also implicated in refractive errors in humans (Figure 6A).

We further explored the pathways enriched for the genes changing in mouse and chick retinas, as well as those overrepresented in human refractive error susceptibility genes (Supplementary Table 12), with a subset highlighted for their relevance in retinal and refractive development (Figure 6B, Supplementary Table 13). Pathways associated with all three analyses included, again, the TGF-beta signaling pathway and circadian entrainment processes. These findings reveal the mechanistic similarities between experimental myopia in mouse and chick retinas and our knowledge of refractive error genetics in humans. They highlight the complexities of myopia as a multifactorial condition influenced by both conserved and species-specific pathways.