Abnormal electrophysiological and molecular phenotypes are present in the retina and the SCNs of Afh mutants

Here we studied how Afh/Afh animals process light, with a particular focus during the transition between light–dark and vice versa.

Therefore, we decided to move retrogradely along the nodes of the RHT (Fig. 2a). Initially, we tested the retinal cone activity of Afh mutants and littermate controls by measuring the flash electroretinogram (fERG) response. The fERG is an electrophysiological measure predictive of the ability of the retina to respond to photoentrainment stimuli [22]. Afh/Afh and+/+mice were exposed to uniform flashes of light of increased luminance generated by a Ganzfeld cupola, and the photoreceptor response was recorded using a gold electrode. A typical fERG waveform consists of ‘a’ and ‘b’ waves, which represent the activity of photoreceptors and bipolar cells, respectively. To measure the cone activity, we removed the rod contribution (i.e., the ‘a wave’) by saturating the background, as previously reported [22]. We exposed mice to the maximum stimulus of 30 cd s/m² of luminance. We observed that the amplitude of the retinal response in Afh mutant mice was significantly decreased (p=0.0124) as compared to that in wild-type mice (Fig. 2b, c), indicating that Afh mutants have reduced responses to light-dependent circadian inputs. The latency of the response was not affected in Afh mutant animals (Fig. 2d).

Open in a separate window

Fig. 2

Altered cone response coupled with a decrease in OPN4 expression in the retina and the SCN of Afh/Afh mice. a Schematic representation of the RHT nodes, starting from the intrinsically photosensitive retinal ganglion cells (ipRGCs) to SCN and hypothalamus. b Flash electroretinogram (fERG) waveform of the cone response in Afh/Afh and+/+animals (n=14 for each genotype) measured in photopic conditions. The saturated negative response and the positive response are represented with letters a and b. c Box plot showing the quantified amplitude values of the b wave between the two genotypes. d Box plot showing the quantified implicit response times after fERG stimulation. e Quantification of western blotting for unglycosylated (53 kDa) and glycosylated (85 kDa) OPN4 isoform protein levels in the two genotypes (+/+n=13, Afh/Afh n=11). OPN4 levels were normalized to β-tubulin levels. f Opn4 expression levels in the SCN in Afh mutants (red) and controls (blue) (n=25+/+, n=24 Afh/Afh). g Relative expression of OPN4 normalized to multiple housekeeping genes at the two different time points, ZT0 and ZT12. The bars represent the average values for three different experiments±SEMs (minimum n=7 for each genotype at each time point). For each box plot (panels C and D), the small square indicates the mean, the central line illustrates the median, and the box limits indicate the 25th and 75th percentiles. The whiskers represent the 5th and 95th percentiles. * p<0.05, ** p<0.01, Student’s t test (c–f), two-way ANOVA plus Bonferroni post-test (g). All image analyses were performed with ImageJ software

Considering that the main functions of the cone cells is associated with visual acuity, and non-visual photic entrainment is only a redundant function of these photoreceptors, it could still be possible that the differences we observed are only incidental to circadian entrainment.

As the photopigment Melanopsin is highly expressed in ipRGCs and plays a crucial role in maintaining circadian entrainment with the external environment by conveying non-visual light signals from the retina to the SCN [23, 24] we investigated OPN4 mRNA as well as the total protein levels in the retinas of Afh/Afh and+/+mice at two different time points in a normal LD cycle, 08.00 (zeitgeber time, ZT0) and 20.00 (ZT12). Since mRNA levels isolated from an heterogenous source as the retina showed high variability (Additional file 4: Figure S3A), we focused our study mainly on protein levels. We first used an antibody able to detect both the unglycosylated (53-kDa) and glycosylated (85-kDa) forms of OPN4 in both wild-type and Afh mice (Additional file 5: Figure S4A). N-glycosylation can be involved in melanopsin maturation and localization. Afh/Afh animals showed significantly lower levels of both 53 kDa (p=0.0046) and 85KDa (p=0.0241) OPN4 expression compared to+/+controls (Fig. 2e). Furthermore, levels of unglycosylated protein were significantly reduced by 30% in Afh/Afh animals compared with+/+controls at ZT0, the beginning of the light phase (+/+1.000±0.09; Afh/Afh 0.64±0.098; p=0.019), but not at ZT12, the end of the light phase (Additional file 5: Figure S4B). Instead, we observed significant reductions in the levels of the OPN4 glycosylated isoform in Afh/Afh mice at ZT12 (+/+1.032±0.121; Afh/Afh 0.641±0.132; p=0.040) but not at ZT0 (+/+1.000±0.112; Afh/Afh 0.812±0.112; p=0.248) (Additional file 5: Figure S4C). We then investigated on the presence of Melanopsin mRNA in the other terminal of the RHT tract, the SCN, due to the paucity of starting material for a biochemical analysis. The SCN tissue was identified using a Six6-specific transcriptional profile (Additional file 6: Figure S5A, B). Opn4 transcript levels were lower in the SCNs of Afh/Afh mice than in those of+/+mice (Fig. 2f:+/+0.85±0.11, Afh/Afh 0.48±0.07, p=0.01). Interestingly, Opn4 levels were low in Afh mutants at ZT0 but returned to wild-type levels by ZT12 (Fig. 2g). Immunohistochemical investigation of SCN slices revealed that OPN4 was expressed in both SCN cells and fibres (Additional file 5: Figure S4D).

