Single-cell in vivo imaging reveals light-initiated circadian oscillators in zebrafish

Circadian clock is a cell-autonomous time-keeping mechanism established gradually during embryonic development. Here we generated a transgenic zebrafish line carrying a destabilized fluorescent protein driven by the promoter of a core clock gene, nr1d1, to report in vivo circadian rhythm at the single-cell level. By time-lapse imaging of this fish line, we observed the sequential initiation of the reporter expression starting at photoreceptors in pineal gland then spreading to cells in other brain regions. Even within pineal gland, we found heterogeneous onset of nr1d1 expression in which each cell undergoes circadian oscillation superimposed over cell-type specific developmental trajectory. Furthermore, we found that single-cell expression of nr1d1 showed synchronous circadian oscillation under light-dark cycle. Remarkably, single-cell oscillations were lost in animals raised under constant darkness while developmental trend still persists. It suggests that light exposure in early clock development initializes cellular clocks rather than synchronizes existing individual oscillators as previously believed.


Introduction
Circadian rhythm evolves to align animal behaviors to periodic daily environmental changes. At the molecular level, vertebrate circadian clock is mainly generated through transcriptional/translational feedback loops of core clock genes [1]. Among them, two transcription factors (TFs), BMAL1 (also known as ARNTL or MOP3) and CLOCK form heterodimers to bind to E-boxes in the promoters and initiate the transcription of their target genes [2][3][4] including Per family genes (Per1, Per2, and Per3) and Cry family genes (Cry1 and Cry2). The activation of these genes results in the formation of PER/CRY complex and thereby inhibits CLOCK/BMAL1 transcriptional activity, forming a negative feedback loop [5]. Nuclear receptor, REV-ERBα (also known as NR1D1), represses the transcription of Bmal1 and itself is under the transcriptional regulation of BMAL1/CLOCK giving rise to second negative feedback loop of circadian clock. The genome-wide regulation by circadian 3 TFs such as BMAL1/CLOCK and REV-ERBα typically leads to thousands of genes showing circadian expression in a given tissue. Although the basic network of core circadian genes is present in almost every cell, many of the circadian-controlled genes are tissue-specific or cell-type specific. Their circadian expression is a result of either tissue-specific binding of circadian TFs [6] or transcriptional cascade from tissuespecific TFs regulated by circadian TFs [7]. At the organismal level, the overt circadian rhythm is governed by an intricate network of circadian oscillators, in which the master pacemaker such as SCN in mammals or pineal gland in zebrafish is believed to play a pivotal role [1,8]. Traditional in vivo or ex vivo studies of cellular circadian clock have relied on luciferase reporter systems driven by core clock gene promoters [9,10]. But luciferase imaging can only achieve single-cell resolution in organotypic slices in culture with high-resolution CCD camera. Transgenic zebrafish lines carrying luciferase reporters driven by core clock gene promoters have been developed in larval zebrafish. But in vivo luciferase imaging of these fish lines lacked single-cell resolution and can only report the population-level circadian rhythm [9,11]. Genetically encoded calcium indicator (GCaMP) has been widely used to monitor in vivo single-cell calcium activity. In comparison, there is still a lack of zebrafish line to report the circadian expression at the single cell level in vivo.
Circadian rhythm has to be established at molecular, cellular, tissue and behavioral levels during animal development. Day-night rhythms in the fetal rat SCN are first detected between E19 and E21 [12]. Rhythmic expression of circadian clock genes is not detected in both in vitro and in vivo mouse embryonic stem cells (ESCs) but only appears when ESCs differentiate into neural stem cells [13,14]. By ex vivo luciferase imaging of mouse fetal SCN, Carmona-Alcocer et al. showed that a few cells in SCN start circadian oscillations on E14.5, widespread synchronized oscillations were formed on E15.5, and then dorsal-ventral phase wave was established at P2 [15].
Zebrafish (Danio rerio) embryos that are in vitro fertilized and transparent provide an accessible model organism for in vivo live imaging, and thereby being widely used in the study of animal development. A functional circadian clock, characterized by free- 4 running activity, rhythmic cell cycle and circadian gene expression, is established after hatching in zebrafish [16]. However, how single-cell circadian clocks were established during embryonic development is still unclear. It is well-documented that early exposure to light-dark (LD) cycle is required for the development of circadian clock in zebrafish larva [9]. But it is still under debate whether the effect of light on clock development is a synchronization of already existing oscillators or an initiation of single-cell clocks. Questions like these can only be addressed by in vivo imaging of single-cell circadian clock reporter.
Here we report a transgenic zebrafish line using destabilized fluorescent protein, Venus-NLS-PEST (VNP), driven by the promoter of a key circadian clock gene, nr1d1.
This system allows us to monitor the development of single-cell circadian rhythm in live zebrafish larva in a cell-type specific manner. We observed that VNP reporter expression undergoes stepwise onset starting at the photoreceptor cells in pineal gland then spread to cells in other brain regions. Using single-cell RNAseq, we characterized the cell types expressing VNP in the whole brain. Within pineal gland, we found that each cell undergoes circadian oscillation superimposed over cell-type specific developmental trajectories. Under LD cycle, cellular expression of VNP shows synchronous circadian oscillation. However, the circadian expression of nr1d1 positive cells was lost while the developmental trend was still present at single cell level in fish raised under constant darkness. Our result suggests that the effect of early exposure of LD cycle on clock ontogeny is an initiation of functional single-cell oscillators rather than a synchronization of already existing oscillators.

