Meta-analysis of transcriptomic changes in optic nerve injury and neurodegenerative models reveals a fundamental response to injury throughout the central nervous system.

Purpose: Injury to the central nervous system (CNS) leads to transcriptional changes that effect tissue function and govern the process of neurodegeneration. Numerous microarray and RNA-Seq studies have been performed to identify these transcriptional changes in the retina following optic nerve injury and elsewhere in the CNS following a variety of insults. We reasoned that conserved transcriptional changes between injury paradigms would be important contributors to the neurodegenerative process. Therefore, we compared the expression results from heterogeneous studies of optic nerve injury and neurodegenerative models.
Methods: Expression data was collected from the Gene Expression Omnibus. A uniform method for normalizing expression data and detecting differentially expressed (DE) genes was used to compare the transcriptomes from models of acute optic nerve injury (AONI), chronic optic nerve injury (CONI) and brain neurodegeneration. DE genes were split into genes that were more or less prevalent in the injured condition than the control condition (enriched and depleted, respectively) and transformed into their human orthologs so that transcriptomes from different species could be compared. Biologic significance of shared genes was assessed by analyzing lists of shared genes for gene ontology (GO) term over-representation and for representation in KEGG pathways.
Results: There was significant overlap of enriched DE genes between transcriptomes of AONI, CONI and neurodegeneration studies even though the overall concordance between datasets was low. The depleted DE genes identified between AONI and CONI models were significantly overlapping, but this significance did not extend to comparisons between optic nerve injury models and neurodegeneration studies. The GO terms overrepresented among the enriched genes shared between AONI, CONI and neurodegeneration studies were related to innate immune processes like the complement system and interferon signaling. KEGG pathway analysis revealed that transcriptional alteration between JAK-STAT, PI3K-AKT and TNF signaling, among others, were conserved between all models that were analyzed. Conclusions: There is a conserved transcriptional response to injury in the CNS. This transcriptional response is driven by the activation of the innate immune system and several regulatory pathways. Understanding the cellular origin of these pathways and the pathological consequences of their activation is essential for understanding and treating neurodegenerative disease.

Introduction

Understanding a tissue’s transcriptional response to injury is a helpful tool to understanding the pathology of any disease. Following optic nerve injury (ONI), transcriptional changes occur that effect the course of injury response and the level to which the tissue is functionally compromised [1,2]. Numerous studies have been performed to quantify transcriptional changes that result from acute or chronic optic nerve injury (AONI and CONI, respectively) in a variety of animal models ranging from zebrafish to non-human primates. AONI models, like optic nerve crush and axotomy, damage the axons of the retinal ganglion cells (RGCs) by pinching or transecting the optic nerve behind the globe of the eye. This axonal damage leads to the rapid degeneration and apoptosis of the RGCs that begins within hours of damage, achieves a maximal rate of cell loss between 5 and 7 days and ends with near complete ablation of the RGC population by 21 days post-injury [3-5].

Among the CONI models used for transcriptome analyses, the most commonly used are models of experimental or inherited glaucoma. Of these, the DBA/2J mouse has been widely studied and characterized. This mouse spontaneously develops elevated intraocular pressure, a risk factor for human glaucoma, as it ages, which leads to the degeneration of the optic nerve and the loss of RGCs starting between 6 and 8 months of age and completing in most animals by 12 months of age [6]. While acute models achieve a highly synchronous loss of RGCs within a narrow time-frame, the DBA/2J chronic model presents with pathology that is more clinically similar to human disease, but that is highly variable from eye to eye and which occurs over a much longer period of time [6].

While optic nerve injury models result in RGC death, many other cell-types are also involved including retinal and optic nerve glia [7]. In analyses of whole tissues, changes in these other cell types are expected to contribute to the overall transcriptomic changes.

Results from microarray and RNA-Seq experiments are often confounding and yield results that are shockingly dissimilar [8]. A meta-analysis comparing the transcriptomic alterations observed in ONI studies is necessary to separate conserved transcriptional changes from artifactual ones that may occur randomly in a study. Biologic processes that are changed in multiple models and multiple species are less likely to be artifacts of a model and more likely to be important global modifiers of disease etiology and progression. For this reason, we sought to compare data from studies of optic nerve injury and brain neurodegeneration that varied in multiple ways including species analyzed, mechanism of injury and method of quantification of transcript abundance.

Using the Gene Expression Omnibus (GEO), transcriptional data was gathered from experiments analyzing transcriptional changes in models of AONI and CONI. We normalized the raw expression data and identified differentially expressed (DE) genes that were either enriched or depleted in the injured/diseased condition relative to the control. To determine if there was a relationship between the expression patterns of each data set, hierarchical clustering of whole transcriptomes was performed on the expression data from mouse experiments that used Affymetrix chips. Lists of DE genes were compared between studies of AONI and CONI. We hypothesized that comparing transcriptional changes between heterogeneous models would reveal genes that were fundamental to the retinal response to injury. This analysis was extended to compare transcriptional changes from ONI models to mouse models of other neurodegenerative conditions. We use Gene Ontology (GO) term over-representation analysis and KEGG pathways to identify biologic processes and pathways that are conserved in the transcriptional changes that follow injury to various areas of the central nervous system (CNS).

Here we report that comparing transcriptional changes resulting from diverse AONI, CONI and neurodegeneration data sets revealed great diversity of non-overlapping DE genes, but within these data sets there was a conserved transcriptional response between models. Lists of enriched DE genes significantly overlapped with one another regardless of species, injury, or tissue. Lists of depleted DE genes from studies examining the same tissue were significantly overlapping but inter-tissue comparisons did not significantly overlap. GO terms associated with commonly occurring enriched genes were often related to the innate immune response, indicating that this response is conserved throughout the CNS response to a host of injuries.

Methods

Summary of optic nerve injury studies

