RNA-seq and GSEA identifies suppression of ligand-gated chloride efflux channels as the major gene pathway contributing to form deprivation myopia

Currently there is no consensus regarding the aetiology of the excessive ocular volume that characterizes high myopia. Thus, we aimed to test whether the gene pathways identified by gene set enrichment analysis of RNA-seq transcriptomics refutes the predictions of the Retinal Ion Driven Efflux (RIDE) hypothesis when applied to the induction of form-deprivation myopia (FDM) and subsequent recovery (post-occluder removal). We found that the induction of profound FDM led to significant suppression in the ligand-gated chloride ion channel transport pathway via suppression of glycine, GABAA and GABAC ionotropic receptors. Post-occluder removal for short term recovery from FDM of 6 h and 24 h, induced significant upregulation of the gene families linked to cone receptor phototransduction, mitochondrial energy, and complement pathways. These findings support a model of form deprivation myopia as a Cl− ion driven adaptive fluid response to the modulation of the visual signal cascade by form deprivation that in turn affects the resultant ionic environment of the outer and inner retinal tissues, axial and vitreal elongation as predicted by the RIDE model. Occluder removal and return to normal light conditions led to return to more normal upregulation of phototransduction, slowed growth rate, refractive recovery and apparent return towards physiological homeostasis.


Scientific Reports
| (2021) 11:5280 | https://doi.org/10.1038/s41598-021-84338-y www.nature.com/scientificreports/ clinical myopia or the source of the rapid changes in axial length in humans [21][22][23][24][25][26][27] as a result of 30 min of prolonged accommodation or water drinking. In animal models of myopia, rapid axial elongation, refractive change and altered gene expression 28 is seen following 6 h of − 10D optical defocus in chicks or within 30mins of removal of form deprivation 29 ; see Fig. 1. One such physiological model of myopia development based on very well established retinal/RPE physiology and extensive literature relating to rapid light/dark induced fluid shifts in the retina [29][30][31][32][33] has been proposed by Crewther 34 as the Retinal Ion Driven Efflux (RIDE) model of myopia. This model postulates that occlusion or acute blur perturbs phototransduction and hence slows the rate of exchange of ions and fluid between photoreceptors, sub-retinal space (SRS) and retinal pigment epithelium (RPE), and consequently attenuates the normal chloride anion driven efflux of fluid from the vitreous across the retina/RPE 30,33,[35][36][37] to choroid. This model is well supported by ultrastructural evidence 38 showing edema in the retina and extensive hyperosmolar stress in the nuclei, mitochondria and basal membranes of the Retinal Pigment Epithelium (RPE) 22,29,39 . Liang et al. 29,39 and Crewther et al. 22 have also employed Scanning Electronmicroscopy with elemental microanalysis to demonstrate increased potassium, sodium and chloride ion concentrations across the posterior eye following form deprivation. Thus we aimed to investigate if the biochemical pathways identified by RNAseq transcriptomics during induction of FDM and short term recovery would refute or add any support for the RIDE model theory.
To date, molecular research into gene models of myopia induction have been dominated by human genome wide association studies (GWAS) 8,[40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57] and limited pathway analysis with subsequent transcriptome analyses of single differentially expressed genes (DEGs) [58][59][60][61][62][63][64][65][66][67][68][69][70][71] or proteins [72][73][74][75][76][77][78][79][80][81][82][83] . The lack of consistent results for the above studies could largely arise due to the stability of non-profound refractions of humans in many of the GWASs Figure 1. Biometric measurements of normal and form deprived (FD) chicken eyes following induction and recovery from form deprivation myopia. Chickens were visually deprived by occlusion of their right eye then given a variable number of hours of normal visual experience (T = 0 h, 6 h, and 24 h). At each timepoint (n = 8), biometric measures were taken including (a) refraction, (b) axial length, (c) vitreous chamber depth (VCD) and (d) anterior chamber (AC) length. Note: VCD is measured from the posterior surface of the lens to the fibre layer of neural retina. AC length is measured from the anterior surface of the lens to the cornea. For AL measurements of specific chicks utilised in the RNA-seq analysis, refer to Supplementary Table 9. and the use of differing species, ocular tissues, durations of form deprivation or hours of optical defocus and relative state of refractive compensation. However, more recent use of comparative pathway analysis techniques, and threshold free pathway analysis of multiple gene sets (ie gene set enrichment analysis (GSEA)) in particular have begun to provide more useful insights highlighting dysregulation in several fundamental biological pathways 28 . In particular, form-deprivation (FD) myopia has now been associated with molecular mechanisms of suppression of neuron structure/growth, signal transduction and inner retinal neurotransmitters associated with Cl − ion transport 63,70 , while transcriptomic changes associated with early optical induction of myopic and hyperopic refractive errors also highlight altered metabolic pathways in response to lens induced myopia (LIM) 69 .
The primary aim of this study was to extend our previous transcriptome microarray study of the differentially affected genes and biochemical pathways during prolonged chick FDM and the early hours of recovery, by using next-generation RNA sequencing technology. This technique possesses greater theoretical sensitivity to further elucidate the molecular pathway underpinning the axial length changes in FDM and during recovery than microarray technology. We have aimed to assess expression of a priori identified gene sets and biological pathways associated with myopia using GSEA rather than individual differentially expressed genes (DEGs) (Reviewed in 28 ). The RIDE model and previous molecular work has suggested that gene families associated with cone phototransduction and outer and inner retinal solute transport and energy metabolism, should characterize pathways likely to be associated with increased ocular growth during the induction of form-deprivation myopia (FDMI) and slower growth rates during refractive recovery (FDMR) following reintroduction of normally modulated light conditions.