Our data suggest an impairment in the expression and modification of both the glycosylated and unglycosylated forms of OPN4 at the transition between the light and dark phase and vice versa. Afh/Afh showed alterations of the fERG activity, in parallel we showed that OPN4 levels were altered also in the SCN of mutant animals at the same shifting moments between light and dark conditions.

Epigenetic phenotypes: Fbxl3 mutation mediates specific DNA methylation changes, including of the Opn4 promoter, in Afh mutant SCNs

Alterations in the OPN4 transcriptional levels in the SCN prompted us to investigate in further details the methylation status of this gene in the SCN. We performed quantitative bisulphite pyrosequencing for Opn4, in Afh/Afh and+/+SCNs collected at both ZT0 and ZT12. An OPN4 CGI was localized in the promoter, a region of a gene body where methylation is associated with repression of transcription [25, 26]. Initially, we showed that a specific region of the Opn4 promoter region was hypermethylated in mutants (Fig. 3a). Afh/Afh mice showed hypermethylation compared to wild-type mice (+/+58.57±1.322; Afh/Afh 64.32±1.322; p=0.0105). The CGI cytosine residues were differentially methylated across the region, especially close to the transcription start sites (TSS) (Fig. 3a). Moreover, methylation level was dynamically altered between ZT0 and ZT12 in the two genotypes, with a significantly higher level at ZT0 in mutant animals compared to controls (p=0.0243).

Open in a separate window

Fig. 3

Discrete and specific differences in the methylation profiles of Afh/Afh SCN vs+/+animals. a Schematic representation of the Opn4 locus highlighting the CGI region (upper panel) at the beginning of the transcription start site (TSS, black arrow), as analyzed by pyrosequencing. Box plot of methylation values for each biological replicate averaged across all of the CpGs in the pyrosequencing amplicon (lower panel); the red lines represent the median, and the whiskers represent the 5th and 95th percentiles. b Pyrosequencing of Opn4 CGI methylation levels in Afh/Afh and+/+mice at the two different time points, ZT0 and ZT12. c Principal component analysis of CGIs methylation vs time of day for+/+and Afh/Afh animals after RRBS sequencing. d Logistic regression of the RRBS data identified 31 differentially methylated CpG islands (CGIs) (with a difference greater than 10%). Quantitation of methylation levels was performed with the SeqMonk software suite considering CGIs with a minimum count of 5 per position and at least 3 observed positions for each CGI across the dataset. The coloured scale on the right indicates the hypermethylated CGIs (in green) and the hypomethylated CGIs (in red). e Relative expression of a subset of genes with differentially methylated CGIs identified by the RRBS analysis at the two different time points, ZT0 and ZT12, normalized to the levels of multiple housekeeping genes. The bars represent the average values of three different experiments±SEM. * p<0.05, ** p<0.01, two-way ANOVA plus Bonferroni post-test (b), Student’s t test (a), two-way ANOVA with Tukey’s post hoc test (e). a–b n=5 for Afh/Afh and n=6 for+/+. c n=6 for each genotype. d n=3 for Afh/Afh and n=3 for+/+at each time point. e Minimum n=7 per genotype at each time point