All expression and platform data used in analysis were retrieved from the GEO database. We sought to include as much data as possible while allowing the data manipulation to remain as uniform and simple as possible. For this reason, the majority of the included data are from Affymetrix array chips. The exceptions are data from a Sentrix chip and from studies that used RNA-Seq. These data were pre-normalized by the study authors and were easily mapped to their gene symbols. Syntactically, we defined a “study” as the collection of expression data published by a group for a particular experiment. Each study contained one or more “datasets” which were comparisons of gene expression between a control condition and an experimental condition (i.e., a time after acute injury). Each study was identified by the name of the first author of the publication attributed to the data in GEO. In all, a total of 14 studies and 55 data sets were used for analysis Appendix 1 summarizes all the data sets used, including citation of the corresponding GEO reference numbers. Individual descriptions of the studies are:

(i) The Agudo study examined retinal gene expression changes in Rattus norvegicus following optic nerve crush or optic nerve axotomy [9]. This study used the Affymetrix Rat Genome 230 2.0 Array and compared expression changes from injured animals to a control group of naïve uninjured animals.

(ii) The Jiang study explored expression changes in Canis lupus familiaris by comparing retinal gene expression between dogs with inherited glaucoma and healthy dogs [10]. This study used the Affymetrix Canine Genome 2.0 Array and dogs of various breeds and ages.

(iii) The Howell study looked at gene expression changes in both the retina and optic nerve during increasing stages of optic nerve degeneration in the DBA/2J mouse (Mus musculus) model of inherited glaucoma [11]. The Affymetrix Mouse Genome 430 2.0 Array was used for this study and separate data sets were generated for the retina and optic nerve. Expression changes in diseased animals were compared with control animals that had the wild-type Gpnmb gene knocked in, which prevented elevated IOP and glaucoma.

(iv) The Lukas study investigated gene expression changes in mouse (C57BL/6) population enriched for retinal ganglion cells 6 h after optic nerve crush [5]. A total of 6000 cells were harvested from the ganglion cell layer by laser captured from frozen sections. RNA from captured cells was amplified as cDNA and screened using the Affymetrix Mouse Genome 430 2.0 Array. Cells from the contralateral eye were used as a control.

(v) The McCurley study explored expression changes in Danio rerio at different times following optic nerve crush [12]. The Affymetrix Zebrafish Genome Array was used and fish that had undergone a sham operation were used as a control.

(vi) The Qu study examined expression changes in the optic nerve of the mouse (C57BL/6) at various times following optic nerve crush [13]. The Affymetrix Mouse Genome 430A 2.0 Array was used and contralateral eyes served as controls.

(vii) The Sharma study looked at expression changes in the retina and optic nerve of the mouse (BALB/c) at various times following optic nerve crush [14]. This study used the Affymetrix Mouse Gene 1.0 ST Array and used the contralateral eye as a control.

(viii) The Steele study compared the retinal expression changes between 8 month old DBA/2J mice and 3 month old DBA/2J mice [15]. They used the Affymetrix Mouse Genome 430 2.0 Array.

(ix) The Templeton study compared the retinal transcriptional response to crush between C57BL/6J and DBA/2J mice [16]. This study used the Sentrix Mouse-6 Expression BeadChip. Retinas from uncrushed animals served as controls.

(x) The Yasuda study used RNA-Seq to quantify gene expression changes in the mouse (C57BL/6) retina two days after optic nerve crush injury [17]. Mice that underwent a sham procedure were used as controls.

Summary of neurodegenerative studies

We collected data from four studies modeling three different neurodegenerative diseases of the CNS.

(xi) The Ferraiuolo study examined expression changes in laser-captured spinal cord neurons in SOD1 G93A C57BL/6J mice of different ages [18]. The SOD1 G93A mouse models Amyotrophic Lateral Sclerosis (ALS) by overexpressing a mutant human SOD1 protein in neurons which leads to age-dependent, progressive phenotypic similarities between the mice and human ALS patients. This study used the Affymetrix Mouse Expression 430A Array and used non-transgenic littermates as a source for control spinal cord neurons.

(xii) The Gjoneska study examined expression changes in mouse hippocampus following accumulation of p25 (C57BL/6J) [19]. Accumulation of p25 leads to a disease that models human Alzheimer disease [20]. This study used RNA-Seq to quantify transcriptional changes and used non-transgenic littermates as controls.

(xiii) The Jonas study measured the expression changes in the mouse (C57BL/6) motor and sensory cortices following myelin oligodendrocyte glycoprotein (amino acids 35–55) induced experimental autoimmune encephalomyelitis (EAE). The EAE model is widely used to study multiple sclerosis [21]. This study used the Affymetrix Mouse Genome 430 2.0 Array and used healthy sensory and motor cortices as controls.

(xiv) The Wakutani study used transgenic TgCRND8 mice (C3H/B6) to explore the transcriptional changes associated with the overexpression of amyloid precursor protein (APP) [22]. The TgCRND8 mouse is a transgenic strain that overexpresses mutant human APP which leads to the accumulation of amyloid-β 40 and 42 as the mouse ages [23]. This study used the Affymetrix Mouse Genome 430 2.0 Array to measure gene expression in the forebrain of TgCRND8 mice at various ages. Non-transgenic mice of the same age were used as controls.

Data normalization and fold change calculation

Affymetrix CEL files were normalized using the RMA normalization method from the “affy” package in R version 2.15.0. For RNA-Seq studies, reads per kilobase of transcript per million mapped reads (RPKM) values provided by the authors were Log2 transformed and used to calculate fold change. For each data set, fold change for each gene was calculated by taking the difference between the Log2 expression value of the injured condition and the Log2 expression value of the uninjured control condition.

Differential expression analysis

We used the R (v2.15.0) limma package to identify DE genes [24]. To be considered differentially expressed, a gene needed to have a p value of less than 0.05, after adjusting for multiple tests using the Benjamini-Hochberg method. We split DE genes into two groups, enriched and depleted, based on a positive or a negative fold change, respectively.

