Identification of novel circadian transcripts in the zebrafish retina

Circadian rhythms are physiological and behavioural changes that recur in about 24 h that allow the organism to adapt to fluctuations in the external environment driven by the day and night cycles imposed by the sun. Empirical observations in natural conditions and elegant experiments in laboratory models established the three tenets of circadian biology: (1) the core clock is internal; (2) the core clock can be entrained to external cues, especially light; and (3) the basic clock mechanism is conserved. The internal nature of the clock became evident when it was seen that the circadian expression of certain genes and the manifestation of circadian behaviour such as sleep and eclosion would continue to show approximately 24 h rhythms irrespective of the external cues. Thus, under constant dark conditions imposed externally, the fruitfly Drosophila melanogaster will continue to display nearly 24 h rhythms, gradually and inexorably slipping from the objective time. However, in the presence of external light cues, the rhythm adjusts to match the cue (Partch et al., 2014).

Transplantation experiments in model organisms (Ralph et al., 1990) and the prevalence of free-run circadian rhythms in certain types of blindness (Sack et al., 1992) have now firmly established that the core clock resides in the suprachiasmatic nucleus (SCN) in the mammalian brain. Many organs and even individual cells have peripheral clocks that are synched to the master clock but when isolated, continue their individual, but eventually dampening, rhythms (Tosini and Menaker, 1996; Moore and Whitmore, 2014). In organisms as distant as the fungus Neurospora crassa, insects like the fruitfly, nocturnal rodents such as mice, and humans, the primary clock mechanism depends on feedback regulation of a set of highly conserved genes (Loudon et al., 2000). The various manifestations of circadian rhythms are brought about by the expression of output genes, under the regulation of the core regulatory transcription factors. The core clock genes and its associated transcriptomes are highly organ specific (Zhang et al., 2014).

The retina is a thin layer of complex tissue along the inner wall of the eye, made up of intricately connected layers of cells with more than 60 functionally distinct neuronal types (Masland, 2012). The rods and cones, with their distinctive shapes, collect the input light and communicate these signals to bipolar cells in the inner nuclear layer (INL). From the INL, the signal is transmitted to the retinal ganglion cells (RGCs) in the ganglion cell layer (GCL) which eventually relay it to the brain via the optic nerve. Furthermore, signals are modulated by amacrine and horizontal cells when they flow from photoreceptors to bipolar cells and from bipolar cells to RGCs, respectively. Since the retina is an extension of the central nervous system, many phenomena relevant to CNS disorders have been studied in the more accessible retinal cells (Athanasiou et al., 2013; Lamb et al., 2007). Retinal cell types are prone to oxidative damage and endoplasmic reticulum (ER) stress arising from the rapid turnover of proteins in addition to light-induced damage to lipids, protein and DNA. This necessitates extensive repair processes and renewal in coordination with circadian fluctuations. Although these repair processes are shared across evolution, the zebrafish (Danio rerio) retina is unique in its ability to regenerate all cell types upon extensive damage, even into adulthood, which is attributed to retinal Muller glial cells (Jorstad et al., 2017).

The Muller glial cell-mediated injury response is marked by the strong upregulation of core clock genes in zebrafish but not in mouse (Sifuentes et al., 2016). In addition to this, many vital functions of the retina are under the influence of diurnal or circadian factors: melatonin production (Huang et al., 2013), photoreceptor disc shedding (Grace et al., 1999), extracellular pH of the eye (Dmitriev and Mangel, 2000) and rod–cone coupling (Ribelayga et al., 2008) to name a few. The segregation of genuine circadian genes from those that merely respond to light and hormonal cues will significantly enhance our knowledge of circadian biology and the role of the retina in providing inputs to the core clock and maintaining a robust peripheral clock. Furthermore, an extensive catalogue of temporally resolved retinal gene expression in zebrafish would be a valuable resource for exploring the gene regulatory networks behind the complex multicellular organization of retinal function.

