NRL-Regulated Transcriptome Dynamics of Developing Rod Photoreceptors

Abstract

Gene regulatory networks (GRNs) guiding differentiation of cell types and cell assemblies in the nervous system are poorly understood because of inherent complexities and interdependence of signaling pathways. Here, we report transcriptome dynamics of differentiating rod photoreceptors in the mammalian retina. Given that the transcription factor NRL determines rod cell fate, we performed expression profiling of developing NRL-positive (rods) and NRL-negative (S-cone-like) mouse photoreceptors. We identified a large-scale, sharp transition in the transcriptome landscape between postnatal days 6 and 10 concordant with rod morphogenesis. Rod-specific temporal DNA methylation corroborated gene expression patterns. De novo assembly and alternative splicing analyses revealed previously unannotated rod-enriched transcripts and the role of NRL in transcript maturation. Furthermore, we defined the relationship of NRL with other transcriptional regulators and downstream cognate effectors. Our studies provide the framework for comprehensive system-level analysis of the GRN underlying the development of a single sensory neuron, the rod photoreceptor.

Keywords: Maf; RNA-seq; basic motif leucine zipper; gene regulation; neuronal development; next generation sequencing; organogenesis; photoreceptor differentiation; retina; transcription.

Figure 1. Study Design and Data Generation

(A) Time course of rod photoreceptor differentiation in mouse retina. Morphologies of developing and mature rods are shown at indicated developmental stages, with major events summarized below. ONBL, outer neuroblastic layer; ONL, outer nuclear layer; OS/IS, outer segment/inner segment. (B) Flowchart of integrated transcriptome and network analysis. WT, wild type; RIN, RNA integration number; IP, immunoprecipitation. (C, D) Comparison of rod transcriptomes generated by directional RNA-seq and exon microarray. (C) Principal component analysis (PCA) of directional RNA-seq (left panel) and microarray (right panel) data. The percentages indicate the amount of variation attributed to each principal component. Small shapes indicate individual samples, and larger shapes show the centroid of each grouping. (D) Correlation between RNA-seq and microarray data. RNA-seq and microarray data were de-duplicated and merged based on gene symbols. The highest FC representative from the RNA-seq transcriptome data was used as the gene representative for RNA-seq. Left Panel: Expression values (in log₂ scale) of individual genes of mature rods (i.e., P28 in RNA-seq and P21 in microarray) were compared with each other using a scatter plot. Log₂ FPKM of one and log₂ signal intensity of seven was used as a cut-off of actively transcribed genes in RNA-seq and microarray, respectively. The lowess regression line is shown in red, and the Pearson’s correlation coefficient is indicated. Right Panel: Differentially expressed genes between newborn (P2) and mature rods (P21 for microarray and P28 for RNA-seq) were identified as having >2 fold change and FDR<0.05 in microarray data or <0.01 in RNA-seq data. The scatter plot shows genes that are significantly differentially expressed genes in both analyses, microarray only, and RNA-seq only (green, purple, and orange dots, respectively). Black dots indicate genes with no significant change in expression in either analysis. FPKM, fragments per kilobase of exon model per millions of reads; FC, fold change; FDR, Benjamini-Hochberg false discovery rate.

Figure 2. Transcriptome Dynamics during Rod Differentiation

(A) Hierarchically-clustered heat map of transcripts (a total of 11,146 transcripts that were expressed ≥5 FPKM at any time point in WT) during rod photoreceptor development (left panel) and correlation coefficient between every possible pair of transcriptomes at various differentiation stages calculated from directional RNA-seq data and presented in a matrix (right panel). A sharp change in transcriptome landscape is apparent between P6 and P10. (B) Dynamic expression patterns of genes specific for photoreceptors (14 rod genes, 11 cone genes and 13 genes expressed in both) and other retinal cell types during rod differentiation. In heat maps, the average expression at each time point is plotted in log₂ scale, and only those transcripts that were expressed at ≥5 in all replicates of at least one time point were included. Color scale is indicated in the bottom.

Figure 3. Transcriptome regulation by NRL

(A) PCA analysis of time series RNA-seq data from Nrl-GFP+ wild type rod photoreceptors and S-cone like photoreceptors from the Nrl-KO retina. A significant change in transcriptome landscape was detected in GFP+ Nrl-KO photoreceptors. (B) Heat map of significantly, differentially expressed protein-coding transcripts. Color scale is indicated. (C) Heat map of differentially expressed genes containing NRL ChIP-seq peaks. An integrated analysis was performed combining differential gene expression data with the NRL ChIP-seq (Hao et al., 2012). The DE transcripts having at least one NRL ChIP-seq peak were considered putative, direct transcriptional targets and were included in the heat map. (D) Identification of transcription regulatory proteins that are putative NRL direct targets. Shown on the top is a genomic locus of one of the differentially expressed genes, Tnfaip3, and NRL ChIP-seq coverage plot.

Figure 4. DNA Methylation of Photoreceptor Genes