Rosetta Stone ortholog table

The Orthology Predictions Search tool on the Human Genome Organization Gene Nomenclature Committee (HGNC) database was used to construct a table that listed each human gene and each gene from different species that is an ortholog of that gene. We used this information to compare genes from mouse, rat, dog and zebrafish that had the same human ortholog(s). This resource was termed the Rosetta Stone ortholog table. Mapping of orthologous genes was conducted using a script written in Python. An example of the process of translation of orthologs from rat and zebrafish data sets to the human designation is given in Figure 1. A URL providing access to the table and analysis script is shown below.

Figure 1

Identifying Orthologous DE Genes Using the Rosetta Stone Ortholog Table. A: A list of DE genes from two rat and one zebrafish data sets being compared. The list has been shortened to facilitate this example of how the orthologs are identified by the common human ortholog. B: The Rosetta Stone Ortholog Table lists all of the orthologs of each human gene for each species in the table. It is important to note that an ortholog from a species may appear multiple times in the Rosetta Stone. If a species has no reported human ortholog, then that index in the table appears as “no orthologue.” C: The lists of DE genes from each data set shown in (A) are now translated into the corresponding human ortholog using an algorithm written in Python and then aligned to show common DE genes in each data set. This translation allows a direct comparison between data sets between species. In this table, nd refers genes that were not detected as DE genes in the data set.

The Rosetta Stone ortholog table mapped genes from dog, mouse, rat and zebrafish that covered 91%, 93%, 88% and 76% of the human genes in the table, respectively. In all comparisons of data sets, regardless if they were cross-species or intraspecies, the genes were first “translated” into the corresponding human orthologs. For dog, mouse, rat and zebrafish respectively, 27%, 34%, 31% and 53% of human genes had multiple orthologous genes in the given species. This potentially inflated numbers of DE genes in the Rosetta Stone relative to the number of DE genes detected for a data set by limma analysis. In figures and tables where the transcriptomes of data sets were compared, multiple orthologs originating from the Rosetta Stone analysis were not pruned from the lists of DE genes for each data set. This was accounted for in Monte Carlo Simulations to determine the statistical significance of overlapping DE genes from compared data sets (see below). Where GO or KEGG analyses were performed, and in the list of highly conserved genes, lists of genes were pruned to account for multiple orthologs. Pruning ensured that the number of human genes submitted for GO or KEGG analysis reflected the number of DE genes that were detected within the data sets.

Pairwise Monte Carlo simulations

To assess if the number of genes shared between two data sets was significantly greater than the expected number of shared genes due to random chance, we performed Monte Carlo simulations [25]. Each simulation compared two data sets, X and Y, for which the number of differentially expressed genes found were DEx and DEy, respectively. DEx genes were randomly selected from the genes present in data set X, and DEy genes were randomly selected among the genes present in data set Y. These genes were then translated into human orthologs using the Rosetta Stone. A total of 10,000 independent simulations were run for each pair of data sets. The proportion of simulations where the number of randomly shared genes was greater than the observed number of shared genes in the real data are reported as the p value (see example in Figure 2).

Figure 2

Monte Carlo analysis of a pair-wise comparison of two data sets. The frequency distribution of outcomes from 10,000 Monte Carlo simulations comparing the overlap of genes between two data sets (Agudo 48 h crush enriched genes and McCurley 4 day crush enriched genes) observed by random chance (red bars). The number of genes shared by the two data sets in the empirical data are actually 79. The p value is the proportion of simulations that produce several randomly shared genes that is greater than, or equal to, the empirically observed number. In this comparison, the number of overlapping genes determined empirically between the two compared data sets is significantly greater than what would be expected by random chance, and we declare this to be a significant overlap for these data. Had the empirical data revealed 56 overlapping genes, then we would have declared no significant overlap of the data sets (p>0.05).

Hierarchical clustering of mouse data set transcriptomes

Hierarchical clustering was performed on mouse data sets that used the Affymetrix Mouse Genome 430 2.0 Array. The ComBat function from the sva R package was used to remove batch effects from the RMA normalized data sets. Limma was used to calculate fold changes for all genes on the microarray. The fold changes for genes found in all data sets, regardless of the significance of the change, were used to calculate Pearson’s correlation between all pairs of data sets. A dendrogram was created using Ward’s method to cluster data sets, with the distance metric calculated as one minus Pearson’s correlation.

GO term and KEGG pathway analysis

Enriched GO terms were identified in gene lists using the Protein Analysis Through Evolutionary Relationships (PANTHER) tool Version 10. PANTHER uses an over-representation test to assess for the over-representation of GO terms associated with a group of genes. We used a p value of 0.05 as a cutoff for statistically over-represented GO terms. We identified biologic pathways association with lists of DE genes using the KEGG Mapper tool. For both analyses, pruned lists of human genes were submitted for analysis.

Scripts used for analysis

The scripts used to normalize data, identify DE genes and run comparative analysis used Python or R and can be found on GitHub. This URL also contains the Rosetta Stone Ortholog table.

Results

Normalization of data and detection of differentially expressed (DE) genes

Before conducting a meta-analysis of independent studies, it was first necessary to transform raw data from these data sets using a common normalization protocol, and then apply a uniform method for determining sets of DE genes. This manipulation of raw data yielded DE gene sets that varied from the DE gene sets reported in the literature using the same data (Appendix 1). Notably, some data sets (i.e., the Sharma study) produced a very small number of DE genes using the limma protocol and comparisons involving these data sets were rarely significant.

Hierarchical clustering of mouse data sets