Material and methods
Animals. Thirty-two male chicks (Leghorn/New Hampshire) were raised with unlimited food and water in a controlled environment on a 12-h light/12-h dark cycle and with the temperature maintained at 30 ± 0.5 °C. Illuminance was maintained at 183 lux during the 12 h day cycle (from 8am to 8 pm) using a 20 W halogen lamp. Twenty-four chicks were monocularly occluded on post-hatch day 5 for 7 days to induce FD myopia. The translucent polystyrene occluders were attached to the periocular feathers of the right eye. On day 12 posthatch, occluders were removed from FD chicks at 10am and animals were either immediately sacrificed, ie 0 h recovery (n = 8), or sacrificed after 6 h (ie at 4 pm; n = 8) or 24 h (n = 8) of the normal light/dark condition (ie 10 h-light/12 h-dark condition and 2 h light to make up 24 h) following the prolonged form-deprivation. Eight separate chicks were included in the analysis as non-occluded controls at the T = 0 h timepoint.
We chose to use a separate batch of control animals rather than use the within subject contralateral fellow eyes. This decision was based on prior evidence showing binocular changes in choroidal blood flow during monocular occlusion 84,85 and similar direction significant changes in refraction and axial length in the fellow eyes 14 . Binocular interaction effects have been suggested to account for the substantial differences in gene expression as observed by He et al. 86 between contralateral eyes and separate controls (see also 75,87,88 , hence we have chosen to use separate control animals for these reasons. All animal work in this study was approved by the La Trobe University Animal Ethics Committee (Approval number AEC 11/68) and is in accordance with the ARVO Guidelines for Use of Animals in Research, Australian NHMRC Animal Ethics requirements and the ARRIVE guidelines 89 .

Ocular biometrics analysis.
Refractive state (dioptres (D)), vitreous chamber depth (VCD in mm), axial length (AL in mm) and anterior chamber (AC in mm) measures were collected from all animals on day 12 post-hatch and after 0 h recovery (n = 8), 6 h recovery (n = 8) and 24 h recovery (n = 8) while animals were lightly anesthetized with an intramuscular injection of a mixture of ketamine (45 mg/kg) and xylazine (4.5 mg/ kg). Refraction in the experimental right eyes were determined by trained ophthalmic practitioners using retinoscopy (Keeler, Vista Diagnostic Instruments) and A-Scan ultrasonography (A-Scan III, TSL; Teknar, Inc. St Louis, USA; 7 MHz probe) was used to measure axial dimensions. Baseline biometric measures were not sought to avoid repeated potential anaesthesia effects on eye growth 67 , RNA quality and sequencing [90][91][92][93] . Chicks were only sedated prior to decapitation. Quantitative data were expressed as Means (+ /− Standard Error). Analysis of Variance (ANOVA) measuring group differences was carried out between same age controls, form-deprived and recovery eyes to determine significant changes in biometric measurements followed by post-hoc tests as required.
RNA isolation and cDNA library preparation. Tissue samples for RNA isolation and sequencing were collected from right eyes only, from four out of eight chicks per time condition. These four chicks were chosen based on comparable axial length measurements (see Supplementary Table 13). Total RNA was isolated from posterior eye retina/RPE/choroid tissue using the Trizol method and the RNeasy mini kit (Qiagen) with the oncolumn DNA digest according to manufacturer's instructions. RNA quality was measured using the Agilent 2100 Bioanalyzer (RNA 6000 Nano Kit; Agilent Technologies, Santa Clara, CA, USA). All samples had a RNA integrity number (RIN) above 8.7. RNA was also quantitated on the Qubit 2.0 Fluorometer (RNA HS kit; Invitrogen, Australia). For library preparation, RNA concentration was calculated using an average measure from the Bioanalyzer and Qubit assays. A total of 2.5 µg of mRNA was purified from total RNA using oligo (dT)-conjugated magnetic beads (Illumina, San Diego, CA, USA). The fragmented mRNA was then subjected to cDNA synthesis using the TruSeq Stranded mRNA kit (Illumina, San Diego, CA, USA) following the manufacturer's low-sample throughput protocol. All cDNA libraries were assessed on an Agilent 2100 Bioanalyzer (DNA 1000 kit; Agilent Technologies, Santa Clara, CA, USA). Size of the final products were approximately 260 bp. cDNA libraries were quantitated using 3 methods: Bioanalyzer (DNA 1000 kit; Agilent Technologies, Santa Clara, CA, USA), Qubit (dsDNA HS assay; Invitrogen) and qPCR (GeneRead Library Quant Array; Qiagen, Germantown, MD, USA). DNA libraries were normalised to 10 nM with Tris-HCl 10 mM, pH 8 104 . Chick Ensemble Gene IDs were mapped to human gene symbol using BioMart (Ensembl release 89; Supplementary Table 11). Where mapping resulted in no known gene, Chick Ensemble Gene IDs were mapped to human gene orthologues (Homo sapiens GRCh38.p7). The conversion of Chick IDs to human gene symbols was done as the MsigDb and associated genesets that are based on human gene annotations. The GSEA technique involves ranking all genes within a sample dataset based on their differential expression between two experimental groups using a basic metric e.g. signal-to-noise ratio (S2N), ratio of average expression from two classes (Ratio), T-test statistic (T-test), or the Pearson correlation coefficient for quantitative studies 101,102,105 . The method then evaluates the general differences in the cumulative distribution in expression of genes in a biological pathway based on a priori knowledge of the gene's biological function (i.e., gene sets from the Molecular Signatures Database (MSigDb)). An enrichment score (ES) for each gene set is calculated. This ES reflects the degree to which a geneset is overrepresented at the top or bottom of a ranked list of genes. The ES is calculated by walking down the ranked list of genes, increasing a running-sum statistic when a gene is in the gene set and decreasing it when it is not. The ES for each pathway reflects the maximum deviation from zero encountered in walking down the list. A normalised enrichment score (NES) is also calculated by GSEA in which differences in pathway size (i.e. geneset size) are considered, allowing for comparisons between pathways within the analysis 101,102 . In this study, the default Signal2Noise metric was used to determine significantly enriched pathways following prolonged occlusion (0 h vs Control) and during the 24 h recovery period (6 h vs 0 h, 24 h vs 6 h, and 24 h vs 0 h). This metric uses the difference of means scaled by the standard deviation 106 . For comparisons purposes, the Pearson's correlation metric was used as recommended for time-series data 106 , to assess changes in gene expression associated with the entire recovery period post occluder removal (i.e., 0 h vs 6 h vs 24 h). The analysis involved 1000 gene set permutations with gene sets limited to 15-500. As recommended by The Broad Institute 106 for exploratory studies, a FDR of 25% was used for all bioinformatic analyses.
To understand which genes contributed to the gene set's enrichment signal, we performed leading-edge subset (LES) analysis. This analysis identifies the core genes by creating a ranked list. The genes located at the top of the ranked list are considered to be upregulated and those at the bottom of the list are downregulated 101 .

Gene validation.
As the primary focus of this study was to uncover the molecular pathway underpinnings of FD myopia, validation of the large number (> 300) of core genes identified from this study using qPCR was not feasible to assess. Instead we aimed to validate the large number of core genes using our previously published microarray study of short term refractive recovery in chick eyes following 10 days of form deprivation 70 . There is also widespread agreement in the literature that highlights good concordance in gene expression between microarray and qPCR methods [107][108][109][110][111][112] . Briefly, in our previous microarray study 70

Results
Ocular biometrics analysis. Refractive state (dioptres (D)) was assessed using standard retinoscopy while vitreous chamber depth (VCD in mm), axial length (AL in mm) and anterior chamber depth (AC in mm) were assessed using A-scan ultrasonography (Fig. 1). Significant myopia (− 18.6D ± 1.44) was achieved after 7 days of occlusion. Six hours of visual experience induced a significant hypermetropic shift of + 7.3D, reducing the mean level of myopia in FD eyes down to -11.3D ± 0.91. Less refractive compensation was seen after 24 h of visual recovery post-occluder removal (− 16.4 ± 1.18D), presumably due to the circadian impact of the regular night period with this measurement being made 2 h after the first 12 h night/dark period (Fig. 1a). Axial length (AL), vitreous chamber depth (VCD) and anterior chamber (AC) showed similar patterns of biometric growth consistent with previously published diurnal rhythm changes [113][114][115] . Prolonged occlusion increased axial length from 9.32 mm for controls to 10.30 mm (Fig. 1b), VCD from 5.65 mm to 6.21 mm (Fig. 1c) and AC (Fig. 1d) from 1.53 to 1.83 mm. In recovering eyes, AL decreased by 0.12 mm at both 6 h and 24 h, however these changes were not significantly different to AL at occluder removal (T = 0 h; p > 0.05). The same trend was also seen in VCD and AC at 6 h with VCD reduced by 0.08 mm and AC by 0.04 mm indicating that the excessive ocular growth response to FD is inhibited by occluder removal and that normal light conditions favour re-emmetropisation.
Differentially expressed genes in prolonged FD and FD recovery. To identify the transcriptomic mechanisms involved in the response to, and recovery from FD, single-gene expression changes were analysed. We first assessed DEG after prolonged occlusion (i.e. 7d of FD with 0 h recovery) and during recovery (6 h vs 0 h, and 6 h vs 24 h) using EdgeR with a moderated t-test and a false discovery rate (FDR) of 5%.
Prolonged occlusion, relative to no-lens controls, induced expression changes in 13 genes with an FDR < 5%. BMP2, ALDH1A2, TNC, SHC4, GSN, SIK1 were down-regulated at 7d of FD relative to controls, while LOC417800, HCK, WNT9A, KIAA1199, CLEC3B, SLCO1C1 and RHOB were upregulated in the FD eye compared to controls (Supplementary Table 1). By comparison, recovery from FD identified 439 transcripts at 6 h compared to 0 h, 386 transcripts between 6 and 24 h, and 216 transcripts at 24 h compared to 0 h (Supplementary  Tables 2-4). To provide a comparison with past transcriptomic studies on chick FD myopia and early recovery 63,70 we analysed all time-points (0 h, 6 h and 24 h) using a moderated F-test which identified a total of 828 transcripts significantly differentially expressed across the 24 h period of normal day/night conditions following 7 days of form deprivation (Supplementary Table 5). The number of overlapping and unique transcripts identified is illustrated in Fig. 2 with the majority of these gene transcripts involved in cellular and metabolic processes (Fig. 3). and recovery from FD, GSEA was employed to identify changes in expression of biological pathways associated with 7 days of FD relative to normal, no lens controls. Prolonged occlusion resulted in suppression of only one pathway, ligand-gated ion channel transport (Table 1; Fig. 4). The leading-edge subset (i.e., core genes) of this pathway implicate suppression of several GABA ionotropic receptors (ie GABA A and GABA A -rho receptors (previously known as the GABA C receptors)) and glycine receptors. Such receptors are associated with Cl − channels 116,117 and would be expected to play a role in transretinal fluid movement from vitreous to choroid 30,31,34 .
The core genes identified in the ligand gated pathway show good concordance with gene expression from our previous microarray study using the same experimental paradigm (Supplementary Table 12).   (Fig. 6). As expected, the data in Fig. 4 suggests that cone receptor phototransduction is still suppressed after only 6 h normal visual experience, i.e. cone phototransduction pathway has not totally recovered its control levels during the first 6 h of visual experience and recovery after occluder removal. An earlier electrophysiological study 38 also found reduced cone sensitivity immediately following prolonged FD which would support such gene pathway suppression in early (6 h) and prolonged recovery (24 h). Indeed, our results comparing gene expression changes between each recovery time-point of 6 h and 24 h also indicate that the pathways underlying cone phototransduction were significantly upregulated when evaluated using the Signal2Noise metric as demonstrated in Figs. 4 and 5 and in Table 2. Interestingly, while transcriptomic changes in mitochondrial metabolism are absent during the first 6 h of recovery, increased mitochondrial metabolism is observed over the next 18 h (Fig. 5; Table 2) after a 1 day/night cycle has been com-  www.nature.com/scientificreports/ pleted. These changes occur in parallel to the increase in expression of the high-energy consuming Na + -Ca 2+ -K + exchanger (NCKX2, also known as SLC24A2), which is characteristic of the dark-adapted retina 118 (Fig. 6).
To determine if the pathways identified during FD recovery using the Signal2Noise metric produced similar pathway findings to the Pearson's correlation metric, the data was then analysed (using GSEA) for pathways correlated to recovery time (0 h, 6 h, 24 h). Twenty-nine pathways were found to be significantly altered during the 24 h recovery period with the top 3 pathways involved in the complement cascade, mitochondrial energy metabolism and retinol metabolism (Supplementary Table 7). Other pathways identified by Pearson's correlation have roles in immunity, lipid metabolism, smooth muscle contraction, apoptosis, protein metabolism, signal transduction, golgi transport, extracellular matrix (ECM), hemostasis and disease. Indeed, the findings from both the Signal2Noise analysis and Pearson's correlation suggests that these metrics can produce comparable results.

Discussion
This study has utilized pathway enrichment analysis (GSEA) to identify biological pathways shown by RNA seq to be involved in profound myopia induced by 7 days of FD and during the first 24 h of early recovery in a normal modulated light environment. Our findings raise intriguing questions regarding the temporal interaction of transcriptome changes and previously described ionic, physiological and morphological changes associated with induction and recovery from FDM 22,29 and theoretical models of myopia.

GABA and glycine mediated chloride ion transport pathways significantly altered in profound FDM.
Consistent with the RIDE model, GSEA identified suppression of the GABA-and glycine-mediated chloride ion transport pathway as the only significantly altered biochemical pathway associated with 7 days of excessive ocular growth, thinning of the retina and choroid and altered ion distribution patterns accompanying profound refractive myopia 29 . The fact that the only pathway showing significant change after 7 days of FDMI was the chloride ion transport pathway is particularly important to understanding development of overly large volume eyes, which is the hallmark sign of myopia. As extensively studied since the early 70s 119 and Reviewed in 120 , chloride channels reside in all plasma membranes and most intracellular organelles and are fundamental to all cell volume regulation, transepithelial transport and regulation of electrical excitability in general. Hence suppression of the Cl-channel pathways across the entire retina/RPE/choroid examined with RNAseq provides an explanation for the ultrastructurally described dehydration seen in FDM tissue in profoundly myopic eyes 22,29,39 . FDM recovery pathway changes -cone phototransduction first. By comparison with the stable conditions during prolonged occlusion, GSEA of RNAseq data showed rapid upregulation of many cellular processes but particularly upregulation of cone phototransduction and BARD1 expression following occluder Table 2. Enriched Pathways Significantly Altered during Recovery from Form-deprivation (FDMR). GSEA was used to identify pathways (KEGG, Reactome, STKE and PID) significantly altered between recovery timepoints 6 h and 24 h, compared to 0 h using the Signal2Noise metric. Gene expression changes at 6 h compared to 0 h were associated with suppression in cone receptor phototransduction. The period between 6 and 24 h post occluder removal was mainly associated with changes in mitochondrial metabolism and other metabolic pathways. Pathways significantly altered after 24 h recovery compared to 0 h involves a range of cellular processes. ES, enrichment score; NES, normalised enrichment score; FDR, false discovery rate. Note: pairwise comparisons with control animals (ie 6 h vs control and 12 h vs control) can be found in Supplementary  Table 6. NOTE: + NES = expression upregulated;-NES = expression downregulated; ≠ present in GSEA using Pearson's correlation metric (Supplementary Table 7); *present in previous data 70 . Upregulation of BARD1 at 6 h recovery. Upregulation of the BARD1 pathway in the first 6 h of normal light conditions highlights that prolonged occlusion, profound myopia and the morphologically thinned choroid is persistently physiologically stressful to cellular function 121,122 . Over the last few years the BARD1 pathway has been reported to play a central role in the control of the cell cycle in response to DNA damage in many diseases and its regulation has been shown to mediate the formation of 'Lys-6′-linked polyubiquitin chains and coordinate a diverse range of cellular pathways such as DNA damage repair, histone ubiquitination and transcriptional regulation to maintain genomic stability [121][122][123] . The presence of BARD1 post profound FDM may be an early indicator of the types of cellular damage observed in human and animal myopia, especially with regard to thinning of the choroid and retina. This may explain why high myopia is a risk factor for secondary ophthalmic diseases later in life.

Implications of suppressed ligand-gated chloride channel pathways for ocular growth.
Significant suppression of the ligand gated chloride channel pathway following prolonged FDM is compatible with the notion that increases in axial length and thinning of the retina and choroid in response to FD occur concurrently with changes in the distribution of physiologically important ions including chloride, sodium, potassium and calcium across the posterior eye 22 . In a review of Ion Channels of the RPE, Wimmers et al. 33 similarly to Gallemore et al. 30 , and Crewther 34 , have highlighted the role of Cl − and K + ions as the key drivers of transepithe- The established role of GABA and glycine in retinal third order neuronal transmission and in Cl − transport 116,117 would also suggest that hydrated Cl − ions could play a particularly important role in axial growth during FD in animals 69,70,83 and human myopia 8 . Indeed suppression of the GABA and glycine ligand-gated chloride iontransport pathway in response to FD (Table 1) is also consistent with previous reports of significantly reduced retinal concentrations of GABA in chick FDM 16,124 and GABA signal following abnormal axial growth in FD guinea pigs 71 . A number of pharmacological studies have also utilized GABA antagonists (e.g. TPMPA) to inhibit the response to FD by inhibiting axial elongation and vitreous chamber depth 16,[125][126][127][128] . Downregulation of GABA and glycine in retinal tissue has also been seen in normally growing chicks over a 48 h period 129 , whereas the same study observed an increase in GABA signaling proteins in the same 6-48 h period when chicks wore negative lens wear and during which refractive compensation was achieved 129 .
Ionotrophic GABA receptor distribution in the posterior eye. Interestingly GABA receptors, that play such an important role in major third order neuronal transmission 130 , are also found in abundance on RPE cells 131,132 suggesting that suppression of such ligand-gated chloride channels in RPE and most retinal cells 133 during FD would be expected to reduce the transretinal fluid efflux towards choroid and result in rapid increase of fluid in the vitreous and axial elongation. Furthermore, as the majority of GABA receptors identified in Table 1 are ionotropic, the role of ion homeostasis, particularly K + and Cl − , in the development of axial elongation and myopia cannot be ignored. This is particularly so given the substantial evidence that fluid flow across the RPE is related to the ionic environment of the retina/RPE/choroid 30 . Spatial distribution of ions using X-ray microanalysis has revealed significant differences in the distribution of ions of phototransduction (Na, K, Cl, Mg and Ca) as well as other physiologically relevant ions (PO 4 3− and SO 4 2− ) known to be important for cellular functioning and structural integrity 22 . Taken together, these findings offer strong support for the role of light, cone phototransduction, RPE mechanisms, ion regulation and neurotransmission in driving myopia development and in no way refute the predictions of the RIDE model 34 .
Metabolic recovery following reintroduction of temporally modulated light. Recovery from refractive myopia and reduced axial, vitreal and anterior chamber growth rates seen during the first 24 h following FD removal appear to be associated with both cone-dominated phototransduction and increased mitochondrial metabolism, especially at 6 h and 24 h post FD. These findings are consistent with our earlier FD microarray study 70 and RNAseq study of refractively compensated negative lens-induced myopia in chick 69,134 . Interestingly, mitochondrial respiratory electron transport chain genes have also been identified in a human study of genetic myopia 135 , highlighting the important contribution of mitochondrial bioenergetics and the phototransduction cascade to ocular function and myopia development. This data is also in line with previous ultrastructural evidence of abnormal photoreceptor elongation, loss of mitochondrial integrity 29,39,136 , and oxidative stress 70 during FD and subsequent reversal of most morphological changes in the 48 h following occluder removal in chick 29,137,138 . Indeed, our earlier microarray study 70 indicated that mitochondrial respiratory complex 1 and 3 genes and mitochondrial reactive oxygen species (mROS) 139 are upregulated during FDMI. The release of mROS has previously been suggested to be a response to cellular stress while coincidently acting as a signalling molecule to facilitate cellular adaptation to this stress 140 . In fact, mROS may upregulate cone photoreceptor pathways and hence neurotransmission in inner retina 141 , chloride transport and homeostasis 142 which are the main processes we have reported previously following termination of occlusion and re-established normal light conditions 22,29 . Thus, the increase in mitochondrial electron transport chain mRNA seen 24 h after occluder removal ( Fig. 5; Table 2) fits with the idea that the dark-adapted (low temporal luminance modulation) retinae require ~ 20% more metabolic activity than the same light-adapted retina 143 . This is not unexpected given that the FD retina has previously been shown to be in a pseudo dark-adapted state 37 and associated with the reaccumulation of K + in the subretinal space and exclusion of fluid under conditions of low temporal luminance modulation 22 . Importantly, the ATPase mechanism of the RPE apical surface is electrogenic and modulates the transepithelial potential which is closely related to the control of fluid flow across the RPE and so would be expected to be upregulated during form deprivation, lower in the first 6 h of normal light and then upregulated again during the ensuing night 34,144 . Immune pathways in recovery. The presence of immune related pathways in the recovery from the physiologically stressful FD condition and the upregulation of the BARD1 pathway suggest a role for the immune system in refractive compensation. Such an association has recently appeared in regards to the relation of clinical blood counts as measures of immune responses and inflammation and high myopia [145][146][147] . The identification of the complement pathway during recovery from FD is consistent with our earlier transcriptomic studies in chick 70,134 and support previous well-described links between mitochondrial bioenergetics and the complement system in the body and in the eye [148][149][150] and in many neurodegenerative diseases such as Alzheimer's disease (AD), Parkinson's disease (PD), and Amyotrophic Lateral Sclerosis (ALS) 151 .
Complement factors have also been associated with previous reports of human and animal refractive errors 70,134,152,153 and should be expected in animal models of FD given the severe choroidal thinning and associated ultrastructural and ionic changes seen in FD eyes post occlusion 22,29 . These changes have also been demonstrated in choriocapillaris fenestration number, choroidal blood vessel permeability in chicks recovering from FDM 154 and seen in oncosis 155,156 . The previously described increase in extravascular space edema 22,29 and lymphatic vessel permeability 154 immediately post-occlusion also supports our results here as well as the recent FD transcriptomics 70 . What was not expected is the involvement of other immune pathways besides complement in FD recovery. www.nature.com/scientificreports/ Limitations. Lastly, while our decision to use retina/RPE/choroid may be considered a limitation, we contend that investigation of the combined changes in expression of genes within the entire biological network of the eye is necessary to understand the photoreceptor induced temporal changes in form deprivation myopia seen in animal models. Indeed, when evaluated experimentally in combination with systematic system neuroscience review techniques, the impact of using inner and outer retina, RPE and choroid in large scale genomic and proteomic studies supports the idea that similar biological mechanisms are associated with FD and lens induced defocus, regardless of the varying combinations of tissues used 28 . Indeed, Fig. 5 illustrates the commonality of the biological processes that underlie the adaptive responses to environmental manipulation of temporal modulation of luminance information across the multiple tissue types of the posterior eye. Thus, we argue that the use of multiple tissues is not an impediment to our GSEA-based interpretation as evidenced by the robustness of our current findings with previous research in human and animals. Furthermore, GSEA has proven a useful tool in identifying significantly enriched biological pathways in large-scale 'omic' studies. However we also acknowledge that large-scale discovery-driven studies may produce false positives. Additionally the limited gene sets and the redundancy in the databases currently available in the Molecular Signature database may affect future interpretations of this data. Future studies may benefit from using a wider range of databases as well as re-analysing data periodically to account for updated gene and pathway information. Additionally, we have not considered the role of circadian rhythms of ocular refraction and gene expression, particularly with regard to the differences in refractive normalization seen at 6 h compared to 24 h recovery from FD. We do not have a circadian matched control for our 6 h recovery group and acknowledge that this may be a limitation in our study. Indeed, 2 circadian-related pathways (RORA Activates Circadian Expression and Circadian Repression of Expression by REV-ERBa) were highly expressed at 6 h recovery compared to the no lens control however this result should be taken with the understanding that the no lens control is not a circadian matched control for the 6 h recovery timepoint. As such, we cannot conclude if circadian rhythms have an influence on refractive development and myopia. However, our no lens controls, 7 days induction, and 24 h recovery groups all experienced the same day/night patterns with tissue collected at approximately 10am. No circadian effects were highlighted suggesting that circadian rhythms effects on gene expression will be minimal and as such, may not be a factor in refractive development and myopia. Further studies should incorporate circadian matched controls to minimise circadian rhythm effects on gene expression.

Conclusions
We conclude that axial myopia is an adaptive response to the environmental perturbation of the visual signal cascade (i.e. Phototransduction). This perturbation in phototransduction consequently affects the ionic environment of the photoreceptors, subretinal space and inner retina as predicted by the RIDE model via GABA and glycine pathway signalling and consequent effects on ligand-gated chloride channels and their role in transretinal fluid flow. Our identification of the different mitochondrial bioenergetic profiles observed between FDMI and FDMR could be considered as a potential molecular hallmark of the myopia condition. Most importantly our findings are also consistent with the predictions of the theoretical implications of the RIDE hypothesis highlighting the importance of light driven osmoadaptive pathways, enhanced mitochondrial bioenergetics and cellular immune responses, during the development of FD myopia and during refractive recovery.

Data availability
The datasets generated during the current study are available in the Gene Expression Omnibus (GEO) repository (https ://www.ncbi.nlm.nih.gov/geo/; Accession # GSE80327).