Screening for in vivo circadian reporters in zebrafish
To monitor circadian rhythm at the single-cell level in live animals, we have screened for in vivo circadian reporters among various combinations of destabilized fluorescent proteins driven by core clock gene promoters in larval zebrafish. We first tested peTurboGFP-dest1 (TGFPD1) encoding a destabilized variant of green fluorescent protein TurboGFP. The plasmids of bmal1a/bmal2/per2/nr1d1:TGFPD1 were constructed by homologous recombination. We observed that none of these plasmids were expressed in zebrafish embryos. We then tested the plasmids containing core clock gene promoters driven Venus-NLS-PEST (VNP), another form of destabilized fluorescent protein. We found that the plasmids of cry2b/bmal1a/per2:VNP were not expressed in zebrafish embryos, while per1a:VNP has too low expression to be used for in vivo imaging. For nr1d1:VNP, it has been reported that zebrafish nr1d1 gene has two promoters, i.e. ZfP1 (distal) and ZfP2 (proximal). ZfP1 is conserved and functionally similar to mammalian Nr1d1 promoter [17]. We found that the fish line containing only the proximal nr1d1 promoter ZfP2 (1.5kb) failed to drive VNP expression. But nr1d1:VNP containing both ZfP1 and ZfP2 (6.2kb) showed robust expression in zebrafish embryos ( Fig. 1a and Supplementary Table 1). Therefore, we chose nr1d1:VNP containing ZfP1 and ZfP2 promoters as the circadian reporter and generated the transgenic fish line Tg(nr1d1:VNP) by the Tol2 system. The 6.2-kb ZfP1 and ZfP2 promoter region to drive VNP expression includes the entire set of known cis-regulatory elements: E-box, RRE, Crx-and Otx5-binding sites (Fig.   1a). We measured the mRNA levels of endogenous nr1d1 and VNP expression in the Tg(nr1d1:VNP) transgenic fish from 3.5 dpf (days post-fertilization) to 7.5 dpf under the 12h/12h LD cycles using real-time PCR. nr1d1 and VNP showed highly correlated expression patterns (Pearson's r = 0.9) indicating that VNP expression faithfully reported the expression of nr1d1 at the mRNA level. In addition, both genes showed higher expression level at the dawn than dusk over days (Fig. 1b), which is consistent with our previous result that nr1d1 shows circadian expression peaking at ZT23 (Zeitgeber time) in larval zebrafish [7]. We examined the spatial distribution of fluorescence-labeled cells in nr1d1:VNP at 7.5 dpf using in vivo two-photon imaging system and found nr1d1:VNP-positive cells in many brain regions including the pineal gland, the optic tectum, and the cerebellum (Fig. 1c, Supplementary Movie 1,2). 6 We next monitored in vivo expression of nr1d1:VNP at the whole-brain scale from 3.5 dpf to 7.5 dpf at ZT0 and ZT12 (Fig. 1d) by time lapse imaging. The expression of nr1d1:VNP was sequentially turned on and underwent gradual increase in distinct brain regions during development, as early as 3.5 dpf in the pineal gland, followed by the optic tectum at 5.5 dpf, and other brain regions such as the cerebellum at later time points (Fig. 1e Supplementary Table 2, Supplementary Movie 1,2). In addition, singlecell tracing revealed that circadian oscillations of nr1d1:VNP expression appear to peak at ZT12 across brain regions (Fig. 1e). The expression pattern of nr1d1:VNP closely resembled that of endogenous nr1d1 reported by Delaunay et al. [18]. So nr1d1:VNP fish can be used as a reporter for both developmental expression and circadian expression of endogenous nr1d1 gene in zebrafish.

Characterization of nr1d1:VNP expressing cells
To identify the cell types expressing nr1d1:VNP in the whole brain, we conducted single cell RNA-seq (scRNA-seq) of ~15,000 cells dissociated from the brain of Tg(nr1d1:VNP) larval fish at 6.5dpf. Among them, 6514 cells were identified with more than 500 genes and were used for the downstream analysis. In total, 26 cells clusters were classified from scRNAseq. They were manually annotated by comparing the marker genes with those from whole-brain single-cell RNA-seq data in adult zebrafish [19] (Fig. 2a). We found that the mRNA of nr1d1:VNP was enriched (adjust p-value < 0.01) in cell clusters of photoreceptors in pineal gland, glutamate neurons in forebrain and hypothalamus, granule cells in cerebellum, dorsal habenula cells as well as some non-neuron cells such as oligodendrocyte, retinal-pigment epithelium-like cell, epidemis, and endothelial cells (Fig. 2b).
We next examined the expression of nr1d1:VNP in pineal gland at the single-cell resolution. Fig. 2c showed an example of 3D reconstruction of nr1d1:VNP signals in the pineal gland. At 4.5 dpf, each pineal gland contains about 60 nr1d1:VNP-positive cells. To better understand the cell type of these cells, we crossed Tg(nr1d1:VNP) with Tg(aanat2:mRFP) [11]. Averaged 3D cell density distribution of nr1d1:VNP-positive 7 cells across six fish (Methods section) revealed that nr1d1:VNP-positive cells were mostly distributed around the lumen of the pineal gland while the cell density was relatively low in the central part of the pineal gland (Fig. 2d). We observed that all nr1d1:VNP signals were overlapped with the aanat2:mRFP signals indicating that all nr1d1:VNP-positive cells were melatonin-synthesizing photoreceptor cells in the pineal gland (Fig. 2c, Supplementary Movie 3). According to our scRNA-seq data, VNP is expressed in both rod photoreceptors and cone photoreceptors (Fig 2e). Indeed, we observed co-labeling of VNP with rod cell marker (xops) and cone cell marker (lws2) in pineal gland by crossing Tg(nr1d1:VNP) fish line with rod fish line Tg(xops:nfsB-mCherry) and cone fish line Tg(lws2:nfsB-mCherry), respectively (Fig 2f). Our scRNA-seq data also suggested that nr1d1:VNP is expressed in proliferative cells.

Developmental dynamics of nr1d1:VNP expression within pineal gland
To examine the development of circadian oscillations more closely in pineal gland, we imaged nr1d1:VNP signals in the pineal gland from 3.5 dpf to 6.5 dpf using two-photon imaging (Fig. 3a，b, Supplementary Movie 4-9, Supplementary Table 5). The majority of cells already showed higher nr1d1:VNP signals in dusk than dawn when the signals first become detectable in the pineal gland ( Supplementary Fig. 2). This pattern does not change significantly during development (ANOVA P = 0.78) (Supplementary Fig.   2). However, we found that VNP positive cells within pineal gland display heterogeneous temporal patterns of nr1d1:VNP expression over development (Fig. 3c).
Some cells showed rapid increase in the baseline level while other cells showed more robust circadian oscillations. We then quantified the circadian and developmental components in each cell by fitting the time-series data with a regression model that combines a step-wise function of time, (B*(-1) (x+1) ) denoting circadian oscillation, with a linear function of time, (Ax + C) denoting the developmental increase (Fig. 3d). After fitting the model, we found that the regression coefficients of developmental effect (A) 8 and circadian oscillation (B) both vary widely across cells but there is a negative correlation between them (Fig. 3e).
To reveal if the heterogeneity of single-cell temporal profiles within pineal gland are due to differences in cell type, we next conducted cell-type specific imaging by crossing nr1d1:VNP with the fish line labeling rod cells in red fluorescent protein Tg(xops:nfsB-mCherry). After imaging nr1d1:VNP signals in the rod and non-rod cells simultaneously in the same fish during development from 3.5 dpf to 6.5 dpf (Supplementary Table 9

Single-cell circadian oscillations in pineal gland
We next monitored the nr1d1:VNP signals in the pineal glands of live zebrafish larvae with higher temporal resolution of every 2 hours at 5 dpf using two-photon imaging Their circadian phases were distributed around ZT11 within a narrow range ( Fig. 4d) with Kuramoto order parameter R=0.86 indicating the phase coherence of single-cell circadian oscillators [21]. Clustering analysis of single-cell VNP traces identified two distinct clusters of cells that can be distinguished by their baseline fluorescence intensity (Fig 4e, f, g). Such difference in baseline level of expression at 5dpf can be explained by the cell-type specific developmental trajectories that we observed above.
We also observed that the cluster with higher baseline level (cluster 1) also showed lower relative amplitude of oscillation (Fig 4h). Similar negative correlation between 9 mean transcriptional activity and relative amplitude of circadian oscillations has been reported in mouse [22]. However, there is no significant difference in circadian phase between two clusters ( Supplementary Fig. 3). In short, circadian clocks in pineal photoreceptor cells are oscillating synchronously under LD in spite of large difference in baseline level.

Light-Dark cycle is essential for the onset of cellular circadian oscillation
It is known that LD cycle is required for the development of circadian clock in zebrafish [9]. However, it is unclear whether the role of LD cycle in clock ontogeny is a synchronization of individual oscillators or an initiation of oscillation. To address this question, we imaged nr1d1:VNP signals in zebrafish larvae raised under the constant dark (DD) from 3.5 dpf to 6.5 dpf using two-photon imaging (Fig. 5a, Supplementary Movie 10-13, Supplementary Table 6). Examination of the expression patterns of all nr1d1:VNP-positive cells from four DD fish showed that although the expression of nr1d1:VNP still increased during development, the circadian expression had been absent in DD cells (Fig. 5b). Indeed, regression analysis showed that the oscillating coefficients of majority of DD cells are not significantly different from zero (Fig. 5c).
However, DD cells have even higher developmental coefficients than LD cells (Fig.   5d). In comparison, when we imaged the fish transferred into DD condition after six LD cycles (LD_DD cells) (Fig 5e，Supplementary Table 7 Table 4). We found that DD cells showed much lower oscillating amplitudes than those in LD cells (Fig. 4i, j). Only 1 out of 142 DD cells showed significant oscillation (p-value<0.05 and relative amplitude>0.1, Fig 5k). Therefore, our result suggests that the cellular circadian oscillations in DD fish were not initiated 10 rather than not synchronized during development.

Discussion
In our study, we constructed a zebrafish nr1d1:VNP fluorescence reporter system and for the first time we can monitor the circadian gene expression at the single-cell resolution in zebrafish larvae. Using this reporter fish, we revealed the interplays among circadian clock, cell-type specific development, and light-dark cycle. The mouse version of nr1d1:VNP was originally developed by Nagoshi et al. in NIH 3T3 cells.
VNP was chosen for its high fluorescence intensity and high folding efficiency [23].
With this reporter, they were able to show that each NIH 3T3 fibroblast cell has a selfsustained circadian oscillator that can be entrained by serum shock. Also using this reporter, two independent studies have demonstrated the tight coupling between circadian rhythm and cell cycle in proliferating NIH 3T3 cells [24,25]. As FUCCI lines to label different stages of cell-cycle are available in zebrafish [26], one can cross nr1d1:VNP fish with cell-cycle reporter lines to examine the relationship between circadian clock and cell cycle in vivo. By crossing fish lines containing cell-type specific fluorescent markers with our nr1d1:VNP fish, one can conduct cell-type specific imaging of circadian rhythm as we have done for rod cells in pineal gland.
Larval zebrafish has also been used in drug screening for chemical compounds affecting circadian rhythm and sleep [27]. Our Nr1d1:VNP fish can be used to investigate the effect of selected drug targets on the synchronization of single-cell oscillators in vivo.
With advanced microfluidic systems to capture and recapture larval zebrafish, one can conduct more intense imaging on the same fish to achieve higher temporal resolution while minimizing the perturbations to the fish [28]. Using newly developed movable platform to track zebrafish larvae and novel volume imaging system, live imaging of cellular circadian rhythm on free-moving zebrafish will become possible [29]. In summary, our single-cell circadian reporter zebrafish line will have broad applications in circadian research.
It was believed that physiological and behavioral rhythms in zebrafish appear during the early stages of zebrafish development [8]. In our study, the expression of nr1d1 appeared gradually in development and step-wise in different brain regions. As the major component of the second negative feedback loop in circadian clock, nr1d1 is directly regulated by BMAL1 and represses the expression of bmal1 [1]. The onset of

Construction of the nr1d1:VNP (Venus-NLS-PEST) Expression Vector
The zebrafish nr1d1:VNP reporter plasmid was constructed as follows. The mouse Nr1d1:VNP reporter plasmid was generously provided by Emi Nagoshi [23] (Department of Genetics and Evolution, Sciences III, University of Geneva). The plasmid carries an intact VNP sequence. VNP (Venus-NLS-PEST) is a nuclear fluorescent protein with a high folding efficiency and short half-life [24]. The zebrafish Nr1d1:VNP plasmid was obtained by homologous recombination with four cassettes.
Targeting cassette was generated by four-step PCR amplification. The first PCR step was performed to amplify the tol2-polyA-Ampicillin gene with an additional 3′ sequence that is homologous to the 5 ′ end region of zebrafish nr1d1 promoter arm of the first PCR was homologous to the 5′ arm of the second PCR, the 3′ arm of the second PCR was homologous to the 5′ arm of the third PCR, the 3′ arm of the third PCR was homologous to the 5′ arm of the fourth PCR, and the 3′ arm of the fourth PCR was homologous to the 5′ arm of the first PCR. The vector construct was verified by sequencing. Fig. 1a shows the schematic of nr1d1:VNP construct design. The structure of 6.2k-bp nr1d1 promoter containing multiple E-box, RRE, and pineal specific elements was also described.

Generation of transgenic Tg(nr1d1:VNP) zebrafish reporter line
Generation of transgenic fish using the Tol2 system was carried out according to the published protocol [33]. Constructs were injected together with capped RNA encoding transposase (10 ng/μl each of DNA and RNA) into fertilized eggs at the one-cell stage.
The injected fish (F0 generation) were raised and screened for integration of the transgene into the germline. Isolation of transgene-positive progeny (F1) was carried out using EGFP imaging using fluorescence stereomicroscpe (Olymous, SZX16). For in vivo imaging, Tg(nr1d1:VNP) zebrafish was crossed with a pigmentation mutant strain (casper mutant), which is transparent in the whole brain except the eyes. Vazyme, China). The pTol-lws2:nfsB-mCherry plasmid and pTol-xops:nfsB-mCherry were co-injected into AB embryos with Tol2 transposase mRNA at one-cell stage. The mCherry signal was screened to identify F1.

In vivo imaging
Fish were randomly selected to be imaged under LD or DD. We were not blinded to LD and DD fish group allocation. For time-lapse imaging in every two hours lasting for one day starting at ZT0, 5.0dpf zebrafish larvae were anesthetized with MS-222 (sigma-Aldrich) (0.01-0.02%) and mounted in low melting-point agarose (1.33%). The fish were maintained in a chamber with circulating egg-water without anesthesia at temperature of 28 ± 0.5 °C. For imaging from 3.5dpf to 6.5dpf every 12 hours at ZT0 15 and ZT12, the zebrafish larvae were anesthetized at each time point with MS-222 (sigma-Aldrich) (0.01-0.02%) and imbedded in low melting-point agarose (0.8-1.0%) before imaging. 3-4 fish were imaged in each imaging session that lasts for half to one hour. Fish were released immediately after the end of each imaging session and in a free-moving state between image sessions. All fish were imaged from dorsal view using two-photon microscope (Olympus) as described previously [34]. Imaging was performed using a 25× (numerical aperture (NA)=1.05) water-immersion objective (Olympus). Excitation was provided by a Ti:Sapphire femtosecond pulsed laser system (Coherent) tuned to 900 nm, which allowed efficient simultaneous excitation of Venus fluorescent protein. Laser power was set to 1.6% for pineal imaging and 2.0% for whole brain imaging. We use two-photon excitation microscopy for two reasons as suggested by Carvalho and Heisenberg [35]. First, the near-infrared wavelength has undetectable effects on fish physiology and behavior. Second, as the two-photon excitation is only achieved near the focal plane, it minimizes photo-bleaching and photo-toxicity. The sample size estimate is based on our previous studies. Two fish under light-dark (LD) condition were imaged for the whole brain imaging. Two fish under LD condition and three fish under dark-dark(DD) condition were imaged for one day every 2 hours. Six fish under LD condition and four fish under DD condition were imaged from 3.5 dpf to 6.5 dpf every 12 hours. Three fish under LD_LD condition and Three fish under LD_DD condition were imaged from 6.0 dpf to 7.5 dpf every 12 hours. Six rod fish were imaged from 3.5 dpf to 6.5 dpf every 12 hours.

Single cell RNA-seq and data analysis
6dpf larval heads were dissected on dissection medium (DMEF/F12 with 2% 100X penicillin-streptomycin) and pineal region were enriched by pipetting the pineal into the tube. Dissociation of the brain cells following the protocol from Miguel A. Lopez-Ramirez et al. [36]. Briefly, Add 300ul papain solution to the dissected tissue at 37 in water heater for 15 minute, pipetting during digest. Then stop the digestion by adding 1.2ml of washing solution. Washing the cells twice using washing solution. In the end, sterilize the cells using 40um pore size filter. Stain cells using trypan blue and count the living cells using a hemocytometer. All the reagents and solutions used were following the materials provided by Miguel A. Lopez-Ramirez et al. [36]. The resulting single cell suspension was promptly loaded on the 10X Chromium system using Chromium Single Cell 3' Reagents v2.
For data analysis, raw sequencing data was converted to matrices of expression counts using the cellranger software provided by 10X Chromium. Reads were aligned to a zebrafish reference transcriptome (ENSEMBL Zv11, release 95). The gene expression matrix were then loaded into Seurat package [37] in R for the following analysis. Cells with less than 500 genes or percentage of mitochondiral genes>0.02 were excluded from the following analysis. The left cells were clustered using the first 50 PCs and with resolution=1, 26 clusters were identified. They were manually annotated by comparing the marker genes with those from whole-brain single-cell RNA-seq data in adult zebrafish [19].

Image analysis
The image data were first converted to 'nrrd' format using a customer Fiji [38] macro code for downstream analysis. Time series of 3D images in 'nrrd' format for each fish were then aligned by CMTK toolkit (parameters: -awr 01 -l fa -g -T 8 -X 10 -C 1 -G 20 -R 2 -E 1 -A '--accuracy 5 --dofs 12' -W '--accuracy 5' for the alignment of pineal gland; -ar 01 -l a -A '--exploration 50 --accuracy 5 --dofs 12' for the alignment of whole brain) to facilitate the single cell tracing. TrackMate [39], an open source Fiji [38] plugin, was applied to perform the single cell nucleus identification and tracing on 3D images over time. Briefly, Laplacian of Gaussian (LoG) filter with a sigma suited to the estimated blob diameter (5.5 micron) was applied to detect the single cell nuclei in each 3D image.
Then HyperStack Displayer was applied to visualize the identified spots, which allows manual editing afterwards. A nearest-neighbor algorithm (parameter: maximal linking distance 5.5 micron) was applied to trace each single cell. After the automatic processing, the identification and tracing of each cell were manually curated using the manual tracking tools of TrackMate [39]. Combining the automated and manual tracking approaches can significantly increase the detection accuracy of single cells. In the end, the cell position and mean image intensities were exported for the following statistical analysis. For the whole brain data, the pineal gland, the optic tectum, and the cerebellum were first manually extracted, then single cells were traced using manual tracking tools of TrackMate [39], and mean intensities of each brain region was used for the following analysis. The cell position and mean intensities for all the fish analyzed were listed in Supplemental Table 2-9.

3D reconstruction of zebrafish pineal gland
To systematically analyze the 3D distribution of the nr1d1:VNP-positive cells, a common 3D pineal gland structure was reconstructed by averaging six zebrafish pineal glands at 5.5dpf. First, we manually reconstructed the 3D pineal structure of six individual 5.5dpf zebrafish pineal glands taking advantage of the clear pineal gland boundary in two-photon imaging. Secondly, all 3D pineal gland structures from different fish were aligned using rigid transformation using CMTK (parameters: -ar 01 -l a -A '--accuracy 5 --dofs 12'). Thirdly, the signals from aligned 3D pineal glands were added and smoothed by Gaussian blur (Sigma = 30) using Fiji package [38]. In the end, a common 3D pineal gland structure was obtained by implementing Li's minimum cross entropy thresholding method using Fiji package [38].
The cells of each pineal gland were then registered to the common 3D pineal gland using the same transform matrix for the alignment of the 3D pineal gland structure. The 'ks' package in R (https://www.R-project.org/) was applied to calculate the cell density distribution. The common zebrafish pineal structure and nr1d1:VNP cell density distribution were visualized by NeuTube software [40] (Fig. 2d).

Statistical analysis
All the statistical analysis was performed using the computing environment R (https://www.R-project.org/). For the time-lapse imaging data from 5.0 dpf to 6.0 dpf every 2 hours, the circadian phase and amplitude of each cell were calculated by fitting the mean intensity of each cell to cosine functions with 24 hours' period and shifting phases as described previously [41]. The fluorescence intensities were log2-tranformed before fitting. Relative amplitude were defined as absolute amplitudes divided by means. R package 'circular' was applied to calculate the average circadian peak between the LD and the DD cells (Fig. 3g, h and 4g, h).     between baseline intensity and relative oscillation amplitude for fish1 (i) and fish2 (j), respectively. *** represents Pearson's correlation p-value <0.001.