To assess the fidelity of data sets generated by different groups, hierarchical clustering was performed using the fold changes of every gene in the transcriptome. To conduct this analysis, a gene must be present in all data sets. Therefore, we limited our analysis to data sets that were created using data from the Affymetrix Mouse Genome 430 2.0 Array. The resulting tree exhibited two main branches (Figure 3). On the upper branch, the optic nerve head data sets tended to cluster together, while the data sets from the retina tended to cluster on the lower branch. Data sets from other regions of the CNS clustered across both branches. Thus, data sets generated by different groups provided enough similarities to accurately cluster according to origin of tissue when assessing the ONH and retina.

Figure 3

Hierarchical clustering of the complete transcriptomes of data sets using the Affymetrix Mouse Genome 430 2.0 Array. The relative level of every gene in the experimental sample was calculated for all the study data sets that had used this array chip. Hierarchical clustering was used to determine how similar the transcriptomes were from samples generated by different groups. Overall, the cluster shows two distinct distant branches, with optic nerve head (ONH) data sets segregating to the upper branch, and retina data sets segregating to the lower branch. The exceptions to this are the “No or Early” retina data sets from the Howell study, which have greater similarity to transcriptomes of the ONH, rather than pathologic retinas. Other data sets from studies of brain-related CNS neurodegeneration appear to distribute between both arms of the cluster.

Comparison of the transcriptional responses of AONI and CONI

The optic nerve crush and axotomy models are commonly used to study degeneration of the optic nerve following acute damage [26]. We tested the hypothesis that the expression changes in AONI models mimic the expression changes observed in studies of CONI. Significant overlap in the transcriptomic changes was observed between AONI studies (Appendix 2) and between CONI studies (Appendix 3). While significant overlap did occur, the numbers of overlapping genes was surprisingly lower than anticipated. For example, when comparing two relatively similar experimental paradigms, the enriched genes in the Agudo 7 day axotomy rat retinal data set to the Templeton 5 day crush DBA/2J mouse retinal data set (Appendix 2), only 10.4% of the DE genes overlapped from the Agudo study and 8.1% of the DE genes overlapped from the Templeton study. Since we performed a common normalization and declaration of DE genes analyses of both these data sets, the lack of greater numbers of overlapping genes potentially underscores technical variations that have originated from each group, or a limitation in the sensitivity of the microarray technology to detect quantitative changes in gene expression, especially when comparing data from different platforms.

For intercomparisons between AONI and CONI data sets, all data sets that contained a non-zero number of DE genes were included. A comparison of enriched genes among retinal data sets revealed significant overlap of DE genes between several AONI and CONI data sets (Figure 4A). Interestingly, acute data sets from rat, mouse and zebrafish significantly overlapped with CONI data sets from dog and mouse. This finding suggests that the retinal response to injury is conserved, not only between acute and chronic models of optic nerve injury, but also evolutionarily from zebrafish to mammals.

Figure 4

Pair-wise comparison of acute (vertical) and chronic (horizontal) optic nerve injury data sets show moderate levels of overlapping DE genes for comparison of (A) enriched and (B) depleted genes. The total number of DE genes identified in the independent analysis is shown for each data set. In each cell, the total number of overlapping genes between two data sets is recorded. Based on Monte Carlo simulations, the cells are color-coded red for a non-significant overlap, or green for a significant overlap.

A comparison of depleted genes among retinal data sets revealed a similar pattern of significance (Figure 4B). The rat, zebrafish and mouse AONI data sets significantly overlapped the CONI data sets from both dog and mouse.

Because degeneration of the optic nerve is critical to the progression of RGC degeneration, we examined how the gene expression changes in the optic nerve following AONI or CONI compared with the expression changes in the retina. Comparing transcriptional changes in the retina and optic nerve head following AONI or CONI revealed a significant overlap of the transcriptomes of several of the retinal data sets with the optic nerve data sets (Appendix 4). This suggests a similar tissue response to acute and chronic optic nerve injury in both the retina and the optic nerve, which could be mediated by cellular elements common to both, such as astrocytes. Evaluation of specific genes and pathways (see below) could help to define this possibility.

The comparison of depleted genes shared between AONI and CONI optic nerve data sets showed significant overlap between acute data sets and the moderate and severe glaucoma data sets, but not between the acute and “no or early” glaucoma data sets. This indicates that downregulation of gene expression occurs early in acute injury paradigms, but at much later stages in the more chronic models, consistent with the idea that chronic injury is less severe.

Gene Ontology (GO) term enrichment of genes common to retinal data sets of AONI and CONI

Next, PANTHER and KEGG were used to characterize the biologic processes and pathways that drove the significant similarity between AONI and CONI transcriptomes. Among retinal enriched genes, there were 47 genes found in at least 2 of the 3 AONI studies and both CONI studies (Appendix 5). The GO terms associated with this list of genes were “Regulation of Cell Migration,” “Intracellular Signal Transduction” and “Regulation of Response to Stimulus” (Table 1). Four genes from the list, RRAS, DUSP16, TNFRSF1A, and FLNB are associated with the MAPK signaling pathway. GNAI2 and CSF1R are part of the RAP1 pathway that regulates cell-cell adhesion and cell motility. CSF1R, EHD2 and HLA-A are associated with endocytosis and ATP6V0A2, CTSS, and LAPTM5 are associated with the lysosomal pathway.

Table 1

Enriched GO terms for DE genes found in both acute and chronic optic nerve injury data sets.

