Molecular Vision 2019; 25:345-358 <http://www.molvis.org/molvis/v25/345>
Received 20 February 2019 | Accepted 03 July 2019 | Published 05 July 2019

RNA sequencing profiling of the retina in C57BL/6J and DBA/2J mice: Enhancing the retinal microarray data sets from GeneNetwork

Jiaxing Wang,1 Eldon E. Geisert,1 Felix L. Struebing1,2,3

1Emory Eye Center, Department of Ophthalmology, Emory University, 1365B Clifton Road NE, Atlanta GA; 2Center for Neuropathology and Prion Research, Ludwig Maximilian University of Munich, Germany; 3Department for Translational Brain Research, German Center for Neurodegenerative Diseases (DZNE), Munich, Germany

Correspondence to: Felix L. Struebing, Center for Neuropathology and Prion Research, Ludwig Maximilian University of Munich, Germany; Feodor-Lynen-Str. 23; 81671 Munich, Germany; Phone: +49 89 2180 78028; FAX: +49 89 2180 78037; email: felix.struebing@med.uni-muenchen.de

Abstract

Purpose: The goal of the present study is to provide an independent assessment of the retinal transcriptome signatures of C57BL/6J (B6) and DBA/2J (D2) mice, and to enhance existing microarray data sets for accurately defining the allelic differences in the BXD recombinant inbred strains.

Methods: Retinas from B6 and D2 mice (three of each) were used for the RNA sequencing (RNA-seq) analysis. Transcriptome features were examined for both strains. Differentially expressed genes between the two strains were identified, and bioinformatic analysis was performed to analyze the transcriptome differences between the B6 and D2 strains, including Gene Ontology (GO) analysis, Phenotype and Reactome enrichment, and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. The RNA-seq data were then directly compared with one of the microarray data sets (Department of Defense [DoD] Retina Normal Affy MoGene 2.0 ST RMA Gene Level Microarray Database) hosted on GeneNetwork.

