Seagrass specimens used for this study are the same as in Ruocco et al. (2020a). In May 2019, entire shoots of Cymodocea nodosa and Zostera marina were collected by snorkeling from a shallow-water (1-2 m depth) mixed meadow near Ilha da Culatra (Faro, Portugal – 36°59′41″ N, 7°50′26″ W) with a significant portion of the sediment. Plants and sediment were then transported to the nearby Ramalhete field station (CCMAR, Centre of Marine Sciences – 37°00′22″ N, 7°58′02″ W) in darkened containers filled with seawater, limiting rhizome and roots breakage. Following removal from the field, shoots were transplanted to six 150-L cylindrical tanks located in an outdoor mesocosm facility exposed to a natural daily light regime. The bottom of each tank was covered with ca. 7 cm washed beach sand. Tanks were supplied with running seawater from the nearby area, previously filtered with sand and UV filters, in an independent open circuit configuration. Seagrass shoots (ca. 20) were randomly assigned to each tank. Plants were left to recover for 24 h following transplantation. At the end of the acclimation period, leaf sections of C. nodosa and Z. marina were collected 6 times during a 24-h cycle, i.e., at midnight (00:00), dawn (04:30), sunrise (06:30), solar noon (13:00), sunset (20:30), and dusk (22:00). All information about sunrise and sunset times were taken from https://www.timeanddate.com/. Irradiance (as Photosynthetically Active Radiation, PAR) and water temperature were monitored continuously in the tanks with a Li-192SA underwater quantum sensor connected to a Li-1400 data logger (Li-Cor), and a HOBO temperature logger, respectively. Daily patterns of PAR and temperature levels are displayed in Supplementary Figure 1. For all analyses, three independent biological replicates were collected for each species and time point (n = 3) and used for both gene expression and carbohydrate analyses. Leaf samples were rapidly cleaned of epiphytes, rinsed with distilled water and immediately submerged in RNAlater^© tissue collection (Ambion, Life Technologies, Waltham, MA, United States). After one night at 4°C, leaf samples were definitely stored at –20°C.

Total RNA from the youngest fully developed leaves of C. nodosa and Z. marina (Ruocco et al., 2019) was extracted with the Aurum™ Total RNA Mini Kit (BIO-RAD, Hercules, CA, United States), following the manufacturer’s protocol. About 5–7 cm-long leaf sections corresponding to about 70–100 mg were ground to a fine powder with a mortar and pestle containing liquid nitrogen. Samples were then homogenized through a Mixer Mill MM300 (QIAGEN, Hilden, Germany) and tungsten carbide beads (3 mm) for 3 min at 20.1 Hz. The quality and purity of the total RNA was checked using NanoDrop (ND-1000 UV-Vis spectrophotometer; NanoDrop Technologies, Wilmington, DE, United States) and 1% agarose gel electrophoresis. RNA was used when Abs260 nm/Abs280 nm and Abs260 nm/Abs230 nm ratios were > 1.8 and 1.8 < x < 2, respectively. The RNA concentration was accurately determined by the Qubit™ RNA BR Assay kit (Thermo Fisher Scientific, Waltham, MA, United States) using the Qubit 2.0 Fluorometer (Thermo Fisher Scientific). The total RNA (500 ng) from each sample was reverse-transcribed into cDNA with the iScript™ cDNA synthesis kit (BIO-RAD), according to the manufacturer’s protocol.

