Retinal genomic fabrics remodeling after optic nerve injury

Glaucoma is a multifactorial neurodegenerative disease, characterized by degeneration of the retinal ganglion cells (RGCs). There has been little progress in developing efficient strategies for neuroprotection in glaucoma. We profiled the retina transcriptome of Lister Hooded rats at 2 weeks after optic nerve crush (ONC) and analyzed the data from the Genomic Fabric Paradigm (GFP) to bring additional insights into the molecular mechanisms of the retinal remodeling after induction of RGC degeneration. GFP considers for the expression of each gene 3 independent characteristics: level, variability and correlation with each other gene. Thus, the 17,657 quantified genes our study generated a total of 155,911,310 values to analyze. This represents 8,830x more data per condition than a traditional transcriptomic analysis. ONC led to a 57% reduction in RGC numbers as detected by retrograde labeling with DiI. We observed a higher Relative Expression Variability after ONC. Gene expression stability was used as a measure of transcription control and disclosed a robust reduction in the number of very stably expressed genes. Predicted Protein-Protein interaction (PPI) analysis with STRING revealed axon and neuron projection as mostly decreased processes, consistent with RGC degeneration. Conversely, immune response PPIs were found among up-regulated genes. Enrichment analysis showed that Complement Cascade and Notch Signaling Pathway, as well as Oxidative Stress and Kit Receptor Pathway were affected after ONC. To expand our studies of altered molecular pathways, we examined the pair-wise coordination of gene expressions within each pathway and within the entire transcriptome using Pearson correlations. ONC increased the number of synergistically coordinated pairs of genes and the number of similar profiles mainly in Complement Cascade and Notch Signaling Pathway. This deep bioinformatic study provides novel insights beyond the regulation of individual gene expression and discloses changes in the control of expression of Complement Cascade and Notch Signaling functional pathways that may be relevant for both RGC degeneration and remodeling of the retinal tissue after ONC.


Introduction
Glaucoma is a multifactorial disease characterized by degeneration of the neurons known as retinal ganglion cells (RGCs), and their long axons which transmit visual information from the retina to the brain [1]. A variety of in vivo animal models have been used to unravel the mechanisms of both RGC degeneration and reorganization of the retina. The animals were subjected to either acute or chronic elevation of intraocular pressure, transection or crush of the optic nerve, as well as the spontaneous glaucoma-like disease (on mice of the DBA/2J strain) [2][3][4]. However, there has been little progress in developing efficient strategies for neuroprotection in glaucoma.
Statement for the Use of Animals in Ophthalmic and Vision Research, and the protocol approved by the Institutional Animal Experimentation Ethics Committee (UFRJ#01200.001568/2013-87).

Optic nerve crush
The rats were divided into two experimental groups: optic nerve crush (n=24, hereafter denoted by ONC) and control, sham operated, (n=19, denoted by CTR). Rats were deeply anesthetized by intraperitoneal injection of a mixture of xylazine (10 mg/kg body weight) and ketamine (60 mg/kg body weight). Under a stereoscopic microscope, the right optic nerve was accessed through an incision made in the superior orbital rim, and the conjunctiva was dissected with forceps towards the back of the eye to expose the retrobulbar portion of the optic nerve. The lesion was done by crushing the optic nerve for 10 s at 3 mm from the optic disc, using a pair of watchmaker's forceps. Sham operated rats underwent the same procedure, except that the forceps were not closed. Animals were euthanized two weeks after the procedure.

Quantification of retinal ganglion cell survival
Cell survival was quantified based on retrograde labeling with the lipophilic tracer DiI (1,1′-dioctadecyl-3,3,3,3′-tetramethylindocarbocyanine perchlorate, Invitrogen) in 11 CTR and 16 ONC retinas. DiI was bilaterally injected into the superior colliculus of one-week old neonate rats, which leads to the labeling of virtually all the RGCs [45]. Seven weeks after DiI injection, rats underwent unilateral ONC. Two weeks after ONC, animals were euthanized by inhalation of carbon dioxide. At this time point, the cells that remained alive were detectable by DiI labeling and quantified. Upon dissection of the eyeballs, retinas in Phosphate-Buffered Saline (PBS, Sigma) were dissected and flattened as wholemounts by making four radial cuts, followed by fixation with 4% paraformaldehyde for 15 min at room temperature, three rounds of washing with PBS, and counterstaining with Sytox Green (ThermoFisher) to visualize cell nuclei. Whole retinas were mounted vitreal side-up on subbed slides, covered with anti-fading mounting media and examined in a confocal epifluorescence microscope (LSM 510 Meta, Zeiss) using a Plan-Neofluar 40x/1.3 objective. Eight photos were taken from each quadrant of the retina, of which 2 were from central retina (~0.9 mm from optic disc), 3 from mid-retina (~2.0 mm from optic disc) and 3 from peripheral retina (~3.7 mm from optic disc), to a total of 32 photos per retina. Cell counts were done in a double-blind manner based on morphology to exclude microglial cells with translocated DiI, and DiI+ RGCs were averaged across evenly distributed fields of each retina. Statistical analysis was done through an unpaired two-tailed t-test using the software GraphPad Prism 5.02 (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004) GraphPad Software, Inc.).