Results: RNA-seq provided an in-depth analysis of the transcriptome of the B6 and D2 retinas with a total of more than 30,000,000 reads per sample. More than 70% of the reads were uniquely mapped, resulting in a total of 18,100 gene counts for all six samples. A total of 1,665 genes were differentially expressed, with 858 of these more highly expressed in the B6 retinas and 807 more highly expressed in the D2 retinas. Several molecular pathways were differentially active between the two strains, including the retinoic acid metabolic process, endoplasmic reticulum lumen, extracellular matrix (ECM) organization, and the PI3K-Akt signaling pathway. The most enriched KEGG pathways were the pentose and glucuronate interconversions pathway, the cytochrome P450 pathway, the protein digestion and absorption pathway, and the ECM-receptor interaction pathway. Each of these pathways had a more than fourfold enrichment. The DoD Normal Retina Microarray Database provided expression profiling for 26,191 annotated transcripts for B6 mouse, D2 mouse, and 53 BXD strains. A total of 13,793 genes in this microarray data set were comparable to the RNA-seq data set. For the B6 and D2 retinas, the RNA-seq data and the microarray data were highly correlated with each other (Pearson's r=0.780 for the B6 mice and 0.784 for D2 mice). These results suggest that the microarray data set can reliably detect differentially expressed genes between the B6 and D2 retinas, with an overall accuracy of 91.1%. Examples of true positive and false positive genes are provided.

Conclusions: Retinal transcriptome features of B6 and D2 mouse strains provide a useful reference for a better understanding of the mouse retina. Generally, the microarray database presented on GeneNetwork shows good agreement with the RNA-seq data, but we note that any allelic difference between B6 and D2 mice should be verified with the latter.

Introduction

The development of high throughput RNA sequencing (RNA-seq) [1,2] is revolutionizing transcriptome analysis. It provides for a direct assessment of gene expression with relatively low background, a potentially unlimited dynamic range, and the ability to detect low abundance transcripts [2]. RNA-seq also provides opportunities to revisit, verify, and improve existing microarray data sets by reinforcing valid findings, and eliminating false positive results that can arise due to hybridization artifacts [3]. Large-scale transcriptome data sets using microarray technology were developed, and have been widely used for decades [4-6]. With the help of the retinal transcriptome microarray databases generated from BXD mouse strains hosted on GeneNetwork [7-11], our group has examined molecular networks contributing to various glaucoma-related phenotypes. We used the BXD strains and tools hosted on GeneNetwork to identify Aldh7a1 as a gene modulating retinal ganglion cell (RGC) death following elevated intraocular pressure (IOP) [12]. We also identified Cdh11 as a gene modulating normal IOP [13]. Finally, we identified Pou6f2 modulating central corneal thickness (CCT) in the mouse as a potential risk of glaucoma in humans [14]. Subsequently, two independent genome-wide association studies in humans identified CDH11 (Gene ID 1009; OMIM 600023) and POU6F2 (Gene ID 11281; OMIM 609062) as genes modulating IOP in humans, and as factors for glaucoma risk [15,16], supporting our findings, and proving the power of this systems biology approach and the microarray data sets hosted on GeneNetwork. The BXD strains represent a recombinant inbred strain set derived from C57BL/6J (B6) and DBA/2J (D2) parents, whose genomic variance forms the foundation of the BXD mice.

In the present study, we generated RNA-seq data from the retinas of B6 and D2 mice. The transcriptome features and differences were investigated. We also compared the RNA-seq data to our microarray data set (Affymetrix MouseGene 2.0 ST Array Gene Level) hosted on GeneNetwork. With the parental RNA-seq data, allelic differences in the BXD recombinant inbred strains that may arise due to hybridization artifacts could be precisely defined. The results suggest that a combination of microarray and RNA-seq data can lead to more valid results, compared to either data set alone.

Methods

Animals

For the RNA-seq analysis, a total of six samples were analyzed, three samples from B6 mice (two males and one female) and three samples from D2 mice (two males and one female). All the animals were between 60 and 100 days of age. Mice were maintained on a 12 h:12 h light-dark cycle in a parasite-free facility with food and water ad libitum. All procedures involving animals were approved by the Animal Care and Use Committee of Emory University, and were in accordance with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research.

RNA extraction

Mice were deeply anesthetized with 15 mg/kg of xylazine and 100 mg/kg of ketamine, and then euthanized with rapid cervical dislocation. Retinas from both eyes were quickly removed under a dissection microscope and directly placed into 160 U/ml Ribolock® (Thermo Scientific, Walton, MA) in Hank’s Balanced Salt Solution (Sigma, St. Louis, MO) on ice. Tissue was stored at −80 °C. RNA was isolated in batches using a Qiacube and the RNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The isolation included on-column DNase1 treatment to remove contaminating genomic DNA. All tissue was harvested between 10 AM and 12 PM to minimize circadian differences in gene expression. RNA integrity (RIN) was assessed on a Bioanalyzer 2100 (Agilent, Santa Clara, CA). The RIN score for all samples ranged from 9.1 to 9.5. RNA was then quantified with spectrophotometry, and 260/280 ratios for all samples were >2.1.

The retinas were dissected from the eye, and during the process, there are traces of RPE still attached to the retina. These traces were small, and to reveal the lack of significant RPE contamination, we examined the expression of RPE genes in the samples. The expression of genes in Figure 1 represents several different ocular cell types, including RPE. These data can be used to determine the makeup of the samples and the low level of RPE contamination.

RNA sequencing and data processing

Six samples of purified retinal RNA underwent library preparation with the Illumina (San Diego, CA) TruSeq Stranded Total RNA kit. Paired-end sequencing (100 bp) was performed by Hudson Alpha (Huntsville, AL). FastQC of the files was performed for quality control. Reads from each replicate were not trimmed before mapping, as this was recently shown to decrease the ability to map, and correlation with microarray data [17]. To avoid mapping issues arising from mouse strain heterogeneity, separate genome indices were built for the B6 and the D2 mouse genome. This was achieved by manually substituting the mm10 B6 reference genome (Ensembl version 86) with D2-specific single nucleotide polymorphisms (SNPs) and indels (from dbSNP 142). The reads were then mapped to each genome using the two-pass mapping mode in RNA-STAR with default parameters. For each sample, strand-specific bam files were sorted and indexed with samtools. Counts per gene were tabulated with the “summarizeOverlaps” internal R function from uniquely mapped reads only. Differential expression analysis was performed using a generalized linear model in edgeR [18] after any genes with fewer than 0.5 counts per million (cpm) in more than three samples were removed. All p values were adjusted using the false discovery rate (FDR). Statistical significance was assumed for genes with FDR<0.001 and a twofold change between strains. For comparison with the microarray data, log2-transformed cpm counts were z-scaled, and multiplied by 2, before a constant of 8 was added to avoid negative values. This is essentially the same normalization procedure that was applied to the microarray data sets hosted on GeneNetwork [6,19,20]. All raw data are deposited in the Gene Expression Omnibus (GEO), accession GSE127942.

Comparison between the RNA-seq and microarray data sets

Correlation tests were performed on genes extracted from the RNA-seq and microarray data sets. Gene symbols were used for the comparisons. In cases where there were duplicate gene symbols represented by different probes in the microarray annotation file, an average value of all probes capturing different parts of the transcripts was calculated to represent the expression level of the gene. The absolute differential expression values for each gene between the two data sets were calculated as the expression level of B6 minus D2 [21]. When the RNA-seq and microarray data sets were compared in detecting differentially expressed genes, the original probe value was used instead of the average value of all probes for a gene.

Bioinformatics

The differentially expressed genes between the B6 and D2 retinas were submitted for enrichment analysis to databases, including Gene Ontology (GO) [22,23], MGI mammalian phenotype ontology [24], Human Phenotype Ontology [25], and Reactome biologic pathways [26] with the R package EnrichR [27,28]. The top hits of the package output ‘combined scores’ and adjusted p values were used to determine significance. The differentially expressed genes were also submitted for Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis [29-31]. Redundancy control was run for each database. For GO analysis, “REVIGO” [32] was used for redundancy reduction, before the results were plotted. For all other enrichment analysis databases, “ReCiPa” [33] was used to reduce redundancy. The redundant pathways were combined into the pathway with a broader spectrum of genes and greater significance.

Results

RNA-seq representation of the retinal transcriptome

Quality evaluation (FastQC) of the RNA-seq files showed that >90% of reads had a Phred quality score of >30 (Appendix 1). Each replicate had >30 million raw reads, of which more than 70% could be uniquely mapped, resulting in a total of 18,100 gene counts for all three samples in the B6 group and three samples for the D2 group. Any genes with fewer than 0.5 cpm in the combined data were removed. Because the retina is a highly heterogeneous tissue composed of various cell types, this cutoff was chosen to include RNA from less common cells in the retina, such as specific RGC subtypes. As the first survey of the data quality, we examined the expression levels of signature genes for different retinal cell types (four of each), including RGCs, amacrine cells, bipolar cells, Müller glia, and photoreceptors [11,34-36] (Figure 1). These genes were compared with marker genes of the cornea and RPE, which should show negligible expression in clean retinal tissues [11]. Commonly used housekeeping genes [37] were also examined as reference control (Figure 1). Most of the retinal cell type-specific signature genes were highly expressed, compared with the housekeeping genes, especially the gene for rhodopsin (Rho) and other photoreceptor markers. The retinas were dissected from the eye, and during the process, there were traces of RPE still attached to the retina. These traces were small, and to reveal the lack of significant RPE contamination, we examined the expression of RPE signature genes in the samples. The expression of genes in Figure 1 represents several different ocular cell types, including RPE. These data are used to determine the makeup of the samples and the low level of RPE contamination.

We also examined two genes with known mutations in the DBA/2J strain contributing to eye phenotypes, tyrosinase-related protein 1 (Tyrp1) and glycoprotein non-metastatic melanoma protein B (Gpnmb) [38,39]. Gpnmb showed a difference in expression between the two strains, while Tyrp1 and the other listed genes were similarly expressed between the B6 and D2 retinas (Figure 1).

Among the housekeeping genes examined, Actb, Gapdh, and Ppia were equally expressed across samples (low standard error of the mean [SEM] levels); however, Rn18s (Encoding 18S rRNA) displayed a surprising variability (Figure 1). This high variability of Rn18s is known to be associated with cell differentiation, forming organs and their maturation [37]. These results suggest that the Actb, Gapdh, and Ppia are better choices as housekeeping genes for the retina.

Transcriptome differences reveal unique features between the B6 and D2 retinas

Differentially expressed genes between the B6 and D2 retinas were extracted from RNA-seq data, and the p values were adjusted for multiple comparisons with the false discovery rate [40]. We found 1,665 genes that were differentially expressed, with 858 of these more highly expressed in the B6 retinas and 807 more highly expressed in the D2 retinas (Appendix 2). The differentially expressed genes were submitted for enrichment analysis, including GO [22,23], Phenotypes in both Mouse [24] and Human [25], and Reactome [26] enrichment using EnrichR [27,28]. The top five pathways for each enrichment are shown in Figure 2. The pathways with statistical significance included the retinoic acid metabolic process and extracellular matrix organization in the GO biologic process (Figure 2A), the endoplasmic reticulum lumen in the GO cellular component (Figure 2C), and the extracellular matrix organization and collagen formation in the Reactome database (Figure 2F). The genes involved in each pathway are summarized in Appendix 3 and Appendix 4.

The differentially expressed genes were also submitted for KEGG analysis, in which the top 30 pathways are listed (Figure 3). The greatest number of differentially expressed genes fell into the metabolic pathway (n=107). The second most involved pathway was the PI3K-Akt signaling pathway (n=41, Appendix 5). There were 19 genes involved in the retinol metabolism pathway, with 16 genes expressed at a higher level in the B6 retinas and only two genes (Rdh9 and Cyp4a12b) showing higher expression in the D2 retinas (Appendix 6). The most enriched pathways were the pentose and glucuronate interconversions, the drug metabolism - cytochrome P450, the ECM-receptor interaction, alpha-linolenic acid metabolism, and protein digestion and absorption, with more than fourfold enrichment for each pathway (Figure 3, Appendix 7- Appendix 8-Appendix 9, Appendix 10, and Appendix 11).

Comparison of RNA-seq and microarray data

RNA sequencing directly quantifies read counts, while microarray technology relies on the binding of a transcript to a predesigned probe sequence. To investigate these technical differences, and understand how they might be reflected in the mouse retinal transcriptome, we contrasted the RNA-seq data with the Department of Defense (DoD) Congressionally Directed Medical Research Programs (CDMRP) Normal Retina Database, which is publicly available on GeneNetwork. This previously established data set uses the Affymetrix MouseGene 2.0 ST Array, which provides expression profiling for 26,191 well-established annotated transcripts, 9,049 non-coding RNAs, and more than 600 microRNAs [8]. Details about the data processing can be found in King et al. [8]. This data set consists of 52 BXD strains, as well as B6, D2, and an F1 cross between B6 and D2 (B6D2F1). A total of 222 microarrays were used to describe the retinal transcriptomes of these combined 55 mouse strains, which were usually represented by four biologic replicates per strain.

A total of 13,793 genes (represented by 16,445 probe sets) in the microarray data set were also identified in the RNA-seq data set (Appendix 12). In the microarray database, some genes were represented by multiple annotated probe sets. When we encountered this situation, the mean value of all probe sets was used to represent the expression level of the gene. The RNA-seq data and the microarray data showed a high correlation; the Pearson’s r value for the B6 retina was 0.780 (Figure 4A) and for the D2 retina was 0.784 (Figure 4B). Thus, there was good agreement between the two data sets, suggesting that gene expression values from either RNA-seq or microarrays are reliable and comparable, except a few genes with either high or low expression.

Identification of differentially expressed genes

The absolute difference in the mean expression of genes between the B6 and D2 retinas from the microarray data set was calculated (see Methods and Appendix 13), and compared to the RNA-seq data, as shown in Figure 5. Pearson’s r value was 0.408 for the differentially expressed values of all genes (p<0.001), indicating moderate agreement between the filtered data sets. When the expression level with differences above 0.5 on a log2 scale (|B6-D2| >0.5) was considered the threshold for the microarray data, 848 genes were identified as differentially expressed, of which 387 (45.6%) were confirmed by the RNA-seq data as true positive, while 461 (54.4%) were identified as false positive (Figure 5). Simultaneously, 15,597 genes did not show differential expression in the microarray data set, of which 14,598 (93.6%) genes were confirmed as true negative, and 999 (6.4%) genes were identified as false negative by the RNA-seq data (Figure 5). Thus, the sensitivity and specificity of detecting differentially expressed genes using the microarray data set were 27.9% and 96.9%, respectively. The precision, i.e., the positive predictive value of the microarray data in detecting differentially expressed genes, was 45.6%, and the negative predictive value was 93.6%. The overall accuracy of the microarray data sets relative to the RNA-seq data was 91.1%.

Although the microarray methodology might overestimate the number of differentially expressed genes, the results suggest that the retinal transcriptome is generally well represented with both technologies. Expression values from microarrays depend on probe hybridization, which can be affected by differences in the genome sequence observed between strains. Thus, genomic variants overlapping the probe targeting area could directly affect the affinity of the probe to the transcript [41,42]. RNA-seq is not dependent upon hybridization of a target to a probe, and thus, this technology can serve to validate the results obtained using a microarray system.

An example for validation of the microarray data by RNA-seq is Stab2. In the microarray data set, clear differential expression of the gene Stab2 was identified between the parental strains, as well as across the BXD strains (Probe_17243763, Max LRS=27). The likelihood ratio statistic (LRS) is used as a measurement of the association or linkage between differences in RNA expression [6], where values exceeding 17 usually approximate statistical significance at the p=0.05 level after correction for multiple testing. The Stab2 probe coverage area contained no SNP, and it was also identified as a differentially expressed gene in the RNA-seq data (Figure 6).

In another case, a clear differential expression of the gene Anapc1 was identified by the microarray data set (Probe_17391461, Max LRS=124.5). However, the RNA-seq data showed no statistically significant differential expression for this gene. Probe verification using the UCSC Genome Browser showed an SNP (rs27412956) in the probe coverage area that may lead to mis-binding toward the D2 sequence. This potentially explains why Anapc1 showed lower expression levels in the D2 mice compared to the B6 mice in the microarray data set (Figure 7), suggesting that an SNP in a probe hybridization locus alters probe binding. Moreover, this mis-binding affected not only the D2 parental strains but also all the BXD strains with D2-derived alleles for this genomic region, resulting in a false-high LRS. In such cases, the RNA-seq data set is helpful in identifying false positive results. Similar false positive cases caused by D2 SNPs were seen in a total of 211 genes (222 probes, Appendix 14).

Discussion

Retinal transcriptome features of B6 and D2 mice

In the present study, we examined the retinal transcriptome of two different mouse strains, C57BL/6J and DBA/2J. Although it is known that there are clear sex differences in gene expression, we chose mixed genders for this study, as this most closely mirrored the sex distribution of BXD mice from the GeneNetwork databases [43]. The B6 mouse is the most widely used mouse strain in research. In many ways, this strain serves as the ultimate inbred mouse strain, not only being the background strain for most transgenic lines but also functioning as the reference genome for the mouse. In contrast, the D2 mouse holds a specific interest for the neuroscience and the vision research community, as this strain is used in behavioral addiction research [44], and develops pigment dispersion glaucoma with onset at approximately 6 months of age [39,45]. When we examined the retinal transcriptome between these two strains using RNA-seq, several differences became apparent. Two mutations are known to be responsible for the DBA/2J eye phenotype: GpnmbR150X, a nonsense mutation, and Tyrp1isa, an allele leading to two amino acid substitutions. However, at the mRNA level, this must not necessarily mean that these two genes are differentially expressed, because we found differential expression only for Gpnmb. Although more research on this topic is warranted, it can be speculated that the nonsense mutation in Gpnbm leads to increased nonsense-mediated decay (NMD), and therefore, early mRNA removal by the NMD machinery [46], which can be seen as different expression between D2 and B6 mice. In contrast, the amino acid substitutions in the Tyrp1 gene do not lead to early transcriptional termination, and therefore, are tolerated at the mRNA level.

In the Gene Ontology analysis for differentially expressed genes, the most significantly enriched pathway was the retinoic acid metabolic process (Figure 2A). The main contributing genes (CYP26A1, CYP26B1, and CYP26C1) come from the CYP26 family, also known as retinoic acid hydroxylase. All three genes are statistically significantly higher expressed in B6 mice than in D2 mice. These genes are also the main contributors to the significant enrichment of the cytochrome P450 pathway and the retinol metabolism pathway in KEGG analysis (Figure 3). Another significantly enriched pathway is the ECM organization, which is also the only significantly enriched pathway of the Reactome analysis (Figure 2A,F). The ECM is known to play multiple roles in the retina, including retinal development [47,48], regulation of retinoid transport, oxygen and nutrient transport [49], vascular homeostasis [50], glia-neuron interactions [51], and even regeneration of RGCs [52]. This difference in the ECM between the B6 and D2 retinas may potentially contribute to the differences in RGC survival [53], and other retinal phenotypes between the two strains, which makes the biologic approach using BXD strains (offspring inbred strains of B6 and D2) more powerful in revealing genetics of retinal diseases, such as glaucoma and age-related macular degeneration.

Pathway analysis using KEGG demonstrated differences in the PI3K-Akt signaling pathway, with 54 genes that were differentially expressed (Figure 3). The PI3K-Akt signaling pathway is an important intracellular signaling pathway for cell cycle regulation. Overactivation of this pathway can reduce apoptosis, allow proliferation, and is associated with specific cancers [54]. This difference in the PI3K-Akt signaling pathway could affect changes in the ability of axons to regenerate. In the retina, Pten silencing is known to increase retinal ganglion cell survival, and promote axonal regeneration after injury [55]. There is a strong difference in the ability of optic nerve regeneration across the BXD strains [56]. Part of this may be due to transcriptome variance between the parental strains (B6 and D2) in this PI3K-Akt signaling pathway, for it is known that AKT is involved in the effects of Pten knockdown [57]. An interesting feature of the Pten-dependent regeneration response is that the parental strains (B6 and D2) are not at the extremes of the phenotype [56]. Thus, axon regeneration displays genetic transgression, and is a complex trait modulated by multiple genomic loci in the BXD strain set.

Furthermore, pathways involved in the innate immune system were also identified in the KEGG pathway analysis (Appendix 15). Many contributing genes were related to the complement system, including C1qb, C1qc, C1ra, C3, and C4b. The retina-intrinsic innate immunity network was reported to be activated by various types of retinal injures, including optic nerve crush [9] and ocular blast injury [7]. The related genes were mostly higher expressed in B6 mice rather than D2 mice (Appendix 7), indicating a greater innate immune response in the retina of B6 mice compared to that of D2 mice.

RNA-seq profiling enhances the transcriptome analysis of microarray data sets from GeneNetwork

The BXD strains are derived from the B6 and D2 progenitor strains, and their genomic variance allows us to map quantitative trait loci (QTLs) with high resolution, and economically, using microarray technology. Although this approach provides large amounts of information, it is not without limitations: Because the core principle behind microarrays is hybridization between two nucleic acid strands, the signal depends on the probe binding to the target sequence. Genomic variants between the B6 and D2 retinas, such as SNPs and indels, may influence hybridization kinetics, and cause incorrect detection of the expression level of genes [3,58]. A sequence that contains a single SNP may show a decreased binding affinity toward the probe, resulting in a less intense fluorescence signal, and thus, a lower readout for the probe capturing that particular transcript [42].

RNA-seq has developed rapidly over the past 10 years, and has now replaced the use of microarrays in transcriptomics [1,2,59,60]. The inherently lower background noise and higher dynamic range make RNA-seq more accurate in detecting low abundance transcripts [2]. In this study, the RNA-seq data from the B6 and D2 mice showed high consistency when compared with the DoD Normal Retina Microarray Data set (Figure 4). Beyond that, a moderate to high correlation of the differentially expressed genes between the two data sets suggested that the microarray data set still performs well at detecting differences between different strains of mice (Figure 5).

Ultimately, the combination of the technologies benefits analysis, especially for the study of recombinant inbred strains: Microarrays are more affordable for high throughput experiments, while RNA-seq analysis of the parental strains can limit the high false positive QTL rate. Note that not all the SNPs have to cause mis-binding between a microarray probe and the cDNA. Microarrays that feature arbitrarily designed mismatch probes still detect approximately 70% of perfect match signals, indicating that hybridization is not a binary measure, compared to RNA-seq read counts [41]. Thus, cross-checking results obtained from microarray data in recombinant inbred populations with RNA-seq data from the respective parental strains refines the analysis. After all, the genome sequences of the BXD strains are inherited from either the B6 or D2 mice. If a differentially expressed gene cannot be confirmed with RNA-seq data, and meanwhile, all the BXD strains with a D2 haplotype at that locus show lower expression, then there is a high chance for this finding to be false positive. If differential expression is not confirmed with the parental RNA-seq data, and the differential expression in BXD strains is not haplotype-dependent, then there are good chances that the finding is valid. These QTLs could be the consequence of inherited modulators in cis or trans (enhancers, suppressors, or other regulatory elements) that affect gene expression. Ultimately, if differential expression is confirmed by the RNA-seq data, then we can be more confident about a real change in gene expression that is allele-based.

Conclusions

Differential expression analysis of B6 and D2 mice using RNA-seq provided a useful reference for better understanding of the biologic differences between these strains. Although the majority of expression differences is mirrored in the microarray data from BXD strains, RNA-seq complements these data sets. Any allelic difference between B6 and D2 mice should be verified with RNA-seq data to avoid false positive QTLs.

Appendix 1. Fastqc plot of retina samples for RNA-seq analysis.

Appendix 2. RNAseq_B6D2-zscaled_Retina

Appendix 3. Enrichlist_DB

Appendix 4. KEGG LIST OF DEGs

Appendix 5. DE genes involvement of the PI3K-Akt signaling pathway.

Appendix 6. DE genes involvement of the Retinol metabolism pathway.

Appendix 7. DE genes involvement of the Pentose and glucuronate interconversions pathway.

Appendix 8. DE genes involvement of the Drug metabolism - cytochrome P450 pathway.

Appendix 9. DE genes involvement of the Alpha-Linolenic acid metabolism pathway.

Appendix 10. DE genes involvement of the ECM-receptor interaction pathway.

Appendix 11. DE genes involvement of the Protein digestion and absorption pathway.

Appendix 12. RNAseq_microarray_comparison_B6D2_Retina

Appendix 13. RNAseq_microarray_one-to-many-comparison_B6D2_Retina

Appendix 14. False positives caused by D2 SNPs

Appendix 15. DE genes involvement of the staphylococcus aureus infection pathway

Acknowledgments

We thank Rebecca King for her technical assistance with the mice and RNA isolation. This work was supported in part by the Vision Core Grant P030EY006360, the Emory Integrated Genomics Core (EIGC) Shared Resource of Winship Cancer Institute of Emory University and NIH/NCI under award number P30CA138292, the Owens Family Glaucoma Fund, and an Unrestricted Grant from Research to Prevent Blindness. FS is supported by the European Union’s Horizon 2020 Research and Innovation Program (grant agreement MSCA-IF-EF-RI 792832).

References

  1. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, Snyder M. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science. 2008; 320:1344-9. [PMID: 18451266]
  2. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009; 10:57-63. [PMID: 19015660]
  3. Pandey AK, Williams RW. Genetics of gene expression in CNS. Int Rev Neurobiol. 2014; 116:195-231. [PMID: 25172476]
  4. Kolesnikov N, Hastings E, Keays M, Melnichuk O, Tang YA, Williams E, Dylag M, Kurbatova N, Brandizi M, Burdett T, Megy K, Pilicheva E, Rustici G, Tikhonov A, Parkinson H, Petryszak R, Sarkans U, Brazma A. ArrayExpress update–simplifying data submissions. Nucleic Acids Res. 2015; 43:D1113-6. [PMID: 25361974]
  5. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Holko M, Yefanov A, Lee H, Zhang N, Robertson CL, Serova N, Davis S, Soboleva A. NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res. 2013; 41:D991-5. [PMID: 23193258]
  6. Mulligan MK, Mozhui K, Prins P, Williams RW. GeneNetwork: A Toolbox for Systems Genetics. Methods Mol Biol. 2017; 1488:75-120. [PMID: 27933521]
  7. Struebing FL, King R, Li Y, Chrenek MA, Lyuboslavsky PN, Sidhu CS, Iuvone PM, Geisert EE. Transcriptional Changes in the Mouse Retina after Ocular Blast Injury: A Role for the Immune System. J Neurotrauma. 2018; 35:118-29. [PMID: 28599600]
  8. King R, Lu L, Williams RW, Geisert EE. Transcriptome networks in the mouse retina: An exon level BXD RI database. Mol Vis. 2015; 21:1235-51. [PMID: 26604663]
  9. Templeton JP, Freeman NE, Nickerson JM, Jablonski MM, Rex TS, Williams RW, Geisert EE. Innate immune network in the retina activated by optic nerve crush. Invest Ophthalmol Vis Sci. 2013; 54:2599-606. [PMID: 23493296]
  10. Freeman NE, Templeton JP, Orr WE, Lu L, Williams RW, Geisert EE. Genetic networks in the mouse retina: growth associated protein 43 and phosphatase tensin homolog network. Mol Vis. 2011; 17:1355-72. [PMID: 21655357]
  11. Geisert EE, Lu L, Freeman-Anderson NE, Templeton JP, Nassr M, Wang X, Gu W, Jiao Y, Williams RW. Gene expression in the mouse eye: an online resource for genetics using 103 strains of mice. Mol Vis. 2009; 15:1730-63. [PMID: 19727342]
  12. Struebing FL, King R, Li Y, Cooke Bailey JN. consortium N, Wiggs JL, Geisert EE. Genomic loci modulating retinal ganglion cell death following elevated IOP in the mouse. Exp Eye Res. 2018; 169:61-7. [PMID: 29421330]
  13. King R, Li Y, Wang J, Struebing FL, Geisert EE. Genomic Locus Modulating IOP in the BXD RI Mouse Strains. G3 (Bethesda). 2018; 8:1571-8. [PMID: 29496776]
  14. King R, Struebing FL, Li Y, Wang J, Koch AA, Cooke Bailey JN, Gharahkhani P. International Glaucoma Genetics C, Consortium N, MacGregor S, Allingham RR, Hauser MA, Wiggs JL, Geisert EE. Genomic locus modulating corneal thickness in the mouse identifies POU6F2 as a potential risk of developing glaucoma. PLoS Genet. 2018; 14:e1007145 [PMID: 29370175]
  15. MacGregor S, Ong JS, An J, Han X, Zhou T, Siggs OM, Law MH, Souzeau E, Sharma S, Lynn DJ, Beesley J, Sheldrick B, Mills RA, Landers J, Ruddle JB, Graham SL, Healey PR, White AJR, Casson RJ, Best S, Grigg JR, Goldberg I, Powell JE, Whiteman DC, Radford-Smith GL, Martin NG, Montgomery GW, Burdon KP, Mackey DA, Gharahkhani P, Craig JE, Hewitt AW. Genome-wide association study of intraocular pressure uncovers new pathways to glaucoma. Nat Genet. 2018; 50:1067-71. [PMID: 30054594]
  16. Khawaja AP, Cooke Bailey JN, Wareham NJ, Scott RA, Simcoe M, Igo RP, , Jr Song YE, Wojciechowski R, Cheng CY, Khaw PT, Pasquale LR, Haines JL, Foster PJ, Wiggs JL, Hammond CJ, Hysi PG, Eye UKB, Vision C, Consortium N. Genome-wide analyses identify 68 new loci associated with intraocular pressure and improve risk prediction for primary open-angle glaucoma. Nat Genet. 2018; 50:778-82. [PMID: 29785010]
  17. Williams CR, Baccarella A, Parrish JZ, Kim CC. Trimming of sequence reads alters RNA-Seq gene expression estimates. BMC Bioinformatics. 2016; 17:103 [PMID: 26911985]
  18. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010; 26:139-40. [PMID: 19910308]
  19. Williams RW, Williams EG. Resources for Systems Genetics. Methods Mol Biol. 2017; 1488:3-29. [PMID: 27933518]
  20. Chesler EJ, Lu L, Shou S, Qu Y, Gu J, Wang J, Hsu HC, Mountz JD, Baldwin NE, Langston MA, Threadgill DW, Manly KF, Williams RW. Complex trait analysis of gene expression uncovers polygenic and pleiotropic networks that modulate nervous system function. Nat Genet. 2005; 37:233-42. [PMID: 15711545]
  21. Bottomly D, Walter NA, Hunter JE, Darakjian P, Kawane S, Buck KJ, Searles RP, Mooney M, McWeeney SK, Hitzemann R. Evaluating gene expression in C57BL/6J and DBA/2J mouse striatum using RNA-Seq and microarrays. PLoS One. 2011; 6:e17820 [PMID: 21455293]
  22. The Gene Ontology C. Expansion of the Gene Ontology knowledgebase and resources. Nucleic Acids Res. 2017; 45D1:D331-8. [PMID: 27899567]
  23. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000; 25:25-9. [PMID: 10802651]
  24. Smith CL, Blake JA, Kadin JA, Richardson JE, Bult CJ, Mouse Genome Database G. Mouse Genome Database (MGD)-2018: knowledgebase for the laboratory mouse. Nucleic Acids Res. 2018; 46D1:D836-42. [PMID: 29092072]
  25. Kohler S, Vasilevsky NA, Engelstad M, Foster E, McMurry J, Ayme S, Baynam G, Bello SM, Boerkoel CF, Boycott KM, Brudno M, Buske OJ, Chinnery PF, Cipriani V, Connell LE, Dawkins HJ, DeMare LE, Devereau AD, de Vries BB, Firth HV, Freson K, Greene D, Hamosh A, Helbig I, Hum C, Jahn JA, James R, Krause R. SJ FL, Lochmuller H, Lyon GJ, Ogishima S, Olry A, Ouwehand WH, Pontikos N, Rath A, Schaefer F, Scott RH, Segal M, Sergouniotis PI, Sever R, Smith CL, Straub V, Thompson R, Turner C, Turro E, Veltman MW, Vulliamy T, Yu J, von Ziegenweidt J, Zankl A, Zuchner S, Zemojtel T, Jacobsen JO, Groza T, Smedley D, Mungall CJ, Haendel M, Robinson PN. The Human Phenotype Ontology in 2017. Nucleic Acids Res. 2017; 45D1:D865-76. [PMID: 27899602]
  26. Fabregat A, Jupe S, Matthews L, Sidiropoulos K, Gillespie M, Garapati P, Haw R, Jassal B, Korninger F, May B, Milacic M, Roca CD, Rothfels K, Sevilla C, Shamovsky V, Shorser S, Varusai T, Viteri G, Weiser J, Wu G, Stein L, Hermjakob H, D’Eustachio P. The Reactome Pathway Knowledgebase. Nucleic Acids Res. 2018; 46D1:D649-55. [PMID: 29145629]
  27. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, Clark NR, Ma’ayan A. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013; 14:128 [PMID: 23586463]
  28. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, McDermott MG, Monteiro CD, Gundersen GW, Ma’ayan A. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016; 44W1:W90–7 [PMID: 27141961]
  29. Kanehisa M, Sato Y, Furumichi M, Morishima K, Tanabe M. New approach for understanding genome variations in KEGG. Nucleic Acids Res. 2019; 47D1:D590-5. [PMID: 30321428]
  30. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017; 45D1:D353-61. [PMID: 27899662]
  31. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28:27-30. [PMID: 10592173]
  32. Supek F, Bosnjak M, Skunca N, Smuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS One. 2011; 6:e21800 [PMID: 21789182]
  33. Vivar JC, Pemu P, McPherson R, Ghosh S. Redundancy control in pathway databases (ReCiPa): an application for improving gene-set enrichment analysis in Omics studies and “Big data” biology. OMICS. 2013; 17:414-22. [PMID: 23758478]
  34. Shekhar K, Lapan SW, Whitney IE, Tran NM, Macosko EZ, Kowalczyk M, Adiconis X, Levin JZ, Nemesh J, Goldman M, McCarroll SA, Cepko CL, Regev A, Sanes JR. Comprehensive Classification of Retinal Bipolar Neurons by Single-Cell Transcriptomics. Cell 2016; 166(5):1308–23 e30.
  35. Welby E, Lakowski J, Di Foggia V, Budinger D, Gonzalez-Cordero A, Lun ATL, Epstein M, Patel A, Cuevas E, Kruczek K, Naeem A, Minneci F, Hubank M, Jones DT, Marioni JC, Ali RR, Sowden JC. Isolation and Comparative Transcriptome Analysis of Human Fetal and iPSC-Derived Cone Photoreceptor Cells. Stem Cell Reports. 2017; 9:1898-915. [PMID: 29153988]
  36. Gamsiz ED, Ouyang Q, Schmidt M, Nagpal S, Morrow EM. Genome-wide transcriptome analysis in murine neural retina using high-throughput RNA sequencing. Genomics. 2012; 99:44-51. [PMID: 22032952]
  37. Kozera B, Rapacz M. Reference genes in real-time PCR. J Appl Genet. 2013; 54:391-406. [PMID: 24078518]
  38. Chang B, Smith RS, Hawes NL, Anderson MG, Zabaleta A, Savinova O, Roderick TH, Heckenlively JR, Davisson MT, John SW. Interacting loci cause severe iris atrophy and glaucoma in DBA/2J mice. Nat Genet. 1999; 21:405-9. [PMID: 10192392]
  39. Anderson MG, Smith RS, Hawes NL, Zabaleta A, Chang B, Wiggs JL, John SW. Mutations in genes encoding melanosomal proteins cause pigmentary glaucoma in DBA/2J mice. Nat Genet. 2002; 30:81-5. [PMID: 11743578]
  40. Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc Natl Acad Sci USA. 2003; 100:9440-5. [PMID: 12883005]
  41. Deng Y, He Z, Van Nostrand JD, Zhou J. Design and analysis of mismatch probes for long oligonucleotide microarrays. BMC Genomics. 2008; 9:491 [PMID: 18928550]
  42. Wang Y, Miao ZH, Pommier Y, Kawasaki ES, Player A. Characterization of mismatch and high-signal intensity probes associated with Affymetrix genechips. Bioinformatics. 2007; 23:2088-95. [PMID: 17553856]
  43. Du M, Mangold CA, Bixler GV, Brucklacher RM, Masser DR, Stout MB, Elliott MH, Freeman WM. Retinal gene expression responses to aging are sexually divergent. Mol Vis. 2017; 23:707-17. [PMID: 29062222]
  44. Mulligan MK, Mozhui K, Pandey AK, Smith ML, Gong S, Ingels J, Miles MF, Lopez MF, Lu L, Williams RW. Genetic divergence in the transcriptional engram of chronic alcohol abuse: A laser-capture RNA-seq study of the mouse mesocorticolimbic system. Alcohol. 2017; 58:61-72. [PMID: 27894806]
  45. John SW, Smith RS, Savinova OV, Hawes NL, Chang B, Turnbull D, Davisson M, Roderick TH, Heckenlively JR. Essential iris atrophy, pigment dispersion, and glaucoma in DBA/2J mice. Invest Ophthalmol Vis Sci. 1998; 39:951-62. [PMID: 9579474]
  46. Baker KE, Parker R. Nonsense-mediated mRNA decay: terminating erroneous gene expression. Curr Opin Cell Biol. 2004; 16:293-9. [PMID: 15145354]
  47. Reinhard J, Joachim SC, Faissner A. Extracellular matrix remodeling during retinal development. Exp Eye Res. 2015; 133:132-40. [PMID: 25017186]
  48. Bishop PN. The role of extracellular matrix in retinal vascular development and preretinal neovascularization. Exp Eye Res. 2015; 133:30-6. [PMID: 25819452]
  49. Ishikawa M, Sawada Y, Yoshitomi T. Structure and function of the interphotoreceptor matrix surrounding retinal photoreceptor cells. Exp Eye Res. 2015; 133:3-18. [PMID: 25819450]
  50. Roy S, Bae E, Amin S, Kim D. Extracellular matrix, gap junctions, and retinal vascular homeostasis in diabetic retinopathy. Exp Eye Res. 2015; 133:58-68. [PMID: 25819455]
  51. Vecino E, Rodriguez FD, Ruzafa N, Pereiro X, Sharma SC. Glia-neuron interactions in the mammalian retina. Prog Retin Eye Res. 2016; 51:1-40. [PMID: 26113209]
  52. Vecino E, Heller JP, Veiga-Crespo P, Martin KR, Fawcett JW. Influence of extracellular matrix components on the expression of integrins and regeneration of adult retinal ganglion cells. PLoS One. 2015; 10:e0125250 [PMID: 26018803]
  53. Templeton JP, Nassr M, Vazquez-Chona F, Freeman-Anderson NE, Orr WE, Williams RW, Geisert EE. Differential response of C57BL/6J mouse and DBA/2J mouse to optic nerve crush. BMC Neurosci. 2009; 10:90 [PMID: 19643015]
  54. Chalhoub N, Baker SJ. PTEN and the PI3-kinase pathway in cancer. Annu Rev Pathol. 2009; 4:127-50. [PMID: 18767981]
  55. Park KK, Liu K, Hu Y, Smith PD, Wang C, Cai B, Xu B, Connolly L, Kramvis I, Sahin M, He Z. Promoting axon regeneration in the adult CNS by modulation of the PTEN/mTOR pathway. Science. 2008; 322:963-6. [PMID: 18988856]
  56. Wang J, Li Y, King R, Struebing FL, Geisert EE. Optic nerve regeneration in the mouse is a complex trait modulated by genetic background. Mol Vis. 2018; 24:174-86. [PMID: 29463955]
  57. Zhang J, Yang D, Huang H, Sun Y, Hu Y. Coordination of Necessary and Permissive Signals by PTEN Inhibition for CNS Axon Regeneration. Front Neurosci. 2018; 12:558 [PMID: 30158848]
  58. Ciobanu DC, Lu L, Mozhui K, Wang X, Jagalur M, Morris JA, Taylor WL, Dietz K, Simon P, Williams RW. Detection, validation, and downstream analysis of allelic variation in gene expression. Genetics. 2010; 184:119-28. [PMID: 19884314]
  59. Taub FE, DeLeo JM, Thompson EB. Sequential comparative hybridizations analyzed by computerized image processing can identify and quantitate regulated RNAs. DNA. 1983; 2:309-27. [PMID: 6198132]
  60. Zhao S, Fung-Leung WP, Bittner A, Ngo K, Liu X. Comparison of RNA-Seq and microarray in transcriptome profiling of activated T cells. PLoS One. 2014; 9:e78644 [PMID: 24454679]