Abstract
Myopia (short-sightedness) affects 1.45 billion people worldwide, many of whom will develop sight-threatening secondary disorders. Myopic eyes are characterized by excessive size while hyperopic (long-sighted) eyes are typically small. The biological and genetic mechanisms underpinning the retina's local control of these growth patterns remain unclear. In the present study, we used RNA sequencing to examine gene expression in the retina/RPE/choroid across 3 days of optically-induced myopia and hyperopia induction in chick. Data were analyzed for differential expression of single genes, and Gene Set Enrichment Analysis (GSEA) was used to identify gene sets correlated with ocular axial length and refraction across lens groups. Like previous studies, we found few single genes that were differentially-expressed in a sign-of-defocus dependent manner (only BMP2 at 1 day). Using GSEA, however, we are the first to show that more subtle shifts in structural, metabolic, and immune pathway expression are correlated with the eye size and refractive changes induced by lens defocus. Our findings link gene expression with the morphological characteristics of refractive error, and suggest that physiological stress arising from metabolic and inflammatory pathway activation could increase the vulnerability of myopic eyes to secondary pathologies.
Introduction
Myopia (short-sightedness) is the most common ocular disorder with rapidly increasing prevalence and severity worldwide (Seet et al., ; Pan et al., ). The increased axial length that characterizes myopic eyes results primarily from expansion of the fluid-filled vitreous chamber (Berman and Michaelson, ), and is accompanied by a decrease in the thickness of the vascular choroid (Hayes et al., ; Wallman et al., 1995; Westbrook et al., 1995) and thinning of fibrous scleral tissue that forms the eye's outer coating (Curtin and Teng, ). These structural changes greatly increase the risk of developing sight-threatening visual pathologies (Saw et al., ; Verhoeven et al., 2015), making the development of treatments to limit excessive ocular growth an important health and socioeconomic priority.
Ocular growth is controlled locally by the retina (Troilo et al., 1987; Wildsoet and Pettigrew, 1988) in a process that involves a complex interaction between an individual's genetic constitution and their environment (Wallman and Winawer, 2004). The role of the visual environment has been investigated using animal models in which rearing with negatively powered defocusing lenses or form deprivation occlusion increases the rate of growth (resulting in a larger myopic eye) while rearing with positively powered defocusing lenses slows growth (resulting in a smaller hyperopic eye; Wallman and Winawer, 2004). Microarray studies examining transcriptome changes in these models have implicated thousands of genes, however relatively few findings have been replicated across studies and it remains unclear which genes are important for controlling ocular growth (Ashby and Feldkaemper, ; Stone and Khurana, ; Stone et al., ).
Many of the transcriptome studies conducted to date have investigated a single ocular growth phenotype (e.g., myopia or hyperopia induction relative to controls; Tkatchenko et al., 2006; Brand et al., ; Mcglinn et al., ; Schippert et al., ). This experimental design makes it challenging to separate the genes involved in growth processes from those responding to the secondary effects of lens-wear (i.e., blur and physiological stress). Experimental designs that concurrently compare myopia and hyperopia induction enable identification of genes with expression profiles that are discriminatory for different ocular growth trajectories. Such genes are more likely to be directly involved in growth processes (thus providing attractive targets for therapy), and may also help to explain why myopic eyes are at a higher risk for secondary pathologies than their hyperopic counterparts. Although several studies have compared transcriptome-wide responses during ocular growth increases and decreases (Shelton et al., ; Ashby and Feldkaemper, ; Stone et al., ), few genes showing sign-of-defocus dependent expression have been identified and it has been suggested that distinctive (rather than bidirectional) genetic responses underlie myopia and hyperopia (Ashby and Feldkaemper, ; Stone et al., ).
Opportunities now exist to improve on the methods used by these past studies. Expression changes in animal models of refractive error have yet to be investigated using next generation RNA sequencing which is more quantitative, sensitive, and reproducible than the microarray technologies used previously (Wang et al., 2009). Additionally, most studies to date have analyzed data at the single gene level or used first generation pathway analyses to associate biological functions with lists of differentially-expressed genes (García-Campos et al., ). Second generation pathway analyses, such as gene set enrichment analysis (GSEA), are able to identify responses that are subtle at the single gene level because they do not require an arbitrary cut-off for differentially-expressed genes (Subramanian et al., ). Although the results of GSEA are more replicable and interpretable than single gene measures (Subramanian et al., ; Manoli et al., ), it has only been employed by one ocular growth study to date (Tkatchenko et al., 2006). Improvements can also be made at the experimental design level where very little is known about longitudinal changes because few studies have analyzed multiple treatment time-points (Brand et al., ; Mcglinn et al., ; Summers Rada and Wiechmann, ; Stone et al., ). Moreover, only one study has directly correlated expression shifts with the ocular growth phenotype of individual animals (Tkatchenko et al., 2006).
Thus, in the present study we used GSEA and a longitudinal design to identify gene expression patterns related to eye size and refraction across three conditions in chick: normal development, optically induced myopia, and optically induced hyperopia. As theories of visually regulated ocular growth hypothesize biological signals that propagate from the retina through to the choroid (e.g., Wallman et al., 1995; Rymer and Wildsoet, ; Crewther et al., ; Feldkaemper and Schaeffel, ), retina/RPE/choroid samples were profiled using RNA sequencing following 1, 2, and 3 days of lens-wear, or no lens rearing. We then used GSEA to identify gene sets correlated with ocular axial length and refraction across lens groups at each time-point.
Materials and methods
Animals and rearing
One hundred male chicks (Leghorn/New Hampshire), obtained from a commercial hatchery, were raised from post-hatch days 0–4 under a 12-h day/night light cycle (beginning at 8 a.m.) in groups of < 25. In the middle of the light cycle on day 5, chicks were randomly assigned to a lens condition (+10 or −10 diopters, or No Lens), and lenses (Polymethyl Methacrylate) attached to Velcro were fixed to the periocular feathers of the right eye. Following a further 1, 2, or 3 days with lenses attached, 10–12 chicks per lens group were anaesthetized (ketamine, 45 mg/kg; xylazine, 4.5 mg/kg i.m.) and right eye refraction and axial dimensions determined by retinoscopy (Keeler, Vista Diagnostic Instruments) and A-Scan ultrasonography (A-Scan III, TSL; Teknar, Inc. St Louis, USA; 7 MHz probe). Chicks were euthanized and their right eyes were enucleated. The retina/RPE/choroid was immediately collected from the posterior eyecup and frozen in liquid nitrogen before being transferred to −80°C. Note that tissue was collected from right eyes only to avoid the confounding influence of right/left eye developmental asymmetries in chick (Rogers and Bolden, ). Similarly, separate control animals were used because monocular treatments in chick can affect blood flow (Shih et al., ; Jin and Stjernschantz, ), refraction and axial length (Wildsoet and Wallman, 1995) in the contralateral eye. All procedures were conducted in accordance with the protocols approved by the La Trobe University Animal Ethics Committee and adhere to the ARVO Statement for the use of Animals in Ophthalmic and Vision Research.
RNA isolation and library construction
Four chicks per lens*time condition were chosen for RNA extraction based on strong and cohesive biometric responses. The selected samples had been collected between 1 and 3 p.m. on days 1–3, and were counterbalanced for the order of collection across lens-groups. Total RNA was isolated from the retina/RPE/choroid using the miRNeasy Mini Kit (Qiagen, Germantown, MD, USA) including DNase digestion. RNA quality and quantity was assessed on the 2100 Bioanalyzer (RNA 6000 Nano Kit; Agilent Technologies, Santa Clara, CA, USA). All samples had an RNA integrity number (RIN) of >8.3. RNA quantity was also assessed on the Qubit 2.0 Fluorometer (RNA-HS assay; Life Technologies, Carlsbad, CA, USA).
Using an average of concentration measures obtained from Qubit and Bioanalyzer assays, 2.5 μg of RNA from each sample was used for library preparation and RNA sequencing. Libraries were prepared using the TruSeq Stranded mRNA LS kit (Illumina, San Diego, CA, USA) with dual indexing according to the manufacturer's instructions. The generated libraries were assessed on the 2100 Bioanalyzer (DNA 1000 kit; Agilent Technologies, Santa Clara, CA, USA) to ensure an average size distribution of approximately 280 bps, then quantified on the Qubit 2.0 Fluorometer (dsDNA HS assay; Life Technologies, Carlsbad, CA, USA) and by qPCR (GeneRead Library Quant Array; Qiagen, Germantown, MD, USA). Libraries were normalized to 10 nM in Tris-HCl (10 nM, pH8.5 with 0.1% Tween 20), pooled, and prepared for cluster generation on the Illumina cBot using the TruSeq SR Cluster Kit V3-cBot (Illumina, San Diego, CA, USA) with denatured template DNA diluted to 7 pM. The flow cell and sequencing reagents (TruSeq SBS Kit V3; Illumina, San Diego, CA, USA) were loaded on the Illumina HiSeq 1500 and a dual-index, single-end, 100 bp sequencing run performed. RNA-Seq data for each sample are available at the NCBI Gene Expression Omnibus under accession number GSE78042 (www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE78042).
Biometric data analysis
Between group statistical comparisons of refraction and axial length were made using Analyses of Variance (ANOVA) with relevant post-hoc tests as required. Both dependent variables were normally distributed (Shapiro Wilks Test p > 0.05 and/or skewness and kurtosis within acceptable range), however refraction data violated the assumption of equal variances (Levene's Test p < 0.05). A conservative significance threshold (α = 0.001) was used for tests of group differences in refractive state to combat any resulting Type I error inflation (Harwell et al., ).
Sequencing data analysis
Pre-processing
Three samples were sequenced in the 1 and 2 day No Lens groups, and 4 samples were sequenced in all other conditions. Read quality was assessed using FastQC (www.bioinformatics.bbsrc.ac.uk/projects/fastqc/), and adapter and low quality (Q-score < 10) sequences removed using CutAdapt (Martin, ) and Trimmomatic (Bolger et al., ). Reads were mapped to the chick genome (GalGal4) using Tophat2 (Kim et al., ) and Bowtie2 (Langmead and Salzberg, ). The number of reads uniquely mapping to each gene was counted using existing gene models with HTSeq (Anders et al., ). Supplementary Table S1 lists the total number of gene counts for each sample, Supplementary Figure S1 shows the variance for each sample group, and Supplementary Figure S2 shows PCA plots for samples at each time-point.
Differential gene expression
Our initial single gene analysis using EdgeR (Robinson et al., ) identified a very large number of differentially-expressed genes following 3 days of positive lens-wear (relative to other conditions). At the suggestion of a reviewer, we subsequently reanalysed our data using the more conservative DESeq2 approach (Love et al., ). We chose to incorporate the latter DESeq2 results in the manuscript as our own comparisons indicated that this method was better able to discriminate genes with previously demonstrated relevance for the treatment factors of interest (myopia and hyperopia induction). Moreover, previous research suggests that DESeq2 is better able to control the Type I error rate (Soneson and Delorenzi, ). The results of the EdgeR analysis are provided in Supplementary Tables S2, S3.
Using DESeq2, we first assessed differential gene expression for each lens condition (myopia or hyperopia induction) relative to age matched no lens controls. Gene counts for all lens conditions were loaded into R (R Core Team, ) separately for each time-point (1, 2, and 3 days). Differentially expressed genes (DEG) were identified using the SARTools (Varet et al., 2016) DESeq2 (Love et al., ) pipeline with default settings and a Benjamini-Hochberg adjusted p-value cut-off of 0.05. DEG in each lens condition were then tested for over-representation of Gene Ontology (GO) level 3–5 terms in ConsensusPathDB (Kamburov et al., ) using all genes measured as the background (FDR q < 0.05). To facilitate interpretation of GO results, the Cytoscape Enrichment Map app (Merico et al., ) was used to cluster over-represented ontologies containing similar genes. Ontologies that clustered together with an overlap co-efficient of >0.55 were collapsed into a single annotation and visualized as chord diagrams using the GOplot package in R (Walter et al., 2015). Results from the original ConsensusPathDB analysis (i.e., before redundant annotations were combined) are provided in Supplementary Tables S4, S5.
Differential gene expression was also assessed across time within each lens group. Gene counts for all time-points were loaded into R (R Core Team, ) separately for each lens condition. Genes differentially expressed in each lens group (negative, positive, and no lens) between 1–2 and 2–3 days were identified. All other processes were conducted as described above for cross-lens comparisons. Hypergeometric tests (p < 0.05) were used to determine whether the overlap in gene findings within and across groups was more than expected as a result of chance (using the total number of genes measured as the reference).
Validation of single-gene findings using previously published data
Rather than performing qPCR validation for a small number of DEG (Hughes, ), we chose to use GSEA to validate our single gene results against a previously published microarray dataset of similar design. This approach was based on the reasoning that the genes classified as significantly up- and down-regulated in our study should be enriched at the top and bottom, respectively, of a microarray dataset testing the same treatment effect (Suárez-Fariñas et al., ).
The most similar available microarray data series was obtained from the GEO Database (GSE24641). In this microarray study, Stone et al. () assessed the effects of 6 h and 3 days of negative and positive lens-wear relative to contralateral control eyes in chick. A comparison of Stone's methods with those of the present study is provided in Supplementary Table S6. The raw microarray CEL files were pre-processed using robust multi-array analysis (RMA), probe sets were median summarized, and log intensity values were imported into the GSEA program. We tested Stone's 6 h and 3 day negative lens data for enrichment of the genes up-regulated and down-regulated in the present study following 1, 2, and 3 days of negative lens-wear using the Signal2Noise metric. We expected that this approach would validate genes that respond robustly to negative lens-wear under the varied conditions encompassed by the two datasets, and help to rule out the influence of small undesirable methodological differences (e.g., contralateral eye effects, exact light intensity, space in the rearing cage, goggle material etc.) likely to be reproduced by a within-lab qPCR validation. The genes differentially-expressed following positive lens-wear were not validated using this approach because GSEA cannot accurately adjust the enrichment statistic for very small gene set sizes (Subramanian et al., ).
Gene set enrichment analysis
GSEA (Subramanian et al., ) was used to analyse the expression of Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathways from the molecular signatures database (mSigDB) (Kanehisa and Goto, ; Liberzon et al., ). The primary analyses were designed to identify KEGG gene sets correlated with ocular axial length and refraction across lens groups at each time-point. Expression values (counts per million) were imported into the GSEA program (Subramanian et al., ). GSEA was conducted with 1000 phenotype permutations using a continuous increasing phenotype label based on the axial length or refraction measure for each sample across lens groups at each time-point (1, 2, and 3 days). Pearson's metric, which uses Pearson's correlation to determine the degree of linear relationship between the gene set and expression profiles, was used for ranking genes. An FDR cut-off of 0.25 was used (as recommended by the Broad Institute GSEA User Guide; http://software.broadinstitute.org/gsea/doc/GSEAUserGuideFrame.html) and gene set sizes limited to 15–500. As with the single gene measures, further GSEA were conducted to assess KEGG pathway expression changes across time within each lens group. A continuous increasing phenotype label was used based on time in hours for each sample. All other processes were as conducted as described above for cross lens comparisons.
Leading edge subset analysis and enrichment map figure generation
Following GSEA, leading edge subset (LES) analysis was used to identify the most relevant genes within each enriched pathway (these LES genes are referred to as “core” pathway genes in the results section). This additional analysis allowed identification of enriched gene sets representing similar biological signals (i.e., gene sets with highly similar core genes; Subramanian et al., ). To cluster these redundant pathways we used the Enrichment Map App (Merico et al., ). For pathways implicated in both refraction and axial length analyses, the LES genes were combined into a single list. A gene set file was then constructed from the LES of each enriched pathway and imported into Cytoscape along with the GSEA results. For each time-point (1, 2, and 3 days), an enrichment map was built using an overlap coefficient cut-off of 0.3. In the resulting network diagrams, each node represents a pathway. Node size indicates the number of core genes in the LES for that pathway, and connections between nodes indicate common LES genes. To further interrogate the basis of pathway correlations with eye size, we also created line graphs showing the mean log2-fold change for each pathway's LES genes across lens groups.
Results
Ocular refraction and axial dimensions
Chicks wearing −10D lenses became myopic and those wearing +10D lenses became hyperopic (Figure 1A). This refractive compensation was accompanied by an increased (−10D) or decreased (+10D) rate of axial growth relative to normally developing eyes (Figure 1B). As expected, there was a strong negative correlation between axial length and ocular refraction (Figure 1C).
Figure 1
Two-way ANOVAs were conducted to compare the effects of lens-wear (+10D, No Lens, −10D) and induction time (1, 2, 3 days) on refraction and axial length. There was a significant main effect of lens-wear [F(2, 90) = 842.96, p < 0.001], but not time [F(2, 90) = 1.24, p = 0.296], on ocular refraction. A significant interaction between lens-wear and time was also observed [F(4, 90) = 24.42, p < 0.001]. Post-hoc tests revealed that the refractive state of all lens groups was significantly different by the earliest 1 day induction time-point (p < 0.001 for all comparisons). In contrast, the main effects of lens-wear [F(2, 91) = 134.52, p < 0.001] and time [F(2, 91) = 24.18, p < 0.001] on ocular axial length were both significant. As expected, an interaction effect was also observed [F(4, 91) = 4.08, p = 0.004]. Post-hoc tests revealed that chicks wearing −10D lenses had longer axial lengths than No Lens controls by 1 day (Tukey HSD p < 0.001), however the axial length difference between +10D and No Lens chicks did not reach significance until 2 days (Tukey HSD p = 0.035). Because positive lens-wear slows axial elongation relative to normal development, the axial length difference between positive and No Lens groups (and thus the absolute effect size) is limited by the rate of growth in the No Lens group. This may explain why the axial length difference between No Lens and positive lens chicks was not significant on day 1.
Genes differentially-expressed between lens groups at each time-point
To provide a basis for comparison with past studies, we first assessed differential gene expression in lens groups relative to no lens controls. In the negative lens-group, 20, 19, and 3 genes were differentially-expressed at 1, 2, and 3 days, respectively (Figures 2A,C). In the positive lens condition, 5, 9, and 2 genes were differentially-expressed at 1, 2, and 3 days, respectively (Figures 2B,D). These single gene results demonstrated good concordance with past microarray studies; many of the DEG have been implicated previously (see Supplementary Table S8 for details), and the DEG following negative lens-wear were enriched at the top and bottom of the most similar available microarray dataset as assessed using GSEA (see Supplementary Figure S3). Notably only one gene, BMP2, was differentially expressed during both myopia and hyperopia induction in a sign-of-defocus dependent manner.
Figure 2
Genes differentially-expressed within lens groups over time
We also analyzed gene expression over time within each lens group (i.e., between days 1–2 and 2–3; see Supplementary Table S9 for detailed results). No genes were differentially-expressed over time during hyperopia induction. This may be because expression measures began at 1 day when substantial refractive compensation had already occurred. By comparison, 32 genes were differentially-expressed during myopia induction (Figures 3B,D), and 57 genes were differentially-expressed during normal development (Figures 3A,C).
Figure 3
Normal developmental expression shifts were not seen in either lens group during the experimental period. There were, however, 13 commonalities between the genes differentially expressed across normal development and the genes differentially expressed when comparing lens and no lens groups at each time-point (Figure 4A). This degree of overlap in differentially-expressed gene findings with normal development was more than expected as a result of chance for both negative (p < 0.001) and positive (p < 0.001) lens groups. For most of the commonalities, lens wear appeared to accelerate the time-course of developmental expression shifts (Figure 4B). Genes showing these accelerated expression patterns in both negative and positive lens-groups were involved in circadian (NFIL3, ARNTL) and phototransduction (BLIMP-1) processes. An additional four phototransduction-related genes displayed accelerated expression patterns in the negative lens-group only (GUCA1A, CNGA3, NET1, CACNA2D4; see also Figure 3C). These findings suggest that lens wear (of both signs) perturbs the timing of developmental changes in circadian and phototransduction gene expression.
Figure 4
Pathways correlated with eye size and refraction
In addition to assessing single gene changes, we used GSEA to identify pathway expression patterns related to eye size and refraction at each time-point across the three growth conditions (normal development, myopia induction, and hyperopia induction). Fourteen pathways were implicated in both refraction and axial length analyses. An additional 12 pathways were significantly associated with either axial length or refraction (Supplementary Table S10). For comparison, we also analyzed pathway expression changes over time within each lens group (Supplementary Table S11; note that because expression data were collected across post-hatch days 6–8 only, interpretation of within group results is limited by the lack of a common “0 h” starting measure on post-hatch day 5).
Structural pathways were positively correlated with eye size and negatively correlated with refraction on day 1 (Figures 5A,B), but not at later time-points. This agrees with biometric measures (Figure 1A) where much of the growth and refractive compensation occurred within the 1st day of lens wear. Several of the identified pathways could be clustered based on common core genes and biological functions (Figure 5A). A cluster of pathways related to extracellular matrix structure was implicated in both axial length and refraction analyses. In addition, tight junction, GPI (glycosylphosphatidylinositol) anchor biosynthesis, and smooth muscle contraction pathways were positively correlated with axial length only. Many of these structural pathways were also up-regulated over time in the positive lens group (following their initial down-regulation; Figure 5C), further emphasizing that expression levels depend on the stage of refractive compensation (i.e., early time-points when refractive compensation is rapid vs. later time-points when refractive compensation is almost complete).
Figure 5
Eight metabolic pathways were positively correlated with eye size and negatively correlated with refraction, and an additional 3 pathways were uniquely associated with either axial length or refraction (Figure 6A). The earliest metabolic processes implicated were fatty acid and sphingolipid metabolism on day 1. These pathways all showed bidirectional expression responses (core genes demonstrated roughly proportional up-regulation in the negative lens group and down-regulation in the positive lens group; Figure 6B). By 2 days expression shifts were evident in pathways downstream of fatty acid metabolism (Figure 6A), however responses diverged across the two lens groups. The citrate (TCA) cycle pathway and a cluster of pathways related to mitochondrial metabolism were strongly up-regulated in the negative lens group, while the butanoate pathway (with core genes primarily related to ketogenesis) was strongly down-regulated in the positive lens group (Figure 6B). The citrate cycle pathway was also down-regulated over time within the No Lens group (Figure 6C). The higher expression of genes from this pathway in the negative lens group at 2 days appeared to result from a failure to follow this normal trajectory of downregulation.
Figure 6
The remaining pathways correlated with eye size and refraction were related to apoptosis and immune processes (Figure 7A). Apoptosis and immune-related pathways were positively correlated with axial length and negatively correlated with refraction on day 1 (Figure 7A). By day 3, a further immune pathway (primary immunodeficiency) was negatively correlated with axial length and positively correlated with refraction (Figures 7A,B). Additionally, the primary immunodeficiency pathway was up-regulated over time in the positive lens group and down-regulated over time in the negative lens group (the only pathway to show such an expression pattern in the present study; Figure 7C). Two pathways from the apoptosis cluster were also up-regulated over time within the positive lens group. These within and across group GSEA results suggest that the profile of immune and apoptosis pathway expression may reverse across the time-course of refractive compensation (particularly for chicks wearing positive lenses).
Figure 7
Discussion
It has previously been suggested that distinct rather than bidirectional genetic mechanisms underlie eye growth increases during negative lens wear and growth slowing during positive lens wear (Ashby and Feldkaemper, ; Stone et al., ). However, we show here that gene expression across a range of structural, metabolic, and immune pathways is correlated with eye size and refractive state during myopia and hyperopia induction in chick. Thus, although responses at the single gene level are primarily distinct, larger transcript networks show extensive subtle bidirectional expression shifts.
The types of pathways implicated, as well as their correlation with eye size and refraction, suggest a close link with the morphological phenotypes of myopia and hyperopia. Expression of extracellular matrix gene sets was positively correlated with eye size and negatively correlated with refraction across lens groups on day 1. Additional tight junction and cell contraction pathways were positively correlated with ocular axial length only. The proteins encoded by these genes play a fundamental role in tissue morphogenesis (Heisenberg and Bellaïche, ) and the osmotic stress response (Brocker et al., ). These findings are consistent with a wide range of studies detailing changes in cell growth (Teakle et al., 1993; Troilo et al., 1996; Beresford et al., ), proliferation (Fischer and Reh, ; Tkatchenko et al., 2006), ionic concentrations (Seko et al., ; Liang et al., ; Crewther et al., ), and ion and water channel expression (Goodyear et al., , ; Zhang et al., 2011) across the posterior eye during myopia induction. Structural changes during refractive error induction presumably alter the availability of metabolites from the choroidal blood supply (Shih et al., ), and the need for the metabolic production of energy and biosynthetic precursors to fuel cell growth (Lunt and Vander Heiden, ). Accordingly, expression of metabolic genes (primarily within fatty acid and mitochondrial metabolism pathways) was also positively correlated with ocular axial length and negatively correlated with refraction at all time-points (although proportional bidirectional shifts across both lens-groups were only seen on days 1 and 3). Finally, concurrent with the changes in fatty acid and structural pathway expression on day 1, apoptosis and immune pathways were positively correlated with eye length and negatively correlated with refraction. Although a role for immune processes in refractive compensation is less established, these findings are consistent with growing evidence in the refractive error field (Lazuk and Slepova, ; Mao et al., ; Mcglinn et al., ; Shelton et al., ; Long et al., ; Gao et al., ), as well as literature linking lipid metabolism and osmotic stress with immune and inflammatory responses (Daynes and Jones, ; Brocker et al., ; Feske et al., ).
Our GSEA results are concordant with the findings of a recent study in mice, where differentially expressed proteins during myopia induction were enriched for cytoskeletal remodeling and cell adhesion processes, as well as unsaturated fatty acid beta-oxidation and oxidative phosphorylation metabolic pathways (Barathi et al., ). Proteomics studies have also implicated APOA1 (Bertrand et al., ), a lipid transport protein in the PPAR signaling cascade. At the transcriptome level, previous studies have linked a range of structural, metabolic, and immune pathways with ocular growth control (Tkatchenko et al., 2006; Brand et al., ; Mcglinn et al., ; Shelton et al., ; Stone et al., ), however our study is the first to identify bidirectional expression responses in these pathways during myopia and hyperopia induction.
Although our pathway findings are primarily novel, at the single gene level our results contain many commonalities with past studies. Here, we replicated past findings suggesting a role BMP and Wnt signaling in co-ordinating structural change (Mcglinn et al., ; Stone et al., ; Zhang et al., 2012; Ma et al., ). Indeed, BMP2 was the only gene to show a sign-of-defocus dependent expression pattern in the present study. Genes mediated by BMP signaling (Nakanishi et al., ; Parisi et al., ; Inai et al., ; Chang et al., ) were also implicated, including HAS2 and PTX3 which were up-regulated during hyperopia induction and CTGF which was down-regulated during myopia induction. These findings are consistent with previous studies linking CTGF expression in retina/RPE (Mcglinn et al., ), and PTX3 and hyaluronan expression in the choroid (Nickla et al., ; Summers Rada et al., ; He et al., ) with periods of altered ocular growth. NPR3, a gene that shows growth-specific expression shifts in the tree shrew sclera (Guo et al., ), was also up-regulated during hyperopia induction. Notably, NPR3 has been associated with fluid transport across the RPE (Mikami et al., ; Dahrouj et al., ), while concurrent up-regulation of the inflammatory biomarker PTX3 and hyaluronan promotes fluid accumulation in the extracellular space of other tissues (Day and De La Motte, ) suggesting that these three transcripts may play a role in choroidal expansion during hyperopia induction.
Genes associated with circadian and phototransduction processes were also differentially expressed during both myopia and hyperopia induction. These findings are similar to those of a past chick microarray study, where several circadian genes were differentially expressed during myopia induction (Stone et al., ). These previous findings contributed to interest in the role of circadian processes as an explanation for the effects of outdoor time on refractive development in children (Stone et al., ). Careful examination of our own data suggests that circadian and phototransduction gene expression changed during the slowing of growth in normally developing eyes. These same genes were differentially expressed when comparing lens wearing and no lens animals because lens wear (of both signs) accelerated the time-course of these developmental expression changes. Thus, the role of circadian genes in refractive compensation appears more complex than previously thought, and may be related to changes in photoreceptor functioning. As noted in the methods (Differential Gene Expression), the relatively small sample sizes used in the present study meant that the outcome of single gene analyses differed depending on the evaluation method used (see also Soneson and Delorenzi, ). In this context our single gene results should broadly be interpreted with caution, however, it should also be noted that the changes to phototransduction and circadian gene expression across all three groups were robust to the analysis approach used. Future studies are now needed to validate and further explore these single gene changes at the mRNA and protein level (e.g., using qPCR and Western blots).
As the ultimate goal of animal studies of ocular growth is to better understand the mechanisms underlying myopia development in human populations, we compared our findings in chick with the results of Genome-Wide Association Studies (GWAS) indexed in the NHGRI catalog (Welter et al., 2014). BMP2 was the only single gene from our study to fall near myopia-associated SNPs (see Supplementary Table S8 for details). At the pathway level, however, the KEGG “Arrhythmogenic Right Ventricular Cardiomyopathy” gene set (which was positively correlated with axial length and negatively correlated with refraction on 1 day) has been previously linked with human myopia (Hysi et al., ). More broadly, genes associated with refractive error in human populations have pleiotropic effects on non-ocular systems where they are strongly associated with systemic phenotypes including trans-fatty acid levels (Hysi et al., ; Hysi, ). That both GWAS and transcriptome methodologies have implicated multiple (but different) genes within these ECM and fatty acid pathways suggests that the co-ordinated actions of broader biological networks are more functionally significant for growth control than single genes.
Although pathway analysis is a valuable tool for investigating ocular growth phenotypes, our novel findings suggest that the outcome of such investigations is strongly influenced by experimental design. Few past transcriptome wide studies have concurrently investigated myopia and hyperopia induction, and none have done so using GSEA. Instead, researchers have used analysis methods that associate networks or pathways with lists of differentially-expressed genes (e.g., Brand et al., ; Shelton et al., ; Stone et al., ). This approach suffers from poor sensitivity, especially when expression changes are subtle (Abatangelo et al., ; Bayerlová et al., ). Because GSEA does not require an arbitrary cut-off for differential gene expression it has a much larger functional range. This greater sensitivity makes GSEA particularly suited to refractive error datasets where expression changes are often modest [for example during hyperopia induction in the present–and previous (Stone et al., )–studies]. However, GSEA is not without limitations; it shows poor specificity in some circumstances (Bayerlová et al., ) and (in the present study) the results are limited to the gene sets available in the KEGG database which has a particular focus on metabolic and signaling pathways (Bauer-Mehren et al., ; García-Campos et al., ). Future studies may benefit from using a wider range of databases and analysis methods. Regarding the latter, methods that incorporate pathway topology are promising but require further benchmarking (Bayerlová et al., ).
Our novel findings may also reflect the tissues profiled; like several previous transcriptome studies (Mcglinn et al., ; Shelton et al., ; Summers Rada and Wiechmann, ; Stone et al., ) we analyzed a combination of posterior ocular tissues. This systems-level approach proved useful in identifying expression shifts likely to be localized (such as sets of photoreceptor-specific genes) through to broader structural and metabolic shifts that presumably affect multiple cell-types. However, it has previously been shown that profiling multiple cell types can obscure localized expression changes, and that some components of the expression response during visually-regulated growth are reversed in sign across different ocular layers (Ashby and Feldkaemper, ; Penha et al., ). Thus, further studies are now needed to localize the expression shifts identified here to individual ocular layers, and to determine whether the bidirectional expression shifts observed play a direct role in mediating ocular growth and refractive change.
In conclusion, our data demonstrate that expression of genes in several structural, metabolic, and immune pathways is correlated with eye size and refraction across a spectrum of ocular growth conditions in chick (normal development, myopia induction, and hyperopia induction). These findings elucidate the transcriptional response underlying the broad morphological changes that occur in the retina, RPE, and choroid during refractive compensation. The involvement of metabolic and immune/apoptosis pathways suggests a further link between structural change and tissue health that may increase the vulnerability of myopic chick eyes to secondary pathologies (see Hayes et al., ; Liang et al., , ; who describe signs of secondary pathologies in chick). Moreover, similar observations in slower primate models (Tkatchenko et al., 2006) and GSEA studies (see above) support the notion that these expression shifts may translate to human myopia where related biological processes [i.e., apoptosis (Xu et al., 1996), mechanical (Saw et al., ; Morgan et al., ) and oxidative stress (Francisco et al., )] have been linked with the development of sight-threatening secondary disorders.
Additional information
Accession code: GSE78042.
Conflict of interest statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Statements
Author contributions
NR, LG, NH, and SC conceived, designed and coordinated the study. NR collected the data. NR and NH conducted the data analysis. NR wrote and revised the manuscript with input from LG, NH, and SC.
Acknowledgments
We would like to thank Melanie Murphy for help with biometric data collection, and David Crewther for valuable comments on the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins.2016.00390
References
1
AbatangeloL.MagliettaR.DistasoA.D'AddabboA.CreanzaT. M.MukherjeeS.et al. (2009). Comparative study of gene set enrichment methods. BMC Bioinformatics10:275. 10.1186/1471-2105-10-275
2
AndersS.PylP. T.HuberW. (2014). HTSeq–A Python framework to work with high-throughput sequencing data. Bioinformatics31, 166–169. 10.1093/bioinformatics/btu638
3
AshbyR. S.FeldkaemperM. P. (2010). Gene expression within the amacrine cell layer of chicks after myopic and hyperopic defocus. Invest. Ophthalmol. Vis. Sci.51, 3726–3735. 10.1167/iovs.09-4615
4
BarathiV. A.ChaurasiaS. S.PoidingerM.KohS. K.TianD.HoC.et al. (2014). Involvement of GABA transporters in atropine-treated myopic retina as revealed by iTRAQ quantitative proteomics. J. Proteome Res.13, 4647–4658. 10.1021/pr500558y
5
Bauer-MehrenA.FurlongL. I.SanzF. (2009). Pathway databases and tools for their exploitation: benefits, current limitations and challenges. Mol. Syst. Biol.5, 290. 10.1038/msb.2009.47
6
BayerlováM.JungK.KramerF.KlemmF.BleckmannA.BeißbarthT. (2015). Comparative study on gene set and pathway topology-based enrichment methods. BMC Bioinformatics16:334. 10.1186/s12859-015-0751-5
7
BeresfordJ. A.CrewtherS. G.CrewtherD. P. (1998). Anatomical correlates of experimentally induced myopia. Aust. N. Z. J. Ophthalmol.26, S84–S87. 10.1111/j.1442-9071.1998.tb01383.x
8
BermanE. R.MichaelsonI. C. (1964). The chemical composition of the human vitreous body as related to age and myopia. Exp. Eye Res.3, 9–15. 10.1016/S0014-4835(64)80003-8
9
BertrandE.FritschC.DietherS.LambrouG.MüellerD.SchaeffelF.et al. (2006). Identification of apolipoprotein A-I as a “STOP” signal for myopia. Mol. Cell. Proteomics5, 2158–2166. 10.1074/mcp.M600073-MCP200
10
BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics30, 2114–2120. 10.1093/bioinformatics/btu170
11
BrandC.SchaeffelF.FeldkaemperM. P. (2007). A microarray analysis of retinal transcripts that are controlled by image contrast in mice. Mol. Vis.13, 920–932.
12
BrockerC.ThompsonD. C.VasiliouV. (2012). The role of hyperosmotic stress in inflammation and disease. Biomol. Concepts3, 345–364. 10.1515/bmc-2012-0001
13
ChangH. M.ChengJ. C.FangL.QiuX.KlausenC.TaylorE. L.et al. (2015). Recombinant BMP4 and BMP7 downregulate pentraxin 3 in human granulosa cells. J. Clin. Endocrinol. Metab.100, E365–E374. 10.1210/jc.2014-2496
14
CrewtherS. G.LiangH.JunghansB. M.CrewtherD. P. (2006). Ionic control of ocular growth and refractive change. Proc. Natl. Acad. Sci. U.S.A.103, 15663–15668. 10.1073/pnas.0607241103
15
CurtinB. J.TengC. C. (1958). Scleral changes in pathological myopia. Trans. Am. Acad. Ophthalmol. Otolaryngol.62, 777.
16
DahroujM.AlsarrafO.LiuY.CrossonC. E.AblonczyZ. (2013). C-type natriuretic peptide protects the retinal pigment epithelium against advanced glycation end product-induced barrier dysfunction. J. Pharmacol. Exp. Ther.344, 96–102. 10.1124/jpet.112.199307
17
DayA. J.De La MotteC. A. (2005). Hyaluronan cross-linking: a protective mechanism in inflammation?Trends Immunol.26, 637–643. 10.1016/j.it.2005.09.009
18
DaynesR. A.JonesD. C. (2002). Emerging roles of PPARs in inflammation and immunity. Nat. Rev. Immunol.2, 748–759. 10.1038/nri912
19
FeldkaemperM.SchaeffelF. (2013). An updated view on the role of dopamine in myopia. Exp. Eye Res.114, 106–119. 10.1016/j.exer.2013.02.007
20
FeskeS.WulffH.SkolnikE. Y. (2015). Ion channels in innate and adaptive immunity. Annu. Rev. Immunol.33, 291–353. 10.1146/annurev-immunol-032414-112212
21
FischerA. J.RehT. A. (2000). Identification of a proliferating marginal zone of retinal progenitors in postnatal chickens. Dev. Biol.220, 197–210. 10.1006/dbio.2000.9640
22
FranciscoB.-M.SalvadorM.AmparoN. (2015). Oxidative stress in myopia. Oxid. Med. Cell. Longev.2015:750637. 10.1155/2015/750637
23
GaoT.-T.LongQ.YangX. (2015). Complement factors C1q, C3 and C5b-9 in the posterior sclera of guinea pigs with negative lens-defocused myopia. Int. J. Ophthalmol.8, 675. 10.3980/j.issn.2222-3959.2015.04.06
24
García-CamposM. A.Espinal-EnríquezJ.Hernandez-LemusE. (2015). Pathway analysis: state of the art. Front. Physiol.6:383. 10.3389/fphys.2015.00383
25
GoodyearM. J.CrewtherS. G.MurphyM. J.GuimmarraL.HaziA.JunghansB. M.et al. (2010). Spatial and temporal dissociation of AQP4 and Kir4.1 expression during induction of refractive errors. Mol. Vis.16, 1610–1619.
26
GoodyearM. J.JunghansB. M.GiummarraL.MurphyM. J.CrewtherD. P.CrewtherS. G. (2008). A role for aquaporin-4 during induction of form deprivation myopia in chick. Mol. Vis.14, 298–307.
27
GuoL.FrostM. R.SiegwartJ. T.Jr.NortonT. T. (2014). Scleral gene expression during recovery from myopia compared with expression during myopia development in tree shrew. Mol. Vis.20, 1643–1659.
28
HarwellM. R.RubinsteinE. N.HayesW. S.OldsC. C. (1992). Summarizing Monte Carlo results in methodological research: the one-and two-factor fixed effects ANOVA cases. J. Educ. Behav. Stat.17, 315–339.
29
HayesB.FitzkeF.HodosW.HoldenA. (1986). A morphological analysis of experimental myopia in young chickens. Invest. Ophthalmol. Vis. Sci.27, 981–991.
30
HeL.FrostM. R.SiegwartJ. T.Jr.NortonT. T. (2014). Gene expression signatures in tree shrew choroid during lens-induced myopia and recovery. Exp. Eye Res.123, 56–71. 10.1016/j.exer.2014.04.005
31
HeisenbergC. P.BellaïcheY. (2013). Forces in tissue morphogenesis and patterning. Cell153, 948–962. 10.1016/j.cell.2013.05.008
32
HughesT. R. (2009). ‘Validation’ in genome-scale research. J. Biol.8:3. 10.1186/jbiol104
33
HysiP. (2015). Commonalities of genetic mechanisms between refractive error and systemic phenotypes: A systematic review of evidence of genetic pleiotropy in genome-wide association studies, in The 15th International Myopia Conference (Wenzhou).
34
HysiP. G.MahrooO. A.CumberlandP.WojciechowskiR.WilliamsK. M.YoungT. L.et al. (2014a). Common mechanisms underlying refractive error identified in functional analysis of gene lists from genome-wide association study results in 2 European British cohorts. JAMA Ophthalmol.132, 50–56. 10.1001/jamaophthalmol.2013.6022
35
HysiP. G.WojciechowskiR.RahiJ. S.HammondC. J. (2014b). Genome-wide association studies of refractive error and myopia, lessons learned, and implications for the future. Invest. Ophthalmol. Vis. Sci.55, 3344–3351. 10.1167/iovs.14-14149
36
InaiK.BurnsideJ. L.HoffmanS.TooleB. P.SugiY. (2013). BMP-2 induces versican and hyaluronan that contribute to post-EMT AV cushion cell migration. PLoS ONE8:e77593. 10.1371/journal.pone.0077593
37
JinN.StjernschantzJ. (2000). Regional blood flow in the myopic chick eye during and after form deprivation: a study with radioactively-labelled microspheres. Exp. Eye Res.71, 233–238. 10.1006/exer.2000.0871
38
KamburovA.WierlingC.LehrachH.HerwigR. (2009). ConsensusPathDB–a database for integrating human functional interaction networks. Nucleic Acids Res.37, D623–D628. 10.1093/nar/gkn698
39
KanehisaM.GotoS. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res.28, 27–30. 10.1093/nar/28.1.27
40
KimD.PerteaG.TrapnellC.PimentelH.KelleyR.SalzbergS. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol.14:R36. 10.1186/gb-2013-14-4-r36
41
LangmeadB.SalzbergS. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods9, 357–359. 10.1038/nmeth.1923
42
LazukA. V.SlepovaO. S. (1994). Study of immune reactions to collagen in patients with myopia. Vestn. Oftalmol.111, 14–16.
43
LiangH.CrewtherD. P.Gillard CrewtherS. G.BarilaA. M. (1995). A role for photoreceptor outer segments in the induction of deprivation myopia. Vis. Res.35, 1217–1225. 10.1016/0042-6989(94)00241-D
44
LiangH.CrewtherS. G.CrewtherD. P.JunghansB. M. (2004). Structural and elemental evidence for edema in the retina, retinal pigment epithelium, and choroid during recovery from experimentally induced myopia. Invest. Ophthalmol. Vis. Sci.45, 2463–2474. 10.1167/iovs.03-1009
45
LiberzonA.SubramanianA.PinchbackR.ThorvaldsdóttirH.TamayoP.MesirovJ. P. (2011). Molecular signatures database (MSigDB) 3.0. Bioinformatics27, 1739–1740. 10.1093/bioinformatics/btr260
46
LongQ.YeJ.LiY.WangS.JiangY. (2013). C-reactive protein and complement components in patients with pathological myopia. Optom. Vis. Sci.90, 501–506. 10.1097/OPX.0b013e31828daa6e
47
LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15:550. 10.1186/s13059-014-0550-8
48
LuntS. Y.Vander HeidenM. G. (2011). Aerobic glycolysis: meeting the metabolic requirements of cell proliferation. Annu. Rev. Cell Dev. Biol.27, 441–464. 10.1146/annurev-cellbio-092910-154237
49
MaM.ZhangZ.DuE.ZhengW.GuQ.XuX.et al. (2014). Wnt signaling in form deprivation myopia of the mice retina. PLoS ONE9:e91086. 10.1371/journal.pone.0091086
50
ManoliT.GretzN.GröneH. J.KenzelmannM.EilsR.BrorsB. (2006). Group testing for pathway analysis improves comparability of different microarray datasets. Bioinformatics22, 2500–2506. 10.1093/bioinformatics/btl424
51
MaoJ.LiuS.WenD.TanX.FuC. (2006). Basic fibroblast growth factor suppresses retinal neuronal apoptosis in form-deprivation myopia in chicks. Curr. Eye Res.31, 983–987. 10.1080/02713680600910510
52
MartinM. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J.17, 10–12. 10.14806/ej.17.1.200
53
McglinnA. M.BaldwinD. A.TobiasJ. W.BudakM. T.KhuranaT. S.StoneR. A. (2007). Form-deprivation myopia in chick induces limited changes in retinal gene expression. Invest. Ophthalmol. Vis. Sci.48, 3430–3436. 10.1167/iovs.06-1538
54
MericoD.IsserlinR.StuekerO.EmiliA.BaderG. D. (2010). Enrichment map: a network-based method for gene-set enrichment visualization and interpretation. PLoS ONE5:e13984. 10.1371/journal.pone.0013984
55
MikamiY.HaraM.YasukuraT.UyamaM.MinatoA.InagakiC. (1995). Atrial natriuretic peptide stimulates Cl- transport in retinal pigment epithelial cells. Curr. Eye Res.14, 391–397. 10.3109/02713689508999937
56
MorganI. G.Ohno-MatsuiK.SawS.-M. (2012). Myopia. Lancet379, 1739–1748. 10.1016/S0140-6736(12)60272-4
57
NakanishiT.KimuraY.TamuraT.IchikawaH.YamaaiY.-I.SugimotoT.et al. (1997). Cloning of a mRNA preferentially expressed in chondrocytes by differential display-PCR from a human chondrocytic cell line that is identical with connective tissue growth factor (CTGF) mRNA. Biochem. Biophys. Res. Commun.234, 206–210. 10.1006/bbrc.1997.6528
58
NicklaD. L.WildsoetC.WallmanJ. (1997). Compensation for spectacle lenses involves changes in proteoglycan synthesis in both the sclera and choroid. Curr. Eye Res.16, 320–326. 10.1076/ceyr.16.4.320.10697
59
PanC. W.RamamurthyD.SawS. M. (2012). Worldwide prevalence and risk factors for myopia. Ophthalmic Physiol. Opt.32, 3–16. 10.1111/j.1475-1313.2011.00884.x
60
ParisiM. S.GazzerroE.RydzielS.CanalisE. (2006). Expression and regulation of CCN genes in murine osteoblasts. Bone38, 671–677. 10.1016/j.bone.2005.10.005
61
PenhaA. M.SchaeffelF.FeldkaemperM. (2011). Insulin, insulin-like growth factor-1, insulin receptor, and insulin-like growth factor-1 receptor expression in the chick eye and their regulation with imposed myopic or hyperopic defocus. Mol. Vis.17, 1436–1448. 10.1167/iovs.09-4615
62
Summers RadaJ. A.WiechmannA. F. (2009). Ocular expression of avian thymic hormone: changes during the recovery from induced myopia. Mol. Vis.15, 778–792.
63
Summers RadaJ. A.WiechmannA. F.HollawayL. R.BaggenstossB. A.WeigelP. H. (2010). Increased hyaluronan synthase-2 mRNA expression and hyaluronan accumulation with choroidal thickening: response during recovery from induced myopia. Invest. Ophthalmol. Vis. Sci.51, 6172–6179. 10.1167/iovs.10-5522
64
R Core Team (2013). R: A Language and Environment for Statistical Computing.Vienna: R Foundation for Statistical Computing.
65
RobinsonM. D.MccarthyD. J.SmythG. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics26, 139–140. 10.1093/bioinformatics/btp616
66
RogersL. J.BoldenS. W. (1991). Light-dependent development and asymmetry of visual projections. Neurosci. Lett.121, 63–67. 10.1016/0304-3940(91)90650-I
67
RymerJ.WildsoetC. F. (2005). The role of the retinal pigment epithelium in eye growth regulation and myopia: a review. Vis. Neurosci.22, 251. 10.1017/S0952523805223015
68
SawS. M.GazzardG.Shih-YenE. C.ChuaW. H. (2005). Myopia and associated pathological complications. Ophthalmic Physiol. Opt.25, 381–391. 10.1111/j.1475-1313.2005.00298.x
69
SchippertR.SchaeffelF.FeldkaemperM. P. (2008). Microarray analysis of retinal gene expression in chicks during imposed myopic defocus. Mol. Vis.14, 1589–1599.
70
SeetB.WongT. Y.TanD. T.SawS. M.BalakrishnanV.LeeL. K.et al. (2001). Myopia in Singapore: taking a public health approach. Br. J. Ophthalmol.85, 521–526. 10.1136/bjo.85.5.521
71
SekoY.ShimokawaH.PangJ.TokoroT. (2000). Disturbance of electrolyte balance in vitreous of chicks with form-deprivation myopia. Jpn. J. Ophthalmol.44, 15–19. 10.1016/S0021-5155(99)00177-X
72
SheltonL.TroiloD.LernerM. R.GusevY.BrackettD. J.RadaJ. S. (2008). Microarray analysis of choroid/RPE gene expression in marmoset eyes undergoing changes in ocular growth and refraction. Mol. Vis.14, 1465–1479.
73
ShihY.-F.FitzgeraldM. E.NortonT. T.GamlinP. D.HodosW.ReinerA. (1993). Reduction in choroidal blood flow occurs in chicks wearing goggles that induce eye growth toward myopia. Curr. Eye Res.12, 219–227. 10.3109/02713689308999467
74
SonesonC.DelorenziM. (2013). A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics14:91. 10.1186/1471-2105-14-91
75
StoneR. A.KhuranaT. S. (2010). Gene profiling in experimental models of eye growth: clues to myopia pathogenesis. Vis. Res.50, 2322–2333. 10.1016/j.visres.2010.03.021
76
StoneR. A.McglinnA. M.BaldwinD. A.TobiasJ. W.IuvoneP. M.KhuranaT. S. (2011). Image defocus and altered retinal gene expression in chick: clues to the pathogenesis of ametropia. Invest. Ophthalmol. Vis. Sci.52, 5765–5777. 10.1167/iovs.10-6727
77
StoneR. A.PardueM. T.IuvoneP. M.KhuranaT. S. (2013). Pharmacology of myopia and potential role for intrinsic retinal circadian rhythms. Exp. Eye Res.114, 35–47. 10.1016/j.exer.2013.01.001
78
Suárez-FariñasM.LowesM. A.ZabaL. C.KruegerJ. G. (2010). Evaluation of the psoriasis transcriptome across different studies by gene set enrichment analysis (GSEA). PLoS ONE5:e10247. 10.1371/journal.pone.0010247
79
SubramanianA.KuehnH.GouldJ.TamayoP.MesirovJ. P. (2007). GSEA-P: a desktop application for Gene Set Enrichment Analysis. Bioinformatics23, 3251–3253. 10.1093/bioinformatics/btm369
80
SubramanianA.TamayoP.MoothaV. K.MukherjeeS.EbertB. L.GilletteM. A.et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U.S.A.102, 15545–15550. 10.1073/pnas.0506580102
81
TeakleE. M.WildsoetC. F.VaneyD. I. (1993). The spatial organization of tyrosine hydroxylase-immunoreactive amacrine cells in the chicken retina and the consequences of myopia. Vis. Res.33, 2383–2396. 10.1016/0042-6989(93)90117-F
82
TkatchenkoA. V.WalshP. A.TkatchenkoT. V.GustincichS.RaviolaE. (2006). Form deprivation modulates retinal neurogenesis in primate experimental myopia. Proc. Natl. Acad. Sci. U.S.A.103, 4681–4686. 10.1073/pnas.0600589103
83
TroiloD.GottliebM. D.WallmanJ. (1987). Visual deprivation causes myopia in chicks with optic nerve section. Curr. Eye Res.6, 993–999. 10.3109/02713688709034870
84
TroiloD.XiongM.CrowleyJ. C.FinlayB. L. (1996). Factors controlling the dendritic arborization of retinal ganglion cells. Vis. Neurosci.13, 721–733. 10.1017/S0952523800008609
85
VaretH.Brillet-GuéguenL.CoppéeJ.-Y.DilliesM.-A. (2016). SARTools: a DESeq2-and edgeR-based R pipeline for comprehensive differential analysis of RNA-Seq data. PLoS ONE11:e0157022. 10.1371/journal.pone.0157022
86
VerhoevenV. J.WongK. T.BuitendijkG. H.HofmanA.VingerlingJ. R.KlaverC. C. (2015). Visual consequences of refractive errors in the general population. Ophthalmology122, 101–109. 10.1016/j.ophtha.2014.07.030
87
WallmanJ.WildsoetC.XuA.GottliebM. D.NicklaD. L.MarranL.et al. (1995). Moving the retina: choroidal modulation of refractive state. Vis. Res.35, 37–50. 10.1016/0042-6989(94)E0049-Q
88
WallmanJ.WinawerJ. (2004). Homeostasis of eye growth and the question of myopia. Neuron43, 447–468. 10.1016/j.neuron.2004.08.008
89
WalterW.Sánchez-CaboF.RicoteM. (2015). GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics31, 2912–2914. 10.1093/bioinformatics/btv300
90
WangZ.GersteinM.SnyderM. (2009). RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet.10, 57–63. 10.1038/nrg2484
91
WelterD.MacarthurJ.MoralesJ.BurdettT.HallP.JunkinsH.et al. (2014). The NHGRI GWAS Catalog, a curated resource of SNP-trait associations. Nucleic Acids Res.42, D1001–D1006. 10.1093/nar/gkt1229
92
WestbrookA. M.CrewtherS. G.BeresfordJ. A.AllenM.KellerI.CrewtherD. P. (1995). Formoguanamine-induced inhibition of deprivation myopia in chick is accompanied by choroidal thinning while retinal function is retained. Vis. Res.35, 2075–2088. 10.1016/0042-6989(94)00282-Q
93
WildsoetC.PettigrewJ. (1988). Experimental myopia and anomalous eye growth-patterns unaffected by optic-nerve section in chickens-evidence for local-control of eye growth. Clin. Vis. Sci.3, 99–107.
94
WildsoetC.WallmanJ. (1995). Choroidal and scleral mechanisms of compensation for spectacle lenses in chicks. Vis. Res.35, 1175–1194. 10.1016/0042-6989(94)00233-C
95
XuG. Z.LiW. W.TsoM. O. (1996). Apoptosis in human retinal degenerations. Trans. Am. Ophthalmol. Soc.94, 411.
96
ZhangH.WongC. L.ShanS. W.LiK. K.ChengA. K.LeeK. L.et al. (2011). Characterisation of Cl-transporter and channels in experimentally induced myopic chick eyes. Clin. Exp. Optom.94, 528–535. 10.1111/j.1444-0938.2011.00611.x
97
ZhangY.LiuY.WildsoetC. F. (2012). Bidirectional, optical sign-dependent regulation of BMP2 gene expression in chick retinal pigment epithelium. Invest. Ophthalmol. Vis. Sci.53, 6072–6080. 10.1167/iovs.12-9917
Summary
Keywords
myopia, hyperopia, refractive error, transcriptome, RNA sequencing, gene expression, chick, retina
Citation
Riddell N, Giummarra L, Hall NE and Crewther SG (2016) Bidirectional Expression of Metabolic, Structural, and Immune Pathways in Early Myopia and Hyperopia. Front. Neurosci. 10:390. doi: 10.3389/fnins.2016.00390
Received
06 June 2016
Accepted
09 August 2016
Published
30 August 2016
Volume
10 - 2016
Edited by
Michael F. Miles, Virginia Commonwealth University, USA
Reviewed by
Claudia Vianna Maurer-Morelli, State University of Campinas, Brazil; Marita Pauline Feldkaemper, Eye Clinic, Germany
Updates
Copyright
© 2016 Riddell, Giummarra, Hall and Crewther.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Sheila G. Crewther s.crewther@latrobe.edu.au
This article was submitted to Neurogenomics, a section of the journal Frontiers in Neuroscience
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.