Primers for target genes (19) involved in key plant metabolic pathways potentially exhibiting circadian changes i.e., sugar metabolism, photosynthesis and light harvesting, photoreception and circadian clock, growth (Harmer et al., 2000), were newly designed with the primer analysis software Primer3 v. 0.4.0 (Untergasser et al., 2012) or selected from previous studies (see Table 1). The design conditions included the primer length (18–23 bp), Tm (∼60 C), GC content (50%), and product size (100–250 bp). To identify target sequences for Genes of Interests (GOIs), the corresponding protein sequences of A. thaliana were downloaded from TAIR¹ and used as queries for BLASTP (Altschul et al., 1990). For Z. marina, sequences were blasted against the genome available in ORCAE² (Olsen et al., 2016). For C. nodosa, potential homologous sequences were identified through sequence similarity search against an available transcriptome (Ruocco et al., 2017) or public databases. The BLASTP analysis was carried out setting a threshold of homology to 1e-10. The top hit for each query was selected and the homology level of each pair of target sequences of Z. marina vs. C. nodosa, was reported (Supplementary Tables 1a,b). To verify the evolutionary relationships among target genes identified in seagrasses in comparison to land plants, homologous sequences of GOIs of other three plant species (Oryza sativa, Zea mays, and Populus trichocarpa), representatives of monocots and eudicots lineages, were downloaded from public databases (i.e., KEGG, Phytozome v13, and UniProt) and included in our analysis (Supplementary Table 2). Multiple sequence alignments were performed using COBALT (Papadopoulos and Agarwala, 2007) with default parameters. A further phylogenetic analysis was carried out on GOIs for which potential orthologues sequences were identified between Z. marina and C. nodosa using the Phylemon2 toolkit (Sánchez et al., 2011; Supplementary Table 3). Poorly aligned sequence regions were removed with trimAl v1.3 (Capella-Gutiérrez et al., 2009) using the strictplus option. For each dataset, a maximum likelihood phylogenetic tree was inferred using PhyML (Guindon et al., 2010) implemented in the online tool Phyml Best AIC Tree v.1.02b. Topology optimization of best model was calculated for each independent dataset and the support of nodes was calculated using the SH-like procedure (Anisimova and Gascuel, 2006). Generated trees were visualized in UGENE v.33.0 (Okonechnikov et al., 2012) and graphically edited in Inkscape (Supplementary Figure 2).

Daily patterns of transcript abundance were obtained via RT-qPCR. All reactions were carried out as outlined in Ruocco et al. (2019). Briefly, each reaction consisted in 5 μl Fast SYBR^® Green Master Mix (Applied Biosystems), 1 μl cDNA (1:5 diluted) template, and 4 μl of 0.7 pmol μl^–1 primers. The thermal profile of the reactions was as follows: 95°C for 20 s, 40 times 95°C for 1 s, and 60°C for 20 s. All RT-qPCR reactions were conducted in triplicate, and each assay included three no-template negative controls. To normalize target gene expression, we used the eukaryotic initiation factor 4A (eIF4A) as a reference gene (RG), which was previously demonstrated to exhibit a stable expression in Z. marina and C. nodosa along a daily cycle (Rasmusson et al., 2017; Ruocco et al., 2020a) and under a range of different conditions (Ransbotyn and Reusch, 2006; Winters et al., 2011; Olivé et al., 2017). RT-qPCR efficiencies for all primer pairs were calculated from the slopes of the standard curves of the threshold cycle (CT) vs. cDNA concentration with the equation E = 10^–1/slope. The relative quantification of the transcript levels was obtained following (Schmittgen and Livak, 2008). In details, the negative differences in the cycles to cross the threshold value between the RG and the respective GOI (-ΔCT) were calculated according to the equation:

Mean –ΔCT values were then calculated for biological replicates at each time point (n = 3) from individual –ΔCT values. Normalized values of diurnal transcript levels of GOIs have been graphically visualized using SigmaPlot v.12.5.

To assess the pattern of sugar accumulation/consumption throughout the day, the total content of non-structural carbohydrates (TNC; soluble sugars and starch) was analyzed in seagrass leaf tissues (n = 3) at five sampling times [i.e., at sunrise (06:30), solar noon (13:00), sunset (20:30), midnight (0:00), dawn (04:30)]. For TNC analysis, leaves (ca. 50 mg) were finely ground in liquid nitrogen, and subsequent extraction was carried out following Ruocco et al. (2020b). Briefly, soluble sugars were solubilized by three sequential extractions with 80% (v/v) ethanol at 80°C for 15 min. After centrifugation (3,000 rpm for 10 min), the ethanol extract was used for the determination of soluble sugars’ content, while the pellet was hydrolysed for starch determination (24 h at RT) with 3 ml NaOH 0.1 M. For both soluble sugars and starch, 3% aqueous phenol (0.25 ml) and 95%–97% H₂SO₄ (2.5 ml) were mixed with 1 ml sample in glass tubes and the solution was allowed to rest for 30 min. Absorbance was then read at 490 and 750 nm with an Agilent 8453 UV–Vis spectrophotometer. TNC content was calculated using sucrose calibration curves (standard sucrose 99%, Biorad).