GO Term	ES	Genes
Regulation of Cell Migration	6.28	RAP2A, RRAS2, MSN, STAT3, CSF1R, SULF1, CORO1A, RRAS, AIF1, EDN2
Intracellular Signal Transduction	4.23	MOV10, RAP2A, RRAS2, RAB20, GNAI2, STAT3, RHOC, CSF1R, PSMB8, RRAS, PRDX4, TYROBP, COL1A2, AIF1, EDN2, TNFRSF1A
Regulation of Response to Stimulus	2.67	LY86, RAP2A, GNAI2, B2M, C1QB, BIRC5, STAT3, RHOC, SESN3, CSF1R, SULF1, PSMB8, SOX2, RRAS, DUSP16, TYROBP, COL1A2, TULP3, AIF1, EDN2, CTSS, TNFRSF1A, HLA-A
Depleted
GO Term	ES	Genes
Neurofilament Bundle Assembly	>100	NEFM, NEFL, NEFH
Response to Amine	33.65	PPP3CA, NME1, ARG1, KCNC1, CALM1
Synapse Organization	15.14	SNCB, RAB3A, SNCG, CACNB4, ATP2B2, CACNG2, NRXN3, NRXN1, UNC13B,
Neurotransmitter Secretion	14.52	CPLX1, RAB3A, GAD1, NRXN3, NRXN1, UNC13B
Regulation of Cation Channel Activity	13.68	FGF12, CACNB4, RASGRF1, CACNG2, NRXN1, CALM1
Modulation of Synaptic Transmission	11.5	GRM1, SLC6A1, RAB3A, CLSTN1, PPP3CA, GNAI2, SNCG, RASGRF1, GRM8, ATP2B2, NRXN1, UNC13B
Axon Development	7.02	CRMP1, EPHA5, RAB3A, NEFM, NEFL, CTNNA2, NEFH, KIF5C, NRXN3, NRXN1
Neuron Projection Morphogenesis	6.11	CRMP1, EPHA5, NRN1, RAB3A, NEFL, CTNNA2, NEFH, KIF5C, NRXN3, NRXN1
Brain Development	6.08	EPHA5, NME1, SLC6A17, ATP2B2, OGDH, KCNC1, NEFL, CTNNA2, NDRG4, PITPNM1, INA, ATP2B2, HPCA, NRXN1, CALM1
Ion Transport	3.98	CPLX1, SLC24A4, KCNA3, SLC6A1, RAB3A, PPP3CA, CHRNA6, GRIA3, SLC6A17, CACNB4, KCNC1, PITPNM1, CLIC4, ATP2B2, CACNG2, GRIA2, UNC13B, CALM1

The analysis was performed with all genes that occurred in 2 of 3 acute optic nerve injury studies and both chronic optic nerve injury studies. Analysis was performed for both (A) enriched genes and (B) depleted genes. The GO terms listed (along with corresponding DE genes identified) were selected based on their Enrichment Score (ES).

Among retinal depleted genes, 74 were found in at least two of five AONI studies and both CONI studies (Appendix 6). The GO terms associated with this list of genes were related to neuronal cellular processes like “Synaptic Organization and Transmission” and “Neuronal Development” (Table 1). The genes in this list fall within KEGG pathways that are neuronally related like glutamatergic synapse and neuroactive ligand-receptor interaction. This result demonstrates that the gene transcripts that are depleted in the retina following optic nerve injury are expressed primarily in neurons.

Comparison of retinal AONI and CONI with neurodegenerative data sets

Having observed a conserved response to injury between different species and injury paradigms in two tissues after optic nerve damage, we tested the hypothesis that many of the transcriptional responses found in AONI and CONI were also conserved among models of neurodegeneration in different regions of the CNS. Transcriptomic changes in AONI and CONI data sets were compared to mouse models of amyotrophic lateral sclerosis, experimental autoimmune encephalomyelitis, and Alzheimer disease (Figure 5). The comparison of enriched genes demonstrated significant overlap of DE genes between both AONI and CONI with all examined models of neurodegeneration (Figure 5A). Taken together, these results suggest a transcriptional response to heterogeneous insults that is conserved across species, insults and spatial location of the insult within the CNS.

Figure 5

Pair-wise comparison of optic nerve injury (vertical) and brain-related CNS neurodegenerative disease (horizontal) data sets show widespread overlap of (A) enriched DE genes, but not (B) depleted DE genes. The total number of DE genes identified in the independent analysis is shown for each data set. In each cell, the total number of overlapping genes between two data sets is recorded. Based on Monte Carlo simulations, the cells are color-coded red for a non-significant overlap, or green for a significant overlap.

In contrast to the comparison of enriched genes, comparing depleted genes between AONI, CONI and neurodegeneration data sets revealed almost no significant overlap of DE genes (Figure 5B). Since the depleted genes in the retina and optic nerve resulted primarily from loss of signal from dying neurons, we hypothesized that the absence of significant depleted DE gene overlap is due to loss of signal from genes that are selectively expressed by specialized neuronal populations in various CNS locations. We tested this hypothesis by looking for highly represented GO terms among pruned lists of DE genes from each neurodegeneration data set. The pruned list of depleted DE genes in the Ferraiuolo 120-day data set was not over-represented for any Biologic Process GO Terms. The list of depleted genes for the Gjoneska 6-week data set was over-represented for the metabolic GO term “Cholesterol Biosynthetic Process” and for neuronal GO terms like “Regulation of Neuronal Synaptic Plasticity,” “Memory” and “Positive Regulation of Neuron Projection Development.” The Jonas Motor Cortex data set’s list of depleted genes was over-represented for the GO term “Retinoid Metabolic Process” and terms associated with cellular migration and morphogenesis. The Jonas Sensory Cortex data set’s list of depleted genes was over-represented for GO Terms related to central nervous system development, chemotaxis and angiogenesis. These results suggest that the neuronal gene expression patterns in different parts of the CNS are variable enough to prevent significant overlap of depleted genes following injury.

