# Remodeling of the cycling transcriptome of the oyster Crassostrea gigas by the harmful algae Alexandrium minutum

Laura Payton, Mickael Perrigault, Claire Hoede, Jean-Charles Massabuau, Mohamedou Sow, Arnaud Huvet, Floriane Boullot, Caroline Fabioux, Helene Hegaret, Damien Tran

<figure><img src="image_1.png" /></figure>

## To cite this version:

Laura Payton, Mickael Perrigault, Claire Hoede, Jean-Charles Massabuau, Mohamedou Sow, et al. "Remodeling of the cycling transcriptome of the oyster Crassostrea gigas by the harmful algae Alexandrium minutum." Scientific Reports, 2017, 7, 10.1038/s41598-017-03797-4. hal-01605998

HAL Id: hal-01605998

https://hal.science/hal-01605998

Submitted on 26 May 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

OPEN

# Remodeling of the cycling transcriptome of the oyster *Crassostrea gigas* by the harmful algae *Alexandrium minutum*

Laura Payton<fnref n="1" />,², Mickael Perrigault<fnref n="1" />,², Claire Hoede³, Jean-Charles Massabau²,⁴,⁵, Mohamedou Sow<fnref n="1" />, Arnaud Huvet⁴, Floriane Boullot⁴,⁵, Caroline Fabioux<fnref n="1" />,⁵, Hélène Hegaret⁴,⁵ & Damien Tran<fnref n="1" />,²

As a marine organism, the oyster *Crassostrea gigas* inhabits a complex biotope governed by interactions between the moon and the sun cycles. We used next-generation sequencing to investigate temporal regulation of oysters under light/dark entrainment and the impact of harmful algal exposure. We found that ≈6% of the gills' transcriptome exhibits circadian expression, characterized by a nocturnal and bimodal pattern. Surprisingly, a higher number of ultradian transcripts were also detected under solely circadian entrainment. The results showed that a bloom of *Alexandrium minutum* generated a remodeling of the bivalve's temporal structure, characterized by a loss of oscillations, a genesis of de novo oscillating transcripts, and a switch in the period of oscillations. These findings provide unprecedented insights into the diurnal landscape of the oyster's transcriptome and pleiotropic remodeling due to toxic algae exposure, revealing the intrinsic plasticity of the cycling transcriptome in oysters.

Received: 10 January 2017

Accepted: 4 May 2017

Published online: 14 June 2017

The physiological and behavioral processes of living organisms oscillate. Temporal organization is a necessary adaptation to cope with an ecosystem governed by various periodic changes<fnref n="1" />. At a cellular level, the circadian clock consists of a self-sustainable auto-regulatory network of transcriptional and translational feedback loops synchronized by environmental zeitgeber, which produce oscillations in the expression of clock-controlled genes (CCGs)². Most studies on the circadian network have used terrestrial organisms, while little is known about these networks in marine organisms, which inhabit complex ecosystems governed by sun-earth-moon trajectories and interactions. The study of clocks in marine species is a relevant issue for understanding how clocks evolved along different evolutionary lineages³⁻¹⁰.

In the phylum of marine molluscan, the Pacific oyster *Crassostrea gigas*, an ecologically and economically important marine species, exhibits *in situ* behavior following not only a circadian rhythm but also ultradian and infradian rhythms¹¹, ¹². An endogenous circadian rhythm has been demonstrated in free running conditions and is characterized as weak, plastic, and dual¹³, ¹⁴. Recently, components of a biological clock in *C. gigas* have been identified¹⁵ without any information on the extent of clock output¹⁵. The first study in marine mollusks showed that 40% of microarray transcripts are susceptible to oscillation in the mussel *Mytilus californianus*³. In mammals, next-generation sequencing has indicated that approximately 50% of the genome has the potential to oscillate; based on a sophisticated regulation of gene expression in each peripheral organ, which is still misunderstood, this estimation has been reduced to 10% in a given tissue¹⁶⁻²². In parallel, alterations in the temporal organization of physiological and metabolic processes have been associated with disorders such as obesity, cardiac disease, cancer or disruption of energy balance²³⁻²⁹.

<footnote n="1">¹University of Bordeaux, EPOC, UMR 5805, F-33120, Arcachon, France. ²CNRS, EPOC, UMR 5805, F-33120, Arcachon, France. ³Plate-forme bio-informatique Genotoul, MIAT, Université de Toulouse, INRA, F-31326, Castanet-Tolosan, France. ⁴Ifremer, Laboratoire des Sciences de l'Environnement Marin (LEMAR), UMR 6539 UBO/CNRS/IRD/IFREMER), CS 10070, F-29280, Plouzané, France. ⁵Laboratoire des Sciences de l'Environnement Marin (LEMAR), Institut Universitaire Européen de la Mer, Université de Bretagne Occidentale, UMR 6539 CNRS/UBO/IRD/IFREMER, F-29280, Plouzané, France. Correspondence and requests for materials should be addressed to D.T. (email: d.tran@epoc.u-bordeaux1.fr)</footnote>
In coastal areas, marine organisms are subjected to various stresses, including proliferations of planktonic microalgae producing toxins, called harmful algal blooms (HABs), reaching up to millions of cells per liter. Among these planktonic algae, the dinoflagellate species *Alexandrium minutum* is known to produce saxitoxin (STX) and its derivatives, called paralytic shellfish toxins (PSTs). *A. minutum* has a worldwide distribution and has been associated with massive deaths of marine animals and human health issues³⁰⁻³². The oyster *C. gigas*, as a filter feeder, is highly susceptible to HAB exposure by feeding on harmful algae³³. Much physiological damage in response to *A. minutum* exposure has already been demonstrated in oysters¹⁴, ³⁴⁻³⁹. In particular, it has been shown that *A. minutum* leads to a loss of daily cycle in the expression of genes involved in detoxification, oxidative stress, valve behavior, and crystalline style length, as implicated in digestion⁴⁰.

In the present study, we used RNA sequencing (RNAseq) to investigate for the first time, to our knowledge, (1) the manner in which the whole transcriptome in gills of *C. gigas* is temporally organized under light/dark entrainment and (2) the impact on this temporal organization of an *A. minutum* exposure mimicking a realistic HAB.

## Results

**A circadian transcriptome with a nocturnal pattern remodeled by *A. minutum* exposure.** To assess the cycling transcriptome of *Crassostrea gigas*, a first group of oysters was entrained to a 9-hr light/15-hr dark (L:D 9:15) environment, without tidal cycles, and fed with the non-toxic alga *Heterocapsa triquetra* (*H.t* condition). To explore the effect of PST ingestion, a second group of oysters was exposed to the PST-producing alga *Alexandrium minutum* (*A.m* condition), mimicking a harmful algal bloom (≈400 cells·mL⁻¹). ELISA assays confirmed a PST bioaccumulation of 95.8 ± 13.5 ng·g⁻¹ eqSTX in the gills in the *A.m* group, whereas the *H.t* oysters accumulated no PST (Supplementary Fig. S1). Gills were sampled every 4 hr over 52 hr, and the cyclic transcriptome was analyzed by RNA sequencing (RNAseq) followed by ARSER⁴¹ to evaluate any significant cyclicity in the 20,846 transcripts. For all cycling detection, a false discovery rate (FDR) of <0.05 was considered significant⁴², ⁴³.

