Transcriptome networks in the mouse retina: An exon level BXD RI database.

PURPOSE: Differences in gene expression provide diverse retina phenotypes and may also contribute to susceptibility to injury and disease. The present study defines the transcriptome of the retina in the BXD RI strain set, using the Affymetrix Mouse Gene 2.0 ST array to investigate all exons of traditional protein coding genes, non-coding RNAs, and microRNAs. These data are presented in a highly interactive database on the GeneNetwork website.
METHODS: In the Normal Retina Database, the mRNA levels of the transcriptome from retinas was quantified using the Affymetrix Mouse Gene 2.0 ST array. This database consists of data from male and female mice. The data set includes a total of 52 BXD RI strains, the parental strains (C57BL/6J and DBA/2J), and a reciprocal cross.
RESULTS: In combination with GeneNetwork, the Department of Defense (DoD) Congressionally Directed Medical Research Programs (CDMRP) Normal Retina Database provides a large resource for mapping, graphing, analyzing, and testing complex genetic networks. Protein-coding and non-coding RNAs can be used to map quantitative trait loci (QTLs) that contribute to expression differences among the BXD strains and to establish links between classical ocular phenotypes associated with differences in the genomic sequence. Using this resource, we extracted transcriptome signatures for retinal cells and defined genetic networks associated with the maintenance of the normal retina. Furthermore, we examined differentially expressed exons within a single gene.
CONCLUSIONS: The high level of variation in mRNA levels found among the BXD RI strains makes it possible to identify expression networks that underline differences in retina structure and function. Ultimately, we will use this database to define changes that occur following blast injury to the retina.

Introduction

Large-scale sequencing initiatives have led to a new era in understanding gene and genome functions [1-5]. There is now an acute need for powerful approaches that integrate and analyze massive proteomics/genomics data sets. In vision research, many single gene variants are known to cause vision loss, including retinitis pigmentosa [6-9], Usher syndrome [10,11], and some forms of glaucoma [12]. However, many ocular diseases have a complex genetic basis with multiple chromosomal loci contributing to differences in the susceptibility and severity of the disease. Two prominent examples are glaucoma [13-15] and age-related macular degeneration [16,17]. In addition, the response of the eye and the retina to trauma is driven by a host of different genes expressed in a large number of different cell types.

Until recently, it was extremely difficult to define the genetic and molecular basis of complex diseases or to adequately monitor the response of the eye and the retina to injury. We used a novel and powerful approach that relies on systems biology and a mouse genetic reference panel, the BXD family of recombinant inbred (RI) strains. This resource is particularly well suited to define complex genetic networks that are also active in human diseases. This approach allows us to not only identify specific gene variants involved in retinal disease and response to injury but also place corresponding molecular changes in a global context in the eye and the retina.

The initial efforts of our group explored the genetic diversity of the BXD family of strains to define the genetic networks active in the eye (see data sets and refs [18] and [19]). In this study, we created a new mouse retinal database that offers a more complete description of the mouse transcriptome. This resource uses the genetic covariance of expression across a panel of 52 BXD strains to identify cellular signatures and genetic networks within the mouse retina. The array we used provides expression profiling at the exon level for 26,191 well-established annotated transcripts, as well as 9,049 non-coding RNAs, including more than 600 microRNAs. Using the bioinformatics tools located on GeneNetwork, we examined the cellular signature of RPE cells. We also analyzed a genetic and molecular network involved in neuronal development and axon growth. In both examples, we highlight the specific benefits of the new database with a special emphasis on microRNAs, non-coding RNAs, and the exon level data available with the Affymetrix MouseGene 2.0 ST array.

Methods

All of the procedures used involving mice were approved by IACUC at the Emory University and adhered to the ARVO Statement for the Use of Animals in Research. The Department of Defense (DoD) Congressionally Directed Medical Research Programs (CDMRP) Normal Retina Database uses the Affymetrix MouseGene 2.0 ST Array (May 15, 2015). Robust multiarray average (RMA) analysis and scaling were conducted by Arthur Centeno. This data set consists of 52 BXD strains, C57BL/6J, DBA/2J, and an F1 cross between C57BL/6J and DBA/2J. A total of 55 strains were quantified. There is a total of 222 microarrays. All data from each microarray used in this data set is publicly available on GeneNetwork.

These are RMA expression data that have been normalized using what we call a 2z+8 scale, but without corrections for batch effects. The data for each strain were computed as the mean of four samples per strain. The expression values on the log2 scale ranged from 3.81 to 14.25 (10.26 units), a nominal range of approximately 1,000-fold. After taking the log2 of the original non-logged expression estimates, we converted the data within an array to a z-score. We then multiplied the z-score by 2. Finally, we added 8 units to ensure that no values were negative. The result was a scale with the mean expression of the probes on the array of 8 units and a standard deviation of 2 units. A twofold difference in expression is equivalent to roughly 1 unit on this scale. The lowest level of expression was 3.81 (Olfr1186) from the DoD CDMRP (the Normal Retina Database uses the Affymetrix MouseGene 2.0 ST Array, May 15, 2015). The highest level of expression was rhodopsin for 17462036 (Rho). The highest single value was 14.25.

