Introduction {#s1}
============

The retinal ganglion cell (RGC) is the final output neuron of the retina, projecting through the optic nerve to the brain, where it targets a number of functionally distinct areas: for visual perception, RGC axons travel to the lateral geniculate nucleus (Chalupa and Günhan, [@B6]); for the regulation of circadian rhythms, they pass through the suprachiasmatic nucleus (Guido et al., [@B19]); for eye movements, a group of RGC axons terminates in the superior colliculus (Triplett et al., [@B68]); and for the pupillary light reflex, RCG axons terminate in the pretectal area (Young and Lund, [@B80]). Each of these areas receives input from distinct subtypes of RGCs with unique morphological and molecular signatures. At the present time, over 30 subtypes of RGCs (Baden et al., [@B2]) are estimated to exist. They all receive inputs from other types of retinal neurons (bipolar cells and amacrine cells), and most of them express similar groups of genes that may serve as general RGC markers (Raymond et al., [@B52]; Rodriguez et al., [@B53]). Identifying gene expression patterns in RGCs and their subtypes is currently an active area of research, as demonstrated by the discovery of new subtypes of ganglion cells based on gene expression (Macosko et al., [@B39]; Sanes and Masland, [@B57]).

The death of RGCs in glaucoma or after injury eventually leads to loss of vision (Templeton et al., [@B66]; Zode et al., [@B83]; Munguba et al., [@B44]; Nuschke et al., [@B46]). However, the susceptibility of RGC subtypes to death differs among the distinct subtypes. Some RGCs are resistant to injury, while others appear to be more sensitive to insult, indicating differential gene expression and response to injury among subtypes (Chang et al., [@B7]; Duan et al., [@B12]; Puyang et al., [@B50]). The present study focuses on transcriptional networks within RGCs of the mouse, using gene expression data measured across 55 strains of recombinant inbred BXD mice (King et al., [@B34]) as well as the bioinformatic tools from GeneNetwork (Williams and Mulligan, [@B76]). The analysis begins by examining genes correlated with two relatively general RGC markers, *Thy1* and *Tubb3*. Each of these markers forms a unique network of genes that appears to function independently across many of the RGC subtypes. These networks are functionally different to the point of having distinct transcription factor binding sites. Subtype-specific networks partially overlapping with the *Thy1*-network are also present. In a meta-analysis of previously published data from a microarray study of a mouse glaucoma model (Howell et al., [@B24]), we examine the differential effects of this disease state on transcriptional networks in RGC subtypes and confirm intrinsically photosensitive RGCs (ipRGCs) and alpha-RGCs as more resistant to cell death (Duan et al., [@B12]).

The systems genetics and bioinformatics approach used in the present study demonstrates how signatures of RGCs and their subtypes can be extracted from a complex neural tissue such as the retina.

Results {#s2}
=======

RGC markers segregate into two major correlation networks
---------------------------------------------------------

The present study examines the correlation of gene expression in the retina across the BXD recombinant inbred strain set to define gene networks active in RGCs. The BXD strain set is derived from two parental strains, the C57BL/6J mouse and the DBA/2J mouse. Natural variation in gene expression across strains can be used to identify co-regulated genes with a similar expression pattern, allowing for the construction of genetic networks (Williams et al., [@B75]; Geisert et al., [@B17]; Templeton et al., [@B65]; Keeley et al., [@B30]).

The data used in this study consist of whole retinal samples collected from 55 BXD strains. They can be found on [www.genenetwork.org](http://www.genenetwork.org) under the identifier "DoD Retina Normal Affy MoGene 2.0 ST (May15) RMA Gene Level". Two features of this dataset enhance the quality of the analysis: The first is that the retina is a tissue that can be consistently isolated with minimal contamination by other tissues. The second is the quality of the RNA with an average RNA Integrity Score of 9.43 and a standard error of 0.037 across the 220 samples isolated for this dataset.

The analysis began with two well-characterized markers for RGCs, *Thy1* and *Tubb3*, as they both exhibited substantial variability in mRNA expression levels across BXD strains (Figure [1](#F1){ref-type="fig"}). Expression of *Thy1* ranged from 10.39 in BXD42 to 11.45 in BXD15 (this data is presented on a log~2~ scale and the difference in expression is equivalent to an over 2-fold change). A similar variability in gene expression was observed for *Tubb3*, with expression levels ranging from 9.68 (BXD6) to 10.74 (BXD2). When correlations for *Thy1* and *Tubb3* were made across all microarray data, both produced a highly correlated group of genes. For *Thy1*, the top 100 correlates had an absolute *r*-value (Pearson) greater than 0.89 (Bonferroni-adjusted *p* \< 1e^−12^) and the top 2000 correlates all had an absolute *r*-value greater than 0.77 (adj. *p* \< 1e^−8^). If we examine the Pearson correlation values above 0.60, which corresponds to an adjusted *p*-value of 0.02, then *Thy1* has a total of 8596 correlates. This tightly correlated list of genes forms a potential network, and indicates that the genes in this network are co-regulated across the BXD strains (King et al., [@B34]). Within the list of the top 2000 *Thy1*-correlates, we found several other well-characterized RGC markers, including *Rbfox3* (producing NeuN, Neuronal Nuclei), *Pou4f1* (producing BRN3A), and *Pou4f2* (producing BRN3B) (Table [1](#T1){ref-type="table"}). Interestingly, the *Thy1* correlate list did not contain other known RGC markers, including *Tubb3*.

![**RNA expression across the BXD RI strain set for *Thy1* (left) and *Tubb3* (right) with the means and standard errors**. There was an approximately 2-fold difference in mean expression levels for both genes. This data is given as raw expression values on a log~2~ scale +8.](fgene-07-00169-g0001){#F1}

###### 

**List of all RGC marker genes used in this manuscript**.

  **Symbol**                                              **RGC Subtype Labeled**                                                                **References**
  ------------------------------------------------------- -------------------------------------------------------------------------------------- ----------------------------------------------------------------------------
  Thy1                                                    pan-RGC                                                                                Barnstable and Drager, [@B4]; Raymond et al., [@B52]
  Rbfox3 (NeuN)                                           pan-RGC                                                                                Wolf et al., [@B77]; Schlamp et al., [@B60]; Kim et al., [@B33]
  Pou4f2 (Brn3b)                                          RGC (about 50--60% of total population)                                                Xiang et al., [@B78]; Erkman et al., [@B13]; Jain et al., [@B27]
                                                          ipRGCs (71% of all Melanopsin-positive cells)                                          
  Pou4f1 (Brn3a)                                          RGC (about 60--70% of total population)                                                Erkman et al., [@B13]; Jain et al., [@B27]
  Tubb3 (class III beta-tubulin)                          pan-RGC                                                                                Mellough et al., [@B41]
  Rbpms (Retina binding protein with multiple splicing)   pan-RGC                                                                                Piri et al., [@B49]; Rodriguez et al., [@B53]
  Nefl (Neurofilament light)                              RGC (\~85% of all RGCs)                                                                Ruiz-Ederra et al., [@B55]
  Chrna6                                                  RGC                                                                                    Mackey et al., [@B38]; Munguba et al., [@B45]
  Slc17a6 (Vglut2)                                        RGC                                                                                    Bai et al., [@B3]; Mimura et al., [@B42]; Wässle et al., [@B72]
  Nrn1                                                    RGCs                                                                                   Picard et al., [@B48]; Sharma et al., [@B63]
  Calb2                                                   pan-RGC (87% of all RGCs) transient OFF α-RGCs (tOFF-αRGCs, Huberman et al., [@B25])   Huberman et al., [@B25]; Mojumder et al., [@B43]; Haverkamp et al., [@B22]
  Sncg (gamma-synuclein)                                  RGC                                                                                    Buckingham et al., [@B5]
  Opn4 (Melanopsin)                                       ipRGC                                                                                  Semo et al., [@B62]
  Jam2                                                    J-RGC (5% of all RGCs)                                                                 Daniele et al., [@B9]; Kim et al., [@B32]
  Spp1 (Osteopontin)                                      alpha-RGC                                                                              Ju et al., [@B28]; Sanes and Masland, [@B57]
  Kcng4                                                   alpha-RGC                                                                              Duan et al., [@B11]; Sanes and Masland, [@B57]
  Cartpt                                                  ooDSGC                                                                                 Adams et al., [@B1]; Kay et al., [@B29]
  Hoxd10                                                  ON-DSGC                                                                                Dhande et al., [@B10]

For *Tubb3*, the top 100 genes had an absolute *r* value greater than 0.71 (adj. *p* \< 1e^−6^), and the top 2000 genes showed values greater than 0.52 (adj. *p* = 0.12) and did not contain *Thy1, Pou4f1, Pou4f2*, or *Rbfox3.* This was also true for the 1387 correlates with Pearson *r* \> 0.60 (equal to a significant Bonferroni-adjusted *p*-value of 0.02). However, other RGC markers were present in the *Tubb3* correlation list, including the newly described *Rbpms* (RNA binding protein with multiple splicing) as well as *Calb2* (Calbindin 2) and *Chrna6* (Cholinergic receptor nicotinic alpha 6).

These two networks are relatively independent, with minimal overlap. When examining the correlates of the *Thy1*-network relative to the *Tubb3*-network, only one gene is shared within their top 100 correlations (1%), 31 genes (1.55%) are present in the top 2000 correlations, and 51 genes (0.05%) are in common with the 9982 genes found in both the *Thy1* (8596 genes) and the *Tubb3* (1386 genes) correlation lists with a Pearson correlation above 0.6.

The basis for the segregation of genes into two distinct networks is illustrated by plotting the correlations for combinations of genes from both networks. In Figure [2](#F2){ref-type="fig"}, the expression of *Thy1, Pou4f1, Tubb3*, and *Rbpms* across BXD strains is displayed in scatterplots. These plots demonstrate the tight correlation for *Thy1-Pou4f1* (Figure [2A](#F2){ref-type="fig"}) and *Tubb3-Rbpms* (Figure [2C](#F2){ref-type="fig"}). They further show the lack of correlation between *Thy1* and *Tubb3* (Figure [2B](#F2){ref-type="fig"}) as well as *Tubb3* and *Pou4f1* (Figure [2D](#F2){ref-type="fig"}). The presence of additional RGC markers in the correlation lists of *Thy1* and *Tubb3* together with the minimal overlap of genes indicated the presence of at least two RGC-specific transcriptional networks in mouse RGCs.

![**Scatterplots illustrating the correlation between major RGC markers**. There was a tight correlation between *Thy1* and *Pou4f1* **(A)**, which are part of the same network. Correlation dropped greatly for *Thy1* and *Tubb3* **(B)**, which are not part of the same network. Similarly, *Tubb3* correlated well with Rbpms **(C)** but not with *Pou4f1* **(D)**. Pearson\'s correlation coefficient is given in each plot on the lower right half. Each dot represents one BXD RI strain, and the confidence interval for the smoothing function (dark gray areas surrounding the blue line) is 0.95.](fgene-07-00169-g0002){#F2}

Since we selected *Thy1* and *Tubb3* to serve as primary RGC markers for this analysis based upon literature evidence, we wondered if an alternative unbiased approach would validate our findings. As an independent examination of these networks, we performed weighted gene correlation network analysis (WGCNA) on the whole microarray dataset (Langfelder and Horvath, [@B36]). This method relies on unsupervised clustering of co-expressed genes across all of the BXD RI strains into so-called modules or eigengenes and thus represents an unbiased approach to test the assumption that *Thy1*- and *Tubb3*- networks exist independently from each other. WGCNA created 18 modules of co-expressed genes. Module \#1 contained 1741 of the top 2000 *Thy1* correlates, whereas only 131 genes from the top 2000 *Tubb3* correlates were assigned to this module. The majority of the *Tubb3* correlates belonged to modules \#4 (*n* = 279), \#5 (*n* = 886), and \#6 (*n* = 327), none of which contained any of the *Thy1* correlates. This minor overlap in module affiliation between genes of both networks reaffirms our finding that the *Thy1*- and *Tubb3*-network function individually.

To define the genomic location of upstream modulators for both the *Thy1*- and the *Tubb3*-network, we examined interval maps for the 7 RGC marker genes found in both correlation lists (*Thy1, Rbfox3, Pou4f2, Pou4f1, Tubb3, Calb2*, and *Rbpms*). When we investigated each marker\'s signature quantitative trait locus (QTL), we found that *Thy1, Rbfox3, Pou4f2*, and *Pou4f1*, had similar interval maps. The same was true for *Tubb3, Calb2*, and *Rbpms*. To further investigate this phenomenon, we plotted multiple interval maps as heat maps for each marker and its 20 highest correlated genes (Figure [3](#F3){ref-type="fig"}). In these heat maps, rows correspond to the QTL curve of a single correlated gene plotted across the entire genome. Two color gradients were used to characterize differential expression between strains: A yellow to red gradient identified a transcript whose expression was higher in strains with a B haplotype (allele origin from C57BL/6J), whereas a green to blue gradient represented a transcript whose expression was higher in strains with the D haplotype (allele origin from DBA/2J). The linkage significance (LRS or LOD-score) increased with color intensity. In other words: A deeply colored vertical line characterized a genomic locus that may contain a regulatory element responsible for the differential expression of these genes.

![**RGC markers segregate into two major correlation networks**. For each of the general RGC markers on the right, the eQTL curve for the 20 highest correlated genes was plotted as a heat map. In these, the LRS/LOD score is given in pseudocolors: A yellow to red gradient identifies a transcript whose expression is higher in strains with a B haplotype at that locus (allele origin from C57BL/6J), whereas a green to blue gradient represents a transcript whose expression is higher in strains with the D haplotype (allele origin from DBA/2J). There are several strong and sharp trans-bands extending across *Thy1, Rbfox3, Pou4f2*, and *Pou4f1*, such as on distal Chromosome 1 or 13. There are also trans-bands extending across *Tubb3, Calb2*, and *Rbpms* on mid Chr. 13 and proximal Chr. 14. No overlap is present between the genes separated by the black line, indicating that RGC markers segregate into two major independently regulated gene networks. The panel on the right lists genes that are found in more than one of the top 20 correlations. For example, *Arhgap44* is present in the 20 highest correlates for both *Thy1* and *Pou4f1*.](fgene-07-00169-g0003){#F3}

These bands are thought to identify a genomic locus modulating the expression of the genes in the network across the BXD RI strains. Likely candidates for modulating gene expression include transcription factors, micro RNAs or long noncoding RNAs (Geisert et al., [@B17]; Templeton et al., [@B67]; Williams and Auwerx, [@B74]). Since the analysis is only correlational in nature, we cannot exclude the possibility that loci for individual networks do not have regulatory roles, as they could just be co-regulated with the other identified networks. When comparing heat maps for both networks, no overlap in patterns was observed. Thus, the RGC marker genes segregated into two independently regulated gene networks. From here on, we will refer to these networks as the *Thy1*-network and the *Tubb3*-network. This is an arbitrary nomenclature, mirroring the most prominent RGC marker for each of both networks. For the *Thy1*-network, the strongest modulatory signature was localized at the distal end of Chromosome 1. For members of the *Tubb3*- network, the most prominent genomic signature was on mid-distal Chromosome 13. Several candidate genes exist in those loci; a comprehensive analysis is attached as Supplemental Material (Supplemental Tables [1](#SM1){ref-type="supplementary-material"}, [2](#SM2){ref-type="supplementary-material"}).

There are three potential scenarios that would allow two specific genetic networks to exist within a single cell population. The first is that the mRNA levels of each marker differ across the entire retina. The second explanation is that their mRNA expression levels differ in specific RGC subtypes, leading to unique expression of each protein within each RGC subtype. The third potential mechanism involves differences in the percentage of each of the RGC subtypes from strain to stain. To further examine this, we immunostained retinas with THY1 and Class III beta tubulin. In these retinal whole mounts, the more than 90% of labeled cells were double-labeled of RGCs were both THY1- and Class III beta tubulin- positive (Figure [4](#F4){ref-type="fig"}), indicating that the two markers are co-localized and expressed in the same cell. Despite their protein co-expression, their mRNA level may differ from cell to cell. When examining the double stained cells, the intensity of each stain varied from ganglion cell to ganglion cell. Some cells had approximately equal labeling for THY1 and Class III beta tubulin, while in other cells THY1 staining was more intense than Class III beta tubulin. In a few cells, the labeling of Class III beta tubulin was more intense. It also appeared that some cells were labeled by only one of the two markers. Thus, even though both markers are co-expressed in many RGCs, their expression levels do not consistently correlate with each other, corroborating that their expression is independently regulated.

![**Maximum intensity projections of confocal z-stacks taken from retinal whole-mounts showing the co-localization of THY1 (A) and Class III beta-tubulin (B) in the retinal ganglion cell layer from a C57BL/6 mouse**. Most cells were double-stained for both RGC markers **(C)**, but some cells only expressed Class III beta-tubulin (arrowhead) or THY1 (arrow). Furthermore, the staining intensity was different across cells, and some cells had large somata and were more intensely stained than others (large arrowhead in "merge"). Staining with the secondary antibodies only did not result in unspecific fluorescence (data not shown). Scale bar in C = 100 μm.](fgene-07-00169-g0004){#F4}

Since multiple subtypes of RGCs are known to exist, one could predict that there are also multiple genetic networks within each of these subtypes. We compiled an extended list of RGC markers from the literature. This list consisted of 17 proteins and their respective genes, 6 of which are known to be relatively specific for a single RGC subtype (Table [1](#T1){ref-type="table"}). Calculating the correlation in expression across the BXD RI strain set and displaying these relations in a network graph revealed that indeed two major network hubs formed around *Thy1* and *Tubb3* based on the highest correlations surrounding these two markers (Figure [5](#F5){ref-type="fig"}). The only significant connection between *Thy1* and *Tubb3* existed through 2 of their respective correlates, *Slc17a6* (VGLUT2), and *Chrna6*. Except for *Hoxd10* (a marker of ON-directionally selective RGCs), all the other RGC markers had a direct connection to either *Thy1* or *Tubb3*.

![**RGC markers segregate into two networks, and two major hubs are formed around *Thy1* and *Tubb3***. The only connection between the *Thy1*- and the *Tubb3*-network is through two of their correlates, *Slc17a6* and *Chrna6*. Pearson\'s correlation coefficient was mapped to line color and thickness (high correlation = thicker and red bars).](fgene-07-00169-g0005){#F5}

Functional differences of RGC networks
--------------------------------------

Following the identification of two gene networks related to normal RGC function, we asked whether both networks had any functional differences.

In order to facilitate this analysis, we grouped genes with a similar heat map pattern belonging to the same network into one "synthetic trait" using principal component analysis. These "synthetic traits" provide for more robust functional network analysis than a single marker gene, and they can be used to generate correlation lists specific to a group of genes by reducing the dimensionality of the data (Yin et al., [@B79]; Vanderlinden et al., [@B70]; Graybeal et al., [@B18]).

The synthetic trait for the *Thy1*-network was constructed using *Thy1* and its direct correlates *Rbfox3, Pou4f1*, and *Pou4f2*, while the *Tubb3*-network was constructed using *Tubb3, Calb2*, and *Rbpms*. When we compared the correlations of these networks using WebGestalt-based gene ontology (GO) analysis (Wang et al., [@B71]), each of the two synthetic traits was enriched in genes with very different molecular and cellular functions. The *Thy1*-network was significantly enriched in genes involved in neuron development, synaptic transmission, cation transmembrane transporter activity and voltage-gated channel activity (adj. *p* \< 0.0000001 for all), strongly indicating that this network was highly neuron-specific (Supplemental Figure [S1](#SM4){ref-type="supplementary-material"}). Some examples of the genes associated with neuronal development and synaptic transmission included: *Tnr, NeuroD2, L1cam, Syn2, Grn4*, and *Gabbr1* (see Supplemental Table [3](#SM3){ref-type="supplementary-material"} for the whole list). A close examination of the highest (*r* \> 0.9) correlates also revealed that at least 6 genes (*Gsk3a, Srgap1, Arhgap44, Ncdn, L1cam*, and *Lrrc4b)* were functionally associated with neurite outgrowth.

The *Tubb3* synthetic network contained genes that were associated with completely different biological functions (Supplemental Figure [S2](#SM5){ref-type="supplementary-material"}). Several members of the tubulin family of proteins were at the top of the correlation list, such as *Tubg2, Tuba4a, Tubg1, Tubb5*, and *Tuba1b*. Gene Ontology analysis revealed enrichment in GO terms "guanyl nucleotide binding" and "protein polymerization" (both with adj. *p* \< 0.001), including the genes *Rab11b, Bab8a, Tufm, Rabb1b, Arf6*, and *Rab4b*. The highest correlates of the *Tubb3* synthetic trait were two genes encoding proteins from the Rab family (*Rab15* and *Rab4b*).

Distinct groups of transcription factor binding sites are associated with the two RGC networks
----------------------------------------------------------------------------------------------

Since the *Thy1*- and the *Tubb3*-network were functionally enriched for distinct GO terms, we asked if this difference would be reflected in their transcription factor binding site (TFBS) distribution. We searched the promoter sequences for genes in the *Thy1*-network using the TRANSFAC FMatch algorithm (Matys et al., [@B40]) for over-represented TFBS in comparison to the *Tubb3*-network and vice versa (Tables [2a](#T2){ref-type="table"},[b](#T3){ref-type="table"}). Genes of the *Thy1*-network were \~4 times more enriched in TFBS for *p53* and *Dec1*, two master regulators of cell cycle progression (Qian et al., [@B51]). The *Thy1*-network was also more enriched in TFBS for PPAR gamma and estrogen receptor alpha as well as effectors SP1 and AP1. In contrast, the promoters of genes in the *Tubb3*-network were significantly enriched in TFBS of developmental origin, such as *Pax6, Six6*, proteins of the *Sox* and *Oct* family, as well as *Pou4f1*.

###### 

**Transcription factor binding site enrichment for genes of the *Thy1*-network**.

  **Transcription Factor**   **Enrichment probability (fold) vs. Tubb3**   **Matched promoters *p*-value**
  -------------------------- --------------------------------------------- ---------------------------------
  P53                        3.6921                                        1.17E-04
  TBP                        4.7405                                        2.47E-03
  PPARgamma:RXR-alpha        1.1771                                        1.07E-02
  PPAR direct repeat         1.7507                                        1.34E-02
  LXR, PXR, CAR, COU         1.9258                                        1.93E-02
  DEC1                       4.4442                                        2.11E-02
  FOXJ1                      4.4442                                        2.11E-02
  AP-1                       1.2523                                        3.47E-02
  SP1                        1.1402                                        4.62E-02
  ER-alpha                   2.4295                                        5.18E-02

###### 

**Transcription factor binding site enrichment for genes of the *Tubb3*-network**.

  **Transcription Factor**   **Enrichment probability (fold) vs. Thy1**   **Matched promoters *p*-value**
  -------------------------- -------------------------------------------- ---------------------------------
  Pax-6                      1.3845                                       2.13E-05
  OTX                        1.8106                                       5.81E-05
  SRY                        1.428                                        3.26E-04
  Oct1                       1.6876                                       2.71E-03
  FOXO1                      1.1574                                       3.08E-03
  Sox1                       3.0377                                       3.10E-03
  Nkx6-2                     1.3845                                       3.16E-03
  Tst-1                      1.2976                                       3.55E-03
  Oct4 (POU5F1)              1.8934                                       3.90E-03
  Foxc1                      1.0931                                       4.14E-03
  Foxm1                      1.0472                                       5.44E-03
  SIX6 secondary motif       1.6129                                       5.65E-03
  NF-AT                      1.2657                                       6.66E-03
  Pitx3                      1.2304                                       7.34E-03
  Brn-2                      1.6214                                       9.59E-03
  POU4F1                     4.219                                        9.66E-03
  c-Myc:Max                  1.8081                                       1.14E-02
  Bach1                      10.1255                                      1.25E-02
  Dlx2                       1.4746                                       1.36E-02
  POU2F1                     1.1806                                       1.45E-02

Genomic regulation of subtype-specific RGC markers
--------------------------------------------------

Since our analysis revealed two distinct molecular networks governing normal RGC function, we hypothesized that subtype-specific RGC markers would contain regulatory signatures from either one or both RGC networks.

Thus, we generated heat maps of the highest correlates of each subtype-specific marker (*Cartpt2, Jam2, Kcng4, Opn4, Spp1*, and *Hoxd10*, see Figure [5](#F5){ref-type="fig"}). This comprehensive analysis revealed that two pairs of genes were regulated in a similar fashion: *Cartpt* and *Jam2* as well as *Kcng4* and *Opn4* showed very similar heat maps to each other. When we compared the heat maps for these 4 subtype-specific markers, we found that some of their bands were at identical loci as bands observed in the *Thy1*-network (Figure [6](#F6){ref-type="fig"}). While *Cartpt/Jam2* had the distinct band from distal Chromosome 1, *Kcng4*/*Opn4* shared a band on distal Chromosome 13. It was also noticed that the allelic distribution of correlates for *Cartpt/Jam2* seen as a blue line on distal Chromosome 1 was diametrical to the original *Thy1*-correlates, where the line was mostly yellow-red (see Figure [3](#F3){ref-type="fig"}). This simply implies that the top 20 correlates of *Cartpt* and *Jam2* are in fact correlated inversely to *Cartpt* or *Jam2* itself, whereas the top 20 correlates of *Thy1* are positively correlated to *Thy1* itself.

![**Heat maps for RGC subtype-specific marker genes**. *Cartpt/Jam2* share the trans-band on distal Chromosome 1 with the *Thy1*-network, whereas *Kcng4/Opn4* share the *Thy1*-network trans-band from Chromosome 13. There is no obvious overlap of trans-bands with either network for *Hoxd10*. For *Spp1*, no obvious trans-bands can be appreciated, suggesting that this gene is not part of a transcriptional network in its normal state.](fgene-07-00169-g0006){#F6}

*Hoxd10* had very distinct trans-QTL bands that were not found elsewhere, and its correlates also did not contain any RGC-enriched genes. Since no connection of *Hoxd10* to either the *Thy1*- or the *Tubb3*-network (see Figure [3](#F3){ref-type="fig"}) was observed, the lack of overlap with RGC-enriched genes casts doubt on the specificity of this marker for the identification of RGC subtypes based on gene expression. The eQTL pattern of *Spp1* was the most complex of all markers, as it lacked any clearly detectable *trans*-band.

Overlap of correlates of subtype-specific markers with the *Thy1*-network and the *Tubb3*-network
-------------------------------------------------------------------------------------------------

When we took the top 2000 correlates (an arbitrary cut-off) of each of the six subtype specific markers (*Cartpt, Jam2, Kcng4, Opn4, Hoxd10*, and *Spp1*) and examined the distributions of these markers across the *Thy1*- and the *Tubb3*-networks, we were able to see the relative interplay between these networks (Table [3](#T4){ref-type="table"}). For *Spp1* and *Hoxd10*, the majority of the overlap was seen with genes from the *Tubb3*-network. The picture was quite different for the remaining subtype-specific markers. For *Cartpt, Opn4, Kcng4*, and *Jam2*, the overlap was dominated by similarities with the *Thy1*-network. This analysis was also conducted using only the genes with a Pearson correlation value above 0.6 (adj. *p* \< 0.02) and it displayed a similar trend (Table [3](#T4){ref-type="table"}). Taken together these data indicate that each of the RGC subtypes has a unique interplay with the *Thy1*- and the *Tubb3*-networks. The dramatic differences in overlap between correlates may reflect the differential expression of the *Thy1*- and *Tubb3*-networks in different RGC subtypes. These relationships can also be seen in the network map (Figure [5](#F5){ref-type="fig"}).

###### 

**Overlap of subtype-specific RGC markers with genes from the *Thy1*- and the *Tubb3*-network**.

  **Subtype Marker**   **Top 2000**                           
  -------------------- ------------------------ ------------- -----------------------------
  Cartpt               1022                     54            1076 (27%)
  Jam2                 382                      25            407 (10.2%)
  Kcng4                326                      241           567 (14.2%)
  Opn4                 713                      185           898 (22.5%)
  Hoxd10               10                       19            29 (0.7%)
  Spp1                 29                       756           1785 (19.6%)
  **Subtype Marker**   **Pearson *r* \> 0.6**                 
                       ***Thy1***               ***Tubb3***   **Total Overlap (of 9982)**
  Cartpt               2482                     0             2482 (24.9%)
  Jam2                 1030                     3             1033 (10.3%)
  Kcng4                886                      207           1093 (10.9%)
  Opn4                 1086                     56            1142 (11.4%)
  Hoxd10               4                        17            21 (0.2%)
  Spp1                 13                       9             22 (0.2%)

For the top table, the top 2000 correlates of both networks served as comparison, whereas a Pearson r cutoff of \>0.6 was chosen for the lower table, corresponding to a Bonferroni-corrected value of p \< 0.02.

Susceptibility of RGC marker genes to glaucomatous nerve damage
---------------------------------------------------------------

The identification of RGC subtypes has raised the question whether or not there is a difference in susceptibility to nerve damage, subsequent cell survival, or axon regeneration. In order to test this hypothesis on a transcriptional level, we used publicly available microarray data generated from DBA/2J glaucomatous eyes (Howell et al., [@B24]) and correlated glaucoma severity score (GSS) with expression levels of RGC markers (Figure [7](#F7){ref-type="fig"}). As expected, expression levels of all general RGC markers decreased as nerve damage increased (implicated by the negative correlation between GSS and RGC markers in Figure [7](#F7){ref-type="fig"}, dashed blue line). This was also the case for two subtype-specific markers, *Jam2* and *Cartpt*. Interestingly, expression levels of markers for alpha-RGCs (*Kcng4* and *Spp1*) and ipRGCs (*Opn4*) did not correlate to GSS (cutoff *r* \< 0.6), suggesting that these two RGC subtypes are differentially susceptible to nerve damage.

![**Correlation network of Glaucoma Severity Score (GSS) to RGC markers during glaucoma progression**. The GSS is a visual grading system identifying axonal damage in the optic nerve and consists of 4 stages: no damage, light, medium, and severe damage. These stages are inversely correlated (dashed blue) to the expression levels of most RGC markers, suggesting that as RGCs die and axonal damage increases, mRNA expression of RGC marker genes decreases (most likely due to decrease in RGC number). The decrease of RGC marker gene expression is strongly correlated across glaucoma stages (red, *r* for all \>0.9). Markers for ipRGCs (*Opn4*) and alpha-RGCs (*Kcng4* and *Spp1*) do not correlate to GSS, possibly indicating their preferential survival.](fgene-07-00169-g0007){#F7}

Discussion {#s3}
==========

Using the DoD CDMRP Retina expression microarray dataset to examine correlates of known RGC markers across the BXD RI strain set, this study revealed two distinct gene networks regulating normal RGC function. Detecting these transcriptional networks with high statistical power required tissue from 55 RI strains with 4 biological replicates per strain. The bioinformatic tools on GeneNetwork allowed us to extract gene networks from a complex structure such as the mouse retina (Geisert et al., [@B17]; Keeley et al., [@B30]). Interestingly, almost all of the general RGC markers described in the literature (see Table [1](#T1){ref-type="table"}) selectively group with one of two genetic networks, the *Thy1*-network or the *Tubb3*-network. Both *Thy1* (Barnstable and Drager, [@B4]) and *Tubb3* (Snow and Robson, [@B64]) are believed to be generalized markers for RGCs. Nonetheless, these markers segregate into two distinct genetic networks. Two scenarios could explain the biological basis of correlated genes that are co-regulated in BXD RI mice: (I) differences in cell number of each RGC subtype across the BXD strains, or (II), differences in gene expression within individual RGC subtypes (Geisert et al., [@B17]; Keeley et al., [@B30]). Currently, we are not able to distinguish between these two potential explanations. Our analysis also revealed that these networks have unique biological functions.

Both the *Thy1*- and the *Tubb3*-networks contained large groups of genes that were consistent with a neuronal phenotype. However, they were each enriched in distinct functional roles. The *Thy1*-network appeared to be functionally involved in neuronal development and maintenance, as well as axon guidance. Some of the highest correlates of this network (*Gsk3a, Srgap1, Arhgap44, Ncdn, L1cam*, and *Lrrc4b)* were all previously reported to affect dendrite or axon pathfinding and outgrowth (Schwaibold and Brandt, [@B61]; Cherry et al., [@B8]; Ip et al., [@B26]; Galic et al., [@B16]). Therefore, this network could be directly involved in RGC process extension, potentially including dendrites and axons projecting to the brain. Furthermore, the network was enriched in transmembrane ion transporters. This included several calcium- and sodium-gated channel subunits, all of which have been previously localized to RGCs (Lipton and Tauck, [@B37]; Farrell et al., [@B15]).

The *Tubb3*-network contained genes that differed in function from those of the *Thy1*-network. The biological processes enriched in this group included protein polymerization and organic substrate metabolic processes. Under molecular processes, the *Tubb3*-network was enriched for GTP binding proteins and GTPase activity, two mechanisms known to be essential for cellular transport involving the cytoskeleton (Roychowdhury and Rasenick, [@B54]; Schappi et al., [@B58]). The highest two correlates of the Tubb3-network were *Rab15* and *Rab4b*, two G-proteins that are known to play an important role in endosome formation and vesicle movement along actin and tubulin networks (Zuk and Elferink, [@B84]; Falk et al., [@B14]). Several other genes in this list were associated with intracellular trafficking, such as sorting nexin 32 (*Snx32*, a molecule that links transport vesicles to dynactin (Wassmer et al., [@B73]), vesicle associated membrane protein 1 (*Vamp1*, a SNARE protein important for linking vesicles to the target membrane; Hasan et al., [@B21]), *Pigu* (a GPI anchor protein; Guo et al., [@B20]), and *Tmem9* (part of the lysosomal membrane; Kveine et al., [@B35]). These genes were highly correlated with and specific to the *Tubb3*-network.

These findings suggest that the *Tubb3*-network is functionally associated with cytoskeleton and vesicle transport in RGCs. Since neurons are known for the significant amount of cytoskeletal proteins necessary to maintain dendritic and axonal integrity (Kevenaar and Hoogenraad, [@B31]), it can be hypothesized that the *Tubb3*-network is involved in maintenance of these structural elements.

Taken together, this analysis indicates that each network has distinct molecular functions and that both of these networks can exist independently within a single RGC.

The functional differences between networks were mirrored by the transcription factor binding site (TFBS) analysis for promoters of their respective genes. Members of the *Thy1*-network showed significant enrichment in TFBS for proteins regulating cell cycle, growth and apoptosis (Tu et al., [@B69]). Genes belonging to the *Tubb3*-network, however, showed an over-representation of TFBS for transcription factors known to establish cell fate decisions during development, particularly in later stages. Proteins such as OTX1, PAX6, SIX6, or POU4F1 belonged to the latter group (Zagozewski et al., [@B81]). *Pou4f1* was not originally found in the correlates of the *Tubb3*-network, which might seem surprising. A possible explanation for this is that the *Thy1*-network, which correlates well with the expression of *Pou4f1* across the BXD strains, acts upstream of the *Tubb3*-network. In a developmental sense, this could suggest that members of the *Thy1*-network establish expression of *Pou4f1*, which can then bind to members of the *Tubb3*-network to influence their transcription. Because *Pou4f1* was originally not found to be part of the *Tubb3*-network but the *Thy1*-network instead, we hypothesize here that basic cell cycle tasks such as survival decisions are managed by members of the *Thy1*-network, while members of the *Tubb3*-network regulate the fine-tuning of RGC development, possibly even influencing subtype commitments.

The identification of neuronal subtypes has gained a considerable amount of attention in neuroscience across a variety of scientific disciplines (Hoshino, [@B23]; Russ and Kaltschmidt, [@B56]). This is especially the case for RGC subtypes (Kay et al., [@B29]; Dhande et al., [@B10]; Sanes and Masland, [@B57]). The establishment of transgenic mouse models has recently transformed the field, and new paradigms are starting to emerge. For example, it is becoming increasingly evident that RGC subtypes seem to be differentially susceptible to nerve injury to the point of having different chances of survival or regeneration (Kay et al., [@B29]; Dhande et al., [@B10]; Pérez de Sevilla Müller et al., [@B47]; Duan et al., [@B12]; Sanes and Masland, [@B57]). In order to develop new therapeutic approaches to nerve regeneration, it may be necessary to tease apart the transcriptional programs of injury-susceptible and non-susceptible RGC subtypes.

When we investigated the heat maps of subtype-specific gene correlates, we found that two pairs---*Kcng4/Opn4*, or *Jam2/Cartpt*---showed similarity to two different heat map bands from the *Thy1*-network. *Kcng4* marks alpha-RGCs, while *Opn4* marks ipRGCs, and these two subtypes were recently found to be most resistant to optic nerve crush (Pérez de Sevilla Müller et al., [@B47]; Duan et al., [@B12]). Since the regulatory pattern of those two gene pairs was very similar based on heat map analysis of their co-varying genes, it would be interesting to systematically investigate if these upstream modulatory loci are responsible for increased neuronal survival or regeneration. Bioinformatic analysis of RGC marker susceptibility to glaucomatous nerve damage supported the notion that alpha-RGCs and ipRGCs were also most resistant to neurodegeneration. The decrease in mRNA expression of alpha-RGC markers (*Kcng4, Spp1*) and ipRCG markers (*Opn4*) did not correlate with GSS (a visual grading scale for optic nerve axon damage), whereas the overwhelming majority of RGC markers did.

Interestingly, the discovery of *Spp1* marking a particular RGC subtype had only been made after crush injury to the optic nerve in a mouse model (Pérez de Sevilla Müller et al., [@B47]; Duan et al., [@B12]). In this study, it was found that alpha-RGCs secreting Osteopontin (the protein made from the gene *Spp1*) were most resistant to optic nerve crush among all RGC subtypes, and that this protective effect was due to Osteopontin. We were not able to identify putative upstream modulatory loci for this gene due to the lack of distinct trans-QTL bands. This suggests that *Spp1* is part of a genetic network that needs to be activated by neuronal injury before it can be co-regulated in a more RGC-specific way. This phenomenon has been observed elsewhere for other genes following retinal injury (Templeton et al., [@B67]).

In summary, we have identified multiple loci modulating RGC function in the BXD mouse strain set, and we have provided *in silico* evidence for the differential susceptibility of RGC subtypes to neurodegeneration and cell death. One caveat of this study is the fact that the performed analysis is only a correlative one. Proving a biological cause will require experimental manipulation of the mentioned genes and examination of the effects *in vivo* or *in vitro*. Nevertheless, our findings enhance the understanding of the RGC\'s normal transcriptome, as they are the first to describe gene regulatory networks for some of their subtypes. They may serve others and us as a reference for future studies on RGC subtype identification and their susceptibility to injury.

Materials and methods {#s4}
=====================

Animals
-------

All of the procedures involving mice were approved by IACUC at Emory University and adhered to the ARVO Statement for the Use of Animals in Research.

Microarray datasets
-------------------

Two microarray datasets were used in this project as part of a comprehensive meta-analysis.

i.  The Department of Defense (DoD) Normal Retina Database (May2015). This is the most comprehensive retina microarray dataset and creation is described in King et al. ([@B34]). The DoD Normal Retina Dataset consists of 222 microarrays from 55 different strains BXD mice, the parental strain C57BL/6J, the parental strain DBA/2J and an F1 cross.

ii. The Howell et al. ([@B24]), DBA/2J Glaucoma Retina M430 2.0 (Sep11) RMA database. This dataset consists of retinal tissue dissected from 40 DBA/2J mouse eyes at 10.5 months of age showing varying levels of RGC damage due to naturally occuring glaucoma in this mouse strain (as graded by visual inspection of the optic nerve). Twenty eyes served as negative controls. The generation of this dataset is described in Howell et al. ([@B24]).

All of the used datasets are publicly accessible through [www.genenetwork.org](http://www.genenetwork.org).

Statistical analysis and plot generation
----------------------------------------

GeneNetwork provided the platform for correlation analysis, principal component generation, and linkage analysis. In general, datasets were queried for gene symbols, downloaded from GeneNetwork, and additional analysis was performed in R whenever necessary. *P*-values mentioned in relation to Pearson\'s coefficient throughout this paper are based on pair-wise comparisons. All *p*-values were Bonferroni-adjusted for 36,012 genes, which is equal to the number of genes captured on the microarray after accounting for replicated and wrongly annotated probes. Plots were generated with R and the dplyr and ggplot2 packages. Cytoscape version 3.2.1 was used to generate gene network graphs.

Gene ontology analysis
----------------------

GO analysis was performed using WebGestalt (Wang et al., [@B71]). The top 500 genes of each network were used to compile a gene list and duplicates were removed. GO term enrichment was calculated using Affymetrix MouseGene 2.0 ST probe set IDs against a background dataset from the same chip. GO trees are appended as expanded view data. All *p*-values presented in the GO analysis are corrected for multiple comparisons using the Bonferroni method.

Candidate gene analysis
-----------------------

A list of candidate genes for putative upstream modulators of both networks was created by extracting genes with cis-QTLs from the loci showing trans-bands. For a cis-QTL considered to be a good candidate gene, it should fulfill three conditions: (i) be expressed above background, (ii) correlate well with the trait used to create the trans-bands, and (iii) have a nucleotide variation that alters either protein structure or other regulatory elements (such as promoters or enhancers). These lists are appended as Supplemental Tables.

Transcription factor binding site analysis
------------------------------------------

Gene lists including the top 500 correlated genes in each network were compiled, duplicates were removed, and uploaded to BIOBASE. The TRANSFAC FMatch algorithm was used to search for overrepresented neuron-specific transcription factor binding sites in the gene lists using the "best supported promoter" model for each gene (Matys et al., [@B40]). As a background dataset for both networks, the opposite network was chosen. Resulting *p*-values were Bonferroni-adjusted. The enrichment window for a gene\'s promoter was chosen to be −5000 to +100 bp of its best supported transcription start site. "Minimize false positives" was selected as cut-off method. The TRANSFAC data version was 2016.2, and only high-quality matrices were used for the analysis.

WGCNA analysis
--------------

Weighted gene network analysis was performed using the WGCNA package in the R environment (Langfelder and Horvath, [@B36]). A soft thresholding power of 10 was used to calculate the adjacency matrix based on the criterion of scale-free topology (Zhang and Horvath, [@B82]). Modules were identified using the following parameters for the blockwiseModules function: minimum module size of 100, merge cut height of 0.25, and maximum block size of 2000. The top 2000 correlates of *Thy1* and *Tubb3* were then each merged with their module affiliation by Affymetrix Probe Set ID and the resulting list was checked for overlaps.

Immunohistochemistry
--------------------

Two 60 day-old Thy1-CFP C57BL/6 transgenic mouse (Jax identifier B6.Cg-Tg(Thy1-CFP)23Jrs/J) were deeply anesthetized with a mixture of 13 mg/kg of xylazine and 87 mg/kg of ketamine and intracardially perfused with 0.9% saline followed by 4% paraformaldehyde. Retinas were dissected and flat-mounted on glass slides following a standard protocol. Retinas were blocked in 4% bovine serum albumin in PBS overnight. One retina was stained with anti-TUJ1 (a gift from Anthony Frankforter, dilution 1:1000) and anti-GFP (Invitrogen, 1:500) antibodies overnight and then labeled with appropriate Alexa Fluor-conjugated secondary antibodies. One retina served as negative secondary antibody control (data not shown). High-resolution Z-stacks were captured on a Nikon confocal microscope with the Nikon C1 software throughout the entire ganglion cell layer only. Z-stacks were collapsed using FiJi (Schindelin et al., [@B59]) and the whole image was adjusted for contrast and brightness using Adobe Photoshop.

Accessibility of data {#s5}
=====================

The data presented in this article is publicly available under [www.genenetwork.org](http://www.genenetwork.org) and can be downloaded under <http://genenetwork.org/webqtl/main.py?FormID=sharinginfo&GN_AccessionId=709>.

Author contributions {#s6}
====================

FS, EG provided the initial design of the study. FS performed experiments and bioinformatic analyses with assistance from EG. FS, EG wrote the paper with input from RW and RL. RW supplied the datasets through the GeneNetwork platform and is responsible for its maintenance. RL has provided intellectual framework for this project as part of an ongoing collaborative effort between his laboratory and the laboratories of EG, RW.

Conflict of interest statement
------------------------------

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This work was supported by DoD CDMRP Grant W81XWH1210255 from the USA Army Medical Research & Materiel Command and the Telemedicine and Advanced Technology (to EG), NIH Grant R01EY017841 (to EG), Vision Core Grant P30EY006360 (to P. Michael Iuvone), and Unrestricted Funds and a Stein Innovation Award from Research to Prevent Blindness (to the Ophthalmology Department at Emory University, the University of Tennessee Health Science Center, and the Bascom Palmer Eye Institute, University of Miami Miller School of Medicine). FS is supported by the institutional training grant T32EY007092-30 (to P. Michael Iuvone). The authors would like to thank Rebecca King for technical assistance with the retinal flat mounts and Dr. John Nickerson for his comments on the manuscript.

Supplementary material {#s7}
======================

The Supplementary Material for this article can be found online at: <http://journal.frontiersin.org/article/10.3389/fgene.2016.00169>

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

###### 

Click here for additional data file.

[^1]: Edited by: Igor Ponomarev, University of Texas at Austin, USA

[^2]: Reviewed by: Claudio V. Mello, Oregon Health & Science University, USA; Michael F. Miles, Virginia Commonwealth University, USA

[^3]: This article was submitted to Neurogenomics, a section of the journal Frontiers in Genetics