In order to have a secondary analysis, we conducted a complete genome wide investigation of methylation in the SCN. We profiled genome-wide DNA methylation in the SCNs of Afh mutants as compared to those of littermate controls through reduced-representation bisulphite sequencing (RRBS). RRBS enriches for CG-rich sequences, such as CpG islands (CGIs) that are associated with the majority of gene promoters and provides highly quantitative data with nucleotide resolution [27]. An average of 767,898 CpG sites with at least 5 reads were detected in each dataset (Additional file 7: Figure S6A and Additional file 8: Tables S3 and Additional file 9: Table S4). Of these CpG sites, an average of 66.3% were located inside CGIs, while 33.65% were located elsewhere in the genome (Additional file 7: Figure S6B and Additional file 8: Table S3). An average of 16,715 CGIs was covered across all datasets, with at least 15,562 CGIs detected in each sample (Additional file 8: Table S3). This coverage represents 67.6% of the CGIs reported for the mouse genome [28]. We used principal component analysis (PCA) to verify the clustering of the different datasets, which resulted in no obvious separation between the two groups according to genotype, (Fig. 3c). Next, we analyzed the data by combining all individual datasets by genotype (+/+n=6 and Afh/Afh n=5), and used logistic regression to highlight differentially methylated CGIs. We looked then for methylation differences between the two genotypes greater than 10%. 3 out of the 31 differentially methylated CGIs identified in Afh/Afh (Sept1, Notch3 and Nrp2) that matched known gene regions (Fig. 3d and Additional file 9: Table S4), survived additional Replicate test for statistical difference.

The role of DNA methylation in gene expression control depends on genomic position. Approximately 50% of CGIs in the mouse genome are associated with TSSs, with the remaining 50% equally distributed between intragenic and intergenic areas [29], although these CGIs can demarcate alternative TSSs. We found that differences in DNA methylation occurred predominantly at intragenic CGIs (e.g., Ttc7, Snx9, Cdk2ap1, Best2, Adam11, Opn4, Brd4 and Sept1), although some occurred at CGIs located in intergenic regions (Additional file 10: Table S5). We then further assessed the mRNA levels of genes identified by RRBS that are involved in cell–cell communication, rhythmic processes and transcriptional regulation, and we observed genotype-dependent increases in transcript levels at ZT0 for all the targets investigated, suggesting light-dependent activation differences of these genes between Afh/Afh and+/+mice (Fig. 3e).

Our results indicate that differential methylation of a specific region of OPN4 promoter in Afh/Afh animals, is not associated with massive genome methylation differences as revealed by RRBS. Nevertheless, ZT0 increased expression of genes identified by RRBS and involved in cell–cell communication and coordination mechanisms, suggests an impairment of these pathways in the SCN, specifically during the dark–light transition.

The expression of Dnmts and Tets is altered across the Afh peripheral terminals of the RHT

Because of the abnormal molecular, electrophysiological and epigenetic dynamics in Afh mutants, we sought to explore regulatory epigenetic processes across the RHT nodes. We focused on the expression of enzymes involved in DNA methylation and demethylation. DNA methylation is determined by different families of proteins [30]; DNA methyltransferase 1 (DNMT1) maintains DNA methylation, DNMT3s are responsible for de novo methylation [31], and the ten-eleven translocation (TET) proteins catalyze oxidative demethylation reactions.

We observed that Dnmt1 expression levels were significantly higher in Afh homozygous mice than in wild-type mice in both the retina (Fig. 4a) and the SCN (Fig. 4b). An Increase in Dnmt1 could be detected in the rest of the hypothalamus, although it did not reach significative difference (p=0.28, Fig. 4c). Dnmt3a and Tet3 mRNA levels were upregulated in the retina (Fig. 4a). In the hypothalamus instead only Tet1 (Fig. 4c). These data indicated a specific dysregulation of methylation-mediated processes at the input and output of the circadian pacemaker (SCN and retina), but not to the rest of the hypothalamus.

Open in a separate window

Fig. 4

Afh/Afh mice show specific methylation enzyme expression anomalies in the different nodes of the retinohypothalamic tract (RHT). a Retinal expression levels of the transcripts for genes involved in DNA methylation in the Afh/Afh (red) and+/+(blue) genotypes; (minimum n=7 for each genotype). b Expression levels of methylation enzymes in the SCNs of Afh mutants (red) and controls (blue) (minimum n=18 for each genotype). c Hypothalamic expression levels of the transcripts of methylation enzymes in Afh mutants (red) and controls (blue) (n=12 for each genotype). d RT-qPCR of methylation-mediating enzymes at the two different time points, ZT0 and ZT12; (minimum n=7 per genotype at each time point). All expression levels are normalized to those of multiple housekeeping genes. All bars represent the average±SEM for at least two different experiments. * p<0.05, **p<0.01, two-way ANOVA plus Bonferroni post-test (panel D), Student’s t test (panels A–C)