(A) Correlation between DNA methylation level and gene expression of rod- and cone-specific genes. Degree of DNA methylation (%) in promoter or gene body region of individual rod- and cone-specific genes (green triangle and blue circle, respectively) are plotted against gene expression level (in log₁₀FPKM). Please note that some promoters and gene bodies were not covered in RRBS analysis and thus not included in the plot. (B) Plotted is correlation between DNA methylation level (%) in promoter or in gene body region and gene expression (in log₁₀FPKM) of non-photoreceptor retinal genes (red triangle) and housekeeping genes (blue circle).

Figure 5. In-silico De Novo Analysis and Validation of Putative Protein Coding Transcripts

(A) Heat map of standardized log₂ FPKM values (with one offset) of 222 previously un-annotated protein-coding transcripts (NPCs) across all time points in rods from wild type and S-cone-like photoreceptors from Nrl-KO retina. Individual replicates are shown to display the degree of variability intrinsic to this type of analysis. The three largest clusters in the dendrogram are highlighted in blue (1), magenta (2), and green (3), respectively. (B) Functional stratification of conserved protein domain identified by HMMER. Identified domains were grouped according to their described function in Pfam database. The relative frequency of occurrence was calculated for domains observed at least 10 times among the significant protein domain hits. (C) qRT-PCR validation of three select NPC transcripts. Relative enrichment of each transcript compared with expression of housekeeping gene Actb was plotted for P2 and P28 wild type and Nrl-KO retinas. Double asterisks indicate statistically significant changes. Agarose gel images are shown below. (D) Schematic representation of a NPC, TCONS_00123129. NRL and CRX ChIP-seq peaks are located in close proximity to the TSS of this transcript and placental mammalian evolutionary conservation peaks (as available from USCS genome browser) overlap with the predicted gene structure, highlighted by a dashed box. Histograms for RNA-seq and NRL- and CRX-ChlP-seq indicate the coverage of all aligned reads across the genomic area. Below the RNA-seq histogram are shown the individual sequence reads. Red or blue blocks indicate discrete reads in antisense or sense orientation, respectively, with thin blue lines marking the splice portions of the reads. (E) In-situ hybridization of the previously un-annotated transcripts TCONS_00044530 and TCONS_00069915 at P6 and P28 in wild type retina. The signal for both transcripts follows the dynamic expression pattern seen in panel A and Figure S3. Blue, DAPI nuclear stain; Red, in situ signal in merged images; White, in situ signal in single channel images. Scale bar, 50 µm. GCL, ganglion cell layer; INL, inner nuclear layer; NBL, neuroblastic layer; ONL, outer nuclear layer.

Figure 6. Alternative Splicing and Promoter Usage Events Identified during Rod Differentiation

(A) Types of alternative splicing events detected and their frequencies in time-wise and group-wise comparisons. Splicing events are color coded to match the relevant graphs. (B) Differential splicing events during maturation of rods and Nrl-KO S-cone like photoreceptors. Shown is the genomic area of Clta gene (upper lane in the left panel), which comprises seven exons and six introns, and RNA-seq coverage plots of P2 and P28 wild type rods in the area of exons four through seven (bottom lanes in the left panel). Solid boxes and lines represent exons and introns, respectively, and untranslated regions are illustrated as thinner boxes. The fifth and sixth exons of Clta gene are differentially spliced and highlighted in blue. Individual spliced reads demonstrate the splice isoform that are generated from Clta at both time points. Translated exons, untranslated exonic regions and introns are indicated as thick boxes, thin boxes and lines, respectively. RNA-seq coverage plots are shown as gray histograms. For Clta gene, individual sequence reads (red and blue boxes for antisense and sense orientation) with spliced area (thin blue lines) are also indicated. (C) Alternative promoter usage of Hcls1. RNA-seq reads show only partial coverage of known exons of Hclsl (exons 18–14), with another exon introduced upstream to exon 8 (red asterisk). Dashed boxes indicate NRL and CRX ChlP-seq peaks. In situ hybridization for a retinal isoform of Hcls1 (probe designed to hybridize to approximately 1000 bases starting in the previously un-annotated exon) shows significant increase in expression from P6 to P28. DAPI nuclear stain is indicated in blue, and in situ signal is indicated in red and in white in merged or in single channel images, respectively. Scale bar, 50 µm. GCL, ganglion cell layer; INL, inner nuclear layer; NBL, neuroblastic layer; ONL, outer nuclear layer.

Figure 7. Gene Regulatory Network Analysis

(A) Regulatory hubs of NRL-centered GRN. (B) Temporal expression profile of secondary hubs. Expression values of four genes in rods (blue line) and Nrl-KO photoreceptors (red dashed line) are plotted in log₂ scale for each time point. When multiple transcripts are expressed, expression level of the most highly expressed transcript is indicated. (C) NRL-centered network and highly represented GO terms among target genes of NRL and secondary hubs. Circles indicate individual nodes (i.e., regulators and target genes) with primary regulatory hub NRL and four other secondary hubs labeled, and lines indicate predicted interaction between regulatory hubs and targets. Some targets are regulated by one regulatory gene, while others have regulatory interaction by two or more hubs. Target genes that belong to highly represented GO terms grouped into C1 through C5 were highlighted in indicated color.