A multivariate analysis was first used to assess the overall signal of all 19 GOIs (based on –ΔCT values). In details, a two-way permutational multivariate analysis of variance (PERMANOVA) based on the Bray-Curtis matrix, was conducted with the Primer 6 v.6.1.12 and PERMANOVA + v.1.0.2 software package. The analysis consisted of two fixed factors: Species (Sp) with two levels (i.e., C. nodosa and Z. marina) and Time (Ti) with six levels (i.e., midnight, dawn, sunrise, solar noon, sunset, and dusk). A principal component analysis (PCA) based on the Bray-Curtis matrix and clustering analysis were also performed on the multivariate gene-expression datasets with the software PAST v.3.03 (Hammer et al., 2001). Subsequently, univariate analyses (one-way ANOVAs and two-way ANOVAs), were used to assess the effects of time and species on daily transcript levels of individual GOIs and TNC levels (soluble sugars and starch). One and two-way ANOVAs were performed using the statistical package STATISTICA (StatSoft, Inc. v. 10, Brookline, MA, United States). Data normality was tested using the Shapiro–Wilk test, and the variance homogeneity was verified using Levene’s test or Cochran test. When assumptions of normality or homoscedasticity were not met, data were Box-Cox transformed. The Student–Newman–Keuls (SNK) post hoc test was used whenever significant differences were detected.

Principal component analysis and clustering analysis were used to separate samples collected at different times of the day and were based on multivariate gene-expression data. In C. nodosa, both approaches identified two main groups: a first including dawn, sunrise and solar noon samples and a second comprising sunset, dusk, and midnight samples (Figure 1A). Transcripts responsible for such separation are highlighted in the PCA biplot and their weightings are outlined in Supplementary Table 4. The most positively correlated transcripts with the PC1 (x-axis) were the circadian clock component LHY, and those involved in sucrose synthesis and transport (SPS1 and SUT2), whereas LHCA1 and FBA1 drove the separation of samples along the PC2 (y-axis). In Z. marina, PCA and clustering analysis identified three main groups: a first including sunrise, solar noon and sunset samples; a second including midnight and dawn samples; and a third comprising only dusk samples (Figure 1B). As evident from the biplot, the most positively correlated transcript with the PC1 was the transcript for the light harvesting protein LHCA1, while circadian clock and photoreception-related transcripts (i.e., LHY, CRY, and NPH1) were negatively correlated with the PC1. On the contrary, transcripts involved in photosynthesis (e.g., psbA, FD, and RBCS) and those involved in sugar metabolism (e.g., GWD and SPS1) were among the most positively correlated with the PC2. Overall, there was a significant difference in the diurnal pattern of accumulation of the targeted transcripts between the two seagrass species, as highlighted by the two-way PERMANOVA (Table 2 and Supplementary Table 5). A significant effect of time (Ti) (P < 0.05), species (Sp) (P < 0.001), as well as a significant Ti×Sp interaction (P < 0.001), was identified (Table 2). In further support of this, differences between species for each GOI based on two-way ANOVAs, are outlined in Supplementary Table 6.

Two-way ANOVAs revealed a significant difference of soluble sugar accumulation between the two species (P < 0.05; Supplementary Table 7). Daily levels of starch accumulation were affected by both species (Sp) and time (Ti) (Supplementary Table 7). In C. nodosa, starch level peaked at sunrise (sunrise ≠ dawn, solar noon, sunset, midnight), in respect to the rest of the day (Figure 6 and Table 5). In Z. marina, starch abundance was almost constant during the day (Figure 6).

Here, we provided a first description of daily regulation of key pathways in two seagrass species, C. nodosa and Z. marina, co-occurring at the same geographic location, thus exposed to identical natural variations in light and temperature. In both species, genes involved in sucrose metabolism (e.g., SPS1) and core circadian clock components (i.e., LHY), in addition to those encoding for light harvesting proteins (i.e., LHCA1) or involved in redox regulation (e.g., FD) exhibited clear diurnal rhythms. Such genes were also those making a greater contribution to sample distribution in the PCA (Figure 1 and Supplementary Table 4).

Despite this common behavior, the variations we observed in their daily transcriptional timing, were species-specific (Table 2 and Supplementary Tables 5, 6), highlighting the importance of intrinsic biological, and likely ecological, attributes in determining the periodicity of functions that go beyond the (common) effects of external stimuli (Müller et al., 2014; Ruocco et al., 2020a). Five GOIs (i.e., SPS1, GWD, FD, LHY and LHCA1) showed a significant Ti×Sp interaction, indicating their daily pattern of expression was substantially different between the two seagrass species (Supplementary Table 6).