To test whether methylation-mediating enzymes follow an expression pattern related to the LD cycle in the SCN, Retina and Hypothalamus, we compared time points at ZT0 and ZT12. We observed that the increased Dnmt1 expression in Afh/Afh SCNs occurred mainly at the onset of light (+/+0.92±0.38; Afh/Afh 2.18±1.35; p=0.0077), but not at the offset of light (Fig. 4d).

No difference in methylation enzymes between ZT0 and ZT12 was detected in the Hypothalamus (Additional file 4: Figure S3B). We observed instead similar result for Dnmt1 and Dnmt3a in the Retina, both presenting an increased expression in Afh/Afh animals at ZT0 (Additional file 4: Figure S3C).

Our findings indicate that enzymes involved in de novo methylation (Dnmt3a), maintenance (Dnmt1) and active demethylation processes (Tet1 and Tet3) are altered throughout the RHT terminals of Afh/Afh mice, suggesting a change in the overall epigenetic landscape of these animals in concomitance with dark–light transition.

Primary neuronal culture experiments

Primary cortical neurons were cultured as previously described [55]. Multi-well plates were coated the day before culture using 0.1 mg/ml poly-D-lysine (Sigma-Aldrich) and incubated overnight in a sterile incubator at 37 °C with 5% CO₂. Embryos were individually dissected at E17–E18 to obtain cortices. 2 ml of 0.125% trypsin (Thermo Fisher Scientific) was added to each cortex and HBSS-diluted 0.25 mg/ml DNAsi (Sigma-Aldrich) was added for 30 min at 37 °C. Trypsin digestion was blocked using 2 ml of Neurobasal medium (Gibco) containing 2% B27 supplement, 1% penicillin/streptomycin, 1% L-glutamine (Life Technologies), and 10% heat-inactivated FBS (Gibco). Cells were centrifuged for 5 min at 1200 rpm and then resuspended by pipetting in 2–3 ml of complete Neurobasal medium plus FBS. Cells debris were removed by centrifuging at 700 rpm for 7 min. Neurons were then resuspended in complete culture medium without FBS, counted with trypan blue dye (Sigma-Aldrich) and then plated at a concentration of 500,000 cells/well. For micro-electrode arrays (MEAs) electrophysiological experiments, neurons were plated at a final concentration of 36–40,000 cells/ml. Cells were plated on to 60-channel 6-well and 60-channel single-well MEAs previously coated with poly-D-lysine and laminin (Sigma-Aldrich) to promote cell adhesion (final density of approximately 1,200 cells/mm²), as previously reported [56, 57]. Medium was changed by half twice per week.

Gene expression analysis in primary synchronized neuronal cultures

Gene expression profiles where computed using the 2^-DDCT method. DCT for a circadian time point was computed as the difference in CT between target gene and GAPDH house-keeping gene. DDCT was computed by subtracting the DCT of a time point with DCT computed at circadian time 6. This procedure resulted in a time series of relative expression values for each gene and sample. We fit the collection of time series for individual genes with a sinusoidal function f(t)=L+A sin(\phi+2\pi t/T) parametrized by translation L, amplitude A, phase \Phi and period T in order to estimate the expression periodicity.

Regression results and statistics were obtained by orthogonal distance regression, using the Python library spicy.odr. Statistics includes mean, STD, confidence interval and significance p values value for each regression parameter.

MEA electrophysiology