In the *H.t* oysters, 1300 (1247 + 53) significant circadian transcripts were identified, corresponding to 6.2% of total transcripts in the *H.t* gills (Fig. 1A–C; Supplementary Table S2). Among these 1300 transcripts, only 53 oscillated in both the *H.t* and *A.m* oysters (Fig. 1A,C). The total number of circadian transcripts in the *A.m* group (630) was two-fold lower than that in the *H.t* oysters (Supplementary Table S2). A large number of transcripts (577) oscillated exclusively in *A.m* (Fig. 1A,D). In *H.t*, the phase lag distribution of circadian transcripts showed a nocturnal and bimodal pattern, with one peak during the beginning of the night (N1, ZT9-ZT19, 40% of circadian transcripts) and another peak at the end of the night (N2, ZT19-ZT24, 47% of circadian transcripts) (Fig. 1E). Finally, only 13% of transcripts peaked during the daytime (D, ZT0-ZT9). In *A.m*, the nocturnal phase peak transcripts decreased to 63% (N1+N2), and diurnal transcripts reached 37% of circadian transcripts (Fig. 1G). An increase in daytime peak expression with *A. minutum* exposure was also observed in the 53 circadian transcripts common in both *H.t* and *A.m* (Fig. 1F). Furthermore, we observed a trend toward a relative decrease in significant transcripts with low amplitude in favor of transcripts oscillating with high amplitude when exposed to *A.m* compared to *H.t* (Fig. 1H–J).

**A nocturnal pattern of circadian core clock gene expression and valve activity behavior modeled with *A. minutum* exposure.** Oyster clock gene expression in the gills of oysters was analyzed by RT-qPCR and exhibited nocturnal phase peaks in *H.t* (Fig. 2A,B). There were significant differences between day and night phase expression for *Clock*, *Cryptochrome 1* (*Cry1*), *Period* and *Timeless* (*Tim1*), but not for *Ror*, *Rev-erb*, *Cryptochrome 2* (*Cry2*) or *Bmal*, peaking just before the night or at the beginning of the day (Fig. 2A, *H.t*, first column of histograms). Fig. 2A and B reveal an effect of *A. minutum* on the clock gene temporal patterns previously observed in *H.t*. First, the overall levels of expression were higher in *A.m* for *Bmal* (*p* = 0.026) and *Tim1* (*p* = 0.038) (Fig. 2A). In addition, none of the 8 clock genes showed differential expression between day and night in *A.m* (Fig. 2A, *A.m*, first column of histograms), leading to a loss of nocturnality. Furthermore, in *A.m*, the amplitude of oscillation decreased for *Bmal*, *Period*, *Tim1* and *Cry2* and increased for *Ror* and *Rev-erb* (Fig. 2A, second column of histograms). A two-way ANOVA confirmed a day/night effect on clock gene expression in *H.t* (*p* < 0.001) but not in *A.m* (*p* = 0.492) and a significant effect of *A. minutum* on gene expression during the night (*p* < 0.001) (Fig. 2C). Finally, significant correlation (*p*-value < 0.0001 for both conditions) between qPCR and RNAseq results on global clock gene expression in *H.t* and *A.m* conditions, comforting results of RNAseq analyses (Supplementary Fig. S2).

Figure 2D shows the mean valve opening duration (VOD) in the *H.t* and *A.m* oysters. Control oysters exhibited cycling and nocturnal valve activity with a rapid increase in VOD at the beginning of the night, corresponding to the N1 phase illustrated in Fig. 1E. This increase in VOD was followed by a decrease to basal levels, observed during the rest of the night and during daytime. The *A.m* oysters exhibited conserved cyclic valve activity, but the mean VOD was significantly increased (*p* < 0.001) and a trend for diurnal activity was observed corresponding to the D phase, shown in Fig. 1G. Figure 2E shows that during the day, the VOD increased from 20.6 ± 1.2% in *H.t* to 38.0 ± 2.6% in *A.m* (*p* < 0.001), whereas during the nighttime, no significant difference in VOD was observed (*p* = 0.305). Finally, the nocturnality observed in *H.t* was disrupted in *A.m*. Indeed, the day/night difference in *H.t* (*p* < 0.001) was abolished in *A.m* (*p* = 0.763).

**Discovery of an ultradian transcriptome remodeled by *A. minutum* exposure.** Analysis of the cyclic transcriptome provided evidence for a large amount of transcripts with ultradian periodicity. Within this range, two major clusters of transcripts appeared in both *H.t* and *A.m* conditions (Fig. 3A): transcripts oscillating with a period in the range of 8–11 hr (Ultradian 1, U1) and transcripts oscillating with a period in the range of 12–16 hr (Ultradian 2, U2). For the U2 transcripts, the peak of period length distribution was shortened in *A.m* oysters. Ultradian transcripts (U1 + U2) represented 18.1% of the transcriptome in *H.t* and 22.4% in *A.m* (Supplementary Table S2). In *H.t*, 2185 (1891 + 294) transcripts were U1 (Fig. 3B–D) whereas in *A.m* oysters,
<figure><img src="image_4.png" /><figcaption>**Figure 1.** A circadian transcriptome with a nocturnal pattern remodeled with *A. minutum* exposure. (A) The number of circadian transcripts. (B–D) Heat maps of expression organized by phase for the 300 most significant transcripts oscillating only in *H.t* (B), common in both (C) and only in *A.m* (D). (E–G) Phase lag, i.e., peak of expression along the cycle, distributions (frequency in number of transcripts) and corresponding pie charts for all oscillating genes in *H.t* (E), *A.m* (G) and both (F). “D”: phase lag distribution during the day (ZT0-ZT9). “N1” (ZT9-ZT19) and “N2” (ZT19-ZT24): phase lag distribution during the first and second part of the night. (H–J) Amplitude distribution (frequency in number of transcripts) for all circadian transcripts in *H.t* (H), *A.m* (J), and both (I), classed by amplitude range.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 2.** A nocturnal pattern of circadian core clock gene expression and valve activity behavior modeled by *A. minutum* exposure. (A) Real-time qPCR analysis of circadian clock gene expression expressed in fold change (FC, mean ± SE, n = 5–6) for *H.t* and *A.m*. Gray areas correspond to nighttime. The *p*-values show different levels of expression over 28 hr between *H.t* and *A.m*, while stars indicate the different levels of expression at each time point. The first histograms show the mean levels of gene expression (FC). Different letters correspond to significant differences (*p* < 0.05). The second histograms show amplitudes of clock gene expression in *H.t* and *A.m*. (B) Heat maps of clock gene expression over 28 hr, organized by phase for *H.t* and *A.m* (RT-qPCR). (C) Two-way ANOVA testing the effect of day/night and *H.t/A.m* conditions on overall clock gene expression. (D) Mean hourly valve opening duration (VOD, %; mean ± SE; n = 16 per condition) of *H.t* (blue) and *A.m* (orange) oysters during the two days of sampling. The *p*-value corresponds to the difference in the mean VOD between the two conditions. (E) Mean VOD during daytime and nighttime for each condition. Different letters indicate significant differences (*p* < 0.05).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Discovery of an ultradian transcriptome remodeled with *A. minutum* exposure. (A) Period distribution of ultradian transcripts in *H.t* and *A.m*. (B) The number of U1 transcripts in *H.t* and *A.m*. (C–E) Heat maps of expression organized by phase for the 300 most significant U1 transcripts oscillating only in *H.t* (C), common in both (D), and only in *A.m* (E). (F) The number of U2 transcripts. (G–I) Heat maps of expression organized by phase for the 300 most significant U2 transcripts oscillating only in *H.t* (G), common in both (H), and only in *A.m* (I). (J) Venn diagrams of circadian (C), U1 and U2 transcripts in *H.t* and *A.m*.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 4.** Ultradian oscillations observed in clock subsystem elements. Heat map expression of orthologous genes implicated in the clock network in mammals and insects, organized by phase of oscillation in *H.t.* "Rhythms", results of cycling analysis (circadian, U1, U2). "Mean level", the mean level of expression (mean ± SE). \*Significant difference at *p* < 0.05.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 5.** Functional implications of cycling transcripts. Enrichment analysis of biological processes of total cycling transcripts (circadian, U1, and U2) in *H.t* (**A**), in *A.m* (**B**), and uniquely in *A.m* (**C**). The ordinate axis shows the GO terms and significance (FDR) of enrichment while the abscissa axis shows the percentage of cycling for each GO term.</figcaption></figure>