In C. nodosa, a clear separation occurred between “light” (i.e., dawn, sunrise and solar noon) and “dark” (i.e., dusk, sunset and midnight) hour responses, where solar noon and midnight represented the most different time points. Yet, similar gene-expression responses were observed during light-to-dark and dark-to-light transition times. On the other hand, in Z. marina, “light-hour responses” were shifted in time, as sunrise, solar noon and dusk clustered together, while dawn-related responses were closer to those occurring at midnight. In general, in Z. marina, larger differences were detected during transition times (dusk vs. dawn), in comparison to C. nodosa (Supplementary Table 5).

Globally, the most frequent time for the max. expression (average transcript level) of GOIs was the middle of the day, followed by midnight/dawn in C. nodosa, while the min. of expression was recorded at sunset. In Z. marina instead, max. average expression occurred later, at the end of the day (i.e., during sunset/dusk), and minimum at dawn. This might suggest a differential diurnal control of growth rates between the two species, which could provide them with physiological and/or ecological benefits (Müller et al., 2014). The regulation of plant growth generally occurs through the control of carbon metabolism and light signaling, both processes being coordinated by distinct components of the circadian clock (Nozue and Maloof, 2006; Graf and Smith, 2011; Dodd et al., 2014). In Arabidopsis, growth rates peak toward the end of the night when plants are not photosynthetically active and have to properly manage metabolic resources until the next morning (Wiese et al., 2007). This regulation is also achieved via LHY/CCA1, which could have a role in ensuring the correct timing of starch breakdown at night, and ELF3 acting to repress growth in the light, thus favoring the accumulation of carbohydrate reserves (Graf et al., 2010; Nusinow et al., 2011). In our experiment, we do not have direct measurements of diurnal growth rates in C. nodosa and Z. marina, and this has never been assessed in seagrasses. However, the shifted peaks of LHY abundance and the slightly different daily patterns of carbohydrate accumulation between species (see below), could point for a different timing of growth based on a different strategy of mobilization and allocation of resources at night, coordinated by the circadian clock (Wiese et al., 2007; Graf and Smith, 2011).

In C. nodosa, low sugar levels occurring by the end of the night, could be responsible for the induction of transcripts involved in sucrose synthesis and transport in the early morning (starting from dawn), such as SPS1 and SUT2, indicating that carbohydrates have declined to critical levels, and plants have become net consumers of carbon (Bläsing et al., 2005). Subsequent to the induction of such genes, leaf soluble sugars peaked in the light, between sunrise and solar noon, demonstrating an activation of the related biochemical pathways. On the contrary, transcripts encoding proteins responsible for sucrose breakdown e.g., SUS, were induced during the night, similar to transcripts involved in starch mobilization. Glucan water dikinase, a key enzyme regulating nocturnal starch breakdown in leaves, and other starch-mobilizing enzymes (e.g., Beta-amylase and Fructose-bisphosphate aldolase) revealed cycling transcription in A. thaliana under a LD cycle and free-running conditions (Harmer et al., 2000; Lu et al., 2005) concomitant to cycling of starch breakdown. In Z. marina, diurnal rhythms of transcripts related to sugar metabolism were shifted in time in respect to what observed in C. nodosa. For instance, SPS1, involved in sucrose synthesis, peaked later in the day (between solar noon and sunset). On the contrary, transcripts encoding enzymes responsible for sugar and starch degradation/mobilization (e.g., GWD) peaked slightly earlier than C. nodosa. Endogenous sugar clock entrainment optimizes plant growth and fitness regulating carbon homeostasis between source and sink tissues (Ohara and Satake, 2017). The observed differences in gene expression could be related to different patterns of photosynthate accumulation/translocation between species: during the day, C. nodosa could quickly use triose-phosphates produced by the Calvin cycle for sucrose synthesis to support leaf maintenance and export everything else to rhizomes, where sugars can be used for rhizome expansion or stored as starch until required; on the other hand, Z. marina seems to invest more in the leaves, accumulating (exporting less) sugars and starch during the day and using starch during the night to support sucrose synthesis and export (Geiger et al., 2000).