Electrical activity of 20 DIV neuronal cultures was recorded on both 6-well and single-well MEAs (Multichannel Systems, MCS; Reutlingen, Germany) consisting of 60 TiN/SiN planar round electrodes (30-μm diameter; 200-μm centre-to-centre inter-electrode distance). The activity of all cultures was recorded using the MEA60 System (MCS). Signals were first amplified 1200x, sampled at 10 kHz, and acquired through the data acquisition card and MC_Rack software (both from MCS). Thermal stress, evaporation and osmolarity were constantly monitored using a controlled thermostat (MCS), a polydimethylsiloxane (PDMS) cap and a custom chamber with controlled atmosphere, as previously reported [56]. The experimental protocol consisted of a control phase lasting 180 min during which recording was performed in Neurobasal complete medium (2% B27, 1% penicillin/streptomycin, 1% L-glutamine). Then, 100 nM dexamethasone (Sigma-Aldrich) was added to the cultures by direct pipetting into the medium. Time 0 was set at 24 h after treatment. First ten minutes of recording were discarded to avoid perturbation of the firing rate. The total number of experiments performed by following the above protocol (i.e., the number of recorded wells) was as follows: 24 h recording without treatment,+/+ n=5, number of embryos=4, and Afh/Afh n=5, number of embryos=5; and 24 h recording with treatment,+/+ n=5, number of embryos=4, and Afh/Afh n=5, number of embryos=5. The data were high-pass-filtered at 300 Hz with the online software MC_Rack (MCS) to selectively consider multiunit activity (MUA) only. Spikes were detected online using a fixed threshold multiple of the standard deviation of the noise (-5σn). The data analysis was performed using offline custom software developed in MATLAB (MathWorks) called SPYCODE that presents a number of different tools suitable for multichannel neural recording [58]. Cortical neuronal cultures presented two patterns of activity: random spikes and bursts. The bursting behaviour was characterized by a highly dense packed of spikes usually occurring simultaneously in many channels. The burst detection method was used as previously described [59]. After the identification of spikes and bursts, we analyzed several parameters describing the electrophysiological patterns, such as the firing rate (the mean number of spikes per second calculated on the active channels). Moreover, we analyzed the number of synchronized bursts, called network bursts (NBs), in all active channels using a custom algorithm [57].

The recent studies have found that the majority of NBs are led by a small group of cells called major burst leaders [60], which are defined as neurons that have a probability of conducting more than 6% of the total number of detected NBs. To quantify how the NBs propagate within the network, we calculated the minimum propagation delay between each detected MBL and the other electrodes, called followers [61]. We considered the difference between the delay of the first follower and that of the last follower (Fig. 3d). Finally, we analyzed the spike train correlation among 60 channels. The cross-correlation function represents the probability of observing a spike in one channel i at time t+τ (τ=3 ms) given that there is a spike in a second channel i+1 at time t. We considered only the channels with a peak of correlation higher than 0.1. To quantify the changes in the synchronicity, we evaluated the coincidence index, which represents the ratio between the cross-correlation area around zero (±3 ms) and the total area. If the data had normality and equal variance, a two-way ANOVA was used with Tukey post-test, while if not, data were ranked and then analyzed by ANOVA with post-test.

Circadian clock simulation

We implemented mathematical simulations of the circadian loop as described in a stochastic model of the circadian clock in neurons [21]. All reaction rates and parameters were set as in [21], with the exception of the CRY1 degradation rate, which was modulated to mimic the lengthening of circadian periods due to the after-hours mutation. The code used for all simulations was written in the C-programming language and consisted of a modified version of the Gillespie algorithm [21] with the addition of a random seed generator used for multiple simulations.

In each run of the algorithm, we simulated a population of 100 neurons for 16 days, and we recorded the PER2 levels with a 5-min resolution. For each cell, we extracted PER2 periodicity by fitting the expression traces with a sinusoidal function using the “optimize.curve_fit” function from the Python package SciPy.

To extract the CRY1 degradation rate corresponding to the after-hours mutation, we performed a number of simulations decreasing the degradation rate starting from the wild-type scenario indicated by Ko et al. [21]. Higher degradation rates resulted in shorter circadian periodicity with an exponentially decaying profile. The in silico results were fitted with an exponential ordinary least squares regression model. The resulting exponential equation was used to find the degradation rate that reproduced the after-hours circadian period (approximately 26.7 h [12]).

Flash electroretinogram (fERG) assay

+/+(n=14) and Afh/Afh (n=14) mice were used. Mice were kept in constant darkness for 1 h before midday (12:12 LD cycle) as previously reported [22] and then anesthetized with urethane by intraperitoneal (IP) injection (Sigma-Aldrich). Mice were gently restrained using a mouth bar on an electrophysiological stage as described previously [62]. After adding a few drops of Tropicamide 1% reagent (VISUfarma) to the eyes to dilate pupils, mice were kept in the dark for 10 min. The recording and ground electrodes were then mounted. Visual stimuli, consisting of uniform flickers of light of different luminance values (0.003, 0.01, 0.1, 1, 10 and 30 cd s/m²), were applied using a scotopic background on a Ganzfeld dome (Bioptica Mangoni). To saturate rod response and measure cone activity, fERG was measured in photopic conditions applying a background of 20 cd s/m². After 5 min of background illumination, a single flash with a luminance of 30 cd s/m² was applied. Retinas were harvested at the end of the experiment. The amplitude and implicit responses of the b wave were measured on the electroretinographs.

Western blotting