After observing a conserved transcriptional response in AONI, CONI and neurodegeneration studies, it was surprising that no enriched genes were DE in all data sets. To identify the most conserved genes across all studies, the inclusion criteria were relaxed to include enriched genes that were DE in two AONI studies, two CONI studies and studies modeling two different neurodegenerative diseases. This list included these 15 genes: C1QB, CD68, CP, CSF1R, CTSS, DECR1, DUSP15, HLA-A, LAPTM5, LY86, MPEG1, MSN, SERPINA3, SESN3, TYROBP (Table 2). This list was too small to identify highly represented GO terms among genes shared between AONI, CONI and other neurodegenerative models. To generate a list of ample size for GO term analysis, we relaxed the inclusion criteria further to identify genes that were present in at least one AONI study, one CONI study and at least 2 different models of neurodegenerative disease. These criteria identified 118 enriched genes (Appendix 7). The GO terms highly-represented by this list of genes were related to immune processes like the complement cascade, interferon signaling and lymphocyte migration (Table 3). Among the numerous KEGG pathways associated with this list of genes were the complement pathway, PI3K-AKT signaling pathway, TNF signaling pathway, the Toll-Like Receptor signaling pathway, the JAK-STAT signaling pathway and the phagosomal and lysosomal pathways. Eight genes in the list are associated with the complement pathway. The complement cascade drives phagocytosis, which is also associated with several genes in the list. FCGR1A, FCGR2B, FCGR3A, C3 and TLR2 are all phagocytosis-promoting receptors present in the list. Furthermore, TUBB and 5 different cathepsins are members of the phagosome pathway that were present in our list. TUBB is a scaffold protein that phagosomes travel along and cathepsins are the hydrolases that have diverse function including the digestion of phagocytosed products following acidification [27,28].

Table 2

Highly represented genes between optic nerve injury (ONI) and CNS neurodegenerative disease (ND) data sets are listed.

Gene	Average FC (Min - Max)	Associated GO Terms
C1QB	3.68 (1.99–12.58)	Complement Activation, Innate Immune Response
CD68	3.55 (1.80–13.15)	Cellular Response to Organic Substance
CP	3.13 (1.36–12.73)	Ferroxidase Activity, Cellular Iron Ion Homeostasis, Aging
CSF1R	2.71 (1.52–10.91)	Inflammatory Response, Cellular Proliferation, Axon Guidance
CTSS	3.85 (1.47–18.95)	Toll-Like Receptor Signaling Pathway, Antigen Processing and Presentation, Extracellular Matrix Dissassembly
DECR1	1.63 (1.18–2.55)	NADPH Binding, Protein Homotetramerization
DUSP15	1.63 (1.18–2.10)	Transforming Growth Factor Beta Receptor Signaling Pathway, Positive Regulation of JNK Cascade
HLA-A	2.85 (1.47–7.67)	Immune Response, Antigen Processing and Presentation of Peptide Antigen via MHC class I
LAPTM5	2.52 (1.54–5.63)	Protein Binding, Transport
LY86	2.78 (1.44–11.51)	Apoptotic Process, Innate Immune Response
MPEG1	3.24 (1.49–15.53)	Integral Component of Membrane
MSN	1.91 (1.21–3.73)	Leukocyte Cell-Cell Adhesion, Leukocyte Migration
SERPINA3	6.19 (1.48–12.92)	Inflammatory Response, Acute Phase Response, Platelet Degranulation
SESN3	1.55 (1.14–2.55)	Cellular Response to Amino Acid Stimulus, Negative Regulation of TORC1 signaling, Regulation of Response to Reactive Oxygen Species
TYROBP	3.00 (1.79–12.05)	Cellular Defense Response, Intracellular Signal Transduction, Innate Immune Response

All the DE genes that are shared between these data sets are enriched. The individual genes, along with the average fold change (FC) among the shared all data sets, and the GO terms associated with that gene are is listed. The inclusion criteria for this list were enrichment of the DE gene in at least two acute ONI, two chronic ONI, and two neurodegeneration data sets.

Table 3

Enriched GO terms for genes occurring in optic nerve injury (ONI) and CNS neurodegenerative disease (ND) data sets are listed.

GO Term	ES	Genes
Positive Regulation of Apoptotic Cell Clearance	95.11	C2, CCL2, C4B, C3
Detection of Molecule of Bacterial Origin	60.53	TREM2, TLR1, C4B, TLR2
Response to Interferon-Alpha	49.93	IFITM3, IFITM2, IFITM1, AXL, OAS1
Positive Regulation of Antigen Processing and Presentation	39.16	TREM2, CD74, PYCARD, SLC11A1
Regulation of Complement Activation	26.01	C2, A2M, C4B, SERPING1, C3
Neuron Projection Regeneration	25.22	APOD, TSPO, GFAP, LAMC1, FAS
Regulation of Lymphocyte Migration	22.49	CXCL10, MSN, APOD, PYCARD, CCL2
Type I Interferon signaling pathway	21.14	IFITM3, IFITM2, IFITM1, IFI35, GBP2, OAS1, HLA-A
Complement Activation, Classical Pathway	20.81	C2, C1QB, C4B, SERPING1, C1QA, C1QC, C3
Interferon-Gamma-Mediated Signaling Pathway	19.02	ICAM1, MT2A, HLA-DQB1, TRIM22, GBP2, OAS1, FCGR1A, HLA-A

To obtain sufficient DE genes to conduct an analysis of common GO terms, the inclusion criteria were DE genes that appeared in at least one acute ONI, one chronic ONI, and two neurodegeneration data sets. The GO terms listed (along with corresponding DE genes identified) were selected based on their Enrichment Score (ES).

Discussion

Inter-study comparison of gene expression data

Considering the low concordance among data sets obtained from similar injury paradigms, we conclude that a sizeable fraction of the data produced by whole transcriptome studies may not be biologically relevant. Similar discordance has been reported in other comparisons of microarray studies [8,29,30]. The cohort of DE genes that are not shared between transcriptomes may result from both biologic and technical variability. Various, subtle environmental stimuli can have an impact on gene expression as can the age, sex and social context of an animal [31,32]. Biologic variability between samples has been demonstrated to be the major contributor to heterogeneity in gene expression studies [33]. These sources of biologic variability are further confounded by technical variability which are explored in detail elsewhere [33,34].

