A defined set of genetic instructions encodes functionality in complex organisms. Delineating these unique genetic signatures is essential to understanding the formation and functionality of specialized tissues. Vision, one of the five central senses of perception, is initiated by the retina and has evolved over time to produce rod and cone photoreceptors that vary in a species-specific manner, and in some cases by geographical region resulting in higher order visual acuity in humans. RNA-sequencing and use of existing and de novo transcriptome assemblies allowed ocular transcriptome mapping from a diverse set of rodent and primate species. Global genomic refinements along with systems-based comparative and co-expression analyses of these transcriptome maps identified gene modules that correlated with specific features of rod versus cone retinal cellular composition. Organization of the ocular transcriptome demonstrated herein defines the molecular basis of photoreceptor architecture and functionality, providing a new paradigm for neurogenetic analyses of the mammalian retina in health and disease.
Comparative analysis of genomes can reveal specific pathways that underlie functional (1) and phenotypic variation (2). Transcriptome studies across distinct species have shown a large degree of conservation, even when comparing primates to amniotes (3). Thus, model organisms such as mice (4) and rats (5) have been used to study complex diseases and identify disease-specific genomic loci. Comprehensive genotype to phenotype maps gleaned from such organisms have been essential in identifying evolutionarily conserved pathways affecting human disease (6). A confounding feature of interspecies genetic analysis, however, is that gene expression differences between species can be selectively neutral, and thus have little functional significance (7). Therefore, systems biological approaches must be used in concert to reveal functionally relevant patterns of transcriptional organization and provide novel insights into human genomic functionality. Such approaches have shed light on critical genetic features of neuronal cell types in the human brain (8) and those pathways affected by disease (9).
Neuronal tissue is well suited for this genomic approach as it evolves significantly slower due to a more fine-tuned expression profile compared to other organs (10,11). Conserved genes can highlight the preservation of core functionality whereas divergent genes and gene families can underlie species-specific functions in those tissues (12). Gene expression changes can also reflect alterations in cellular composition and connectivity, which are required for precise functioning in neuronal tissues. An excellent example is the neural retina, composed of diverse cell types that support different ratios of rod and cone photoreceptors that vary in populations with species as well as geographic in primates (13). Recent application of RNA-sequencing (RNA-Seq) technology (14) has begun to uncover complex genetic features of the retina (4,15–17), supporting tissues (18), and specific retinal cell populations (19,20). With dysfunction of the photoreceptor and adjacent supporting cell layers in the retina as the driving force for most retinal degenerative diseases (21) more complete transcriptome analyses of the retina, in both health and disease states is a high priority. A new approach for understanding visual function at the genetic level is essential, because genetic interrogation of human tissue is limited. In addition, the primary rodent models studied, Rattus norvegicus (rat) and Mus musculus (mouse) display a rod-rich retina. Although higher order mammals like humans only have approximately 5% cone photoreceptors in the retina, a central region termed the macula, most specifically the fovea, is rich in cone photoreceptors (22) and is responsible for high-resolution vision. Thus, restricting genetic interrogation to these two rodent models cannot reveal the more complex retinal disease aetiologies afflicting the central cone-rich region of primates (23). Less studied rodents, Arvicanthis niloticus (Nile rat) and Ictidomys tridecemlineatus (ground squirrel), are diurnal cone-dominant rodent species that reflect features of human vision with respect to electroretinogram readings (24), neuronal connections (25), and genomic karyotype (26), but are limited due to lack of comprehensive genomic annotation. Taken together, these four rodent species feature a spectrum of rod versus cone photoreceptor ratios and behavioural patterning that can be used to delineate genetic features that shape the complex human retina.
In this work, RNA-Seq of ocular tissue from Rattus norvegicus, Mus musculus, Arvicanthis niloticus, Ictidomys tridecemlineatus rodent species allowed transcriptome analyses across diverse physiologic retinal architectures that support varying photoreceptor populations. Comprehensive analysis of gene co-expression networks across rodent species revealed enrichment of gene signatures in specific retinal neuroanatomical structures. Ocular tissue from higher order mammals, Macaca fascicularis and Homo sapiens, were also incorporated into the study to elucidate regional specific transcriptome features that support higher order visual function. This approach separated conserved genetic elements essential for ocular function from those divergent features that have shaped higher order primate vision. This work elucidates novel genetic features of the retina and specific photoreceptor subtypes needed to identify pathologic changes affecting visual function and guide future therapies for such retinopathies.
Physiologic consequences of photoreceptor gradation
The vertebrate eye evolved over the course of 600 million years to produce the retinal architecture of modern day mammals (27). Visual function is mediated by a complex interplay between the photoreceptor cells and the neighbouring retinal pigmented epithelium (RPE) to drive phototransduction (28) and the retinoid cycle (29) (Fig. 1) . Reliance on the universal visual chromophore, 11-cis-retinal, by both rod and cone photoreceptors indicates that distinct metabolites in the retina must exist to support the difference in the photoreceptor response kinetics (30). Nocturnal species, such as Rattus norvegicus and Mus musculus, rely on scoptic vision and thus exhibit a retina with an abundance of rod photoreceptors, whereas diurnal species, such as Arvicantis niloticus and Ictidomys tridecemlineatus, depend more on photopic vision and thus feature increasing cone dominance in the retina (Fig. 1). Non-polar retinoids extracted and analyzed from eyes revealed that enrichment of particular retinal metabolites was related to the photoreceptor gradation across rodent species (Fig. 1). In nocturnal species, there was enrichment of retinyl-aldehydes, reflecting the high concentration of visual pigment in rod cells, whereas in diurnal species, there was an increase in retinyl-esters, proposed to be preferentially used in vertebrate cone vision (31). These differing retinal metabolic patterns indicate that distinct genetic features exist to mediate vision in unique rod versus cone-dominant retinas. To an extent, mouse models can be used to delineate the genetic drivers of distinct photoreceptor populations (23). Understanding the physiologic consequences of photoreceptor gradation relates to the retinal architecture in primates, wherein the peripheral retinal architecture is rod dominant and the central, macular region is cone dominant (22). Therefore, we used RNA-Seq to analyze the genetic signatures unique to varying photoreceptor populations in the physiologic retinal environment across these rodent species.
RNA-seq and genomic refinement for global ocular transcriptome mapping
To understand the genetic signatures unique to distinct rod versus cone-dominant photoreceptor populations, we conducted RNA-Seq of rodent species (Rattus norvegicus, Mus musculus, Arvicanthis niloticus, Ictidomys tridecemlineatus) as well as higher order primates (Macaca fascicularis and Homo sapiens). Prepared libraries of poly-adenylated RNA from whole eye and retinal tissue were sequenced with the Illumina platform, yielding ∼1.96 billion reads (32). These were further refined manually with RefSeq gene models (RNor, GRCm38, Mmul, GRCh37) (33), most notably for Macaca fasicularis. Over 50% of such gene models were missing in the rhesus assembly (34), even essential visual genes such as rod opsin (Rho). Missing homologous genes were identified in each assembly and reads were mapped to nucleotide sequences from the RefSeq gene models. The chromosomal architecture of Mus musculus was then used as a scaffold to visualize coverage of the 16,930 genes detected (≥1 RPKM (35)) across rodent and primate species (Fig. 2A). This global approach provided sensitive assessments of expression levels, highlighting the gradation across species of key genes such as short-wavelength cone opsin, Opn1sw, and rod opsin, Rho (Fig. 2B). The landscape of gene expression across rodent and primate eyes revealed that 71% of reads mapped to 12,029 homologous genes (Fig. 2C), emphasizing that conservation of core visual function across species restricted dramatic transcriptome divergence. Divergent genetic signatures were further investigated to understand possible patterning aspects that have permitted more fine-tuned visual function in primates.
Divergent expression profiles implicate adaptive genetic signatures in the retina
A quantitative map of gene expression across the typical model rodent and primate retinas is illustrated in Fig. 3A, highlighting the genetic signatures unique to a particular species. Breakdown of the categorical enrichment profile of the 2,088 divergent genes between rodents and primates demonstrated species-specific genetic trends (Fig. 3B). Interestingly, both rodents and primates had a large percentage of unique transcripts with no annotated function. The most markedly distinct enrichment profile was seen in transcription/translation genes, with primates featuring enrichment relative to rodents, even when compared to all orthologous transcripts. Further analysis of this transcriptional enrichment pattern provided examples of how the differential expression of genes can underlie unique retinal features driving adaptive aspects of mammalian vision. The Rax family of homeobox genes are essential for vertebrate retinal development (36), wherein our analysis revealed that Rax2, a gene involved in the neural retina development driven by BMP signalling (37), had expression restricted to primates. Most categories of genes, however, failed to show any marked differences in enrichment between rodents and primates. Therefore, to identify those genetic elements unique to a particular photoreceptor composition, more homogeneous photoreceptor populations would be required to distinguish factors unique to rod and cone photoreceptor ratios that shape the unique architecture of the primate retina.
De novo gene assembly of non-model rodent species
We reasoned that the non-model rodent species Arvicanthis niloticus and Ictidomys tridecemlineatus would provide unique information regarding physiologic retinal environments abundant in cone photoreceptors. De novo assembly was carried out and used as a scaffold for quantitative RNA-seq transcriptome analysis (38). The assembled Arvicanthis niloticus eye transcriptome consisted of 124,439 transcripts totaling 100,656,031 bases, with an N50 value of 1,595 bp for all transcripts. Assembled Ictidomys tridecemlineatus transcriptomes consisted of 165,655 transcripts totaling 103,326,677 bp for the eye and 214,069 transcripts totaling 135,830,374 bp for the retina, with N50 values of 2,205 bp and 1,977 bp, for the eye and retina, respectively (Fig. 2A, we identified expression of 14,799 genes out of the 16,930 genes detected across species with existing genomic assemblies, representing 87% coverage of the de novo assembled transcripts. These novel ocular transcriptome assemblies allowed us to expand from single gene comparisons to a more systematic approach to identify biologically relevant patterns unique to rod- and cone-dominant retinas.
Gene co-expression network analysis of the rodent eye delineates modules of functionally related genes
Functional organization of the ocular transcriptome was investigated by applying weighted gene co-expression network analysis to four physiologically diverse rodent species. As shown in Fig. 1, Rattus norvegicus, Mus musculus, Arvicanthis niloticus, and Ictidomys tridecemlineatus exhibited distinct structural patterning of the retina from rod-dominant to cone-dominant, observed on the genomic level as well as in terms of photopigment expression (Fig 4A). These modules were then interrogated to reveal functional relationships among genes. Analysis of the genes in Fig. 4A revealed that interspersed between known rod and cone photoreceptor genes, such as structural rod disc genes (
Module residence permitted gene function annotation, as genes in modules have a strong likelihood of participating in specific functions represented by their particular modules (Fig. 4B and C). Modules were also examined for enrichment for known retinal disease-causing genes (RetNet,
Geographic expression profiles reveal genetic factors driving higher order vision in the primate retina
In addition to whole eye and retinal tissue from Macaca fascicularis, three biological replicates of peripheral retina and four biological replicates of central retina tissue were isolated and sequenced for gene expression analysis. Of the 11,060 genes identified in Macaca fascicularis that were expressed at a level of ≥1 RPKM, 807 were differentially expressed by at least 2 fold (P value ≤0.05) between peripheral and central retina samples (Fig. 5A). Similar to the rodent co-expression analysis, known rod and cone photoreceptor genes were found enriched in the peripheral and central retina (Fig. 5B). In the peripheral retina, where metabolically active rod photoreceptors are abundant, there was a predictable enrichment of genes involved in cellular metabolism. There was also an enrichment of binding/cell adhesion genes in the peripheral retina, as rods have a strong structural role in maintaining retinal integrity. In the central retina, where there is a specialized cone-rich architecture for high-resolution vision, there was an enrichment of genes related to signal transduction, transcription and transport functional categories. Interrogation of these functional categories revealed overrepresentation of key gene families such as cytoskeletal protocadherin proteins, developmental BMP/WNT signalling proteins, glutathione peroxidase metabolic proteins, orphan G protein-coupled receptors involved in signal transduction, and homeobox and Kruppel-like transcription factors. The differential expression of these gene families highlighted the functional relevance of particular pathways in the maintenance of the specialized primate retina.
Visual function has been tailored to fit the evolving needs of various species. Rodents such as Rattus norvegicus and Mus musculus have developed rod-dominant retinas for improved scotopic vision for their nocturnal lifestyle. In contrast, diurnal patterning has led Arvicanthis niloticus and Ictidomys tridecemlineatus rodent species to develop cone-dominant retinas that permit improved photopic vision. Higher order primates have further developed a central cone-rich macula that enables high-resolution vision. At the genomic level, Mus musculus and Homo sapiens differ in ∼2.5% of the bases (39) whereas Macaca fascicularis and Homo sapiens differ in only ∼1.2% of bases aligned between species (40). RNA expression level differences seemingly drive the observed phenotypic diversity (41). Thus analytic approaches must expand from single gene comparisons to systems level approaches to reveal the functional relevance of novel genes as drivers of this phenotypic patterning (42). Our work provides a comprehensive analysis of visual adaptation in rod and cone-enriched architectures across species and even specific geographical portions of the primate retina. More importantly, comparative network analyses identified gene modules that correlate with specific neuroanatomical retinal structures. Furthermore, novel genes pertaining to specific photoreceptor cell types in the retina and their relationships to human visual function and associated pathology were elucidated.
Starting from a phenotypic characterization of rodent species with varying rod versus cone photoreceptor ratios and behavioural patterning (Fig. 1), it became evident that key genetic drivers must exist to account for phenotypic diversity in such closely related species. Native and genetically modified mouse models have been employed to approximate various human retinal disease conditions (17), limiting its relevance. This problem prompted our examination of diverse photoreceptor (rod and cone) architectures to delineate physiologically relevant signatures. The power of this current study lies in the examination of the ocular transcriptome in aggregate across species. Furthermore, in studying whole ocular tissue rather than isolated cell populations we could reduce confounding variables that result from loss in transcript integrity from disruptions of cellular connectivity that occur in preparation of single cell sequencing techniques (43).
Initial analysis of ocular transcriptome reads from species with existing genomes revealed global gene expression profiles, highlighting those essential conserved genes and those with unique signatures that have diverged (Fig. 2). Most importantly, it reinforced the need to use non-model rodent species to better understand physiological consequences of cone enrichment. Examination of genetic signatures between species (Fig. 3) revealed gene families such as GPCR kinases, represented by Grk1 and Grk7, with the former expressed across all species, whereas the latter was restricted to primates and cone dominant rodents, such as Ictidomys tridecemlineatus (44). Restricting genetic analysis to traditional Rattus norvegicus and Mus musculus species would not allow one to understand the specialized cone rich macula of the human retina, prompting the inclusion of non-model cone-dominant rodent species in the genetic analysis. In doing so, this work led to the sequencing and first de novo assemblies of the ocular transcriptomes from Arvicanthis niloticus and Ictidomys tridecemlineatus species. Incorporating these novel genetic assemblies with gene co-expression analyses was essential to reveal their functional significance in the unique cellular composition of the primate retina.
In large datasets as involved in this work, examining distinct modules of gene expression changes across species (45) that are enriched in certain functional categories can highlight phenotypically relevant networks. Co-expression analyses identified genes that corresponded to retinal environments enriched for specific photoreceptor populations (Fig. 4), providing a new methodology for annotating gene function in the retina. Moreover, this enabled the generation of hypotheses related to gene function by examining the module membership. Two such examples are highlighted here. One relates to the cone photoreceptor enrichment of Zic3, a zinc finger transcription factor involved in a variety of developmental processes thought to function in retinal progenitor cell differentiation and proliferation. Zic3 was more highly expressed in cone-dominant rodents and specifically in neuronal tissues across organs (Fig. 4), consistent with previous work (46). The second was Fam169a, a gene of unknown function that was enriched in rod-dominant rodents, but across organs, its expression was restricted to the eye, specifically the retina (Fig. 4). The power of modular analysis is that it can identify genes with a strong likelihood of participating in specific functions represented by each module. In the case of Fam169a, its module is enriched in cell adhesion, metabolism, and transcription related genes as well as in known rod photoreceptor retinopathy genes. This approach can similarly be expanded to other unknown genes to elucidate how they relate to health and ultimately influence disease in the retina. Localization techniques can further elucidate the specific neuroanatomical role of these novel genes.
These results were leveraged against data obtained from Macaca fascicularis in which geographical transcriptome dynamics were examined by comparing reads from central versus peripheral retina samples representing cone and rod enriched areas, respectively (Fig. 5). Comparison of gene expression revealed a predominant up-regulation of gene expression in the central region compared to the periphery. This speaks to the increased complexity of the central macular region as it has been shown that increased expression levels can correspond to higher order modelling when proceeding from non-human primate brains to human brains (47). Key gene families related to metabolism, signal transduction and transcription were enriched in a geographic manner in the primate retina consistent with their function in maintaining the specialized architecture in the human retina (Fig. 6). The neural retina consists of post-mitotic photoreceptor cells that rely on the neighbouring RPE cells for a circadian coordinated mechanism of maintenance (48). These highly metabolic cells are subjected to considerable oxidative stress generated by exposure to light and high oxygen tension. Our work reinforces the idea that key metabolic gene products are enriched in particular retinal environments to mediate this homeostatic process. Glutathione-mediated oxidative stress protection by the GPX family of proteins (Gpx1, Gpx3, and Gpx8) seems essential to the peripheral retina, displaying rod-photoreceptor enrichment, and has been implicated to have disrupted expression in age-related retinal phenotypes (4). The differential expression and localization of different metabolic gene products highlights the complex retinal architecture required for normal function and how mutations in gene products can lead to specific retinal phenotypes such as Stargardt disease. Genes linked to this disease, Abca4 and Elovl4, were both found in this study to be enriched in photoreceptor clusters, especially the rod photoreceptor layer and the peripheral retina in Macaca fascicularis. Increased understanding this disease at a molecular level has allowed systems pharmacological approaches to devise new treatment paradigms (49).
Another avenue of discovery of this work has been elucidating the complex factors that dictate the maintenance of post-mitotic photoreceptor cells. Whereas progress has been made in understanding the transcriptional network involved in rod photoreceptor specification and maintenance, much less is known about how cone photoreceptors are committed and maintained. Understanding retina cell specification, especially those factors that drive cone cell development and survival, is crucial in developing future regenerative therapies for the retina (50). Recent work has shown that a multifunctional BMP, TGF-β and WNT network of genes can drive cone photoreceptor development from human embryonic stem cells (51). Supporting this point, the central primate retina has a reduced expression of Bmp4 compared to the peripheral retina with an increase of inhibitory factors in this multifunctional network such as Wif1 in cone clusters and enhancement of activating factors like Wnt4 and Lrp6 (52) in rod clusters. Even in signal transduction gene products, we see that orphan GPCRs are integrated into this hierarchical network, as Gprc5b mediates neurogenesis via a beta-catenin pathway (53) and is enriched in cone populations. There is an observed cone enrichment pattern of the BMP interacting FGF network (54) (Fgf1, Fgf2, Fgf9, Fgf17) and even a novel downstream gene of Fgf2 such as Ccdc115 (55). Moreover, in the central retina, we noted enrichment of Kruppel-like factors, such as Klf4 that have been implicated in stem cell renewal (56) and are directly influenced by the BMP network (57). Such findings suggest that the central cone-rich macula of primates can have stem cell functionality that allows their life-long survival and maintenance. Taken together these analyses highlight a set of key transcriptional frameworks that are instrumental in the function and maintenance of the adult retina. With technologies emerging that can generate three-dimensional retinal tissue platforms (58), this could shape the future of stem cell based therapeutics as these tissue models are used to restore vision (59).
In summary, this work represents a comprehensive analysis of ocular transcriptome mapping to unveil adaptive elements underlying rod versus cone photoreceptor specialization in mammals. Moreover, the use of gene co-expression relationships in rod and cone photoreceptor-dominant retinas revealed known genes along with novel genetic elements that are crucial for future study, especially in determining changes caused by retinal diseases (60). Interestingly, we found enrichment of known disease causing retinal genes residing in rod clusters in the co-expression analysis and the primate peripheral retina that likely stem from the reliance on ocular genomes of rod dominant species. The elucidation of transcriptome maps from ocular tissue of cone-dominant mammals in this work could facilitate the identification of novel disease genes related to cone retinopathies. Such associations can enable identification of clinically relevant pathways suitable for therapeutic intervention that have previously been unavailable. With the CRISPR-CAS system used to study the effects of human disease through genetic changes in mice (61), enabling the study of multiple gene mutations (62) and extending to gene editing in primates (63), improved understanding of human disease should result as candidate genes are increasingly identified. This work offers a novel outlook on the structural and functional organization of the mammalian retina, providing a paradigm shift as to how to decipher the genotype to phenotype relationship to better predict outcomes and design therapeutic interventions through a systems-based approach (49).
Materials and Methods
C57BL/6J mice (Mus musculus) were purchased from the Jackson Laboratory (Bar Harbor, ME, USA). Nrl-deficient mice (64) with a B6 background were kindly provided by Dr. Anand Swaroop (National Eye Institute, Bethesda, MD, USA). Rhodopsin knockout mice (65) with a B6 background were obtained from Dr. Janis Lem (Tufts University, Boston, MA, USA). All mice used were 4 weeks of age. Four-week-old Long-Evans rats (Rattus norvegicus) were purchased from Harlan Laboratories (Madison, WI, USA). Nile rats (Arvicanthis niloticus) at 6 weeks of age, were kindly provided by Dr. Laura Smale (Michigan State University, East Lansing, MI, USA). Wild caught ground squirrels (Ictidomys tridecemlineatus) were purchased from TLS Research (Bloomingdale, IL, USA) and also kindly provided by Dr. Wei Li (National Institutes of Health, Bethesda, MD, USA). Mice, rats, Nile rats, and ground squirrels were all housed in the Case Western Reserve University (CWRU) animal facility where they were maintained on standard chow diets in a 12h light (∼10 lux)/12h dark cycle. Enucleated macaque (Macaca fascicularis) eyes and retinal tissue were obtained from Ricerca Biosciences (Painesville, OH, USA). All animal procedures and experiments were performed in accordance with U.S. animal protection laws and were approved by the CWRU Animal Care Committees and conformed to both the recommendations of the American Veterinary Medical Association Panel on Euthanasia and the Association for Research in Vision and Ophthalmology.
Clinical evaluation of the human patient from whom the retinal tissue was obtained was carried out at the Cleveland Clinic Cole Eye Institute (Cleveland, OH, USA). This research conformed to the tenets of the Declaration of Helsinki.
HPLC analyses of retinoids from the eyes of rat, mouse, nile rat, and ground squirrel species
After dark-adaptation for 3–12 h, animals were euthanized and their eyes were enucleated and flash frozen in liquid nitrogen. To extract retinoids, each eye was transferred to a glass homogenizer in 1 ml of 40 mM NH2OH in a 1:1 solution of PBS:methanol. Each eye was thoroughly homogenized and transferred to a 15 ml Falcon tube, and then kept at room temperature for 20 mins to allow derivatization of retinyl-aldehydes into their retinyl oxime metabolites. Next, 4 ml of hexanes were added to each tube, vigorously shaken to facilitate extraction of retinoids, and then centrifuged at 5000g for 10 min, after which the upper organic phase was collected in a glass test tube. The organic solvent was dried in a rotary SpeedVac, and the residual retinoids were re-dissolved in 250 ml of hexanes, vortexed, and then transferred to standard 2 ml vials with 0.4 ml inserts for HPLC analysis. These analyses were performed with a normal phase column (Agilent Zorbax Sil 5 µm 4.6 x 250 mm) (66). For retinoid separation, a linear gradient of 0.5% ethyl acetate in hexane was applied over 15 min followed by 35 min of 10% ethyl acetate in hexane at a constant flow rate of 1.4 ml/min. Retinyl-esters and retinyl-oxime metabolites were detected at wavelengths of 325 nm and 360 nm, respectively. The column was regenerated by a 10 min equilibration with 0.5% ethyl acetate in hexane after each run.
Rodent species (Rattus norvegicus, Mus musculus, Arvicanthis niloticus, and Ictidomys tridecemlineatus) were euthanized by cervical dislocation and eyes and retinal samples were dissected out and immediately placed in RNALater stabilization reagent (Qiagen, Valencia, CA, USA). For geographic retinal tissue collection from macaques, each retina was cut into quadrants, and the central region containing the macula was carefully dissected out with scissors to separate it from peripheral retinal tissue before being placed in RNALater. The human retinal sample consisted of a tumour-free hemi-retina that was immediately placed in RNALater from a donor patient requiring enucleation for a large ocular melanoma.
RNA-Seq library preparation. Library preparation for Illumina RNA-Seq (Illumina Inc., San Diego, CA, USA) was carried out as previously described (17). The library preparation for isolated eye and retinal tissue in RNALater, obtained from human and macaque species, was performed similarly to that used for rodent tissues except for the larger scale homogenization required. Each rodent species-specific library preparation was pooled from multiple RNA samples: five mice for Mus musculus, five rats for Rattus norvegicus, five Nile rats for Arvicanthis niloticus, and four ground squirrels for Ictidomys tridecemlineatus. Human and macaque samples were each prepared from one eye. Biological replicates were prepared for both Mus musculus and Macaca fascicularis eye and retina samples.
RNA-seq, read mapping and determination of reads per kilobase per million reads (RPKM)
Libraries of Rattus norvegicus, Mus musculus, and Homo sapiens were sequenced by single-end sequencing, whereas prepared libraries of Arvicanthis niloticus, Ictidomys tridecemlineatus, and Macaca fascicularis were sequenced by paired-end sequencing with the Illumina Genome Analyzer IIx or HiScan SQ. Data for Rattus norvegicus, Mus musculus, Macaca fascicularis and Homo sapiens were processed and aligned with the University of California–Santa Cruz (UCSC; Santa Cruz, CA, USA) rat (rn4), mouse (mm9), rhesus (rheMac3), and human (hg19) genome assemblies and transcript annotations, respectively, using the Genomic Short-read Nucleotide Alignment Program (GSNAP), manual extraction of uniquely-mapped reads, HTseq for raw read counts of genes, and manual calculation of RPKM (35) statistics by gene. Details of sequencing and read statistics are provided in
Transcriptomes of Arvicanthis niloticus and Ictidomys tridecemlineatus were assembled de novo. The processed reads were assembled using velvet (67), oases, and Trinity (68) with multiple k-mers ranging from 21 to 67 with an increment of 2 whereas reads assembled with ABYSS and trans-ABYSS (69–71) were done with multiple k-mers ranging from 40 to 78 with an increment of 2. The assemblies were combined and the duplicated transcripts were collapsed together to remove redundancy. Transcriptome annotation was done by first searching against the SWISS-PROT database using the BLASTx algorithm (72) and an E value cutoff of 1e-5. Transcripts with the best matches to proteins from viruses and bacteria were considered as potential contaminants and thus were removed from further analysis. The resulting filtered BLASTx results were then used for Gene Ontology (73) analysis with the blast2GO program (74). Protein domains were predicted for all such transcripts by searching the PFAM database (75) with the HMMER program (76). The expression value for each transcript was calculated by aligning the original Illumina reads to the transcriptome using the Bowtie program (77), with all mapping reads reported. The expression value for each transcript was calculated as RPKM (35) using a customized Perl script. Details of these sequencing and read statistics are described in
The processed and raw fastq files for Mus musculus (accession numbers GSE38359, GSE29752, GSE48974, GSE84927), Rattus norvegicus (GSE84932), Arvicanthis niloticus (GSE84928), Ictidomys tridecemlineatus (GSE84931), Macaca fascicularis (GSE84929), and Homo sapiens (GSE84930) have been deposited in the NCBI GEO database and are available for public use.
Analysis of RNA-seq data
Concordance of R2 ≥ 0.95 has been demonstrated between biological replicates (4) and was noted in this study from biological replicates isolated from pooled or single tissue samples. Orthologous gene relationships across species were extracted with the NCBI biosystems database (78). Genes were then categorized for their biological and functional roles with Gene Ontology using AmiGO software (79). From the total list, an initial set of 4049 genes that displayed differential expression across the rodent species of Rattus norvegicus, Mus musculus, Arvicanthis niloticus, and Ictidomys tridecemlineatus, were considered for gene co-expression network analysis. Co-expression networks were constructed from RNA-Seq data by calculating Pearson correlations across these rodent species. Concordance of gene expression with the Pearson correlation was such that the gene expression pattern was weighted against the pigment composition in rodent species. A threshold (R2 ≥0.5) was set to identify those with moderate to strong correlations. Genes grouped with coherent expression profiles that exhibited topographical overlap were assembled into discreet modules by average linkage hierarchical clustering (80,81), with the gene modules corresponding to branches of a dendrogram. For each module, a heat map of transcript expression was produced with rows corresponding to genes and columns corresponding to the rodent eye sample. These modules were then examined to identify both the global and individual networks to which they corresponded (82).
Histological and IHC procedures were performed as previously described (83). Cross sections of mouse eyecups were incubated with primary antibodies against Rho and Opn1mw. Signals were detected with Cy3- or Alexa488-conjugated secondary antibodies. Nuclear staining was achieved with 4’,6-diamidino-2-phenylindole (DAPI). A Zeiss LSM 510 inverted Laser Scanning Confocal Microscope was used to image the sections and the analysis was accomplished with ImageJ software.
Semi-quantitative real-time PCR
Total RNA was purified by using the RNeasy Mini Kit with on column DNase treatment (Qiagen, Valencia, CA). Selected genes (
Experimental results were analyzed by an independent two sample t test. P values of 0.05 or less were considered statistically significant. Data presented graphically in figures are shown as means ± standard deviations.
We thank Dr. Leslie T. Webster Jr. for helpful comments on this manuscript. We also thank Dr. Akiko Maeda (CWRU) for primate tissue isolation, the laboratories of Dr. Laura Smale (MSU) for Nile rat animals to establish a breeding colony at CWRU and Dr. Wei Li (NIH) for providing ground squirrel animal tissue for retinoid studies. Ms. Simone Edelheit in the CWRU Genomics Core performed the sequencing.
Conflict of interest statement. None declared.
This work was supported by funding from the National Institutes of Health EY021126 and U01 EY025451, the Arnold and Mabel Beckman Foundation, and the Canadian Institute for Advanced Research (CIFAR). K.P. is the John H. Hord Professor of Pharmacology and Distinguished University Professor.