Cases used to generate this data set

Almost all animals were young adults between 60 and 100 days of age. We measured expression in conventional inbred strains, BXD recombinant inbred (RI) strains, and reciprocal F1s between C57BL/6J and DBA/2J.

BXD strains

The first 32 of the strains were from the Taylor series of BXD strains generated at the Jackson Laboratory (Bar Harbor, ME) by Benjamin A. Taylor. BXD1 through BXD32 were started in the late 1970s, whereas BXD33 through 42 were started in the 1990s. BXD43 and higher were bred by Lu Lu, Jeremy Peirce, Lee M. Silver, and Robert W. Williams starting in 1997 using B6D2 generation 10 advanced intercross progeny. This modified breeding protocol doubles the number of recombinations per BXD strain and improves the mapping resolution [20]. All of the Taylor series of BXD strains and many of the new BXD strains are available from The Jackson Laboratory. Several strains were specifically excluded from the data set. For BXD43 and higher, the DBA/2J parent carried the Tyrp1b mutation and the GpnmbR150X mutation; these two mutations produce pigment dispersion glaucoma. Mice that carried these two mutations were not included in the data set: BXD53, BXD55, BXD62, BXD66, BXD68, BXD74, BXD77, BXD81, BXD88, BXD89, BXD95, and BXD98. In addition, BXD24 was omitted, since it developed a spontaneous mutation, rd16 (Cep290) that resulted in retinal degeneration and was renamed BXD24b/TyJ [21]. Several additional strains were excluded due to abnormally high Gfap levels observed in the Full HEI Retina (April 2010) data set: BXD32, BXD49, BXD70, BXD83, and BXD89.

Tissue preparation protocol