2452 (2158 + 294) transcripts were U1 (Fig. 3B,D,E), representing 10.5% and 11.8% of the transcriptome in *H.t* and *A.m*, respectively (Supplementary Table S2). Note, only 294 U1 transcripts were common to *H.t* and *A.m* (Fig. 3B,D).

Subsequently, in *H.t*, we observed 1576 (1348 + 228) U2 transcripts (Fig. 3F–H), whereas in *A.m*, 2207 (1979 + 228) transcripts were U2 (Fig. 3F,H,I), representing 7.6% and 10.6% of the transcriptome, respectively (Supplementary Table S2). Finally, only 228 U2 transcripts were common to *H.t* and *A.m* (Fig. 3F,H). A gene-to-gene comparison of cycling transcripts in *H.t* and *A.m* showed that *A. minutum* exposure led to a switch in the period range of many transcripts (Fig. 3J). Indeed, 13% and 14% of circadian transcripts in *H.t* became U1 and U2 in *A.m*. By contrast, 14% of U1 in *H.t* became U2 in *A.m*, while 15% of U2 in *H.t* became U1 in *A.m*. When exposed to *A. minutum*, the switch from an ultradian period to a circadian one was scarce, 4% for both U1 and U2 transcripts.

**Ultradian oscillations observed in circadian clock network elements.** A cycling transcriptome is the result of sophisticated regulation achieved by an integrated and complex network of core clock genes, kinases, deacetylases, and transcription factors (TFs)², ⁴⁴, ⁴⁵. In *C. gigas* (Fig. 4), we identified 19 orthologs of vertebrate and insect genes associated with the circadian network²¹, ⁴⁴, ⁴⁶⁻⁵⁵. Circadian oscillations were observed for 2 of these orthologs in *H.t*: the TFs *AP-1* and *CTCF*. Interestingly, in *H.t*, 3 transcripts were ultradian: the TF *HLF isoform 1* and the kinases *JNK-3* and *JNK interacting protein 4*. We noted a remodeling of clock elements in *A.m*, with an increase in ultradian oscillations. Indeed, in *A.m*, 9 transcripts were ultradian while only 1 transcript, *HLF isoform 1*, was circadian. In addition to remodeling the rhythmic status in *A.m*, the mean expression levels of 7 of the 19 clock network element transcripts were modified, among which 5 were increased and 2 were decreased.