Furthermore, the deviation in the number of DE genes we detect from the number of DE genes reported by the authors of the various studies demonstrates that the statistical method chosen to analyze data has a dramatic impact on the number of DE genes detected for the same raw data. Using a uniform approach to analyze heterogeneous studies was necessary to compare the results from different studies. Importantly, this method of normalization and DE gene detection did not work well for all the data we attempted to analyze, which illustrates one difficulty in comparing the data from gene expression studies with different study designs, controls and numbers of microarrays.

Although the transcriptional response to acute injury is variable among different strains of the same species [16,35], we were able to observe a conserved transcriptional response between acute and chronic optic nerve injuries and between any two species examined. Since both AONI and CONI insult RGC axons in the optic nerve [26], we speculate that this conserved response is the result of shared damage signals from axons.

Common pathways in response to neuronal injury

Analysis of KEGG pathways implicates RAS signaling in the retinal response to both the AONI and CONI paradigms. RAS signaling is known to cause increased phagocytosis, MAPK signaling and reorganization of the actin cytoskeleton in the cell. These RAS signaling targets were represented among the GO terms identified and in the KEGG pathways, which suggests that RAS signaling controls these cellular functions after CNS injury.

The RAP1 GTPase signaling pathway was common to both AONI and CONI. In RGCs, RAP1 is important for transducing neurotrophic signals, particularly NGF, from the distal axon to the cell body via signaling endosomes [36]. This signaling induces sustained MAPK activation in PC12 cells and is important in neuronal migration [36,37]. Further work is needed to characterize the role of RAP1 in retinal response to injury, but the presence of pathway members in both optic nerve injury paradigms hints at its importance.

The PI3K-AKT pathway was common to ONI and neurodegenerative paradigms. The cellular processes affected by PI3K-AKT activation are diverse and overlap with several conserved biologic processes we identified with GO and KEGG including cellular migration, proliferation, TNF signaling and autophagy [38,39]. In neurons, PI3K-AKT acts in a neuroprotective fashion by inhibiting degenerative MAPK signaling immediately following injury, suggesting that modulation may be therapeutically beneficial in a host of neurodegenerative diseases [40,41]. It is difficult to speculate on the implication of PI3K-AKT involvement in the CNS injury response because all detected enriched genes associated with this pathway were either cell surface receptors or extra-cellular matrix proteins which act upstream of PI3K. Further analysis is needed to identify which cells activate PI3K-AKT following injury and to tease out the implication of this mechanism on pathological outcome.

The activity of the JNK and p38 mediated MAPK pathway in neurodegeneration is well characterized. Numerous groups have described the role of JNK activation in the degeneration and apoptosis of RGCs [42-44]. JNK activation and subsequent nuclear accumulation of c-JUN lead to neuronal apoptosis [1]. Furthermore, a recent study identified a pivotal role for MAPK signaling in axon degeneration, which precedes apoptosis and must be prevented to maintain functional neurons [41]. Our findings demonstrate that these pathways are activated in ONI and neurodegeneration data sets, indicating that therapeutic intervention into this pathway may be beneficial for a host of human neurologic diseases. As expected, the depleted genes shared by AONI and CONI models were related to neuronal cellular processes, indicating that they resulted from the degeneration and death of RGCs. The identification of pathologically relevant depleted genes requires a more granular analysis focusing on early time points following injury.

Innate immune response pathways common to injury in the CNS

The observation that the transcriptional response to nerve injury is conserved throughout the CNS is evidence that these tissues have a common response to various injury stimuli. The magnitude and persistence of this response depends on the nature of the injury, but the activated biologic pathways are redundant for all tested CNS insults. Analysis of the common genes using Gene Ontology and KEGG pathways implicates the innate immune system as a major contributor to this similarity. The observation of several complement cascade related genes being shared between optic nerve injury and neurodegenerative studies agrees with previous studies demonstrating induction of the cascade in neurodegeneration and illustrates that this pathway is regulated transcriptionally [11,18,45,46], although the effect on disease progression remains controversial and may be disease or component specific [47,48].

Our findings demonstrate that transcriptional activation of the tumor necrosis factor (TNF) pathway is common to CNS injuries. TNF-α and its receptor are both more prevalent following injury and this has been shown to contribute to retinal ganglion cell death in several injury paradigms and in human glaucomatous tissue [49,50]. The TNF pathway genes shared between ONI and neurodegeneration data sets were downstream effectors that influenced cell adhesion, leukocyte recruitment and intracellular signaling (Bcl3 and Socs3).

TLR signaling has emerged as an important mechanism for glia to react to injury in the CNS [51-55]. TLR receptors on glia respond to damage associated molecular patterns that can come from external stimuli, like bacterial LPS, or from neuronally released signals (i.e., HSPs), which cause a myriad of glial responses that generally exacerbate the pathologic outcome [56,57]. TLR knockout animals and those treated with TLR antagonists have abrogated neuronal cell loss following injury relative to wild-type or untreated mice [58-60]. Multiple TLRs contribute to neuronal cell loss following artificial injury and it is unclear which are truly important for human neurodegenerative disease. It is likely that there is redundant TLR signaling following CNS injury.

The GO term “Leukocyte migration” was highly represented in the list of enriched genes common to ONI and neurodegeneration data sets, suggesting that leukocytic infiltration may be a common phenomenon among CNS diseases. Circulating leukocytes can be recruited into CNS tissues with glia activated by TNF-α or the TLR-4 ligand LPS [61]. The specific role of infiltrating immune cells is varied and controversial, but their presence in a host of neurodegenerative models implicates them as a factor in the progression of CNS disease [62-64].