We compared the expression profiles of core clock genes in zebrafish, mouse and baboon to establish that zebrafish recapitulates the salient features of the diurnal behaviour of the baboon. Since temporal gene expression patterns of the zebrafish retinal tissue were not available, we generated RNA-seq data and carried out extensive annotation of the retinal transcriptome of the zebrafish, during the light and dark phases of a day. Then, we validated the expression of selected novel transcripts and studied their response to light and dopamine. Lastly, we used free run experiments under constant dark conditions, to unambiguously assign a circadian role to a subset of the novel transcripts. This study not only provides a valuable resource for understanding the zebrafish retina but also uncovers novel genes that encode proteins and non-coding RNAs whose functional role in circadian gene expression is currently unknown.

Assam wild-type (ASWT) zebrafish (Danio rerio) strain was used for all the experiments. Adult fish used were six months to one year old. Maintenance and experiments were strictly within the guidelines provided by Institutional Animal Ethics Committee (IAEC) of the CSIR Institute of Genomics and Integrative Biology, India with timed feeding and uniform temperature of 28±2°C. Standard practice for use of zebrafish was followed, as recommended by Westerfield (2000). Zebrafish were anaesthetized using Tricaine (Sigma) for whole eye cup dissection. Dim red light was employed for eye cup dissection at dark time points.

From zebrafish whole eye, RNA was isolated using RNeasy kit (Qiagen). Approximately 1 μg total RNA was used for library preparation using Truseq stranded RNA library (Illumina). Ribosomal RNA (rRNA) in each sample was depleted as per manufacturer's instructions (Illumina) followed by fragmentation. Complementary DNAs (cDNAs) were prepared after rRNA depletion. Then, sequencing adapters were ligated and the enriched libraries were subjected to 101×2 paired end sequencing using Hiseq 2500 platform following the regular protocol (Illumina). The FASTQ sequencing reads were adapter-trimmed along with a minimum length cut-off of 50 bases using Prinseq-lite. The reads were aligned to zebrafish genome assembly (Zv9) using TopHat (v.2.0.11) followed by reference-based assembly using Cufflinks (v.2.2.1). Then differentially expressing transcripts were identified using Cuffdiff (v.2.2.1).

Tissue or whole animal RNA was isolated using RNeasy kit (Qiagen). cDNA was prepared using a random nanomer-based reverse transcriptase core kit (Eurogentec). qRT-PCR was carried out using SYBR Green master mix (Applied Biosystems) for detection in Light cycler LC 480 (Roche). All primers used for qRT-PCR are given in Table S2.

Depletion of tyrosine hydroxylase positive (Th⁺) amacrine cells in one eye was achieved by 2 µl intraocular injection of a 1:1 mixture of 6-hydroxydopamine (6-OHDA) and pargyline (Sigma) for two consecutive days at zeitgeber time (ZT) 4. Corneal incision and intraocular injection were followed as described in Li and Dowling (2000). One week after injection, ablation of Th⁺ amacrine cells from the retina was analysed by immunolabelling and qRT-PCR against tyrosine hydroxylase protein and mRNA, respectively.

Whole retinas were dissected out of the eye cup. Specimens were fixed in 4% paraformaldehyde and 5% sucrose overnight followed by 0.1 mol l⁻¹ glycine treatment for 10 min at room temperature. For antigen retrieval, tissue was placed in sodium citrate buffer (10 mmol l⁻¹ sodium citrate, 0.05% Tween-20; pH 6) and incubated in boiling water for 10 min then cooled to room temperature followed by washing with 0.5% Triton X-100 in PBS. Samples were blocked for 2 h (5% BSA,1% DMSO, 1% Tween-20, 0.1% sodium azide in PBS) and incubated overnight at 4°C in primary antibody (1:500, tyrosine hydroxylase, Sigma-Aldrich, SAB2701683) in blocking buffer with 2% BSA. Secondary antibody (1:500, FITC conjugated, Invitrogen) was then added for 2 h at room temperature. Tissue was mounted and imaged using a Leica TCS SP8 confocal microscope.

Whole-mount in situ hybridization (WISH) was performed by following the standard zebrafish protocol (Thisse and Thisse, 2008). Adult eye cryosections were subjected to an extra step of acetylation treatment and slides with sections were treated with 0.25% acidic anhydride in 0.1 mol l⁻¹ ethanolamine (pH 8) to prevent non-specific binding of probe on slides.