**Functional implications of cycling transcripts.** Overall, 5061 and 5289 transcripts oscillated (circadian and ultradian), respectively, in *H.t* and *A.m* (24.3% and 25.4% of the transcriptome, respectively (Supplementary Table S2). An enrichment analysis of the biological processes for all oscillating genes in *H.t* showed an enrichment of DNA replication initiation and DNA duplex unwinding gene ontology (GO) terms (Fig. 5A), accounting for 0.3% of total cyclic transcripts in *H.t*. The same analysis conducted with *A.m* oysters showed an enrichment accounting for 2.8% of cyclic transcripts (Fig. 5B). Pseudouridine synthesis, mRNA and ncRNA (non-coding
RNA) processing, nucleotide-excision repair and catabolic processes of proteins and cellular macromolecules were enriched in the *A.m* cycling transcriptome. Finally, Fig. 5C shows that the transcripts oscillating exclusively in *A.m* enrichment in RNA processing account for 1.3%. For further analysis, Supplementary Table S3 provides the gene ID, a gene description, and their cyclic status under the other condition (*H.t* or *A.m*) for the 30 most significant circadian, U1 and U2 transcripts under *H.t* and in *A.m* conditions.

**A. minutum** effect on gene expression level throughout the whole transcriptome.
Finally, we analyzed the impact of *A. minutum* exposure on the level of gene expression in the gills of *C. gigas*. Differential expression analysis was conducted using 13 time samples of each condition as replicates (summarized in Fig. 6A). At a FDR of <0.01, 740 transcripts were differentially expressed (DE); 89.5% of transcripts were upregulated and 10.5% were downregulated in *A.m* (Fig. 6A,B). Among these 740 transcripts, we conducted GO enrichment analyses to assess their functional implications (Fig. 6C). We found that in enriched biological processes, most of the transcripts were upregulated, especially those involved in the oxidative-reduction process, transmembrane transport and ATP catabolic process. Figure 6D shows a pie chart detailing the percentage of rhythmic status for the 740 DE transcripts (left pie chart) and for the whole transcriptome (right pie chart). We showed that 58% of the differentially expressed (DE) transcripts were affected by remodeling of their cyclic status compared to 39% in the whole transcriptome. This result highlighted the overrepresentation of transcripts with remodeling cycling status in differentially expressed transcripts. Specifically, transcripts with a gain of rhythmicity tend to be over-represented in DE transcripts, affecting 35% of the transcripts, compared to 18% in the whole transcriptome. Furthermore, Supplementary Fig. S3 showed the rhythmic status of the DE transcripts associated with all the significantly enriched GO terms. We noted that all the enriched GO terms were associated with remodeling of the cycling transcriptome. Finally, for further analysis, Supplementary Table S4 provides lists of the 25 most significant up- and downregulated transcripts, with their log fold change (FC) value and rhythmic status.

Discussion

Little is known about the strategy of molecular biological rhythms in marine organisms. Our work provides sig-
nificant insights for the first time, to our knowledge, into the temporal regulation of the gill transcriptome of the
oyster *Crassostrea gigas*. Such gills are of importance as one of the primary tissues implicated in the mechanisms
for adapting to environmental changes, such as exposure to a bloom of the PST-producer *A. minutum*. We pro-
vide evidence for the complexity and plasticity of the cycling transcriptome when oysters face a disruption in their
environment. We showed that for conditions *H.t* and *A.m*, at least 42% of the different transcripts were able to
exhibit rhythmicity (Supplementary Table S2). Strikingly, heat maps of the whole transcriptomes (Fig. 7) suggest
wide cyclic gene expression, strongly implying a more extensive temporal organization than specifically indicated
in the present work. The extensive deep sequencing effort performed in this study provide sufficient data to vis-
ualize the global oscillation of the transcriptome. Indeed, read depth sequencing is a key factor that determines
the accuracy of measurements made by high-throughput sequencing, specifically for low-amplitude cyclic tran-
scripts¹⁹. For example, in fruit flies¹⁹, characterized by a transcriptome of about 50 000 transcripts, 18 million
reads per sample were estimated to be sufficient to characterize most of the circadian transcripts. According to
these results, the sequencing of the 20 846 transcripts in the gills performed in this study, with a read depth of
24.5–33 million reads per sample (Supplementary Table S1), allowed to visualize possible cyclic transcripts even if
the statistical analysis was not always significant. To our knowledge, the whole transcriptome organized according
to the phase of the periodicity is not shown in most of published studies, likely due to a lack of sequencing effort
in term of read depth analysis.

The results of the transcriptome analysis revealed circadian oscillations of 6.2% of the gill transcriptome (Fig. S1A–C). This proportion of oscillating transcripts is consistent with circadian transcript levels of 3–10% reported in specific mammalian tissues¹⁶⁻²². In the gills of the mussel *Mytilus californianus*, microarray analysis identified that ≈23% of the transcripts were circadian³. Circadian transcripts in our analysis were characterized by a nocturnal and bimodal pattern (Fig. 1E). Nocturnal peaks of expression for circadian genes in the gills were correlated with the expression of clock genes and the nocturnal behavior of the oysters (Fig. 2). The nocturnal pattern of behavior and core clock genes expression was in accordance to the time of year (winter) when was performed the experiment. Indeed, previous studies demonstrated the dualism of the circadian rhythm in *C. gigas*. Oyster behavior and clock gene expression profiles were nocturnal during autumn and winter¹³, ⁴⁰, and diurnal during spring and summer¹³, ¹⁵. This present study, performed in winter, provided an extensive characterization of nocturnal pattern of gene expression. Although the maximum VOD occurs during the first part of the night (N1), most circadian transcripts peaked during the last part of the night (N2), when the VOD returned to its basal value (Figs 1E and 2D). These findings clearly suggested that the circadian genes in N1 and N2 are not involved in similar physiological processes. In N1, the oysters remained open most of the time, suggesting intense activity related to the surrounding environment, such as gill filtration, feeding, or respiration. By contrast, during N2, the oysters were closed most of the time, and other molecular mechanisms took over. Surprisingly, in addition to circadian transcripts, we also observed two larger groups of ultradian transcripts in *H.t*: 10.5% oscillating under a period range of 8–11 hr (U1) and 7.6% under a period range of 12–16 hr (U2). In the crustacean *Eurydice pulchra*, it has been shown that ultradian rhythms, ≈12.4-hr cycles, were driven by a dedicated circatidal pacemaker that is independent from the circadian system¹⁰. In *M. californianus* exposed to simulated daily and tidal entrainments, ≈2% of transcripts oscillated at a periodicity of 10–14 hr³. However, no endogenous circatidal rhythm has been demonstrated in *C. gigas*¹⁴, and oysters in this experiment were not exposed to ultradian zeitgeber, such as tidal cues. Several studies have identified clock transcripts with an ultradian period of 12 hr in the mouse liver¹⁸, ⁴⁹, ⁵⁶⁻⁵⁸. A system involving a layered circadian transcriptional cascade has been proposed, which could lead to the generation of ultradian genes by an interaction with TFs, which is in turn associated with the induction of clock genes and CCGs⁴⁵. In this study, we identified oyster orthologs of the mammalian and insect circadian subsystems, such
A

Differential expression analysis.

<table><thead><tr><th rowspan="2"></th><th colspan="2">All logFC</th><th colspan="2">logFC > |0.5849|</th><th colspan="2">logFC > |1|</th></tr><tr><th>Down</th><th>Up</th><th>Down</th><th>Up</th><th>Down</th><th>Up</th></tr></thead><tbody><tr><td>FDR < 0.001</td><td>464<br>9 %</td><td>91 %</td><td>463<br>9 %</td><td>91 %</td><td>245<br>10 %</td><td>90 %</td></tr><tr><td>FDR < 0.01</td><td>740<br>10.5 %</td><td>89.5 %</td><td>700<br>11 %</td><td>89 %</td><td>298<br>14 %</td><td>86 %</td></tr><tr><td>FDR < 0.05</td><td>1152<br>17 %</td><td>83 %</td><td>941<br>17 %</td><td>83 %</td><td>380<br>24.5 %</td><td>75.5 %</td></tr></tbody></table>

C

<figure><img src="image_13.png" /></figure>

D

<figure><img src="image_14.png" /><figcaption>**Figure 6.** Effect of *A. minutum* on gene expression level throughout the whole transcriptome. (A) Number of differentially expressed (DE) transcripts in *A.m* for different significances (FDR) and fold changes (log FC). (B) Fold change distribution. The black bars represent equally expressed transcripts, while the orange bars represent DE transcripts. The white dotted line shows the limit between up- and downregulated transcripts. (C) Enrichment analysis of biological processes of DE transcripts with an FDR of <0.01. The ordinate axis shows the GO term and significance of enrichment (FDR) while the abscissa axis shows the percentage of DE transcripts enriched for each GO term. Positive and negative values indicate up- and downregulated transcripts, respectively. (D) The right pie chart shows the cycling transcripts' status for the DE transcripts. The left pie chart shows the cycling transcripts' status for the whole transcriptome in *A.m*.</figcaption></figure>

as TF *HLF isoform 1*, kinase *JNK 3*, and *JNK interacting protein 4*, which oscillated with an ultradian periodicity, suggesting that the oyster circadian network was involved in ultradian rhythmicity (Fig. 4, H.t). The unusual period of ultradian transcripts (U1 and U2) may be explained by the 9 L: 15D regime, mimicking the winter photoperiod observed in the field at the time of the experiment. Indeed, both periods of oscillation, U1 and U2, were
<figure><img src="image_15.png" /><figcaption>**Figure 7.** Remodeling of cycling transcriptome. Heat map expression of the whole transcriptome (20,846 transcripts), organized by significance, phase, and period in *H.t* and *A.m* respectively. Significant circadian, U1, and U2 transcripts are highlighted by a line in each condition. A summary of remodeling processes due to *A. minutum* exposure is presented on the right.</figcaption></figure>

close to the unequal durations of light and night (9 L: 15 D) (Fig. 3A). This regime differed from the classical 12 L: 12 D cycle used in most chronobiological studies.

HABs are currently increasing due to concomitant factors such as global warming, world shipping, and agriculture, and exert a significant effect in terms of ecological and physiological consequences on marine organisms such as filter-feeding bivalves³³. In oysters, much damage has been shown in terms of genotoxicity, lipid metabolism impairment, immunotoxicity, increase of mucus in gills, myopathies, and modifications of filtration and valve behavior¹⁴, ³⁴⁻³⁹, ⁵⁹, ⁶⁰. Here, we discovered an important remodeling of the cycling transcriptome. Compared to *H.t*, *A. minutum* exposure led to (1) a loss of 68% of previous cyclic transcripts; (2) a gain of 69% of *de novo* cyclic transcripts; and (3) a switch of periodicity in 21% of cyclic transcripts, as shown in Supplementary Table S2. At the circadian level, exposure to *A. minutum* led to a ≈50% decrease in the number of cyclic transcripts, with a modification in phase lag and a loss of cyclic transcripts with low amplitudes (Fig. 1, Supplementary Table S2). The modification of phase lag expression was associated with a loss in nocturnal and bimodal pattern, leading to the absence of a clear pattern of daily expression (Fig. 1E–G). Similar trends were observed in the expression of clock genes and valve opening behavior (Fig. 2). As explain earlier, diurnal and nocturnal behavioral patterns alternate in *C. gigas* throughout the year in a predictable way¹³. This dualism could be an adaptive response to variations in energy needs during the year¹³. Consequently, the disruption of the nocturnal circadian pattern observed in this study could lead to a loss of fitness of the oyster within its biotope. Note, in *A.m*, we observed an increase in *Target of EGR1* expression (Fig. 4). As *EGR1* is involved in the amplitude of circadian transcripts in the mouse⁵¹, *Target of EGR1* could be implicated in the increase of circadian transcripts with high amplitudes in *A.m* (Fig. 11). Remodeling of circadian gene profiles (loss and gain of oscillation, as well as modification of phase lag expression) has been previously demonstrated in the mouse. Evidence for a modulatory role of the transcriptional factor *KLF15* in the diurnal pattern of the circadian genes in the mouse heart has been shown⁶¹. Indeed, in the heart of KLF15 knockout mice, authors observed 805 circadian transcripts compared to 1335 in control heart, with only 332 common circadian transcripts. Authors demonstrated that KLF15 was involved in the generation of transcripts oscillations and biphasic patterns as well as in the inhibition of oscillations of others transcripts and the quiescent expression of other genes associated to the circadian subsystem⁶¹. Similar weak number of commonly detected circadian transcripts was observed in liver of mice fed with high fat diet compared to normal diet¹⁷. This reprogramming involved both an impairment of CLOCK-BMAL chromatin recruitment and a cyclic activation of surrogate pathway through the transcriptional regulator PPARγ¹⁷.

At the ultradian level, we observed a 24% increase in ultradian transcripts in *A.m* compared to *H.t* (Supplementary Table S2). An analysis of oyster orthologs of TFs and kinases involved in the circadian network of mammals and insects also showed an increase in ultradian oscillations in *A.m* (Fig. 4). Interestingly, the circadian TFs *HLF* and *DBP2* are known to control the expression of enzymes and regulators involved in detoxification and
drug metabolism in mice⁶²,⁶³. In our study, *HLF isoform 2*, *HLF isoform 3*, and *DBP2*, which were arrhythmic in *H.t*, exhibited an ultradian oscillation in *A.m*. Thus, the increase of ultradian transcripts in *A.m* could be linked to the implementation of detoxification mechanisms in response to PSTs bioaccumulation, through the TF *HLFs* and *DBP2*. The TF ortholog of *AP-1* was observed to switch from circadian to U1 oscillation in *A.m*. This result suggested an acceleration of cell proliferation, apoptosis, and abnormal cell cycles, i.e. cellular processes in which *AP-1* is a regulator⁶⁴. The kinase *GSK-3β* and casein kinase 2α, which are also arrhythmic in *H.t*, gained ultradian rhythmicity. Both kinases are known to act on TIM or BMAL1 phosphorylation and could induce a shortened period of cyclic processes, offering interesting leads to explain the decrease of U2 period length observed in *A.m*⁴⁷, ⁶⁵, ⁶⁶.

In addition to the remodeling of the cycling transcriptome, *A. minutum* exposure led to a change in the expression level of 740 transcripts (FDR of <0.01), with 89.5% exhibiting overexpression (Fig. 6), suggesting physiological reactions against this harmful alga. Fig. 6D revealed that 58% of differentially expressed transcripts appeared to be affected by cyclic remodeling. Moreover, all enriched biological processes (GO terms) were impacted by a remodeling of the cyclic transcriptome (Supplementary Fig. S3), suggesting a pleiotropic impact of *A. minutum*.

Mechanisms associated to the remodeling of cyclic transcriptome of *C. gigas* during *A. minutum* exposure remained unknown. *Alexandrium minutum* is known to produce paralytic shellfish toxins (such as STX) but also bioactive extracellular compounds that can exhibit toxic activities by direct contact with gills⁶⁷. An effect of *A. minutum* directly on core clock genes or on the clock subsystem genes, through genotoxic impacts could be suggested, as previously observed¹⁴. Moreover, previous studies on chronotoxic effects in mammals and mollusk gastropods such as *Aplysia* have shown that the potential action of clock neurons was blocked by TTX, a homologue of STX toxin. The master clock kept time, but did not send the circadian message to other regions⁶⁸, ⁶⁹. Another study reported that TTX damped and desynchronized the suprachiasmatic nucleus cellular rhythms in *Per1* expression⁷⁰.

To conclude, this original study of the temporal transcriptome of *C. gigas* revealed highly cyclic gene expression and highlighted the complexity of clock mechanisms, with an intriguing predominance of ultradian genes in a circadian context. Moreover, the results showed pleiotropic impacts of *A. minutum* exposure on the oyster gill transcriptome. The transcriptional remodeling associated with exposure to harmful algae revealed plasticity in the temporal organization of gene expression. To date, the consequences of this remodeling plasticity remain unknown. This response could be an adaptive reaction to the toxicity of the harmful algae. It may also reveal a consequent chronotoxicity effect leading to a disruption in the temporal structure and an uncoupling of the transcriptional rhythmicity associated with the physiological needs. In any case, this effect led to a loss of fitness of the oyster within its environment.

**Experimental procedures.** *Animals and experimental treatment.* Diploid Pacific *C. gigas* oysters (1.5 years old, n = 156) from the Bay of Arcachon (France) (72.8 ± 0.5 mm shell length; 40.6 ± 0.4 mm shell width; 38.5 ± 0.8 g total fresh weight; mean ± SE) were divided among 6 tanks isolated from external vibration using an anti-vibrating bench and continuously supplied with homogenized seawater of constant composition (T °C = 15.0 ± 0.2 °C; pH = 8.0 ± 0.1; salinity = 35.8 ± 0.01‰). Experiments were conducted in December to avoid gametogenesis. Oysters were acclimated for 7 days to 9 hr light-15 hr dark conditions mimicking the natural light regime in winter time fed with the non-harmful alga *Heterocapsa triquetra* of pre-exposure (continuously supplied with algae at a flow rate of 100 mL.min⁻¹ at a concentration of 1000 cells.mL⁻¹; 208.3 ± 44.5 cells.mL⁻¹ measured in *H.t* tanks; 208.7 ± 38.1 cells.mL⁻¹ measured in the *A.m* tanks, mean ± SE). Then, the *H.t* treatment group (3 replicate tanks) remained exposed to *H. triquetra* whereas the *A.m* treatment group (3 replicate tanks) was exposed to the harmful alga *Alexandrium minutum* for 6 days (continuously supplied with algae at a flow rate of 100 mL.min⁻¹ at a concentration of 1000 cells.mL⁻¹; 159.1 ± 15.7 cells.mL⁻¹ measured in *H.t* tanks; 430.9 ± 160.8 cells.mL⁻¹ measured in the *A.m* tanks, mean ± SE). During the two last days of exposure, gill tissue (n = 6 oysters per condition) was collected every 4 hr for 52 hr, under natural light during light phases and under dim red light during dark phases, for a total of 13 sampling times (day 1: Zeitgeber time ZT 4.5, ZT 8.5, ZT 12.5, ZT 16.5, ZT 20.5; day 2: ZT 0.5, ZT 4.5, ZT 8.5, ZT 12.5, ZT 16.5, ZT 20.5; day 3: ZT 0.5; ZT 4.5). Zeitgeber time 0 (ZT0) corresponds to the light switch-on.

**Algae.** The paralytic shellfish toxin-producing dinoflagellate *A. minutum* (strain AM89BM, Halim) and the non-toxic dinoflagellate *H. triquetra* (strain HT99PZ-Ehrenberg, 1840) were grown in 10 L and 80 L phytoreactors using f/2 medium⁷¹. During the exponential growth phase, *A. minutum* produced 1.3 ± 0.1 pg saxitoxin per cell⁷². *Heterocapsa triquetra* was chosen as control for several reasons. As *A. minutum*, *H. triquetra* belonged to the dinoflagellates and both species exhibited similar characteristics of shape and size with 19–28 µm and 23–29 µm cell sizes for *H. triquetra* and *A. minutum*, respectively. Algae concentrations were assessed by a Beckman Coulter Z2 (Beckman Coulter Inc., USA).

**Paralytic Shellfish Toxins (PST).** Measurement of PST content in gill tissues of animals used for RNA-seq analysis was done during the 52 hr of sampling (n = 72). PST (saxitoxin and its derivatives) were assessed by ELISA (Abraxis, Novakits, France) according to manufacturer's instructions. Briefly, individual samples were crushed on ice, centrifuged, and diluted. Standard solutions and samples were added in duplicate (50 µl) into wells of an ELISA microtiter plate coated with a secondary sheep anti-rabbit antibody. Then, the saxitoxin horseradish peroxidase (saxitoxin-HRP) conjugate solution and the primary antibody solution (rabbit anti-saxitoxin) were added successively. Plates were incubated at room temperature and washed four times using washing solution. Substrate solution containing tetramethyl-benzidine was added to each well and plates were incubated. Finally, reactions were stopped by a solution containing sulfuric acid and absorbance was read at 450 nm. Color intensity was inversely proportional to the concentration of PSTs in the samples.
**RNA sequencing (RNAseq).** Total RNA from gills was extracted from samples using TRI® Reagent (Invitrogen). After purity (Nanodrop) and quality (Agilent 210 BioAnalyzer, RIN value > 8.5) checks, RNA samples were normalized and pooled (n = 6) in equal quantities of total RNA, according to time and condition. Messenger RNA purifications were performed to generate 26 libraries from 4 µg of total RNA using the Illumina TruSeq Strand mRNA Sample Prep kit according to the manufacturer's protocol. Sequencing of 100 bp paired-end fragments was performed on four pooled libraries, each one on two lanes, on the same flowcell, using an Illumina Hiseq 2500 at the genomic platform Génopole Toulouse/Midi-Pyrénées (INRA Auzeville, France).

**RNAseq data analysis.** Adapter and base quality trimming (Phred score < 20) was performed with Trim Galore (trim_galore_v0.3.7). The read depth was fixed at 47.0 million reads/sample (down-sampling normalization, Supplementary Table S1) for cycling transcriptome analysis, according to¹⁹. Relative log expression normalization (RLE)⁷³ was applied instead of down-sampling for differential expression analysis (edgeR version 3.8.6⁷⁴). Sequenced reads were aligned by STAR (STAR_2.4.2a⁷⁵) to both the genome and mitochondrial genome of *Crassostrea gigas*⁷⁶, Ensemble 1.30; and NCBI reference NC_001276) and assigned by FeatureCount (1.4.5 version⁷⁷) to annotations of the genome in public databases. Expression levels of genes were presented as cpm (counts per million). Filtering of transcripts with 0 cpm in the whole set of libraries was applied. Indeed, considering the homogenized read depth of samples (Supplementary Table S2), 0 cpm value for a transcript that have at least 1 cpm in at least one library was considered as a low expression and not a missing data. Finally, a total of 20,846 and 20,796 transcripts were obtained for analysis of cycling transcripts and differentially expressed transcripts, respectively.