Microarray gene expression
Animals were euthanized by inhalation of carbon dioxide. Each retina was freshly dissected with clean instruments and immediately frozen in liquid nitrogen. For each experimental group (CTR or ONC), four separated retinas were profiled individually. We applied an optimized protocol [46] for RNA extraction (Qiagen RNeasy mini-kit), reverse transcription (adding fluorescent tags) and hybridization onto Agilent G2519F 60mer two-colour gene expression rat 4×44k arrays in the "multiple yellow" design that provides maximum flexibility in comparing the conditions and 100% usage of the resources [47]. RNA concentrations before and after reverse transcription were estimated with a NanoDrop ND2000 Spectrophotometer, and purity was assessed with Agilent RNA 6000 Nano kit in an Agilent 2100 Bioanalyzer. RIN (RNA Integrity Number) of at least 8.0 was considered as satisfactory. RINs of the selected ONC samples were: 8.20, 8.40, 9.10, 8.50, while those for the CTR samples were: 8.40, 8.20, 8.50, 8.60. In each array, 825 ng of Cy3 (green, g) or Cy5 (red, r)-labeled RNA from each retina were co-hybridized for 17 h at 65°C with that of the same experimental group in the combinations: CTR1(g)CTR2(r), CTR3(g)CTR4(r), ONC1(g)ONC2(r), ONC3(g)ONC4(r). Chips were scanned with an Agilent G2539A dual laser scanner at 5 μm pixel size/20-bit, and raw data were collected with Agilent Feature Extraction software v. 11.1.1.

Data processing
We disregarded all spots showing signs of local corruption or with foreground fluorescence less than twice their background in any one of the 8 profiled samples. Data were normalized by an in-house developed, iterative method that alternates intra-and inter-array normalization to the median of the background-subtracted fluorescence of the valid spots, until the overall maximum error of estimate reaches less than 5% [48]. Normalized expression levels were organized into redundancy groups composed of all spots probing the same gene and represented by the weighted average of the values of individual spots.
The REV (REV = median of the Bonferroni like-corrected chi-square interval estimate of the pooled coefficient of variation) was taken as a statistical estimate of the expression variability of one gene among biological replicas [49].  The genes were then ordered according to decreasing variability, such that the first percentile (GES < 1) contained the most unstably expressed and the 100th percentile (GES > 99) the most stably expressed genes [50].
A gene was scored as significantly regulated in ONC with respect to CTR if the absolute fold-change |x| exceeded the cut-off (CUT) computed for that gene (Eq. 1), AND the p value of the heteroscedastic t-test for the equal expressions was < 0.05. Thus, we replaced an arbitrary uniform absolute fold-change cut-off (e.g. 1.5x) for all genes with a value that accounts for the combined contributions of the technical noise and biological expression variability of each gene independently [51]. This composite criterion, in which the cut-offs were computed with Bonferroni-type corrections applied to the redundancy groups [52][53], eliminates most of the false positives without increasing the number of false negatives (Eq. 2).