The publicly available GEO dataset GSE106121 was used in our analysis. Custom code provided online (https://bitbucket.org/Bastiaanspanjaard/linnaeus) was used in single-cell transcriptome analysis using Seurat (Macosko et al., 2015; Butler et al., 2018).

Circadian oscillations in gene expression have been documented in different tissues and sub-regions of the brain (Zhang et al., 2014; Mure et al., 2018) in a variety of organisms (Boyle et al., 2017; Kuintzle et al., 2017; Hatori et al., 2014). We compared the reported gene expression profile of 5 days post-fertilization (dpf) zebrafish embryos (Li et al., 2013) and liver of the adult fish (Boyle et al., 2017) with recently reported data from baboons (Mure et al., 2018) and mice (Zhang et al., 2014). By comparing these published datasets, generated under free run conditions (dark/dark) in zebrafish and mouse and diurnal conditions for baboon, we found that many key circadian genes showed similar peak expression times in zebrafish and baboon, but not in the nocturnal rodent model (Fig. 1A). We reasoned that diurnal and circadian data can be compared since core clock genes are expected to oscillate irrespective of external cues. Mure et al. (2018) showed that overall transcriptional activity was muted in the baboon liver early in the night. In close agreement, zebrafish liver also shows a biphasic pattern of transcriptional activity (Fig. 1B). Besides these similarities to the circadian regulation of diurnal mammals, zebrafish offer additional advantages such as conserved retinal architecture, rapid development and large clutches of transparent embryos. Despite extensive work on the zebrafish retina (Gestri et al., 2012; Link and Collery, 2015; Avanesov and Malicki, 2010) and circadian behaviour (Cahill, 2002; Hirayama et al., 2005), there is no catalogue of temporally resolved gene expression patterns in the zebrafish retina.

We collected whole eyecups from adult zebrafish exposed to alternating light–dark cycles. Total RNA isolated from the whole eyecup was collected at two time points corresponding to the peak of dopamine and melatonin expression: ZT4 and ZT16 (Fig. 2A); ZT0 represents the switch from dark to light and ZT14 represents the switch from light to dark conditions. RNA-seq libraries constructed from the pooled total RNA of four eyecups at each time point were subjected to Illumina short read sequencing. Zebrafish reference guided transcriptome assembly was performed using TopHat and Cufflinks (Trapnell et al., 2012), and novel transcripts with no previous annotations were identified (Fig. S1). Expression values were assigned to both novel transcripts and the transcripts already annotated in the zebrafish genome version (Zv9). Transcripts were found to be more or less uniformly distributed across the genome. At the extremes, chromosome 5 has the highest and chromosome 24 has the lowest density of the transcripts, respectively (Fig. 2B). In total, 39% (ZT4) and 31% (ZT16) of known genes were identified in our study, suggesting that the zebrafish eye expresses a substantial diversity of genes. This is in agreement with previous reports of mouse and human retinal gene expression profiles (Storch et al., 2007; Pinelli et al., 2016).

To segregate genes that show ultradian fluctuations, we compared the gene expression profiles of zebrafish eyes at ZT4 and ZT16, i.e. 4 h after light (ZT4) and 2 h after onset of darkness (ZT16). The transcripts with significant variation at the two time points were identified using Cuffdiff. After imposing an arbitrary cut-off of 5 FPKM (fragments per kilobase of transcript per million mapped reads) either at ZT4 or ZT16, 120 transcripts were found to be differentially expressed, with 87 showing higher expression during light conditions, at ZT4; while the remaining 37 transcripts were induced during darkness, at ZT16 (Fig. 2C).Ten transcripts amongst this set of diurnal eye transcripts are known to code for clock components (Table S1). Notably, per1b, per3, per2, cry1a, cry1b and nr1d2b showed higher expression during the day, in agreement with previous reports in zebrafish (Huang et al., 2018; Cahill, 2002). Having established that the set of 120 differentially expressed genes show the expected diurnal regulation, we looked more closely at 21 novel transcripts. Seven novel transcripts showed low sequence complexity, making them unsuitable for further validation. The remaining 14 transcripts, named pc1 to pc14 were selected for validation by qRT-PCR. Despite repeated attempts with different primers, pc4 and pc9 could not be amplified consistently. Seven out of 12 tested transcripts (pc1, pc2, pc3, pc5, pc6, pc7 and pc12) showed mutually corroborating differential expression in RNA-seq and qRT-PCR. However, two transcripts (pc13 and pc14) showed discordance between qRT-PCR and RNA-seq (Fig. 3A).

The RNA-seq experiments and qRT-PCR experiments sampled only two time points, one each under light and dark conditions. To establish the diurnal rhythmicity of the selected transcripts, we repeated the qRT-PCR experiments using samples collected at 4 h intervals over three consecutive days with light and dark cycle entrainment (Refinetti et al., 2007). Most of the novel transcripts peaked at ZT4, although pc12, like the positive control gene arylalkylamine N-acetyltransferase-2 (aanat-2) (Gothilf et al., 1999; Falcón et al., 2001; Appelbaum et al., 2006), peaked at ZT16. This is in agreement with the pattern observed in the previous experiments since pc12 was strongly downregulated at ZT4 and upregulated at ZT16 (Fig. 3B).

Next, we separated the retina comprising the thin inner layer of cells and in this fraction enriched for retinal cells and studied expression of the novel genes. All but pc14 showed enrichment in the retina. These genes showed a range of retinal enrichment, some (pc1, pc6 and pc8) even higher than that of aanat-2 (Fig. 4A). In zebrafish, aanat-1 and aanat-2 are both expressed in the retina, but only aanat-2 shows persistent circadian expression in constant darkness (Gothilf et al., 1999).

Next, we carried out qRT-PCR experiments at various developmental stages (28 hpf, 48 hpf, 5 dpf) (Fig. 4B) to identify the peak expression times of these transcripts during zebrafish retinal development. Behavioural responses to visual stimuli suggest that the zebrafish eye starts detecting light between 2.5 and 3.5 dpf. It is well established that the zebrafish optic cups are clearly discernible in embryos at 24 hpf (Easter and Nicola, 1996) while the period between 24 hpf and 36 hpf is marked by rapid cell division and the formation of post-mitotic neurons. By about 50 hpf, the rod cells start expressing opsin genes and by 5 dpf, all five types of photoreceptor cells: rods, short single cones, long single cones, and short and long members of double cone pairs are morphologically distinct. All the novel retinal transcripts identified here, without exception, showed highest expression at 5 dpf (Fig. 4B), the stage by which neurogenesis in the zebrafish retina is complete (Nawrocki et al., 1985).

With the rapid advances in transcriptome sequencing, it is now clear that protein-coding sequences account for only a small fraction of the genes in mammalian genomes. Although to a lesser extent, non-coding RNAs are also abundant in lower organisms where they carry out various functional roles in gene regulation. Ribo-seq profiling identifies the transcripts found in association with ribosomes, thus providing an experimentally verified resource of likely protein-coding transcript (Ingolia et al., 2011). We used publicly available Ribo-seq data (Chew et al., 2013) to classify the novel transcripts identified in this study into potential non-coding RNA and protein-coding transcripts. The chromosome coordinates used to query the GWIPs-viz browser (Michel et al., 2014) that allows visualization of the Ribo-seq data are listed in Table S4. The transcripts pc1 and pc6 were detected by RNA-seq but not in the Ribo-seq, suggesting that they are non-coding transcripts (Fig. 4C,D). The RNA-seq/Ribo-seq data was inconclusive for pc3, pc4 and pc14, due to poor coverage in these experiments. Nine of the 14 transcripts (pc2, pc5 and pc7–pc13) were detected in Ribo-seq experiments, suggesting a protein-coding role. We translated the RNAs and identified predicted proteins arising from these transcripts. The largest novel protein, potentially coded from these transcripts was 1980 amino acids long and showed significant similarity to ryanodine receptor 2 (RYR2) at the protein level. In rodents, RYR2 is known to regulate circadian output from the SCN neuron (Aguilar-Roblero et al., 2007). These observations further reinforced that the newly identified transcripts play roles in circadian gene expression.

As mentioned above, one of the unique advantages of the zebrafish retina is the layered arrangement of cells, which facilitates visualization of cell types and spatial expression patterns by a combination of in situ hybridization and nuclear staining in tissue sections. We carried out whole-mount in situ hybridization (WISH) of the 5 dpf embryo eye to detect the expression of the novel transcripts using labelled antisense riboprobes. Most of the transcripts (pc1, pc5–pc8, pc10, pc13) showed a clear expression pattern in the eye and brain, while no signal was observed with the use of respective sense probes (Fig. 5A). Closer examination of the expression patterns in the retinal layers using cryo-sections showed that all but pc13 showed INL enrichment. pc13, a potential protein-coding transcript was alone expressed in the GCL (Fig. 5A). WISH on 5 dpf embryos, and adult eye cryosections also revealed the diurnal rhythm in the expression of transcripts pc1 and pc10. Both these transcripts showed a rapid clearance of the transcript at ZT16 (dark) following high expression in the INL at ZT4 (light) (Fig. 5B,C).

Spanjaard et al. (2018) recently reported the single cell expression profiling of zebrafish embryos at 5 dpf. We used these data to identify the cell type of the differentially expressed transcripts identified in our study (Fig. S2). Although several novel transcripts were not detected due to the limited coverage of single cell transcriptome dataset, expression of pc2 and pc10 could be unambiguously assigned to retinal R1, R2 and R3 clusters of cells (Fig. 6A,B). In agreement with our in situ hybridization results, these cells are located in the INL of the retina. Furthermore, pc10-expressing cells are also positive for cabp5b, a marker of bipolar cells (Glasauer and Neuhauss, 2016). In future, this expression overlap could be used for cell type-specific functional studies. Next, we used the single cell transcriptomics data to identify the cells of the retina where there is the highest convergence of expression of core clock genes. As shown in Fig. 6C, the majority of the core clock genes show the highest expression in the R3 cell cluster, which also express pc2 and pc10.

Endogenously produced dopamine and light are the two major factors known to regulate retinal diurnal and circadian physiology (Witkovsky, 2004). Hence, we decided to investigate the role of dopamine and light on the differential expression of novel transcripts. Dopamine-producing amacrine cells are present in the INL layer, the main source of dopamine in the retina, can be selectively ablated using a combination of 6-OHDA (6-hydroxydopamine) and pargyline (Li and Dowling, 2000). We used intraocular injections of a mixture of 6-OHDA and pargyline (Fig. 7B) and confirmed ablation by immunohistochemical staining for Th1 (Fig. 7C), a specific marker of dopaminergic amacrine cells. The bilateral uninjected control eye and the injected eye were compared to study the effect of loss of dopaminergic neurons in the presence and absence of light. Fish were deprived of light by extending the dark period for an additional 4 h after the control fishes were transferred to light in the standard day–light cycle (Fig. 7A). The positive control, per2 was upregulated in the presence of light, as reported previously (Besharse et al., 2004). This regulation was abolished in the absence of amacrine cells but could be restored by external intraocular administration of dopamine, clearly showing the combined and mutually dependent effects of light and dopamine on per2 mRNA expression (Fig. 7D). We then extended the analysis to all 14 novel transcripts identified in our transcriptomics study and validated by qRT-PCR. None of the transcripts was affected by light or dopamine, suggesting that the circadian fluctuations seen in our experiments are driven by intrinsic regulation, independent of external cues (Fig. 7E,F).

Next, we compared the expression pattern of these novel transcripts under constant darkness to verify that even without external light, the transcripts show rhythmic expression (Fig. 8A). Total RNA isolated from whole eyecups from fish maintained under constant dark conditions (Koike et al., 2012) (D/D) or entrained (L/D) condition were isolated and subjected to qRT-PCR analysis for known circadian genes (per2, aanat-1 and aanat-2) and novel transcripts identified in this study. per2 was rhythmic only in response to light, whereas aanat-2 continues to show rhythmic expression even in constant darkness (Fig. 8B). Amongst the novel transcripts, eight (pc1, pc3, pc5–pc8, pc10 and pc11) peaked early in the day irrespective of the light cue while pc13 and pc12 peaked during the early night. The latter two transcripts were arrhythmic under constant darkness (Fig. 8C).

The retina is an extension of the central nervous system, composed of a complex network of neurons and glia. By their very function, the cells of the retina undergo rapid, dynamic changes during the course of the day and night cycles. The ability of these cells to anticipate and respond rhythmically, yet tweaking the response according to seasonal changes in day–night cycles, is a critical factor for the survival of the organism. Unlike rodent models, the diurnal nature of zebrafish behaviour makes it an excellent model for human circadian biology. Furthermore, our comparative assessment of core circadian gene peak expression between zebrafish and baboon placed zebrafish more closely to baboon.

Although the zebrafish retinal core clock component shows marked similarities to the mammalian retinal clock, it also differs in specific features. The genome duplication events in the evolutionary lineage of the zebrafish genome suggest that many circadian genes may have duplicated copies that are free to accumulate mutations and develop modified activity. This is even true for many zebrafish core circadian genes which exist as multiple paralogues (Liu et al., 2015).

To the best of our knowledge, this is the first study to report the temporally resolved, circadian gene expression of adult zebrafish whole eye especially highlighting novel retinal transcripts. The sequencing of the human genome led to the surprising insight that only 2% contains protein-coding genes while a significant, previously unknown fraction gives rise to transcripts that do not code for proteins. Currently, the core clock from insects and mammals has predominantly protein components, perhaps since the mechanistic studies on the circadian clock pre-dates the discovery of non-coding RNAs as a significant layer of gene regulation. It is intuitively appealing to speculate that the unbiased discovery of novel circadian transcripts, several of which do not have a strong coding potential may pave the way for a deeper understanding of the regulation of clock components. The novel transcripts reported here include two non-coding RNAs (pc1 and pc6), with absolutely no protein-coding potential or association with ribosomes in any tissue. In light of recent evidence that many authentic long non-coding RNAs are associated with the ribosome (Carlevaro-Fita et al., 2016), it is premature to assume that all the ribosome-associated novel transcripts reported here give rise to proteins. Zebrafish is amenable to knockdown gene experiments to functionally characterize these transcripts in future for their role in circadian gene regulation and retinal physiology.

Besides these known genes, our analysis also revealed 18,000 apparently novel transcripts. It has been shown in several studies that the increased sequencing depth in NGS studies can add a large number of fragmented, low expression transcripts that may be artifacts of the assembly process, and in any case are not amenable to further studies. Although the primary clock mechanism, limited to a handful of genes, is extensively studied, in recent years, a series of transcriptomics studies have shown that an overwhelming number of genes are regulated by the circadian clock (Hughes et al., 2009). The novel players in circadian biology reported here may provide insights into the distinctive mechanisms in operation in the peripheral clocks. The retina is particularly interesting in this context since its input is necessary for SCN circadian rhythm (Paul et al., 2009).

Recent advancement in the single cell transcriptome rendered high confidence cell type catalogues on many organs, including retinal cells in zebrafish. Other than obtaining cell type identity for differentially expressed transcripts, our analysis also indicates that the majority of the novel differentially expressed transcripts from zebrafish eyecups were expressed in the inner nuclear layer. Initially, it was thought that self-sustained circadian clocks of the photoreceptors in the retina drive rhythmic melatonin release. Ruan et al. (2006) showed that the INL and GCL express the core clock genes and in vitro cultured mouse retina continued to show sustained oscillations after several weeks, in spite of the degradation of the photoreceptors. In agreement with the model proposed by Ruan et al., we find that majority of core clock genes in zebrafish retina were expressed in the inner nuclear layer based on the single cell transcriptome. Furthermore, species-specific circadian transcripts may provide an explanation for adaptations unique to the species and its habitat that tweak the basic clock mechanism to produce species-specific variations, as in the case of blind cavefish (Cavallari et al., 2011).

In summary, we present the circadian transcriptome of zebrafish retina and identify several novel clock-regulated genes. Spatially, many of these RNAs are expressed in the INL, which has earlier been shown to harbour cell-autonomous clocks in mice. The temporally resolved expression of these novel transcripts show that at least some of these transcripts show sustained oscillations even in the absence of external cues.

The authors acknowledge Mr Asgar Hussain Ansari for his critical inputs in RNA-seq analysis and R-packages used and Rijith Jayarajan for assistance in sequencing. We also acknowledge CSIR-IGIB zebrafish facility, imaging facility and high computing facility for the infrastructure provided.