All the above highlights the importance that sugar levels could have for diurnal gene regulation in seagrasses (Contento et al., 2004; Bläsing et al., 2005). Yet, the rate of starch degradation at night, which is known to be under control of the circadian clock (Müller et al., 2014), could also be of fundamental importance for the carbon economy of seagrass species with large effects on their short and long-term productivity and deserve further investigation (Sulpice et al., 2009; Graf and Smith, 2011).

In C. nodosa, transcripts involved in photosynthetic carbon fixation, as well as light harvesting proteins, were triggered by light, and co-ordinately peaked in the middle of the day, thus supporting the synthesis and export of sucrose needed for plant metabolism and storage. This pattern of accumulation of photosynthesis-related transcripts is very similar to the one previously detected in shallow-growing P. oceanica plants (Procaccini et al., 2017) and matches what is observed in Arabidopsis under constant light conditions, where photosynthesis genes peaked near the middle of the subjective day (Harmer et al., 2000). However, the diurnal pattern of accumulation of such transcripts was one of the most striking differences observed between the two analyzed species. Indeed, in Z. marina, transcripts for photosynthesis-related genes were slightly induced at sunrise, and then there was a down-regulation around noon followed by a further increase between sunset and dusk. We do not have a conclusive explanation for this pattern, however photo-physiological data collected in the same experiment indicated that one of the most noticeable difference between C. nodosa and Z. marina was their susceptibility to high irradiance levels at noon (Hofman, 2019). Specifically, Z. marina was much less affected by high light than C. nodosa, whose effective quantum yield (ΔF/Fm′) dropped significantly more at noon and exhibited a much higher non-photochemical quenching (NPQ) (Hofman, 2019). Besides, the max. abundance level of AOX1A, which is involved in energy dissipation mechanisms via alternative mitochondrial respiration, also differed across species, peaking at dusk in C. nodosa and at noon in Z. marina (ns). This highlights that the two species could adopt different strategies for dissipating excess of energy throughout the day, both during photosynthesis and respiration processes, and this could be somehow linked to the differential transcriptional timing we observed (e.g., for psbA) (Phee et al., 2004). However, studies at the protein level should be conducted, as transcriptome and proteome could have very different scales of regulation during the day (Graf et al., 2017).

As discussed above, diurnal rhythms observed in C. nodosa and Z. marina could be regulated by sugar availability and the circadian clock. At this regard, LHY could play a key role, as its diurnal changes are highly significant in both species. Our data confirmed that the peak of abundance of this transcript occurs at the beginning of the day (dawn/sunrise) in C. nodosa, as observed in many terrestrial species (McWatters and Devlin, 2011; Haydon et al., 2013b; Dodd et al., 2015), while in Z. marina its peak was slightly earlier (midnight/dawn). Around dawn, the phase of the circadian oscillator is indeed adjusted in response to low light intensity detected by photoreceptors as the sun rises, entraining primary metabolism, before the so-called “metabolic dawn” that occurs later due to the increasing concentration of sugars (Haydon et al., 2013a,b; Dodd et al., 2015).

The anticipation in timing of the peak in expression of LHY in Z. marina before dawn was quite unexpected. Differences in diurnal rhythms of circadian clocks components (i.e., LHY/CCA1 and TOC1) are known to occur across ecotypes (Slotte et al., 2007). Shifts in phase of LHY/CCA1 have been tracked in A. thaliana mutants exposed to artificial light conditions (Mizoguchi et al., 2002), while changes in the transcriptional peaks of several core clock components were observed in A. thaliana plants growing under different photoperiods (Flis et al., 2016). Besides the photoperiod cues, dynamic adjustment of phase and period of the plant core clock are induced in response to sugars and several biotic and abiotic stimuli (Webb et al., 2019). Further studies under controlled conditions should be performed to investigate which are the cues that are more important for controlling the circadian responses of the different seagrass species.

Transcripts encoding blue photoreceptors CRY1 and the NPH1 significantly oscillated during the day in Z. marina, where they had a very similar expression pattern, with max. expression levels at the beginning of the light phase (i.e., dawn). This is similar to terrestrial model land plants, where diurnal rhythms of mRNA levels of photoreceptors (both cryptochromes and phytochromes) exhibit their max. in the light phase, anticipating lights-on and lights-off signals (Tóth et al., 2001), and inducing change of gene expression in both nuclear and plastid transcripts (Wang et al., 2014).