Analysis of predicted protein-protein interactions
Interactions among proteins encoded by differentially expressed genes were assessed using STRING (https://string-db.org/, [54]), which provides uniquely comprehensive coverage and ease of access to information of both experimental and predicted interaction. A protein-protein interaction network was constructed, in which the interactions of ONC vs. CTR differentially expressed genes (DEGs) were mapped to STRING, based on information including the sequence characters and structures. Each protein-protein interaction stored in the STRING database receives a confidence score between zero and one. These scores indicate the estimated probability that an interaction is biologically significant, specific and reproducible. To highlight the most biologically relevant interactions, we used a confidence score cut-off of 0.4.

Enrichment analyses
PathVisio3 (www.pathvisio.org, [55]) was used to identify pathways significantly altered, and the Gene Ontology Knowledgebase (www.geneontology.org) was used to identify main Biological Processes, Molecular Functions and Cellular Components affected by ONC. The pathways were ranked based on a standardized difference indicator (Z score). Pathways with Z score > 2 and p value < 0.05 were accepted as significantly altered. Graphical representation of the molecular pathways was based on the Kyoto Encyclopedia for Genes and Genomes (KEGG, www.genome.jp/kegg/, Kanehisa Laboratories, Japan; [56]).
The EnrichR platform (http://amp.pharm.mssm.edu/Enrichr/, [57]) contains a collection of diverse gene set libraries available for analysis and was used to disclose disease-related perturbations by differentially expressed genes, according to the Gene Expression Omnibus (GEO) database.
Panther (Protein ANalysis THrough Evolutionary Relationships) classification system (http://www.pantherdb.org/) version 15.0 was used to analyze enrichment of genes within a specific subset of genes affected by ONC.

Coordination of expression
We used the variation in gene expression among multiple biological replicas to calculate the pair-wise Pearson correlation coefficients p between the levels of expression of gene pairs. Two genes were scored as: synergistically expressed if their expression levels had a positive covariance within biological replicas; antagonistically expressed when they manifested opposite tendencies (i.e., negative covariance) or independently expressed, when their transcription levels were not correlated (close to zero covariance). In the case of four biological replicas, the (p < 0.05) significant synergistically expressed genes have ρ > 0.90, antagonistically expressed genes have ρ < −0.90; and independently expressed ones have | ρ | < 0.05. The set of correlation coefficients between the expression level of a particular gene and of each other gene within the biological replicas forms the coordination profile of that gene. The profiles of two genes can be similar, opposite or neutral. Two genes were considered as having similar coordination profiles when there is significant overlap of their synergistically, antagonistically and independently expressed partners; genes were considered as having opposite coordination profiles when most of the synergistically expressed partners of one gene were antagonistically expressed partners for the other, and most of the independently expressed partners of one gene were synergistically or antagonistically expressed for the other [52][53].

Quantitative Real-Time PCR
To validate selected genes, a new set of animals (n = 4 per group) was used for qRT-PCR. Total RNA (2 μg) isolated from retinas using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), was reversely transcribed into cDNA with the Su-perScript™ III Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA). qRT-PCR was done with SYBR Green Mix The PCR program was: denaturation at 95°C for 5 s, annealing at 60°C for 30 s, elongation at 68°C for 20 s for 40 cycles. Each sample was applied in technical duplicates, and fold changes were calculated using the 2-(ΔΔCT) method of relative quantification. For each gene tested, all experimental and CTR retinas were assessed in the same experiment. The data were normalized to the geometric mean of the two housekeeper genes and expressed as fold change. Statistical analysis was done using GraphPad Prism Software (Mean ± SEM) using unpaired two-tailed Student's t-test. Differences were accepted as significant when p < 0.05.

Histopathological observations
At 14 days after ONC, there was a robust loss of RGCs retrogradely labeled with DiI, as expected. CTR retinas displayed an average of 1.780±134 cells/mm2, whereas ONC retinas showed a reduction of approximately 60% of RGCs, to 768±65 cells/mm2 (Figure 1A-C). analysis showed that the majority (47%) are expressed in RGC and 12% were described to be decreased in other retinal microarray studies of RGC degeneration. (e) Differentially expressed genes after ONC were plotted as log2 values against their p values in log10. RGC-specific genes which were previously shown to be altered by ONC were highlighted in green. ***: p<0.001, unpaired Student's T test.

Transcriptomic alterations in rat retina after ONC
The experimental details, as well as both raw and normalized expression data from the microarrays were deposited and are publicly available at Gene Expression Omnibus (GEO -accession number GSE133563). Among 17,657 unigenes quantified in our microarrays, expression of 193 (1.1 %) was significantly altered in ONC samples with respect to CTR, of which 127 (65.8 %) were up-and 66 (34.2%) were down-regulated. Table 1 presents the most down-regulated 50 genes and Table 2 the most up-regulated 50 genes.
Quantitative RT-PCR was used to validate two of the highest up-regulated genes (Cd74 and C3) and three markers of RGCs (Tubb3, Nefm, Nell2) that were found to be down-regulated in ONC retinas as compared to CTR ( Table 3). Cd74 and C3 were found to be up-regulated by qPCR by 28.56-and 3.24-fold, respectively. Tubb3, Nefm and Nell2 were downregulated by -3.93-, -19.27-and -4.15-fold in ONC retinas when compared to CTR ones, consistent with the results of the arrays ( Table 3).

Protein-protein interaction encoded by the most down-regulated genes
We used the STRING platform to predict protein-protein interaction (PPI) networks. Among the down-regulated genes, only 22 showed a network of 21 interactions, with a medium confidence score of 0.4. Enrichment analysis using the Gene Ontology database indicated that the main Biological Processes altered by ONC were Intermediate Filament Cytoskeleton Organization, Intermediate Filament bundle Assembly and Response to Acrylamide, all of which include neurofilaments (Figure 2A), a major component of the optic nerve axons. Among Molecular Functions, ONC affected networks of Protein Binding, Binding and Protein Heterodimerization Activity ( Figure 2B). Finally, among cellular components, ONC affected Cell Projection, Axon and Neuron Projection ( Figure 2C). All such functional interactions are consistent with retrograde degeneration as a consequence of ONC.

Figure 2. Protein-protein interaction (PPI) network composed with down-regulated genes after the ONC.
Among 66 down-regulated genes, 21 interactions between 22 genes were found using STRING software, with a confidence cut-off of 0.4. Nodes labeled with the encoding gene symbol indicate proteins and the lines represent the corresponding interactions. The confidence score of each interaction is mapped to the line thickness (the thicker the line, the more evidence to support the interaction). The network was then enriched according to Gene Ontology database. The categories of Gene Ontology are depicted: (a) Biological Processes; (b) Molecular Function and (c) Cellular Components and the 3 most significantly enriched categories were used to color the nodes of the interaction networks.

Protein-protein interaction encoded by the most up-regulated genes
Consistent with previous studies [15-18; 28], the group of 50 most up-regulated genes ( Table 2) includes inflammatory response-related genes, such as Cd74 molecule, major complex, class II invariant chain (25.6x), H-2 class II histocompatibility antigen gamma chain (MHC class II-associated invariant chain) (Ia antigen-associated invariant chain) (Ii) (CD74 antigen) (11.6x), follistatin-like 3 (5.3x), serpin peptidase inhibitor, clade G, member 1 (3.5x) and Eph receptor A2 (3.3x). Similar to the down-regulated set of genes, the EnrichR platform indicated glaucomatous/neurodegenerative transcriptomic profile as the main systemic syndrome or disease (Table S2). Further, PPI analysis showed that out of the up-regulated genes, 64 had at least one interaction, that form a network of 129 interactions with a medium confidence score of 0.4 (Figure 3). Followed by enrichment analysis, the main Biological Processes identified were: Immune System, Regulation of Cell Proliferation and Innate Immune Response ( Figure 3A). Among Molecular Functions, Binding, Protein Binding and Ion Binding categories presented the most pronounced interactions, with 33, 25 and 24 altered genes, respectively ( Figure 3B). As for Cellular Components, interactions were found in Extracellular Space, Extracellular Region and Extracellular Region Part categories ( Figure 3C). Both Molecular Function and Cellular Components are consistent with immune system signaling, the main Biological Process identified among the up-regulated genes.

ONC alters the stability of gene expression
In addition to the average expression level, the REV among biological replicates provided an important parameter of comparison between CTR and ONC retinas. We found that the REV profile of genes after ONC was very similar to, but with a slight shift towards higher variability than CTR retinas (Supplementary Figure S1A). REV values were then used to estimate the GES of individual genes in each condition. We found a decrease in the number of very stably expressed genes (REV < 10), from 232 in CTR retinas to 48 after ONC. In order to analyze whether the GES values of individual genes were similar between the experimental groups, the GES value for each spot from CTR retinas was plotted against the GES value for the same spot from ONC retinas (Supplementary Figure S1B). Only a few genes, delimited by the red dotted circles at Supplementary Figure S1B, displayed clearly distinct GES values between the two conditions. Panther was used to analyze enrichment of the genes within the red dotted circles (Supplementary Figure  S1B) and indicated that most of these genes are related to eye structure and development ( Table 5). This result is in agreement with the structural changes in the retina after ONC as expressed by RGC degeneration (Figure 1). Particularly, crystallin genes underwent a dramatic change in their GES, from stable expression in CTR sample (GES > 60) to a very unstable expression (GES < 6), indicating active participation in the retinal changes at 14 days after ONC.

Pathway enrichment analysis in ONC
We used the PathVisio3 platform to identify which biological pathways were the most affected among genes upregulated by ONC. "Complement Activation, Classical Pathway" and "Complement and Coagulation Cascades" pathways appeared with the highest Z scores, followed by "Delta-Notch Signaling Pathway" and "Notch Signaling Pathways", all related with immune response ( Table 6). Additionally, "Oxidative Stress" and "Kit Receptor Signaling Pathway" were found significantly altered in ONC retinas when compared with controls. In addition, many among the top 50 up-regulated genes belong to the Complement Cascade and Notch Signaling Pathways ( Table 6). One representant from "Complement Cascade" (C3), "Oxidative Stress" (Cyba) and "Kit Receptor Signaling" (Fyn) pathways were further validated by RT-qPCR ( Table 3). We plotted the results obtained with PathVisio3 into pathway templates from either the Kyoto Encyclopedia of Genes and Genomes (KEGG) or WikiPathway to visualize significantly altered genes, both up-and down-regulated, in all pathways. We found more significantly up-regulated genes (red), with higher foldchange within the "Complement Cascade Pathway", "Notch-" and "Delta-Notch" Signaling Pathways (Figures 4-6) than within Oxidative Stress and Kit Receptor Signaling Pathway (Figures 7-8).
In the "Complement Cascade Pathway", we found no down-regulated genes and five significantly up-regulated genes (red) belonging to both "Classical" and to the "Alternative Pathways": C3 (2.17x), Cfb (2.26x), C1s (2.43x), C4 (2.41x) and Serping1 (3.49x) (Figure 4). Two other genes belonging to the "Membrane Attack Complex" were also upregulated: vitronectin (1.83x) and clusterin (1.71x) (Figure 4A-B). Next, we examined the number of coordination interactions (synergistic and antagonistic) between gene pairs with significant pairwise Pearson correlation coefficients from the "Complement Cascade Pathway". This analysis allows an estimate of the degree of network interlinkage within the Pathway. Expression coordination of two genes may be an indication that the encoded proteins are linked in one or more functional pathways, and therefore their expression levels should be in a particular proportion [60]. Among the CTR retinas we found 29 positive coordinations ("synergisms") represented by magenta lines in Figure 4A. After the ONC, the number of synergisms raised to 36, and 2 independent coordinations were found (black dotted lines, Figure  4B). For all pathways examined in PathVisio3, no negative coordination ("antagonism") was found, while an increase in positive coordinations was observed following ONC.
We also examined the coordination profiles of all genes from the "Complement Cascade Pathway" against all other genes in our microarray datasets in each experimental condition. Such analysis identifies gene pairs with either similar or opposite coordination profiles, based on the overlap of their synergistically, antagonistically or independently expressed partners. We found that the synergistic coordination among the "Complement Cascade Pathway" genes increased from 21% in CTR retina to 36% after ONC, indicating a stronger alignment of the genes to the network, whereas there were no opposite coordination profiles. As an example, the coordination profiles of the pairs C3-C1qb, Serping1-C1qb and Serping1-Vtn are shown in Figure 4C-E. In CTR retinas, the coordination profiles of the paired genes were neutral (black dots), whereas ONC retinas showed a robust increase in similarity (oranges dots), with R 2 values of at least 0.9. The "Delta-Notch" and "Notch Signaling Pathways" were also significantly altered by ONC, with five genes upregulated when compared to CTR retinas: Cntf (1.55x), Notch1 (1.94x), Fhl1 (1.53x), Smad1 (2.12x) and Dtx2 (1.51x) ( Table  6). The WikiPathway template was used to visualize the altered expression in Notch1 and Deltex2 genes ( Figure 5A-B) and the number of coordination interactions between pairs of genes by Pearson's correlations was represented by magenta (synergism) or black dotted (independent coordination) arrows. The number of synergistically expressed genepairs increased from 20 in the CTR condition to 52 after ONC (Figure 5A-B). No independently expressed genes were observed in the ONC dataset. To further assess whether genes in the "Notch Signaling Pathway" are committed to this network, we examined the coordination profiles of genes from this pathway against all other genes in our microarray datasets in each experimental condition. We found an increase of similar coordination from 22% in CTR retinas to 28% after ONC. As an example, two neutral profiles (black dots) in the CTR retinas for the Dll1-Notch1, Dvl3-Dtx3 and Notch4-Notch1 pairs (R 2 =0.08, 0.18 and 0.63, respectively), changed to similar coordination (orange dots) (R 2 >0.9) after ONC (Figure 5C-E). These results indicate that ONC induces strong gene networking in the "Notch Signaling Pathway". We performed similar analysis of Pearson's correlations for genes present in the Delta-Notch Signaling Pathway ( Figure  S2). Using KEGG templates we showed that ONC led to up-regulation of four transcripts (Notch1, Cntf, Smad1 and Fhl1) and no significant down-regulated genes ( Figure S2). Whereas the CTR group had 38 synergistic correlations, ONC led to an increase in the number of such interactions to 76 ( Figure S2). As previously stated, another two pathways were altered by ONC with significant Z scores as indicated by Path-Visio analyses: "Oxidative Stress" and "Kit Receptor Pathway". Two genes were significantly up-regulated in the "Oxidative Stress Pathway" (Nfix and Cyba) and no down-regulated genes were altered by ONC (Figure S3). In the CTR retinas, Nfix showed synergistic interactions with Xdh and Maoa, which were abrogated in the ONC, where Nfix had synergism only with Cyba ( Figure S3A, Table 3). The overall number of synergisms in this Pathway was increased from 57 in the CTR to 82 in the ONC groups.
Finally, we analyzed the Kit Receptor Signaling Pathway ( Figure S4A). Cblb, Fyn and Tec were up-regulated by ONC and no significantly down-regulated genes were verified in this pathway ( Figure S4B, Table 3). The number of synergistic correlations was also increased by ONC in this Pathway, 43 in CTR and 53 in ONC (Figure S4). We chose one representative gene from each pathway (Complement, Oxidative Stress and Kit Receptor Signaling) to perform a RT-qPCR validation in a new set of retinas ( Table 3) and found that C3, Cyba and Fyn were consistently up-regulated in rat retinas after ONC.   In this study, we conducted a transcriptomic analysis of rat retina at 14 days after ONC, a time-point when the population of RGCs was reduced to only 40% of the original population, as compared to controls. We identified differentially expressed genes, at a total of 127 genes significantly up-regulated, and 66 down-regulated. We examined changes in GES and found a reduction in the number of very stably expressed genes from 336 to 124. Using bioinformatic tools to propose protein interaction (STRING) and main altered molecular pathways (PathVisio3), we found that the Complement Cascade and Notch Signaling Pathway were the main affected pathways. Although protein-protein interaction determined with STRING (based on data mining in genome wide datasets) may be a good predictor of real interactions (Szklarczyk et al., 2019), we deepened our analyses with additional bioinformatics and mathematical modeling of gene networks of ONC retinas, which has been shown by our group to be a consistent method for such predictions. Coordination analysis for both pathways showed an increase in the number of synergistically expressed gene pairs and of genes with similar coordination profile.
Our analysis was carried out on whole-retina extracts, and the transcriptome profiles analyzed represent the adaptations of the retinal tissue to the RGC degeneration promoted by the ONC. Although the observed regulation is triggered by the ONC, which leads to RGC death, changes at 14 days after injury likely involve the activation of the retinal glia (.i.e. astrocytes, Müller cells, and microglial cells), and also other retinal neurons besides RGCs, as signs of retinal remodeling. These changes may impact potential treatment strategies to ameliorate the secondary degeneration associated with CNS insults [41].
Unique among transcriptomic studies of the retina is that beyond the regulation of gene expression, we also disclosed REV, changes in GES and coordination following ONC. GES score ranks the genes according to the strength of the cellular homeostatic mechanisms to keep their expression variability within narrow limits. GES, that includes the most stably expressed genes in the 100th percentile, actually identifies the genes most critical for the survival and phenotypic expression of the tissue. By contrast, the most unstably expressed genes (GES in the first percentile) most likely empower the cell to respond/adapt to the environmental fluctuations [52][53]. Moreover, the expression coordination analysis refines the functional pathways, by the identification of gene pairs whose expression fluctuations among biological replicas are positively (synergistically), negatively (antagonistically) or neutrally (independently) correlated. In previous papers [50,[52][53], we have speculated that genes whose encoded proteins are linked in functional pathways should coordinate their expression to optimize a kind of "transcriptomic stoichiometry". Changes in the coordination profiles indicate remodeling of the functional pathways [44; 49-53; 59-61].
Expression of individual genes depends on local conditions that, although similar, are not identical among biological replicas. We assume that expression of key genes is kept by the cellular homeostatic mechanisms within narrow intervals while that of non-key genes is less restrained to readily adapt to environmental changes. In our results we found a higher variability in the REV profile of genes with ONC ( Figure S1A), as expected in pathological conditions, suggesting that CTR mechanisms are also affected.
Thus, by analyzing GES changes following the ONC, we found clues about retina cells changing priorities in controlling the expression level of certain genes. Our results identified, at 14 days after ONC, that only a few genes presented an altered profile in their GES, with enrichment in genes of eye structure ( Table 7). The crystallin genes were found to have relevant changes in their GES, thus acquiring a more unstable expression profile. Crystallins have chaperone-like properties involved in an increase of cellular resistance to stress-induced apoptosis [62]. Up-regulation and down-regulation of crystallin expression has been reported upon different cellular stresses, and those alterations are viewed as a cellular response against environmental and metabolic insults [63][64][65]. Crystallins have been related with both RGCs [66] and amacrine cells [67], supporting its active participation in the remodeling of retinal tissue after injury, possibly through a decrease in resistance to stress, as well as leading to tissue stabilization. Curiously, two other genes presented a drastic change in GES: Adad2 (GES (CTR) = 79 to GES (ONC) = 5), a gene related with RNA editing [68], and RGD1308160 (GES (CTR) = 61 to GES (ONC) = 10), a gene that, according to BlastP, has a loose homology with the C. elegans gene smc-4, involved in chromosome stabilization [69]. These data suggest that both such mechanisms are very active at 14 days after ONC.
Typical RGC-expressed genes, underrepresented in ONC retinas compared to CTR ( Figure 1D-E), reflect RGC death after optic nerve damage. Among the 66 down-regulated genes, 20% have been described as expressed by ganglion cells [5-6; 21]. Enrichment analysis of the down-regulated genes showed cytoskeletal organization as the main biological process, with genes of intermediary filaments and of proteins that interact with them. Neurofilaments are particularly abundant in axons, and essential for their growth, maintenance of caliber, transmission of electrical impulses, velocity of impulse conduction, as well as regulation of enzymes function and the structure of linker proteins [70][71][72].
Additional biological roles have been attributed to neurofilaments such as synaptic plasticity. Marked decrease in the expression of the neurofilament genes was first shown in brain tissue from Alzheimer's disease patients as compared to CTRs [73]. It is generally accepted that CNS disorders are divided between early and late changes, associated with dysfunction of neuronal activity caused by perturbations of synapses [74]. Neural reorganization after stroke is initiated as cellular reactions to degeneration. While neurons die in an ischemic region, axons and synapses degenerate in further brain regions, promoting regenerative responses that lead the growth of new connections among surviving neurons [75]. RGC dendritic arbors also showed significant reduction after 2 weeks of elevation of intraocular pressure, a common model of chronic glaucoma [76]. In a primate glaucoma model, a correlation was established between abnormalities of parasol RGC dendrite morphology and function [77]. Morphological changes of RGC dendritic arbors have been detected before axon thinning or soma shrinkage, suggesting that dendritic abnormalities may precede degeneration of other ganglion cell domains. Alterations of neurofilament expression in various neurological disorders not only underlie axonopathy, but also synaptic remodeling because of its vital roles in synaptic function [78][79]. Since at 14 days after ONC the retinal tissue still harbors 40% of the RGCs, the decreased expression of genes related with cytoskeletal organization, as intermediary filaments, may be representative of retinal remodeling, not restricted to the injured axon but also at the level of dendrites. This was further reinforced by the predicted protein-protein interactions (STRING) analysis that showed "Axon" and "neuron projection" cellular components were significantly altered. PPI along with Pearson correlation analyses performed herein take advantage of the concept that expression of individual genes is linked to each other so that the protein amounts would respect the "stoichiometry" of the biochemical reactions [80].
Many components of the innate immune system are expressed intrinsically by retinal cells (Retinal Database, GeneNetwork.org), and innate immune networks are activated by various types of insults to the retina, such as ONC, partial transection of the optic nerve, DBA2J mouse models of glaucoma, among others [6; 10; 15; 21; 81-83]. In the present study, enrichment analysis disclosed the Complement Cascade as the most prominent biological process among differentially expressed genes after ONC (Table 8), which was further confirmed by predicted PPI as "innate immune response" and "immune system process". Complement activation was shown to be up-regulated in the retina as early as 2 days after ONC [84]. Such activation, specifically including C1q, C3, and Cfi [20], suggests that the innate immune system plays an important role in retinal immune response through glia cells, and specifically in the response of the retina to injuries. It has been suggested that, at an early stage of a mouse glaucoma model, the Complement Cascade mediates synapse loss, in particular C3 [85][86][87][88]. At late-stage glaucoma, the Complement Cascade is thought to play its traditional role of opsonizing and removing the cellular debris from widespread RGC death [89]. Recently, it was described that C1q marks a subset of RGC in the embryonic retina and that knockout of C3 receptor, which is only found in microglia, resulted in increased RGC numbers, indicating a direct relationship with microglia-mediated RGC elimination [90].
Our data showed that ONC leads to an increase in synergistic coordination, most evident in the Classic Pathway of Complement Cascade, which indicates a stronger gene network in this pathway. It is particularly interesting that the serping1 gene showed increased synergistic coordination at 14 days after ONC (Figure 4). Serping1 is a C1 inhibitor known to block the initiation of the Complement Cascade through the classical pathway, therefore limiting inflammation. It is indirectly involved in the production of bradykinin, a peptide that promotes inflammation by increasing the permeability of blood vessel walls and also inhibits leukocyte recruitment into ischemic cardiac tissue, an effect independent of protease inhibition [91]. Therefore, the up-regulation of serping1 in ONC retinas may help CTR the inflammatory response after injury, both by modulation of the Complement Cascade, as well as preventing infiltration of peripheral inflammatory cells through the blood-retina barrier.
Our enrichment analysis also showed up-regulation of the Notch Signaling Pathway at 14 days after ONC. Notch signaling is known to play important roles during retinal development, in both RGC differentiation [92] and proliferation of Müller glia in acutely damaged chick retina [93]. Notch is expressed in most Müller glia cells at low levels in undamaged retina, and blockade of Notch activity prior to damage is protective to amacrine, bipolar and ganglion cells [94]. Also, overexpression of Notch1 has been observed in several neurodegenerative diseases [95][96]. Besides Notch, the deltex2 gene was also up-regulated (Figure 5). A recent study showed that Dtx together with its interacting partner, Hrp48, down-regulated Notch signaling and induced cell death during Drosophila development [97]. At 14 days after ONC, the up-regulation of both Notch and deltex2 suggests that Notch signaling may act upon retinal Müller cells via a non-canonical pathway, thereby positively influencing RGC degeneration. Notch signaling mediates many different intercellular communication events, influencing retinal Muller cells, ganglion cells and microglial cells. It is known that Notch-1 signaling in neurons during ischemic conditions may enhance apoptotic signaling cascades, while activation of notch-1 in microglial cells and leukocytes may exacerbate inflammatory processes that contribute to neuronal death [98]. Also, several studies found that Notch signaling serves important functions in the regulation of neurite outgrowth and maintenance. Both Sestan et al., [99] and Berezovska et al., [100] found that Notch signaling influence dendritic morphology: while activation inhibits neurite outgrowth or causes their retraction, inhibition promotes neurite extension. These data, together with our results suggest that Notch signaling may confer a positive effect on RGC degeneration through an exacerbation of inflammatory processes as well as promoting retinal tissue remodeling, by altering dendritic morphology.
Oxidative Stress Pathway was also found to be enriched in our analyses, and this showed to occur due to Cyba and Nfix up-regulation in ONC retinas. In the CNS, normal NAD(P)H oxidase function appears to be required for processes including neuronal signaling, but overproduction of ROS contributes to neurotoxicity and neurodegeneration. Conversely, RGCs, were found to express NAD(P)H oxidase subunits under physiologic conditions and after ischemia [101]. NOX proteins are homologs of NAD(P)H oxidases that generate reactive oxygen species [102]. NOX proteins interacts with and is stabilized by Cyba (or p22phox) in intracellular vesicles, which, when activated by cytosolic p47phox, form a protein complex that transports electrons from cytoplasmic NADPH to extracellular oxygen to generate superoxide (O2−). Cyba was also found to be increased in a rat model of diabetic retinopathy [103] and in retinal pigmented epithelial cells under Angiotensin-II-induced oxidative stress [104].
NFI transcription factors a/b/x were shown to be temporally expressed during retinal development by retinal progenitor cells, bipolar cells and Müller glia [105]. Moreover, Nfix is known to have a role on CNS development and Nfixnull mice display major brain malformations in the cortex, cerebellum and hippocampus [106, for reviews]. In the skeletal muscle, Nfix can also regulate myogenesis [107,108], however deletion of Nfix leads to decreased muscle regeneration and has been implicated in the pathogenesis of muscular dystrophies [109]. The knock-down of Nfix in a myoblast cell line was protective against oxidative stress [110], thus indicating that this transcription factor may have a versatile role in health and disease. In the context of retinal degeneration, the up-regulation of Cyba and Nfix transcripts reinforces the hypothesis that oxidative stress takes place and accelerates RGC death in ONC models [111].
Oxidative stress and inflammation are closely related pathophysiological processes, one of which can be easily induced by another. Thus, both processes are simultaneously found in many chronical pathological conditions, such as glaucoma [112]. The immune activity is a necessary intrinsic response to promote tissue cleaning, healing, and regeneration process after injury [113]. However, unbalanced or sustained immune response, turns the beneficial potential into potentiation of the injury process. The mechanical injury together with hypoxia, both promoted by the ONC could be the triggers for the activation of inflammation, oxidative stress and metabolism alterations. Notch-1, c-Kit and complement signaling are known as linkers between inflammation, metabolism, oxidative stress and dendritic plasticity, acting into RGC and Müller cells [114][115][116][117][118][119]. Our results support the hypothesis that not only the events intrinsic to RGCs, but also environmental signals from other cells, such as Müller cells, are critical for cell death stimuli, and RGC-glia interactions are critically important for different aspects of RGC neurodegeneration. Improvement of the current understanding of the role of RGC-glia interaction, mitochondria, and the immune system in neurodegeneration will potentiate the development of innovative treatment methods for neuroprotection.

Conclusions
In summary, our data indicate that ONC induces a robust, long-lasting alteration of at least two main pathways (Complement Cascade and Notch Signaling Pathway), including a significant increase in the expression and coordination of inflammation-related genes. The coordination analyses of genes within each pathway and with the entire transcriptome suggest that such biological pathways may act after ONC as major players in changes of gene expression that lead to both degeneration and remodeling of retinal tissue. Further functional studies on the effects of these genes are warranted to clarify their roles in molecular mechanisms within the retinal tissue after damage to the optic nerve.

Institutional Review Board Statement:
Experiments were performed on eight weeks old female Lister Hooded rats (n=43), housed in plastic cages, on a 12-hour light-dark cycle, with water and food ad libitum. All experiments were carried out in accord with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research, and the protocol approved by the Institutional Animal Experimentation Ethics Committee (UFRJ#01200.001568/2013-87). Figure S4. ONC affects the Kit Receptor Signaling Pathway. This pathway was obtained from the WikiPathway platform and used as a template to highlight the effect of the ONC in the retinas. Magenta lanes represent synergistic correlations between pairs of genes. ONC induced an increase in the number of synergistic correlations. Yellow boxes indicate genes that showed no alteration in ONC versus CTR retinas, whereas red-filed boxes represent gene up-regulation (>1.5-fold change, p <0.05).