The list of enriched genes that was conserved between AONI and CONI contained numerous genes related to the lysosomal and endocytic pathways. We also observed several autophagy related genes in the list of enriched genes found in AONI, CONI and neurodegeneration data sets. Since endocytosis and autophagy are activated in ONI models following injury and are impaired in other neurodegenerative “protein accumulation” diseases, these two pathways are vital for neuron health and the ability of neurons to respond to stress [65,66]. It remains unclear if the activation of autophagy in models of ONI is detrimental to the survival of retinal ganglion cells [67,68]. In a protein accumulation model of neurodegeneration, the activation of autophagy increased neuron survival and decreased the number of tau positive cells, suggesting a therapeutic benefit of autophagy activation [69].

Together, these results demonstrate that the innate immune response is activated following neuronal injury. This activation can exacerbate the extent of tissue damage or aid in recovery in different contexts [70]. These findings indicate that this response is significantly similar between optic nerve injuries, autoimmune neuronal injuries and in injuries caused by the overexpression of neurotoxic peptides.

A small, highly conserved group of genes that respond to neuronal injury

We identified 15 genes in two AONI data sets, two CONI data sets and two different models of neurodegeneration. Of the 15 highly conserved genes, C1QB, CD68, CSF1R, CTSS, HLA-A, LY86, MPEG1 and TYROBP are well characterized as members of the innate or adaptive immune response following CNS injury, or are expressed by monocytes, macrophages or microglia.

Another conserved gene, ceruloplasmin (CP), is a ferroxidase that is expressed by astrocytes and Müller cells in the retina and may have antioxidant properties [71,72]. Increased expression of CP has been noted in the rodent retina following optic nerve crush and in human and murine glaucomatous retinas [72,73]. Iron accumulation resulting from deficient CP is a hallmark of neurodegenerative disease [74]. Furthermore, ceruloplasmin has been implicated in controlling the production of nitric oxide by glia following stimulation in the CNS through signaling mechanisms that require further study [75].

The enzyme 2,4-dienoyl-CoA reductase, DECR1, is responsible for the beta oxidation of unsaturated fatty acids in the mitochondria [76]. Activation of this metabolic pathway is downregulated in cancer and may be controlled by AKT [77,78]. A proteomic screen of brain tissue from patients with atypical frontotemporal lobar degeneration identified DECR1 as differentially expressed, but further work is necessary to characterize the role of DECR1 in the healthy CNS and in neurodegeneration [79].

The ortholog mapping of dual specificity phosphatase, Dusp genes, between species is difficult because there are multiple Dusp genes in each species and they do not map to a human gene with the same gene symbol. Therefore, we report DUSP15 as a highly-enriched gene, but the DE genes that were detected in the various data sets were not named Dusp15, even though they are orthologous to human DUSP15 (see methods). From this, we postulate that Dusp genes are transcriptionally enriched following neuronal injury and are potentially important for understanding and treating neurodegenerative disease.

Lysosomal-associated protein multispanning transmembrane 5, LAPTM5, is a lysosomal protein whose function is incompletely understood. LAPTM5 is expressed in immune cells and contributes to the NFκB and MAPK signaling of immune cells following exposure to TNF [80,81]. These factors suggest LAPTM5 may fit well into the narrative of innate immune activation discussed here, but more work is needed to characterize the role of LAPTM5 in neurodegeneration.

Moesin (MSN) is a member of the Ezrin-Radixin-Moesin (ERM) family of proteins which link actin filaments and microtubules to the plasma membrane and play a role in extracellular signal transduction [82]. MSN is important for the maintenance and alteration of cell shape as exemplified by a study in Drosophila, which demonstrated that MSN antagonizes rhodopsin (RHO) activity to maintain cell morphology [83]. MSN can also bind to microtubules to assist in the mitotic process [84]. In the context of neurodegeneration, MSN has been characterized as a substrate of the Parkinson Disease related gene, LRRK2, whose dysregulation leads to irregularities in the cellular cytoskeleton [85]. The ERM family member, Ezrin, has also been identified as a contributor to the progression of Huntington’s disease [86]. Therefore, it appears that cytoskeletal dynamics are disrupted in neurodegeneration, but the impact of ERM family proteins remains to be explored.

Serine protease inhibitor A3, SERPINA3, is a protein that blocks the protease activity of several serine proteases, including Cathepsin G [87]. Although it is normally localized to the extracellular space, SERPINA3 has been seen in the nucleus of cancerous cells. Furthermore, SERPINA3 can bind to DNA, perhaps to prevent DNA polymerization following DNA damage [87]. Importantly, SERPINA3 is the only Serpin gene expressed in astrocytes [87]. SERPINA3 is found in the senile plaques associated with Alzheimer disease, but the implication of its presence in plaques is still being worked out. Finally, other in silico analyses have identified SERPINA3 as a common entity in neurodegenerative disease [88-90].

Sestrin 3, SESN3, is the third and least characterized member of the sestrin gene family. Sestrin genes activate autophagy and prevent the accumulation of reactive oxygen species, although they may accomplish this through indirect means [91]. Transcriptional regulation of sestrin genes is tied to p53 and forkhead transcription factors [91,92]. Interestingly, AKT and RAS activation have been shown to suppress SESN3 expression, which conflicts with the results observed here. It is possible that the expression of these genes comes from different cells types or that transcriptional regulation in the damaged CNS is controlled by some novel mechanism. The role of sestrin genes, particularly Sestrin 3, in neurodegenerative conditions requires further study.

The findings presented here demonstrate a common transcriptional response to heterogeneous injury stimuli. The similarity in the responses of different stimuli result from the innate immune system activation in all conditions. The Rap1 pathway and the presence of cathepsins are understudied aspects of retinal response to ONI that require further study. CP, DECR1, Dusp genes, LAPTM5, MSN, SERPINA3 and SESN3 were identified as understudied genes that were conserved across numerous neurodegenerative paradigms. Understanding the function of these conserved genes and pathways is critical for understanding and treating neurodegeneration.