**Analysis of transcripts.** Sequencing outputs and bioinformatics treatments were carefully analyzed and checked several times to ensure the reliability of data supporting the discussion. Cycling transcripts were detected using ARSER⁴¹ with MetaCycle (version 1.1.0⁷⁸ on R⁷⁹ (32-bit, version 3.2.2)). ARSER algorithm was developed and applied to high-throughput time-series analysis and particularly adapted to experimental designs with a timeframe of 4 h over 52 hr⁴¹, ⁷⁸, ⁸⁰. Post-hocs were corrected for multiple comparisons to remove false discoveries⁴², ⁴³. For all cycling detection, a false discovery rate (FDR) of < 0.05 was considered significant. Note that in 24 hr light/dark entrainment, the external circadian cues may directly generate 24 hr rhythms, in an addition of CCGs group. Amplitudes of oscillating transcripts were calculated as ((high value - low value)/low value)⁵⁷. All these genes with a period length between 20 hr and 28 hr comprise circadian transcripts⁴⁵. Differentially expressed transcripts were analyzed with edgeR (edgeR version 3.8.6), considering the 13 sample times (day 1 ZT 0.5 to day 3 ZT 0.5) of each algal treatment as replicates; the results were corrected for multiple tests using the FDR⁸¹.

**Gene ontology analysis.** Gene descriptions were obtained using the *C. gigas* reference genome⁷⁶ on Biomart, Ensembl Metazoa (http://metazoa.ensembl.org/biomart/martview/). Enrichment analyses were performed using the Fisher test on Blast2GO software (version 2.8), selecting only specific biological processes.

**cDNA synthesis and real-time quantitative PCR (RT-qPCR).** The same samples used for RNA sequencing were analyzed by performing RT-qPCR at 7 sampling times (day 2: ZT 0.5, ZT 4.5, ZT 8.5, ZT 12.5, ZT 16.5, ZT 20.5 and day 3: ZT 0.5). Total RNA from individual samples of gills was extracted as previously described and reverse transcribed using oligo dT17 and Moloney murine leukemia virus reverse transcriptase (M-MLV; Promega). RT-qPCR reactions were performed in replicate (n = 6) on individual samples using Brillant III Ultra-Fast SYBR® Green QPCR Master Mix kit (Agilent)¹⁵. Primer sets of circadian clock genes were designed¹⁵ and are listed in Supplementary Table S5. EF1 (Elongation factor 1) was selected as the housekeeping gene using geNorm, BestKeeper and Normfinder (Wang *et al.* 2012) and relative expression levels were calculated using the comparative Ct method (2⁻ΔΔCt method⁸²). Validation of the RNAseq analysis was done by RT-qPCR assessment of clock genes. Normalization of gene expression was done for both analyses as follow: (Xᵢ - min) / (max - min), where Xᵢ was the gene expression value, min and max were the minimum and maximum values of the gene expression in a given condition and for a specific clock gene (Supplementary Fig. S2).

**Behavioral monitoring.** The valve activity of *C. gigas* (n = 32) was studied using two HFNI (high-frequency non-invasive) valvometers⁴⁰, ⁸³. Lightweight electromagnets (0.1 g) were glued on each valve of each animal. These electrodes were connected to the valvometer by flexible wires which allowed the oysters to move their valves without constraint. The sampling frequency for each individual was 0.2 Hz. Data were processed using LabView 8.0 software (National Instruments). For the VOD, the mean hourly valve opening duration was expressed as the percentage of time that oysters spent with their valves open: 100% indicates that the valves were open for the entire hour and 0% indicates that the valves were closed for the hour. The 48-hr mean hourly opening durations for 16 oysters per condition was used in the analyses.

**Statistical analysis for RT-qPCR and behavioral analysis.** Results are expressed as the means ± the SE. Treatment differences were determined using analysis of variance (ANOVA) after checking assumptions (normality and homoscedasticity of the error term). When assumptions were not met, the non-parametric Kruskall-Wallis test was performed. If the null hypothesis was rejected, the Student-Newman-Keuls method was applied to determine significant differences between conditions. For all statistical results, a probability of p < 0.05 was considered significant. Statistical analyses were performed using Sigma Stat software (Version 13.0, SYSTAT, Chicago, USA).

**Data access.** RNAseq data are available in the ArrayExpress database (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-5345.
## References

1. Yerushalmi, S. & Green, R. M. Evidence for the adaptive significance of circadian rhythms. *Ecol. Lett.* **12**, 970–981 (2009).

2. Zhang, E. E. & Kay, S. A. Clocks not winding down: unravelling circadian networks. *Nat. Rev. Mol. Cell Biol.* **11**, 764–776 (2010).

3. Connor, K. M. & Gracey, A. Y. Circadian cycles are the dominant transcriptional rhythm in the intertidal mussel *Mytilus californianus*. *Proc. Natl. Acad. Sci.* **108**, 16110–16115 (2011).

4. Hur, S.-P. et al. Fish sleeping under sandy bottom: interplay of melatonin and clock genes. *Gen. Comp. Endocrinol.* **177**, 37–45 (2012).

5. Oren, M. et al. Profiling molecular and behavioral circadian rhythms in the non-symbiotic sea anemone *Nematostella vectensis*. *Sci. Rep.* **5**, 11418 (2015).

6. Teschke, M., Wendt, S., Kawaguchi, S., Kramer, A. & Meyer, B. A circadian clock in antarctic krill: an endogenous timing system governs metabolic output rhythms in the euphausid species *Euphausia superba*. *PLoS ONE* **6**, e26090 (2011).

7. Tessmar-Raible, K., Raible, F. & Arboleda, E. Another place, another timer: Marine species and the rhythms of life. *BioEssays* **33**, 165–172 (2011).

8. Wilcockson, D. & Zhang, L. Circatalid clocks. *Curr. Biol.* **18**, R753–R755 (2008).

9. Zantke, J. et al. Circadian and circalunar clock interactions in a marine annelid. *Cell Rep.* **5**, 99–113 (2013).

10. Zhang, L. et al. Dissociation of circadian and circatalid timekeeping in the marine crustacean *Eurydice pulchra*. *Curr. Biol.* **23**, 1863–1873 (2013).

11. Bernard, I. et al. *In situ* spawning in a marine broadcast spawner, the Pacific oyster *Crassostrea gigas*: timing and environmental triggers. *Limnol. Oceanogr.* **61**, 635–647 (2016).

12. Tran, D. et al. Field chronobiology of a molluscan bivalve: how the moon and sun cycles interact to drive oyster activity rhythms. *Chronobiol. Int.* **28**, 307–317 (2011).

13. Mat, A. M., Massabuau, J.-C., Ciret, P. & Tran, D. Evidence for a plastic dual circadian rhythm in the oyster *Crassostrea gigas*. *Chronobiol. Int.* **29**, 857–867 (2012).

14. Mat, A. M., Haberkorn, H., Bourdineaud, J.-P., Massabuau, J.-C. & Tran, D. Genetic and genotoxic impacts in the oyster *Crassostrea gigas* exposed to the harmful alga *Alexandrium minutum*. *Aquat. Toxicol.* **140–141**, 458–465 (2013).

15. Perrigault, M. & Tran, D. Identification of the molecular clockwork of the oyster *Crassostrea gigas*. *PLoS ONE* In press (2017).

16. Akhtar, R. A. et al. Circadian cycling of the mouse liver transcriptome, as revealed by cDNA microarray, is driven by the suprachiasmatic nucleus. *Curr. Biol.* **12**, 540–550 (2002).

17. Eckel-Mahan, K. L. et al. Reprogramming of the circadian clock by nutritional challenge. *Cell* **155**, 1464–1478 (2013).

18. Hughes, M. E. et al. Harmonics of circadian gene transcription in mammals. *PLOS Genet* **5**, e1000442 (2009).

19. Li, J., Grant, G. R., Hogenesch, J. B. & Hughes, M. E. In *Methods in Enzymology* (ed. Sehgal, A.) **551**, 349–367 (Academic Press, 2015).

20. Panda, S. et al. Coordinated transcription of key pathways in the mouse by the circadian clock. *Cell* **109**, 307–320 (2002).

21. Ueda, H. R. et al. System-level identification of transcriptional circuits underlying mammalian circadian clocks. *Nat. Genet.* **37**, 187–192 (2005).

22. Zhang, R., Lahens, N. F., Ballance, H. I., Hughes, M. E. & Hogenesch, J. B. A circadian gene expression atlas in mammals: Implications for biology and medicine. *Proc. Natl. Acad. Sci.* **111**, 16219–16224 (2014).

23. Antunes, L. C., Levandovski, R., Dantas, G., Caumo, W. & Hidalgo, M. P. Obesity and shift work: chronobiological aspects. *Nutr. Res. Rev.* **23**, 155–168 (2010).

24. Doi, M. et al. Salt-sensitive hypertension in circadian clock-deficient *Cry*-null mice involves dysregulated adrenal *Hsd3b6*. *Nat. Med.* **16**, 67–74 (2010).

25. Fonken, L. K. et al. Light at night increases body mass by shifting the time of food intake. *Proc. Natl. Acad. Sci.* **107**, 18664–18669 (2010).

26. Knutsson, A. Health disorders of shift workers. *Occup. Med.* **53**, 103–108 (2003).

27. Sharifian, A., Farahani, S., Pasalar, P., Gharavi, M. & Aminian, O. Shift work as an oxidative stressor. *J. Circadian Rhythms* **3**, 15 (2005).

28. Suwazono, Y. et al. A longitudinal study on the effect of shift work on weight gain in male japanese workers. *Obesity* **16**, 1887–1893 (2008).

29. Viitasalo, K., Puttonen, S., Kuosma, E., Lindström, J. & Härmä, M. Shift rotation and age - interactions with sleep-wakefulness and inflammation. *Ergonomics* **58**, 65–74 (2015).

30. Anderson, D. M., Cembella, A. D. & Hallegraeff, G. M. Progress in understanding harmful algal blooms: paradigm shifts and new technologies for research, monitoring, and management. *Annu. Rev. Mar. Sci.* **4**, 143–176 (2012).

31. Bricelj, V. M. & Shumway, S. E. Paralytic shellfish toxins in bivalve molluscs: occurrence, transfer kinetics, and biotransformation. *Rev. Fish. Sci.* **6**, 315–383 (1998).

32. Hallegraeff, G. M., Anderson, D. M., Cembella, A. D. & Enevoldsen, H. O. *Manual on harmful marine microalgae*. (Unesco, 2003).

33. Food and Agriculture Organization of the United Nations (FAO) & World Health Organization (WHO). Technical paper on Toxicity Equivalency Factors for Marine Biotoxins Associated with Bivalve Molluscs. (2016).

34. Galimany, E., Sunila, I., Hégaret, H., Ramón, M. & Wikfors, G. H. Pathology and immune response of the blue mussel (*Mytilus edulis* L.) after an exposure to the harmful dinoflagellate *Prorocentrum minimum*. *Harmful Algae* **7**, 630–638 (2008).

35. Haberkorn, H. et al. Effects of *Alexandrium minutum* exposure on nutrition-related processes and reproductive output in oysters *Crassostrea gigas*. *Harmful Algae* **9**, 427–439 (2010).

36. Haberkorn, H. et al. Effects of *Alexandrium minutum* exposure upon physiological and hematological variables of diploid and triploid oysters, *Crassostrea gigas*. *Aquat. Toxicol.* **97**, 96–108 (2010).

37. Hégaret, H., Wikfors, G. H. & Shumway, S. E. Diverse feeding responses of five species of bivalve mollusc when exposed to three species of harmful algae. *J. Shellfish Res.* **26**, 549–559 (2007).

38. Mello, D. F., Silva, P. M., da Barracco, M. A., Soudant, P. & Hégaret, H. Effects of the dinoflagellate *Alexandrium minutum* and its toxin (saxitoxin) on the functional activity and gene expression of *Crassostrea gigas* hemocytes. *Harmful Algae* **26**, 45–51 (2013).

39. Tran, D., Haberkorn, H., Soudant, P., Ciret, P. & Massabuau, J.-C. Behavioral responses of *Crassostrea gigas* exposed to the harmful algae *Alexandrium minutum*. *Aquaculture* **298**, 338–345 (2010).

40. Tran, D. et al. The toxic dinoflagellate *Alexandrium minutum* disrupts daily rhythmic activities at gene transcription, physiological and behavioral levels in the oyster *Crassostrea gigas*. *Aquat. Toxicol.* **158**, 41–49 (2015).

41. Yang, R. & Su, Z. Analyzing circadian expression data by harmonic regression based on autoregressive spectral estimation. *Bioinformatics* **26**, i168–i174 (2010).

42. Storey, J. D. & Tibshirani, R. Statistical significance for genomewide studies. *Proc. Natl. Acad. Sci.* **100**, 9440–9445 (2003).

43. Storey, J. D., Xiao, W., Leek, J. T., Tompkins, R. G. & Davis, R. W. Significance analysis of time course microarray experiments. *Proc. Natl. Acad. Sci.* **102**, 12837–12842 (2005).

44. Korenčič, A. et al. Timing of circadian genes in mammalian tissues. *Sci. Rep.* **4** (2014).

45. Westermark, P. O. & Herzel, H. Mechanism for 12 hr rhythm generation by the circadian clock. *Cell Rep.* **3**, 1228–1238 (2013).

46. Duffield, G. E. et al. A role for *Id2* in regulating photic entrainment of the mammalian circadian system. *Curr. Biol.* **19**, 297–304 (2009).

47. Lin, J.-M. et al. A role for *casein kinase 2α* in the Drosophila circadian clock. *Nature* **420**, 816–820 (2002).
48. Martinek, S., Inonog, S., Manoukian, A. S. & Young, M. W. A role for the segment polarity gene shaggy / GSK-3 in the drosophila circadian clock. *Cell* **105**, 769–779 (2001).

49. Ramsey, K. M. et al. Circadian clock feedback cycle through NAMPT-mediated NAD+ biosynthesis. *Science* **324**, 651–654 (2009).

50. Takeuchi, J., Shannon, W., Aronin, N. & Schwartz, W. J. Compositional changes of AP-1 DNA-binding proteins are regulated by light in a mammalian circadian clock. *Neuron* **11**, 825–836 (1993).

51. Tao, W. et al. EGR1 regulates hepatic clock gene amplitude by activating *Per1* transcription. *Sci. Rep.* **5**, 15212 (2015).

52. Vielhaber, E., Eide, E., Rivers, A., Gao, Z.-H. & Virshup, D. M. Nuclear entry of the circadian regulator mPER1 is controlled by mammalian casein kinase I varepsilon1on. *Mol. Cell. Biol.* **20**, 4888–4899 (2000).

53. Xi, Y. & Chen, D. Partitioning the circadian clock. *Science* **345**, 1122–1123 (2014).

54. Yoshitane, H. et al. JNK regulates the photic response of the mammalian circadian clock. *EMBO Rep.* **13**, 455–461 (2012).

55. Zhao, H. et al. PARP1- and CTCF-mediated interactions between active and repressed chromatin at the lamina promote oscillating transcription. *Mol. Cell* **59**, 984–997 (2015).

56. Cretenet, G., Le Clech, M. & Gachon, F. Circadian clock-coordinated 12 hr period rhythmic activation of the IRE1α pathway controls lipid metabolism in mouse liver. *Cell Metab.* **11**, 47–57 (2010).

57. Hughes, M. E., Grant, G. R., Paquin, C., Qian, J. & Nitabach, M. N. Deep sequencing the circadian and diurnal transcriptome of drosophila brain. *Genome Res.* **22**, 1266–1281 (2012).

58. Vollmers, C. et al. Time of feeding and the intrinsic circadian clock drive rhythms in hepatic gene expression. *Proc. Natl. Acad. Sci.* **106**, 21453–21458 (2009).

59. Wildish, D., Lassus, P., Martin, J., Saulnier, A. & Bardouil, M. Effect of the PSP-causing dinoflagellate, *Alexandrium* sp. on the initial feeding response of *Crassostrea gigas*. *Aquat. Living Resour.* **11**, 35–43 (1998).

60. Bardouil, M., Bohec, M., Bougrier, S., Lassus, P. & Truquet, P. Feeding responses of *Crassostrea gigas* (Thunberg) to inclusion of different proportions of toxic dinoflagellates in their diet. *Oceanol. Acta* **19**, 177–182 (1996).

61. Zhang, L. et al. KLF15 establishes the landscape of diurnal expression in the heart. *Cell Rep.* **13**, 2368–2375 (2015).

62. Claudel, T., Cretenet, G., Saumet, A. & Gachon, F. Crosstalk between xenobiotics metabolism and circadian clock. *FEBS Lett.* **581**, 3626–3633 (2007).

63. Gachon, F., Olela, F. F., Schaad, O., Descombes, P. & Schibler, U. The circadian PAR-domain basic leucine zipper transcription factors DBP, TEF, and HLF modulate basal and inducible xenobiotic detoxification. *Cell Metab.* **4**, 25–36 (2006).

64. Fu, L. & Lee, C. C. The circadian clock: pacemaker and tumour suppressor. *Nat. Rev. Cancer* **3**, 350–361 (2003).

65. Itaka, C., Miyazaki, K., Akaike, T. & Ishida, N. A role for glycogen synthase kinase-3β in the mammalian circadian clock. *J. Biol. Chem.* **280**, 29397–29402 (2005).

66. Sahar, S., Zocchi, L., Kinoshita, C., Borrelli, E. & Sassone-Corsi, P. Regulation of BMAL1 protein stability and circadian function by GSK3β-mediated phosphorylation. *PLoS ONE* **5**, e8561 (2010).

67. Borcier, E. et al. Effects of bioactive extracellular compounds and paralytic shellfish toxins produced by *Alexandrium minutum* on growth and behaviour of juvenile great scallops *Pecten maximus*. *Aquat. Toxicol.* **184**, 142–154 (2017).

68. Eskin, A. Neurophysiological mechanisms involved in photo-entrainment of the circadian rhythm from the *Aplysia* eye. *J. Neurobiol.* **8**, 273–299 (1977).

69. Schwartz, W. J., Gross, R. A. & Morton, M. T. The suprachiasmatic nuclei contain a tetrodotoxin-resistant circadian pacemaker. *Proc. Natl. Acad. Sci.* **84**, 1694–1698 (1987).

70. Yamaguchi, S. et al. Synchronization of cellular clocks in the suprachiasmatic nucleus. *Science* **302**, 1408–1412 (2003).

71. Guillard, R. R. L. In *Culture of Marine Invertebrate Animals* (eds. Smith, W. L. & Chanley, M. H.) 29–60 (Springer US, 1975).

72. Ravn, H., Anthoni, U., Christophersen, C., Nielsen, P. H. & Oshima, Y. Standardized extraction method for paralytic shellfish toxins in phytoplankton. *J. Appl. Phycol.* **7**, 589–594 (1995).

73. Anders, S. & Huber, W. Differential expression analysis for sequence count data. *Genome Biol.* **11**, R106 (2010).

74. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. EdgeR: a Bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics* **26**, 139–140 (2010).

75. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. *Bioinformatics* **29**, 15–21 (2013).

76. Zhang, G. et al. The oyster genome reveals stress adaptation and complexity of shell formation. *Nature* **490**, 49–54 (2012).

77. Liao, Y., Smyth, G. K. & Shi, W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. *Bioinformatics* **30**, 923–930 (2014).

78. Wu, G., Anafi, R. C., Hughes, M. E., Kornacker, K. & Hogenesch, J. B. *MetaCycle: an integrated R package to evaluate periodicity in large scale data*. (2016).

79. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/ (2013).

80. Wu, G. et al. Evaluation of five methods for genome-wide circadian gene identification. *J. Biol. Rhythms* **29**, 231–242 (2014).

81. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. *J. R. Stat. Soc. Ser. B Methodol* **57**, 289–300 (1995).

82. Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2^−ΔΔCT^ method. *Methods* **25**, 402–408 (2001).

83. Tran, D., Ciret, P., Ciutat, A., Durrieu, G. & Massabauaj, J.-C. Estimation of potential and limits of bivalve closure response to detect contaminants: Application to cadmium. *Environ. Toxicol. Chem.* **22**, 914–920 (2003).

## Acknowledgements

This work was supported by the financial support of the Nation Research Agency (ANR), ACCUTOX project 13-CESA-0019 (2013–2017) and a ministerial PhD scholarship for L. Payton. We thank O. Bouchez and A. Chaubet (Génopole Toulouse/Midi-Pyrénées, INRA Auzeville, France) for sequencing and discussions. We would thank I. Gonzalez, G. Lefort (Genotoul, INRA Auzeville, France) and Anthony Bertucci (UMR EPOC, Arcachon, France) for their advices in bioinformatics and biostatistics treatment. Authors also thank Dr. K.Flynn for English corrections.

## Author Contributions

Conceptualization, L.P., D.T.; Investigation, L.P., D.T., M.P., A.H., F.B., C.F., H.H., J.C.M.; Samples treatment, L.P., M.P.; Bioinformatics treatment, C.H.; Biostatistics treatment, L.P., M.S., C.H.; Funding acquisition, D.T., J.C.M.; Visualization, L.P., D.T., Writing—Original Draft, L.P.; Writing correction, D.T., M.P., J.C.M., M.S., A.H., C.H., H.H., F.B.; Supervision, D.T.

## Additional Information

Supplementary information accompanies this paper at doi:10.1038/s41598-017-03797-4

Competing Interests: The authors declare that they have no competing interests.
**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_16.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre-
ative Commons license, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons license and your intended use is not per-
mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the
copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017