Two groups of+/+(n=13) and Afh/Afh (n=11) animals were sacrificed by dislocation for biochemical analysis at ZT0 (n=7 and n=6, respectively) and ZT12 (n=6 and n=5, respectively). For protein extraction from the retina, the eyecups were removed, and the retinas were dissected in ice-cold PBS and snap frozen on dry ice. Retinas were lysed by sonication in 150 mM NaCl, 10 mM Tris pH 7.4, 1 mM EGTA, 0.5% Triton X-100 (Sigma-Aldrich), with protease and phosphatase inhibitor cocktail (Roche) and incubated for 30 min on ice. Samples were centrifuged at 20,800×g for 20 min at 4 °C. Supernatants were collected, and protein concentrations were determined using a BCA kit (Pierce) following the manufacturer’s instructions. Equal amounts of proteins were run and separated by 10% sodium dodecyl sulphate polyacrylamide gel electrophoresis (SDS-PAGE) as previously described [63]. Proteins were transferred overnight at 4 °C onto nitrocellulose membranes (GE Healthcare). Membranes were blocked with 5% semi-skimmed milk in PBS (pH 7.4) containing 0.05% Tween 20 (Sigma-Aldrich). Primary antibody incubation was performed overnight at 4 °C using an anti-Melanopsin antibody (rabbit anti-Melanopsin, 1:500, #PA1-780, Thermo Fisher Scientific) in blocking buffer. Membranes were then rinsed at least three times in Tris-buffered saline (pH 7.4) with 0.05% Tween 20 (TBST) and then incubated for 1 h at room temperature with horseradish peroxidase (HRP)-conjugated secondary antibodies (Immunopure peroxidase-conjugated, Thermo Fisher Scientific). Proteins were visualized using Immobilon Western Chemiluminescence Kit (Millipore). Chemiluminescence signals were acquired using an Image Quant LAS 4000 Mini apparatus (GE Healthcare) and densitometric analysis was performed using NIH ImageJ Software [64]. Protein levels are expressed as ratios with respect to β-tubulin levels. The glycosylated OPN4 form was normalized to the level of the unglycosylated form.

RT-qPCR analysis

Transcriptional analysis of the indicated genes was performed using another independent set of 13-week-old+/+and Afh/Afh animals. For the genotype analysis, a cohort of animals of each genotype was sacrificed throughout the circadian day. For time point investigation, one set of mice was culled at ZT0 (at the beginning of the light phase, immediately after the light onset), and one was culled at ZT12 (at the beginning of the dark phase, immediately after the light offset). All tissues were rapidly dissected using a tissue puncher on an ice-cold surface and frozen on dry ice. For the SCNs, experiments were run to produce data from n=24+/+and n=23 Afh/Afh mice for genotype experiments and for n=10 per genotype for the two time points. The prefrontal cortex and the rest of the hypothalamus were harvested as reported above (n=12 for each region and genotype). The eyes were enucleated from 8 mice per genotype, and retinas were rapidly dissected and placed in liquid nitrogen. For the SCNs, total RNA was extracted from snap-frozen tissue using a RNeasy Tissue Micro Kit (Qiagen) and QIAzol reagent (Qiagen) in combination with a Tissue-Lyser apparatus (Qiagen) following the manufacturer’s instructions. RNA samples were quantified with an ND1000 NanoDrop spectrophotometer (Thermo Fisher Scientific). RNA from the retina, prefrontal cortex and hypothalamus was purified and treated as described previously [65]. Reverse transcription of approximately 800 ng of RNA was performed using an ImpromII Reverse Transcription Kit (Promega) according to the manufacturer’s instructions. RT-qPCR was performed using an ABI Prism 7900 RT-qPCR machine (Applied Biosystems) and SYBR Green master mix (Qiagen). The reactions were carried out in triplicate. The data analysis was performed as previously described with minor modifications [65]. All samples were normalized on a panel of three housekeeping genes: Gapdh, ß-actin and Hprt1. The expression levels relative to these housekeeping genes were determined by calculation of the ΔCt, and the data are expressed as 2^−ΔΔCt, where ΔΔCt is the difference between the+/+ZT0 cohort and the other experimental cohorts for the time point experiments and between the+/+pooled samples and the other samples for the genotype analyses. For in vitro experiments, neurons were washed three times with ice-cold phosphate-buffered saline (PBS) solution and lysed with 350 μl of RNeasy Lysis Buffer (Qiagen) or 300 μl of TRIzol (Life Technologies). Cell lysates were collected and pipetted several times through a 22-gauge needle into 1.5-ml microcentrifuge tubes to better lyse the cells, and the lysates were stored at − 80 °C until RNA extraction. RNA was isolated using a RNeasy Micro Kit (Qiagen) according to the manufacturer’s protocol. cDNA synthesis was performed with an ImpromII Reverse Transcription Kit (Promega) according to the manufacturer’s specifications. Real-time PCR was performed and analyzed as above with at least duplicate wells for each sample at each time point. All the primers used in real-time PCR experiments both ex vivo and in vitro are reported in Additional file 12: Table S6.

Reduced-representation bisulphite sequencing (RRBS)

Thirteen-week-old female+/+and Afh/Afh mice were euthanized by dislocation. At least two for each genotype were euthanized at zeitgeber time (ZT) 0 (at the beginning of the light phase, immediately after the light onset) and at ZT12 (at the beginning of the dark phase, immediately after the light offset). SCNs were rapidly dissected using a tissue puncher on an ice-cold surface and frozen in dry ice. SCN tissue punches were validated for region specificity (Additional file 6: Figure S5) using the specific marker gene Six6, which is selectively expressed in the SCNs of adult rodents [66]. The dissected SCNs were used for DNA genomic methylation screening by RRBS and for pyrosequencing assays. Genomic DNA was extracted and purified using a QIAamp DNA Micro kit (Qiagen) following the manufacturer’s instructions. Purified DNA was quantified using an ND1000 NanoDrop spectrophotometer from Thermo Fisher Scientific, and 100 ng was used for RRBS library generation as previously reported with minor modifications [27, 67, 68]. The bisulphite-converted DNA was indexed using KAPA U polymerase (KAPA Biosystems) and RRBS Multiplex TAG primers (10 μM, Sigma-Aldrich). The processed DNA was purified using solid-phase reversible immobilization (SPRI) beads (Agencourt) following the manufacturer’s specifications. A second round of amplification (15 cycles) was performed using KAPA U uracil stalling-free polymerase (KAPA Biosystems). An additional SPRI purification step was performed, and the sample libraries were screened with a Bionalyzer (Agilent Technologies) and then sequenced on an Illumina HiSeq 2500 platform (Illumina). All the primers and probes used for RRBS sample generation are shown in Additional file 13: Table S2. Because of the low complexity at the start of each sequence (MspI fragments), dark cycles were performed (the first 4 bases of each sequence were not recorded). Sequence alignment and methylation calling of the RRBS datasets were performed using Bismark software [69]. CpGs with read depths<5 were discarded. For every analysis, all informative CpGs were used. Mapping efficiency across the two genotypes was high (64.03%±0.38 and 60.87%±0.57 for Afh/Afh and+/+mice, respectively; see Table S3. Moreover, the CpG methylation was 23.3%±0.65 in Afh/Afh mice and 22.3% (±0.65) in+/+mice, confirming coverage levels previously reported at CpG dinucleotides in the mammalian genome [25, 26, 70, 71]. Indeed, the CHG and CHH methylation levels (where H is A, C or T) were~1% in both groups (Additional file 8: Table S3).

To score CpG island (CGI) methylation, cutoffs were applied. The methylation levels were determined for CGIs with information on≥10% of their total CpGs (with a minimum of 3 CpGs) and by averaging individual cytosine methylation levels.

Logistic regression was further applied to identify differentially methylated CGIs using a p<0.05 after correction for multiple comparisons with the Benjamini–Hochberg procedure. A cutoff of a minimum methylation difference of 10% between groups was also used. In addition, a replicate test (t test/ANOVA) (p<0.01) was applied to identify the strongest candidates from the hits obtained by logistic regression.

The dataset analysis was based on the Grcm38 build of the mouse genome and was performed using the SeqMonk software suite. Promoter CGIs were defined as overlapping an annotated transcription start site (TSS), using the University of California, Santa Cruz (UCSC) or Ensemble databases. Intragenic CGIs were defined as overlapping an annotated gene without its TSS. Promoters were defined as the region 2 kb upstream of annotated TSSs. Among the 31 differentially methylated CGIs in Afh/Afh mice, six mapped to intergenic regions and were annotated as ‘null’.

Pyrosequencing

Snap-frozen SCN tissue was treated as previously reported, and 100 ng of genomic DNA was bisulphite converted using an Imprint DNA Modification kit (Sigma-Aldrich, USA) following the manufacturer’s instructions. Pyrosequencing analysis was conducted by NXT-DX (Belgium), as previously reported [72]. The sequences of the primers used in this study are as follows: Opn4F, TTAGTGTGGTTGTTGAGTTG, biotinylation modified; Opn4R, AAAACTTTAAAAATATTCCTATCAC; and Opn4S2, AAAATATTCCTATCACTC.

Immunofluorescence of SCN slices

Wild-type and mutant mice were sacrificed by cardiac perfusion. Animals were anaesthetized using IP injection of 20% urethane solution and were then transcardially perfused with 4% freshly prepared PFA (Sigma-Aldrich) in 0.1 M phosphate buffer (PB; pH 7.4). Brains were extracted, placed in a 2% PFA solution for approximately 2 h and then transferred to a 30% sucrose solution overnight to allow cryopreservation of the brain structures and morphology. Brains were rinsed several times in 0.1 M PBS solution and then cut coronally (50-μm thick slices) using a Microm KS34 freezing microtome (Thermo Fisher Scientific). Free-floating slices were collected in serial order throughout the entire hypothalamic region, considering that the mouse SCNs are between − 0.22 and − 0.82 μm from the Bregma.

Slices were permeabilized and blocked using 3% normal goat serum (NGS) and then 0.3% Triton X-100 in PBS for 2 h. Slices were incubated with OPN4 antibodies overnight at 4 °C (PA1-780, Thermo Scientific, 1/100 primary, 1/500 secondary). Cell nuclei were detected using DAPI diluted 1:300 in PBS for 10 min at room temperature. Images were acquired using an Eclipse Ti A1 confocal inverted microscope (Nikon, Japan). Acquisitions were automatically performed using the motorized Z-stack function of the A1 microscope.

Bioinformatics analysis of transcription factor binding sites (TFBSs)

Gene names of targets of interest from the RRBS screening and of enzymes involved in methylation pathways were imported into the Genomatix software program as previously described [73, 74]. Promoter sequences of the genes were then identified and retrieved using LASAGNA 2.0 software [75]. The promoter regions were further analyzed for common TFBSs using the MatInspector suite of Genomatix [73].

Chromatin immunoprecipitation (ChIP)

Primary cortical neurons from E17–E18+/+embryos were cultured in Neurobasal complete medium supplemented with L-glutamine, penicillin/streptomycin and B27 supplement. Cells were harvested and processed for ChIP assays using a LowCell ChIP kit (Diagenode) according to the manufacturer’s instructions. Chromatin shearing was performed using Bioruptor Plus (Diagenode) for 10 cycles of 20 s on/20 s off. Sheared chromatin equivalent to 100,000 cells was used for immunoprecipitation using 8 µg of anti-Rev-ERBα antibody (13,418 Cell Signaling), and normal rabbit IgG (Sigma-Aldrich) was used as a negative control. The isolated DNA samples were analyzed by RT-qPCR using promoter primers for the following epigenetic targets: dnmt3a (first pair, Fw AACGGTGTCCTTGTCCTC and Rv ATTTCTGCCACCCATAGTCT; second pair, Fw AGGTCTAGTGCCCGTCTG and Rv TGAAGAGGTGGAAGGTTGAAC), dnmt3b (Fw GAGGAACCCAGGTAGTTG and Rv TTCTGCTTCCTGCTTTCA), Tet1 (Fw GCTATTGTTATTTTAGACCCCAAA and Rv ATCTTCCTTTTGAGGAGAATCTG), and Bmal1 (Fw AGCGGATTGGTCGGAAAGT and Rv ACCTCCGTCCCTGACCTACT). The data were normalized to the input and plotted as the % of recovery.

Small interfering RNA (siRNA) experiments

Neurons were obtained from+/+embryos at E17-E18. Cortical neurons at 6 DIV were transfected with Silencer Select Pre-designed siRNA against Rev-ERB (siRNA ID#: s117137) or with control siRNA (Life Technologies) using Lipofectamine RNAiMAX Transfection Reagent (Life Technologies). Neurons were seeded in 6-well plates at a density of 600,000 cells/well. siRNAs were used at 10 μM concentrations according to the manufacturer specifications, as previously reported [55]. Neurons were harvested at 8 DIV and analyzed for gene expression of epigenetic targets.

Drug treatments

Neurons were collected from wild-type littermate embryos at E17–E18 and plated in 12-well plates at a final density of 500,000 cells per well. Neurons were treated with 10 μM GSK 4112 (Sigma-Aldrich) or DMSO in normal Neurobasal complete medium as a control for 24 h. The drug was then removed, and the cells were harvested to investigate the gene expression profiles of DNA methylation genes.

Not applicable.