Optic Nerve Crush Activates the Genomic Fabrics of the Complement Cascade and Delta-Notch Signaling

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 applied systems biology approaches to better understand the molecular mechanisms related with the retinal remodeling after induction of RGC degeneration. 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. Enrichment analysis showed that Complement cascade and Notch signaling pathway were the main affected pathways after ONC. To expand our studies of these two pathways, we examined the coordination of gene expressions within each pathway and with the entire transcriptome. ONC increased the number of synergistically coordinated pairs of genes and the number of similar profiles. This study provided novel findings beyond the regulation of individual gene expression and disclosed changes in the control of expression by Complement cascade and Notch signaling functional pathways important for both RGC degeneration and remodeling of the retinal tissue after ONC.


Animal handling
27 eight weeks old female Lister Hooded rats were 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 (#01200.001568/2013-87).

Optic nerve crush
The rats were divided into two experimental groups: optic nerve crush (n=16, hereafter denoted by ONC) and control, sham operated, (n=11, 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.

Microarray gene expression
We applied an optimized protocol [47] 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 [48]. 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. For each experimental group (CTR or ONC), four separated retinas were profiled individually. RINs of the selected OPC 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. With each array, samples of 825 ng of Cy3 labeled RNA from the retina of a rat subjected to ONC and 825 ng Cy5-labeled RNA from the retina of a CTR rat were co-hybridized for 17 h at 65°C. 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% [49]. 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 Relative Expression Variability (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 [50].
The genes were then ordered according to decreasing variability, such that the first percentile (or gene expression stability score, GES < 1) contained the most unstably expressed and the 100th percentile (GES > 99) the most stably expressed genes [51].
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 [52]. This composite criterion, in which the cut-offs were computed with Bonferroni-type corrections applied to the redundancy groups [53][54], 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/, [55]), 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 cutoff of 0.4.

Analyses of enrichment
Pathvisio3 (www.pathvisio.org, [56]) 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 Optic Nerve Crush. 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; [57]).
The EnrichR platform (http://amp.pharm.mssm.edu/Enrichr/, [58]) 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 GEO (Gene Expression Omnibus) database.

Coordination of expression
We used the variation in gene expression among multiple biological replicas to calculate the pair-wise Pearson correlation coefficients ρ 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 [53][54].

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 SuperScript™ III Reverse Transcriptase (Invitrogen, 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 control 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
As expected, ONC induced robust loss of retinal ganglion cells (RGCs). At 14 days after the procedure, control retinas displayed an average of 1.780±134 cells/mm 2 , retrogradely labeled with DiI, whereas ONC retinas showed a reduction of approximately 60% of RGCs, to 768±65 cells/mm 2 (Figure 1a  Pubmed 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.

Transcriptomic alterations in rat retina after optic nerve crush
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 (GEOaccession number GSE133563). Among 17,657 unigenes quantified in our microarrays, expression of 1707 (9.7%%) was significantly altered in ONC samples with respect to CTR, of which 1167 (6.6%) were up-and 540 (3.1%) were down-regulated. Table S1 in the Appendix presents the most down-regulated 50 genes and Table S2 the most up-regulated 50 genes. Quantitative RT-PCR validated two of the highest up-regulated genes (Cd74 and C3), and their respective fold-change increase in expression, compared with control retinas, were consistent with the results of the arrays (Table S3).

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 (Figure2(a)), a major component of the optic nerve axons. Among Molecular Functions, ONC affected networks of Protein Binding, Binding and Protein Heterodimerization Activity (Figure 2

Protein-protein interaction encoded by the most up-regulated genes
Consistent with previous studies [15-18; 29], the group of 50 most up-regulated genes (Table S2) include 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 S5). 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 43

ONC alters the stability of gene expression
In addition to the average expression level, the relative expression variability (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 ( Figure S1(a) in the Appendix). REV values were then used to estimate the stability of expression (GES) of individual genes in each condition. We found a decrease in the number of very stably expressed genes (REV < 10), from 232 in control retinas to 48 after optic nerve crush. In order to analyze whether the GES values of individual genes were similar between the experimental groups, the GES value for each spot from control retinas was plotted against the GES value for the same spot from ONC retinas ( Figure S1(b)). Only a few genes, delimited by the red dotted circles at Figure S1(b), displayed clearly distinct GES values between the two conditions. Panther was used to analyze enrichment of the genes within the red dotted circles ( Figure S1(b)) and indicated that most of these genes are related to eye structure and development (Table S8). This result is in agreement with the structural changes in the retina after ONC as expressed by RGC degeneration (Figure 1).

Complement cascade and Delta/Notch Pathways are targeted by ONC
We used the PathVisio3 platform to identify which biological pathways were the most affected among genes up-regulated by ONC. "Complement Cascade Pathway" appeared twice and had the highest Z scores, as well as "Notch-" and "Delta-Notch Signaling Pathways", all related with immune response (Table S6). In addition, many among the top 50 up-regulated genes belong to the Complement cascade and Notch signaling pathways (Table S6). 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 fold-change within the "Complement Cascade Pathway", "Notch-" and "Delta-Notch" signaling pathways (Figure 4, Figure 5 and Figure 6) than within Oxidative Stress and Kit Receptor Signaling Pathway (Figures 7  and 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 6). Two other genes belonging to the "Membrane Attack Complex" were also up-regulated: vitronectin (1.83x) and clusterin (1.71x) ( Figure 6A-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 control retinas we found 29 positive coordinations ("synergisms") represented by magenta lines in Figure 6A. After the ONC, the number of synergisms raised to 36, and 2 independent coordinations were found (black dotted lines, Figure 6B). For all pathways examined in PathVisio3, no negative coordination ("antagonism") was found, while an increase in positive coordinations was observed following ONC ( Figure S1).
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 control 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 6C-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 R2 values of at least 0.9. The "Delta-Notch" and "Notch Signaling Pathways" were also significantly altered by ONC, with five genes up-regulated when compared to control 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 7A-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 gene-pairs increased from 20 in the control condition to 52 after ONC ( Figure 7A-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 control retinas for the Dll1-Notch1, Dvl3-Dtx3 and Notch4-Notch1 pairs (R2=0.08, 0.18 and 0.63, respectively), changed to similar coordination (orange dots) (R2>0.9) after ONC ( Figure 5C-E). These results indicate that ONC induces strong gene networking in the "Notch Signaling pathway". down-regulated. We examined changes in gene expression stability and found a reduction in the number of very stably expressed genes from 336 to 124. Using bioinformatic tools to examine interacting proteins (STRING) and main altered pathways (Pathvisio3), we found that the Complement cascade and Notch signaling pathway were the main affected pathways. 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 relative late stages, at 14 days after injury, likely involves the activation of the retinal glia (.i.e. astrocytes, Muller cells, and microglial cells), and also other retinal neurons besides RGC, as signs of retinal remodeling. These changes may impact potential treatment strategies to ameliorate the secondary degeneration associated with CNS insults [42].
Unique among transcriptomic studies of the retina is that beyond the regulation of gene expression, we also disclosed relative expression variability, changes in gene expression stability and coordination following ONC. Gene Expression Stability score (GES) 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 cell. By contrast, the most unstably expressed genes (GES in the first percentile) most likely empower the cell to respond/adapt to the environmental fluctuations [53][54]. 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 [51,[53][54], 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 [45; 50-54; 60-62].
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 6A), as expected in pathological conditions, suggesting that control mechanisms are also affected.
Thus, by analyzing GES changes following the crush of the optic nerve, 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 (Supplementary Table 1). Particularly, crystallin genes underwent a dramatic change in their GES, from stable expression in control sample (GES > 60) to a very unstable expression (GES < 6), indicating active participation in the retinal changes at 14 days after ONC. The crystallin genes have chaperone-like properties involved in an increase of cellular resistance to stress-induced apoptosis [63]. 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 [64][65][66]. Crystallins have been related with both retinal ganglion cells [67] and amacrine cells [68], 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 [69], 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 [70]. 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 control ( Figure  3A), reflect RGC death after optic nerve damage. Among the 66 downregulated 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 [71][72][73].
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 controls [74]. 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 [75]. 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 [76]. Retinal ganglion cell dendritic arbors also showed significant reduction after 2 weeks of elevation of intraocular pressure, a common model of chronic glaucoma [77]. In a primate glaucoma model, a correlation was established between abnormalities of parasol RGC dendrite morphology and function [78]. 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 [79][80]. Since at 14 days after ONC the retinal tissue still harbors 40% of the RGC, 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.
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 6)) and was 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 the immune-privileged environment of the retina 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]. To date, inflammation has been considered secondary to the tissue damage, and is now thought to be part of a protective response of the immune system. It is also believed that an intermediate stage exists between basal homeostatic conditions and true inflammation [91].
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 7). 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 [92]. Therefore, the up-regulation of serping1 in ONC retinas may help control 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 retinal ganglion cell differentiation [93] and proliferation of Müller glia in acutely damaged chick retina [94]. 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 [95]. Also, overexpression of notch1 has been observed in several neurodegenerative diseases [96][97]. Besides notch, the deltex2 gene was also up regulated (Figure 8). A recent study showed that Dtx together with its interacting partner, Hrp48, down-regulated Notch signaling and induced cell death during Drosophila development [98]. 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.

Conclusions
In summary, our data indicate that ONC induces a robust, long-lasting alteration of two main pathways, 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.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Figure S1: ONC affect the expression variability; Table S1: Previously described RGC markers downregulated in ONC vs CTR, Table  S2: Enrichment analysis with genes had a dramatic change in their expression stability after ONC.