The mice were killed by rapid cervical dislocation. The retinas were removed immediately by placing the globe under pressure and cutting the cornea. The lens burst out of the opening in the cornea followed by the retina. In this process, no specific procedures were used to include or exclude the RPE. The retina was placed immediately in 1 ml of 160 U/ml RiboLock (Thermo Scientific Waltham, MA) for 1 min at room temperature. The retina was then transferred to Hank’s Balanced Salt solution with RiboLock in 50 µl RiboLock (Thermo Scientific, RiboLock RNase #EO0381 40 U/µl 2500U) and stored in −80 °C. The RNA was isolated using a QiaCube (Hilden, Germany) and the in-column DNase procedure. All RNA samples were checked for quality before the microarrays were run. The samples were analyzed using the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA). The RNA integrity values ranged from 7.0 to 10. Our goal was to obtain data from independent biologic sample pools for both sexes of each BXD strain. The four batches of arrays included in this final data set collectively represented a reasonably well-balanced sample of males and females.

Affymetrix gene array

In the present study, we used the Affymetrix Mouse Gene 2.0 ST Array. The array was designed with a median of 22 unique probes per transcript. Each probe is 25 bases in length. The arrays provide comprehensive transcriptome coverage with more than 30,000 coding and non-coding transcripts. In addition, there is coverage for more than 600 microRNAs. For some arrays, the RNA was pooled from two retinas, and other arrays were run on a single retina. Dr. XiangDi Wang (University of Tennessee, Health Science Center) was involved in the retinal extractions and isolation of RNA. The Affymetrix arrays were run by two research core laboratories: the Molecular Resource Center at UTHSC (Dr. William Taylor, director) and the Integrated Genomics Core at Emory University by Robert B. Isett (Dr. Michael E. Zwick, director). In a separate set of experiments, we tested a set of arrays from C57BL/6J retinas run at each facility to determine if there were batch effects or other confounding differences in the results. We did not detect any significant difference in the arrays run at UTHSC or at Emory University. Thus, we included both sets of data in the analysis.

Results

The DoD CDMRP Retina Database presents the retinal transcriptome profiles of 52 BXD RI strains on a highly interactive website, GeneNetwork. There are two separate presentations of the microarray data. The first is at the gene level (DoD CDMRP Retina Affymetrix Mouse Gene 2.0 ST (May 15, 2015) RMA Gene-Level Database), and the same data is presented at the exon level (DoD CDMRP Retina Affymetrix Mouse Gene 2.0 ST (May 15, 2015) RMA Exon-Level Database). For analyzing these data sets, a suite of bioinformatics tools is integrated in the GeneNetwork website. These tools identify genes that vary across the BXD RI strains, construct genetic networks that control the development of the mouse retina, and identify the genomic loci that underlie complex traits in the retina. In this paper, we present these two new data sets and illustrate their use with three examples. The first was to identify the genetic signatures of the RPE. The second identified genes that are differentially spliced between the C57BL/6J retina and the DBA/2J retina. In the third example, we looked at the genetic network associated with roundabout homolog 2 (Robo2) gene and the modulation of axonal growth.

Cellular signature of the RPE in the DoD CDMRP Retina Database

The DoD CDMRP Retina Database has a unique signature for RPE cells. When looking at the expression of the RPE marker Rpe65, there was an almost biphasic distribution of expression (Figure 1). Many of the strains expressed low levels of Rpe65 (approximately 7 units on our scale) while other strains had high levels of expression ranging from two- to eightfold higher (8 to 11 units). We believe that this difference is due to the time of day that the retinas were removed from the eye. Many of the retinas were isolated at the University of Tennessee Health Science Center where isolation started at approximately 10:00 a.m. and lights on in the animal colony occurred at 6:00 a.m. These retinas were isolated approximately 4 h after the lights came on. At Emory University, the retinas were isolated starting at 9:00 a.m., and lights on occurred at 7:00 a.m. These later sets of retinas were removed at approximately 2 h after the lights came on in the animal colony. The mean expression level of Rpe65 in the samples isolated at the University of Tennessee was 8.3, and the mean for the samples isolated at Emory University was 10.2, roughly a fourfold difference in expression. This difference in expression between the samples isolated at the different locations was significant using a Student t test (p<0.01). If we examine the correlation between these differences in Rpe65 expression and the genes’ associated circadian rhythm, several genes correlate with the expression of Rpe65 across the data set, including cytochrome 2 (Cry2, r = –0.75) and period homolog 2 (Per2, r = –0.65).

Figure 1

The expression at the gene level of Rpe65 across the BXD strains in the DoD CDMRP Normal Retina Database. The expression levels of Rpe65 are shown for many of the BXD strains as the mean expression and the standard error of the mean. The individual strain identifications are shown along the bottom, and the scale is log2. Notice the low levels of Rpe65 in some stains (DBA/2J, BXD5 BXD12, BXD34, BXD40, BXD48a, BXD60, BXD69, BXD100, BXD101, and BXD102) and the eightfold higher levels of expression in other strains (BXD16, BXD31, BXD42, BXD43, BXD50, BXD56, BXD75, and BXD85). Most of the high expressing strains were isolated at 2 h after light on and the low expressing strains had retinas isolated at least 4 h after light on.

When this data set is used, it is important to remember the difference in Rpe65 expression due to the time of day the retinas were isolated. However, if one is interested in defining a molecular signature of RPE cells and if the level of Rpe65 is directly associated with the number of cells that adhere to the retina, then we can use the level of Rpe65 in each strain to define a set of genes that covary with Rpe65 across the BXD strain set. When we examined the data set for genes with a similar expression pattern across the BXD strains, a list of genes uniquely expressed in RPE was observed (Table 1). This cellular signature represents genes that are expressed within the RPE, including retinal G protein coupled receptor (Rgr), lecithin-retinol acyltransferase (Lrat), retinol dehydrogenase 5 (Rdh5), transferrin (Trf), and retinal pigment epithelium derived rhodopsin homolog (Rrh). This signature can also be thought of as the result of genetic networks that drive gene expression within a given cell type.

Table 1

The list of genes with similar expression patterns across the BXD strains as RPE65.

Symbol	Description	Location (Chr: Mb)	Mean Expr	Sample r
Rpe65	retinal pigment epithelium 65	Chr3: 159.262145	8.84	1.00
Rgr	retinal G protein coupled receptor	Chr14: 37.850676	10.66	0.98
Pon1	paraoxonase 1	Chr6: 5.118090	7.61	0.98
Ttr	transthyretin	Chr18: 20.823751	10.86	0.95
Ermn	ermin, ERM-like protein	Chr2: 57.897524	7.17	0.94
Lrat	lecithin-retinol acyltransferase	Chr3: 82.696501	8.19	0.93
Rdh5	retinol dehydrogenase 5	Chr10: 128.350646	8.76	0.92
Slc6a20a	solute carrier family 6	Chr9: 123.545240	8.36	0.92
Trf	transferrin	Chr9: 103.106331	10.66	0.89
Slc26a7	solute carrier family 26, member 7	Chr4: 14.429577	7.98	0.89
Car12	carbonic anyhydrase 12	Chr9: 66.561493	7.90	0.89
Clstn1	calsyntenin 1	Chr4: 148.960577	12.74	-0.88
Camk2b	calcium	Chr11: 5.869645	8.98	-0.88
Pkm	pyruvate kinase, muscle	Chr9: 59.504175	13.46	-0.88
Thbs1	thrombospondin 1	Chr2: 117.937612	8.37	0.88
Snora17	small nucleolar RNA, H/ACA box 17	Chr2: 26.494759	6.27	0.87
Abl1	c-abl oncogene 1, non-receptor tyrosine kinase	Chr2: 31.543896	8.80	-0.87
Itgb8	integrin beta 8	Chr12: 120.396495	9.95	0.87
Rrh	retinal pigment epithelium derived rhodopsin homolog	Chr3: 129.507326	10.40	0.87
Tiam1	T-lymphoma invasion and metastasis-inducing	Chr16: 89.787356	9.16	-0.87
Chchd3	coiled-coil-helix-coiled-coil-helix domain containing 3	Chr6: 32.740976	10.72	0.87
Gm5567	predicted gene 5567	Chr6: 40.060252	9.84	-0.87
Hdac5	histone deacetylase 5	Chr11: 102.086139	10.80	-0.87
Cntn2	contactin 2	Chr1: 134.406002	8.84	-0.87
Olfr875	olfactory receptor 875	Chr9: 37.580222	7.60	0.87
Gm19522	predicted gene, 19522	Chr16: 42.884483	7.57	0.87
Gm19272	predicted gene, 19272	Chr6: 7.251908	6.40	-0.87
Sez6l2	seizure related 6 homolog like 2	Chr7: 134.094049	10.14	-0.86
Nrxn2	neurexin II	Chr19: 6.418731	10.14	-0.86
Olfr726	olfactory receptor 726	Chr14: 50.703389	7.14	0.86
Sv2a	synaptic vesicle glycoprotein 2 a	Chr3: 95.985074	12.28	-0.86
Tmem161b	transmembrane protein 161B	Chr13: 84.361901	10.91	0.86
Pdcd5	programmed cell death 5	Chr1: 191.101187	11.26	0.86
Fap	fibroblast activation protein	Chr2: 62.339000	7.19	0.86
Cacna1g	calcium channel, voltage-dependent,	Chr11: 94.269705	9.95	-0.86
Rapgef1	Rap guanine nucleotide exchange factor (GEF) 1	Chr2: 29.475240	10.50	-0.86
Nsmce2	non-SMC element 2 homolog	Chr15: 59.205753	10.91	0.86
Plxna1	plexin A1	Chr6: 89.266307	10.09	-0.86
Ube2b	ubiquitin-conjugating enzyme E2B	Chr11: 51.798648	11.28	0.86
Aacs	acetoacetyl-CoA synthetase	Chr5: 125.956184	9.52	-0.86
Acsf3	acyl-CoA synthetase family member 3	Chr8: 125.299405	9.51	-0.86
Slc17a7	solute carrier family 17	Chr7: 52.419291	13.14	-0.86
Srebf2	sterol regulatory element binding factor 2	Chr15: 81.977696	11.31	-0.86
Epb4.9	erythrocyte protein band 4.9	Chr14: 71.001070	9.98	-0.86
Vapa	vesicle-associated membrane protein, associated protein A	Chr17: 65.929392	12.32	0.85
Slc12a5	solute carrier family 12, member 5	Chr2: 164.786302	12.12	-0.85
Card9	caspase recruitment domain family, member 9	Chr2: 26.207696	7.65	-0.85
Cfh	complement component factor h	Chr1: 141.982432	8.59	0.85
Dagla	diacylglycerol lipase, alpha	Chr19: 10.319755	9.57	-0.85
Pde4dip	phosphodiesterase 4D interacting protein	Chr3: 97.493751	9.55	-0.85
Ank1	ankyrin 1, erythroid	Chr8: 24.085316	9.54	-0.85
Zcchc14	zinc finger, CCHC domain containing 14	Chr8: 124.122603	10.39	-0.85
Ssr3	signal sequence receptor, gamma	Chr3: 65.186870	6.09	0.85
Mterfd1	MTERF domain containing 1	Chr13: 67.007904	9.55	0.85
G3bp2	GTPase activating protein (SH3 domain) binding protein 2	Chr5: 92.052145	8.76	0.85
Apba1	amyloid beta (A4) precursor protein binding	Chr19: 23.833366	10.55	-0.85
Acot13	acyl-CoA thioesterase 13	Chr13: 24.909817	10.63	0.85
Rbl1	retinoblastoma-like 1 (p107)	Chr2: 156.971629	9.07	0.85
Cpxm2	carboxypeptidase X 2 (M14 family)	Chr7: 139.234493	8.07	0.85
Mir467h	microRNA 467h	Chr9: 115.291078	6.80	0.85
Ret	ret proto-oncogene	Chr6: 118.101766	10.23	-0.85
Mospd1	motile sperm domain containing 1	ChrX: 50.698185	10.71	0.85
Lrp3	low density lipoprotein receptor-related protein 3	Chr7: 35.984852	9.10	-0.85

The gene symbol and name are listed, along with the chromosomal location, mean expression and correlation to Rpe65.

With the MouseGene 2.0 ST Affymetrix chip, we found not only protein-coding genes that correlate with Rpe65 but also microRNAs and non-coding RNAs. When we examined the top 500 correlates of Rpe65 (all of which have a correlation higher than 0.8 with Rpe65), five microRNAs were present: Mir98, Mir666, Mir449a, Mir301b, and Mir28b (Table 2). Using the bioinformatics tools on TargetScan (Targetscan.org) [22-24], we predicted targets for each microRNA from the top 500 correlates of Rpe65. One microRNA, Mir666, did not appear on the TargetScan website. The remaining four microRNAs appeared on the website. When the microRNAs were scanned for targets, Mir98 had 29 targets in the RPE signature, Mir449a had 14 targets, Mir301b had 13 targets, and Mir28b had one target. This type of analysis may be one approach to constructing and understanding microarray networks within a specific cell type such as the mouse RPE.

Table 2

The predicted targets from the top 500 correlates of Rpe65 are listed in columns below for each of the five microRNAs found to correlate themselves with Rpe65.

Rpe65	Mir98	Mir449a	mir301b	mir28b
Camk2b	–	–	–	Camk2b
Atp2b2	–	–	Atp2b2	–
Cfl2	–	–	Cfl2	–
Dlc1	–	–	Dlc1	–
Eif2c1	Eif2c1	–	Eif2c1	–
Nptx1	Nptx1	Nptx1	Nptx1	–
Pitpnm2	–	–	Pitpnm2	–
Ppp6r1	–	–	Ppp6r1	–
Psap	–	–	Psap	–
Slc17a7	–	–	Slc17a7	–
Snx2	–	–	Snx2	–
Sub1	–	–	Sub1	–
Tet3	–	–	Tet3	–
Zbtb4	–	–	Zbtb4	–
Zcchc14	–	–	Zcchc14	–
2610507B11Rik	2610507B11Rik	2610507B11Rik	–	–
Abr	Abr	Abr	–	–
Ahsa2	Ahsa2	Ahsa2	–	–
Cacna2d2	Cacna2d2	Cacna2d2	–	–
Cntn2	Cntn2	Cntn2	–	–
Dcaf7	Dcaf7	Dcaf7	–	–
E2f5	E2f5	E2f5	–	–
Fbxo10	Fbxo10	Fbxo10	–	–
Mgat5b	Mgat5b	Mgat5b	–	–
Ndst1	Ndst1	Ndst1	–	–
Nrxn2	Nrxn2	Nrxn2	–	–
Pvrl1	Pvrl1	Pvrl1	–	–
Ret	Ret	Ret	–	–
Slc6a1	Slc6a1	Slc6a1	–	–
Tfdp2	Tfdp2	Tfdp2	–	–
Trim67	Trim67	Trim67	–	–
Usp31	Usp31	Usp31	–	–
Agap1	Agap1	–	–	–
Apba1	Apba1	–	–	–
Bsn	Bsn	–	–	–
Fbxl14	Fbxl14	–	–	–
Insr	Insr	–	–	–
Kcnc1	Kcnc1	–	–	–
Nsmce2	Nsmce2	–	–	–
Slc7a14	Slc7a14	–	–	–
Srebf2	Srebf2	–	–	–
Syt7	Syt7	–	–	–
Thbs1	Thbs1	–	–	–

For many Affymetrix probes, there is minimal annotation. For example, within the top 100 correlates of Rpe65, 17 Affymetrix probes are present that have minimal annotation furnished by Affymetrix. For two of these probes, Affy_17203447 and Affy_17204181, there is no annotation from Affymetrix, including the sequence of the probe itself. Thus, for these two probes, no further analysis can be conducted until Affymetrix furnishes sequence data. For the remaining 15 Affymetrix probes, the expression within the retinal transcriptome is low, ranging from 4.6 to 8.8. The sequences of these probes can be related to the mouse genome using the Verify Tool on the probe’s Trait Data and Analysis page on GeneNetwork. The three probes with expression levels higher than 8 all aligned with the mouse genome. Affy_17241598 aligns to a sequence on chromosome 10 that is a predicted protein with no further annotation. Affy_17414264 aligns with a sequence on chromosome 4 that is a non-protein coding gene or gene fragment. Affy_17527409 aligns with a sequence on chromosome 9 that is a predicted protein. When the role of these transcripts in the network associated with Rpe65 is considered, these data are far from informative. Unfortunately, for many of the probes on the Affymetrix Mouse Gene 2.0 ST Array, this is the current state of annotation. We are beginning to improve the annotation, and the identification of probes associated with microRNAs is due to the efforts of members of the Rob Williams group. With time, we believe that the annotation will improve allowing investigators to include these probes in functional genetic networks.

Analysis of differentially spliced genes using the exon dataset

One of the extended features of the Affymetrix MouseGene 2.0 ST Array is its extensive coverage of gene expression at the exon level, and these data are presented in the DoD CDMRP Retina Affymetrix Mouse Gene 2.0 ST (May 15) RMA Exon-Level Database. At the present time, we do not have specific bioinformatic tools that fully investigate the exon-level data set. However, this database can be used to identify genes that are differentially spliced in the DBA/2J mouse and the C57BL/6J mouse. If an exon is expressed in one strain of mouse and not the other strain, the exon will have a large and significant likelihood ratio statistic (LRS) score across the BXD RI strain set. Basically, that individual exon will function like a Mendelian trait being either highly expressed or expressed at a low level. Therefore, to begin the analysis, we identified the exons with LRS scores higher than 60. We identified 2,314 exons, and the highest LRS score was 250. Then we reasoned that if an exon had a significant LRS score and at the gene level there was not a high LRS score, then the selected exon(s) was behaving differently from the other exons within the gene. Of the 2,314 exons with an LRS score above 60, 1,569 exons were part of a gene that did not have a high LRS score. An extensive evaluation of all these exons is beyond the scope of the present paper. Therefore, the top ten exons with LRS scores ranging from 165 to 202 were selected for further analysis. These exons were from ten genes: Cyb5r3, Hmgn2, Kif22, Col18a1, Uba2, Wdtc1, Haus5, Sdc2, Poc5, and Cntn1. In every case, at least one exon was differentially expressed between the C57BL/6J mouse and the DBA/2J mouse. To illustrate the differential splicing seen in these genes, we examined Col18a1 in depth. In the exon data set, there are 26 separate probes for the exons in the Col18a1 gene (see Table 3). Most of the probes have an associated LRS score ranging from 7 to 19. However, one probe (17242250) had an LRS score of 176.3. When we examined the distribution of the expression of this exon across the BXD RI strains, we saw that it is highly expressed in the C57BL/6J mouse retina relative to the DBA/2J retina (Figure 2). In other exon probes for Col18a1, there is a similar level of expression between the C57BL/6J mouse retina and the DBA/2J mouse retina, as well as the remaining BXD strains (Figure 3). It is possible that these differences in probe binding could be due to sequence variants between the C57BL/6J mouse and the DBA/2J. To test this hypothesis, we examined the sequence of the exons with high LRS scores to define the sequence differences between the two strains. Of the ten exons examined, five had single nucleotide polymorphisms (SNPs) within the region recognized by the Affymetrix probe. For these five exons, the difference in expression could be explained by differences in probe binding and not by differential splicing. For the remaining five exons, including that of Col18a1, the sequence in the C57BL/6J mouse was identical to that of the DBA/2J mouse. This type of analysis appears to identify differentially spliced genes in the retina of the BXD RI strain set. Our group is in the process of developing bioinformatic tools to take full advantage of the data from the Affymetrix exon chips. In the near future, this type of analysis may be as simple as a single query on the Trait Data and Analysis page of GeneNetwork.

Table 3

The probes for each of the exons of Col18a1 are presented.

Probe ID	Symbol	Description	Location (Chr, Mb)	Mean Expr	Max LRS
17,242,233	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.515054	7.706218147	12.8
17,242,235	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.516921	7.824909037	12.5
17,242,236	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.517574	7.407200033	11.2
17,242,238	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.521455	9.306981902	18.4
17,242,239	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.521899	8.679945391	11.9
17,242,242	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.523117	6.338236349	12
17,242,246	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.529745	6.939890905	14.2
17,242,248	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.531176	5.754218162	13.7
17,242,250	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.532565	9.958654568	176.3
17,242,251	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.532696	9.252999965	12.1
17,242,254	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.534773	7.823781802	19.3
17,242,255	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.536490	6.731690944	15.9
17,242,256	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.537103	8.810418181	14.5
17,242,259	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.540499	5.251781802	12.4
17,242,260	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.540746	6.876309109	10.9
17,242,261	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.541023	6.412472682	12
17,242,264	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.543517	9.105672802	14.4
17,242,266	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.548057	9.839090885	7
17,242,267	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.550127	10.48203631	11.3
17,242,269	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.552081	6.064109126	9
17,242,270	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.559097	9.295490941	9.2
17,242,271	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.575355	8.399236367	17.5
17,242,272	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.575935	8.588909123	17.4
17,242,273	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.576553	7.142545466	12.7
17,242,274	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.629063	8.688581805	8.5
17,242,275	Col18a1	collagen, type XVIII, alpha 1	Chr10: 76.629242	8.288054553	10.3

Notice that most of the probes have LRS ranging from 7 to 19. However one probe (17242250) has an LRS of 176.3. This latter exon is differentially spliced.

Figure 2

Expression of Col18a1 (Exon 14) is illustrated across the BXD strains in the DoD CDMRP Normal Retina Exon Level DataSet. The expression levels of Col18a1 at exon 14 are shown for many of the BXD strains as the mean expression and the standard error of the mean. The individual strain identifications are shown along the bottom, and the scale is log2. This difference in expression reflects the differential splicing of this exon.

Figure 3

Expression of Col18a1 (Exon 5) is illustrated across the BXD strains in the DoD CDMRP Normal Retina Exon Level Database. The expression levels of Col18a1 are shown for many of the BXD strains as the mean expression and the standard error of the mean. The individual strain identifications are shown along the bottom, and the scale is log2. Thus, this exon is not differentially spliced across the BXD RI strains.

Example of a functional network in the DoD CDMRP Retina Database

To illustrate the features of the new DoD CDMRP Retina Database, we chose a specific gene, Robo2 (roundabout homolog 2) and used it to demonstrate the analytical powers of the database and the bioinformatics tools associated with GeneNetwork. Robo2 is highly expressed in the retina with a mean value of 10.7 across the BXD strain set. The expression within individual strains varies from a low of 10.2 to a high of 11.1. This is a log2 scale and represents approximately a twofold difference in expression (Figure 4). When we examined the database for genes with a similar pattern of expression across the BXD strain set, we found a group of genes that are highly correlated with the expression pattern of Robo2 (Table 4). One example is the third correlate on the list, Ncam2 (Figure 5) with a value of 0.926. Even the 100th correlate on the list (Git1) has a high correlation (r=0.873) with Robo2 (see Supplemental Table 1).

Figure 4

Expression of Robo2 across the BXD strains in the DoD CDMRP Normal Retina Database. The expression levels of Robo2 are shown for many of the BXD strains as the mean expression and the standard error of the mean. The individual strain identifications are shown along the bottom of the plot, and the scale is log2. This variability from strain to strain indicates that the gene is differentially regulated by multiple genomic elements.

Table 4

A list of genes that are highly correlated to Robo2.

Symbol	Description	Location (Chr: Mb)	Mean Expr	Sample r
Robo2	roundabout homolog 2 (Drosophila)	Chr16: 73.892551	10.72	1.00
Cask	calcium	ChrX: 13.094206	10.91	0.94
Ncam2	neural cell adhesion molecule 2	Chr16: 81.200942	9.95	0.93
Gria3	glutamate receptor, ionotropic, AMPA3	ChrX: 38.754305	10.25	0.92
Lphn3	latrophilin 3	Chr5: 81.450227	10.25	0.92
Cnksr2	connector enhancer of kinase suppressor of Ras 2	ChrX: 154.259368	10.37	0.92
Clstn2	calsyntenin 2	Chr9: 97.344814	10.00	0.92
Gnaq	guanine nucleotide binding protein	Chr19: 16.207321	10.63	0.92
Dnajc6	DnaJ (Hsp40) homolog, subfamily C, member 6	Chr4: 101.169253	11.97	0.92
Slc8a1	solute carrier family 8 (sodium	Chr17: 81.772445	9.45	0.91
Dpysl2	dihydropyrimidinase-like 2	Chr14: 67.421701	13.25	0.91
Plch1	phospholipase C, eta 1	Chr3: 63.500156	10.96	0.91
Gria1	glutamate receptor, ionotropic, AMPA1 (alpha 1)	Chr11: 56.824889	9.77	0.91
Fam165b	family with sequence similarity 165, member B	Chr16: 92.301531	10.04	−0.91
Ncam1	neural cell adhesion molecule 1	Chr9: 49.310243	12.95	0.90
B4galt6	UDP-Gal:betaGlcNAc beta 1,4-galactosyltransferase	Chr18: 20.843100	10.76	0.90
Odz1	odd Oz	ChrX: 40.677132	8.81	0.90
Sort1	sortilin 1	Chr3: 108.087009	11.52	0.90
Ubqln2	ubiquilin 2	ChrX: 149.932775	11.29	0.90
Mtmr9	myotubularin related protein 9	Chr14: 64.142447	11.71	0.90
Gabbr2	gamma-aminobutyric acid (GABA) B receptor, 2	Chr4: 46.675190	11.04	0.90

The gene symbol, name of the gene, chromosome location, mean expression and correlation to Robo2 are listed across the top of the table.

Figure 5

The Pearson correlation between Robo2 and Ncam2. Ncam2 was the second highest correlate to Robo2 in the Department of Defense (DoD) Congressionally Directed Medical Research Programs (CDMRP) Normal Retina Database. These data indicate that these two genes are co-regulated with across the BXD strain set. When one gene is high in a strain the other gene is also expressed at a high level.

To define the regions of the genome that modulate the expression of Robo2, we plotted a genome-wide scan for Robo2 (Figure 6). This plot defines regions of the genome that correlate with the level of Robo2 expression, a quantitative trait locus (QTL). In this interval map, there is one significant QTL on chromosome 16 (notice the peak reaches the red line on the scan, p=0.05), and there are two suggestive peaks on chromosome 1 and chromosome 17 (above the gray line, p=0.63). The expression of Robo2 is modulated by genomic elements on chromosome 16. Two types of elements could affect the expression of Robo2: a cis-QTL or a gene with a nonsynonymous SNP. When we examined the significant QTL on chromosome 16 (21–27 Mb), we found there were no significant cis-QTLs at the gene level. With the DoD CDMRP Retina Database, it is now possible to look at the individual probes in exons and introns. When we checked the DoD CDMRP Retina Exon Level Database, we found one probe (Affy_17329472) that lies within the Leprel1 gene. When we checked the location of the probe with the Verify function on GeneNetwork, the probe lies in an intron and may be a non-coding RNA. However, when we examined the RNA-seq data from GeneNetwork, it appeared that this probe was detected in an RNA-seq analysis of the hippocampus and thus may be part of Leprel1 gene itself. Nonetheless, this probe marks a candidate for modulating the expression of Robo2. The second approach was to examine this region for nonsynonymous SNPs. Using the SNP browser in GeneNetwork, we looked at chromosome 16 (21–27 Mb) and found four known genes with nonsynonymous SNPs: Kng2, Kng1, BC106179, and Masp1. This analysis provided us with five candidates for modulating the expression of Robo2.

Figure 6

Genome-wide Interval Map of Robo2. This genome-wide graph displays the quantitative trait loci (QTL) distribution across the Department of Defense (DoD) Congressionally Directed Medical Research Programs (CDMRP) Normal Retina Database. The x-axis plots the locations of the QTLs that control the transcript expression. Positions are measured in megabases from chromosome 1 to chromosome X (1–2,600 Mb). The y-axis plots the “likelihood ratio statistic” (LRS). The significant levels of individual QTLs are color-coded. The red line represents a genome-wide significance level (p=0.05), and the gray line is suggestive. Notice the significant QTL on chromosome 16.

To determine whether this highly correlated set of genes in the Robo2 network have functional relationship(s), we used WebGestalt (WEB-based GEne SeT AnaLysis Toolkit) to examine the top 500 correlates of Robo2 to determine if there were specific functional transcript enrichments. The list of the top 500 correlates of Robo2 was enriched for several biologic processes (nervous system development, synaptic transmission, and neuron differentiations); molecular functions (enzyme binding, post synaptic density [PDZ] domain binding, inorganic cation transmembrane transporter, and metal ion transmembrane transporter activity); and cellular components (cell projection part, neuron projection, intracellular part, and axon genes). This type of analysis plays a critical role in many genetic networks, defining the functional role of the network. In this specific case, the analysis indicates that the Robo2 network is involved in axonal growth and neuronal development.

Discussion

This article announces the release of two new BXD retina databases on GeneNetwork. The first is at the gene level (DoD CDMRP Retina Affymetrix Mouse Gene 2.0 ST (May 15) RMA Gene-Level Database). The second data set is an exon-level analysis of the data presented in the first data set (DoD CDMRP Retina Affymetrix Mouse Gene 2.0 ST (May 15) RMA Exon-Level Database). Here, we emphasize some of the special aspects of these two data sets, including exon-level analysis and the inclusion of microRNAs and many non-coding RNAs. To illustrate many of these new features, we presented three approaches for analysis using the data sets.

The first was an examination of a cell signature within the data set. Within the DoD CDMRP Retina Database, there is a pronounced RPE signature. Some strains demonstrate low levels of expression of RPE65 while other strains have more than 16-fold higher levels of expression. This difference could not be due to differences in expression within the RPE. We believe that this difference is due to differences in the time of day when the retinas were isolated. The retinal samples were collected at two locations. At the University of Tennessee, the samples were usually isolated starting at 10:00 a.m., and lights were turned on in the animal colony at 6:00 a.m. Thus, the retinas were isolated at least 4 h after the lights came on. At Emory University, the retinas were isolated starting at 9:00 a.m., and lights on occurred at 7:00 a.m. (2 h after the lights came on in the animal colony). These differences in the time of day the retinas were isolated may be related to the number of RPE cells that adhered to the retinal samples [25].

Several bioinformatics tools are available to the vision research community. These tools include the NEI Bank project, which provides transcriptome profiling of the tissues of the eye, including mouse and human [26]. The Cepko group at Harvard, Boston, MA has provided the mouse retina serial analysis of gene expression (SAGE) library that includes gene expression of the embryonic and postnatal retina [27,28]. Daiger and his group at the University of Texas Health Science Center, Houston TX have lists of mapped loci and cloned genes associated with inherited retinal disease on the RetNet website. The Gene Expression Nervous System Atlas (GENSAT) now has a section devoted to the Retina Project [29] at GENSAT. The cell-specific labeling in the retina for different genes was illustrated using BAC transgenic mice [30]. The pattern of labeling in the retina defines the retinal cell types that express specific genes. This cellular localization aids in defining the localization of genetic networks in the retina. Finally, we posted the data from the study of glaucoma by Howell et al. [31] on the GeneNetwork website under the BXD eye database. These data are helpful in understanding the role of specific genetic networks in glaucoma (for example, see Templeton et al. [32]).

In conclusion, the DoD CDMRP Retina Database offered on GeneNetwork is a new resource for the vision community, in the ever-expanding variety of bioinformatics tools available. Previously, we offered several BXD microarray databases on GeneNetwork to the vision science community: the transcriptome of the whole eye (Eye M430v2 (September 2008) RMA Database) described in detail by Geisert et al. [18], a normal retina database (Normal Retina (April 2010) RankInv Database) described in detail by Freeman et al. [19], and the retina 2 days after optic nerve crush (ONC Retina (April 2012) RankInv Database) described by Templeton et al. [33]. This new data set offers a unique look at expression at the exon level. In addition, many non-protein coding transcripts are represented in the data set. The bioinformatics tools offered on GeneNetwork and these new databases are a unique resource for the vision research community.