Daily Regulation of Key Metabolic Pathways in Two Seagrasses Under Natural Light Conditions

The circadian clock is an endogenous time-keeping mechanism that enables organisms to adapt to external environmental cycles. It produces rhythms of plant metabolism and physiology, and interacts with signaling pathways controlling daily and seasonal environmental responses through gene expression regulation. Downstream metabolic outputs, such as photosynthesis and sugar metabolism, besides being affected by the clock, can also contribute to the circadian timing itself. In marine plants, studies of circadian rhythms are still way behind in respect to terrestrial species, which strongly limits the understanding of how they coordinate their physiology and energetic metabolism with environmental signals at sea. Here, we provided a first description of daily timing of key core clock components and clock output pathways in two seagrass species, Cymodocea nodosa and Zostera marina (order Alismatales), co-occurring at the same geographic location, thus exposed to identical natural variations in photoperiod. Large differences were observed between species in the daily timing of accumulation of transcripts related to key metabolic pathways, such as photosynthesis and sucrose synthesis/transport, highlighting the importance of intrinsic biological, and likely ecological attributes of the species in determining the periodicity of functions. The two species exhibited a differential sensitivity to light-to-dark and dark-to-light transition times and could adopt different growth timing based on a differential strategy of resource allocation and mobilization throughout the day, possibly coordinated by the circadian clock. This behavior could potentially derive from divergent evolutionary adaptations of the species to their bio-geographical range of distributions.


INTRODUCTION
Earth's rotation causes repetitive changes of day and night that are reflected in diurnal cycles of light and temperature. Most living organisms experience these changes and thus face the challenge of coordinating their lives with such rhythms. This was the primary driver for the emergence and evolution of endogenous clocks that can be set by the rising or the setting of the sun (McWatters and Devlin, 2011). In plants, circadian clocks regulate many aspects of physiology, growth, development and reproduction (Harmer et al., 2000;Harmer, 2009). These produce daily rhythms of metabolism, and interact with signaling pathways controlling daily environmental responses (Dodd et al., 2015;Kim et al., 2017). Ultimately, circadian clocks allow plants to anticipate daily and seasonal changes in the environment, conferring them an adaptive advantage (Ouyang et al., 1998).
Over the last decade, a combination of empirical research and mathematical modeling has led to the identification of several clock components in terrestrial model species and an understanding of their regulative roles within plants' biochemical pathways (Webb and Satake, 2015). Three interconnected modules mainly compose the plants' circadian-clock system: the receptors that perceive environmental signals and provide the inputs, the core oscillator that generates the rhythm itself, and the outputs that determine physiological and metabolic rhythms (McWatters and Devlin, 2011). In Arabidopsis thaliana, the clockwork is regulated by at least three tightly interlocked feedback loops, where clock proteins regulate their own transcription, either directly or indirectly (Barak et al., 2000;Harmer, 2009). These loops of gene expression patterns are synchronized by post-translational modifications of their own proteins.
Different photic cues (e.g., day-length, intensity, spectral composition of sun light and solar position) may contribute to the entrainment of circadian clock of plants (Millar, 2004). These signals are perceived via three regulatory families of photoreceptors including phytochromes (PHY), which are sensitive to red/far-red, and cryptochromes (CRY) and phototropins (PHOT), which adsorb in the UV-blue light (Millar, 2003), although phototropins seems not to have a direct role in light input into the circadian clock (Jones, 2009;Litthauer et al., 2016). The ZEITLUPE (ZTL) family that belongs to the LOV domain photoreceptors is the third family of blue light photoreceptors in plants, with a central role in circadian clock and photoperiodic flowering (Song et al., 2014).
Genome-wide studies in terrestrial plants have revealed that a large portion (up to 89%) of the expressed transcripts show a rhythmic expression under diverse circadian, photocycle and thermocycle growth conditions (e.g., Harmer et al., 2000;Michael et al., 2008;Zdepski et al., 2008;Filichkin et al., 2011;Jończyk et al., 2011). This is especially true for transcripts involved in main physiological pathways, such as sugar synthesis/utilization, which imply a major role of the clock system in coordinating key aspects of the plant metabolism (Harmer et al., 2000;Covington et al., 2008;Haydon et al., 2013a). Several sugar-responsive transcripts show rhythmic changes that match the diurnal variation in sugar levels (Bläsing et al., 2005). On the other hand, sugar levels themselves can entrain circadian rhythms by regulating the expression of clock-components (Haydon et al., 2013b). In general, the circadian oscillator is coupled with carbohydrate biochemistry and there is a tight connection between nuclear-encoded circadian-clock components and processes occurring in the chloroplasts (Dodd et al., 2015). For instance at dawn the phase of the circadian oscillator responds to the low light intensity detected by photoreceptors. Then, a second so-called "metabolic dawn" occurs in response to the first accumulation of sugars due to the sensitivity of clock components such as LHY/CCA1 (Haydon et al., 2013b;Dodd et al., 2015).
The daily dynamic of starch production and consumption is also known to be regulated by the circadian oscillator with feedbacks from carbon availability (Haydon et al., 2013b;Dodd et al., 2015;Seki et al., 2017;Flis et al., 2019). The molecular mechanisms that coordinate the turnover of starch with the external photoperiod are currently unknown, although two different models have been proposed, mainly differing in the way the temporal and metabolic cues are integrated. In the first, it is assumed that the abundance of starch is measured (Scialdone et al., 2013), while in the other, it is assumed that sucrose is measured (Feugier and Satake, 2013;Webb and Satake, 2015;Seki et al., 2017). In the starch sensing models, the circadian clock is a passive timer used to measure the time of day (Scialdone et al., 2013). In the second, sucrose feedbacks to the circadian clock to dynamically regulate the phase of the circadian oscillator (Feugier and Satake, 2013;Seki et al., 2017).
The structure and function of the circadian system might be widely conserved across higher plant species (Serikawa et al., 2008;Filichkin et al., 2011); however, the regulation of clock-output traits can be species-specific such that predictions from model species as Arabidopsis are not always possible. For example, recent studies have highlighted fundamental differences in the circadian behavior of traits like growth rate in dicotyledonous vs. monocotyledons species (Müller et al., 2014). Such differences question if the clock oscillator has the same importance to regulate physiology, metabolism and growth across species (Müller et al., 2014). In addition, several studies have examined the potential ecological implications of endogenous differences in plant photo-perception and timing systems (Resco et al., 2009) on questions such as biogeographical distribution of species and range-shifts under climate changes (Huffeldt, 2020), and on the relevance that circadian regulation of physiology has on the interactions at ecosystem level (Resco de Dios et al., 2016;Lu et al., 2017;Hubbard et al., 2018) and biosphere-atmosphere regulations (Resco de Dios, 2013).
If the transfer of knowledge regarding circadian-regulated patterns is already complex across species inhabiting similar habitats, this is even more difficult in species experiencing completely different habitats, such as seagrasses that live submerged in marine environments. Seagrasses are a polyphyletic assemblage of ∼70 species, which evolved from land ancestors and returned to the sea in the Cretaceous period (Les et al., 1997). Seagrasses are distributed along temperate and tropical coastlines worldwide (Short et al., 2007) where they form marine forests that are among the most productive and valuable costal ecosystems (Ruiz-Frau et al., 2017) and contribute to climate change mitigation (Marbà et al., 2015). The current distribution of seagrass species across latitudes reflects their specific requirements in terms of light and temperature for growth (Short et al., 2007) and phenology (Blok et al., 2018;Yue et al., 2020), which could be influenced by circadian and/or photoperiodic responses.
Some of the biggest challenges imposed by the marine environment in respect to the terrestrial one, are the differences in light intensity and quality (Kirk, 2010), as well as the lower CO2 availability (Larkum et al., 2006(Larkum et al., , 2018. The exposure to such a different habitat could have influenced seagrasses' circadian clock, in terms of composition of gene networks and their subsequent regulation (Olsen et al., 2016). A first evidence of rhythmic diel oscillations of transcripts related to photosynthesis and respiration in seagrasses has come from preliminary reports in Posidonia oceanica and Zostera marina Rasmusson et al., 2017). However, studies of circadian rhythms in marine plants are still way behind in comparison to terrestrial species, notably for the lack of genomes available for most species and the inapplicability of fundamental molecular genetics approaches. This strongly limits the possibility to investigate the interconnection between energy metabolism and clock systems in marine plants, as well as how circadian and photoperiodic responses may influence their geographical distribution, productivity and survival under climate change.
Here, we want to provide a comprehensive overview of daily transcriptional timing of putative core circadian clock components and metabolic pathways such as sugar production/utilization, photoreception, photosynthesis, and cellular respiration, as well as patterns of carbohydrate accumulation, in two seagrass species, Cymodocea nodosa and Zostera marina (order Alismatales). These species are characterized by distinct biogeographical distributional ranges: C. nodosa typically spreads in warmer and temperate waters (from Mauritania to the Canary archipelago and Mediterranean Sea), while Z. marina is the dominant seagrass species in colder water of North hemisphere (Short et al., 2007). We explored daily responses of the species when co-occurring at the same geographic location (Ria Formosa lagoon, Faro, Portugal), thus being submitted to identical light and temperature cycles. This area represents the northern and southern distribution edge in the Atlantic region of C. nodosa and Z. marina, respectively. Our hypothesis is that possible inter-specific variations in daily metabolic rhythms, could reflect e.g., evolutionary adaptations of the species to divergent latitudinal distribution related to their optimal photic and thermal requirements.
Our study was conducted under natural light conditions, thus it aims at uncovering the rhythms of key pathways in the two species under naturally fluctuating conditions.

Seagrass Sampling and Experimental Design
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.

RNA Extraction and Target Genes Selection
Total RNA from the youngest fully developed leaves of C. nodosa and Z. marina (Ruocco et al., 2019) was extracted with the Aurum TM 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 TM 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 TM 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  Gene acronym, protein name, species, primer sequences, percent efficiency (E), correlation coefficient (R 2 ), references, and sequence ID are given. Contigs ID in C. nodosa refers to the transcriptome by Ruocco et al. (2017) or GenBank acc. no, in Z. marina sequences ID refer to those deposited in ORCAE (https://bioinformatics.psb.ugent. be/orcae/) (Olsen et al., 2016). *The same primer couple in C. nodosa and Z. marina. §Primers designed on orthologues sequences in C. nodosa and Z. marina.
A. thaliana were downloaded from TAIR 1 and used as queries for BLASTP (Altschul et al., 1990). For Z. marina, sequences were blasted against the genome available in ORCAE 2 (Olsen et al., 2016). For C. nodosa, potential homologous sequences were identified through sequence similarity search against an available transcriptome  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  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 Gene-Expression Analysis
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 R 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  -CT values. Normalized values of diurnal transcript levels of GOIs have been graphically visualized using SigmaPlot v.12.5.

Total Non-structural Carbohydrate Determination
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 2 SO 4 (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).

Statistical Analysis
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.

Daily Cycle of Gene Expression in Cymodocea nodosa and Zostera marina
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 Significant effects (P < 0.05) of the time of the day on transcript levels are highlighted in bold. P < 0.1 are underlined. For SNK pairwise tests, P < 0.1 are indicated in brackets.
Frontiers in Ecology and Evolution | www.frontiersin.org 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 photoreceptionrelated 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.

Photosynthesis and Light Harvesting
Transcripts encoding for main subunits of Photosystem I (PSI) and II (PSII) (psbA and psbD), RuBisCO small subunit (RBCS) and ferredoxin (FD) in C. nodosa exhibited a very coordinate expression pattern, with max. abundance between solar noon and sunset, and a subsequent slow-down during the rest of the day (Figure 2). A similar trend was observed for transcripts encoding for light harvesting proteins, which exhibited max. abundance at noon and min. from midnight to sunrise (Figure 2). Among photosynthesis-related genes, psbD and FD were the ones significantly affected by the daily time according to the oneway ANOVA (psbD: sunset = dusk; FD: solar noon = dusk) ( Table 3). Expression levels of LHCA1 was also significantly different between many time points during the day ( Table 3).
In Z. marina, photosynthesis-related genes, after a first peak  at sunrise, decreased their abundance until solar noon, and then their level started to increase again until dusk (Figure 2). However, this pattern was only statistically significant for FD (Table 4). LHCA1 peaked between sunrise and sunset ( Table 4).

Sucrose and Starch Metabolism
In C. nodosa, transcripts encoding for proteins involved in sucrose synthesis and transport (SUT2 and SPS1) peaked at the beginning and during the light phase (i.e., at sunrise and solar noon), then their abundance decreased throughout the day, with a min. of expression at midnight (Figure 3). On the contrary, transcripts with a role in sucrose and starch metabolism globally tend to increase their expression levels at the end of the day, peaking during the night phase (i.e., from dusk to dawn) (Figure 3). In details, transcript abundance of SUT2 significantly varied between sunrise and midnight (P < 0.05), while SPS1 exhibited significant variations between sunrise and midnight (P < 0.05), and solar noon and midnight (P < 0.05) ( Table 3). Sucrose synthase (SUS) displayed significant differences between solar noon and midnight (P < 0.05) ( Table 3), while the other GOIs within this category did not exhibited significant changes.
In Z. marina, SPS1 peaked later in the day (i.e., at sunset), while SUT2 showed no differences during the day (Figure 3 and Table 4). Genes related to sucrose and starch metabolism did not exhibited a very coordinate expression pattern, although their levels of abundance were also generally higher during dark hours. GWD exhibited significant changes at dawn, in respect to other time points (Figure 3 and Table 4).

Circadian Clock and Photoreception
In C. nodosa, the circadian clock component LHY displayed a peak of abundance in the early morning (i.e., from dawn to sunrise) and a min. of expression between sunset and midnight (Figure 4 and Table 3). ZTL, as well as bluelight photoreceptors (NPH1 and CRY1) showed no significant expression differences among the selected time points (Figure 4). In Z. marina, the circadian-clock component LHY reached its max. level slightly earlier than C. nodosa, between midnight and dawn (Figure 4 and Table 4). Blue-light photoreceptors exhibited a similar pattern, as they showed max. expression levels between midnight and dawn (Figure 4). According to oneway ANOVA, levels of LHY differed significantly at midnight FIGURE 4 | Daily transcript levels (as normalized -CT data) of GOIs with a function in circadian clock and photoreception in C. nodosa (upper panels) and Z. marina (lower panels). Significant daily changes (P < 0.05) for each GOI based on SNK tests following one-way ANOVAs (Tables 3, 4), are indicated with different letters. Letters are coloured as the respective GOI. Values are mean ± SE for n = 3. and dawn in respect to the rest of the day (Table 4), while NPH1 and CRY1 varied significantly their expression between solar noon and dawn (P < 0.05 for NPH1 and P = 0.06 for CRY1) ( Table 4).

Mitochondrial Respiration and Growth
Regarding respiration-related genes, the transcript for Cytochrome C (COX5B) showed no significant expression changes among time points in both species. The transcript encoding for AOX1A instead exhibited significant oscillation during the day in C. nodosa, with max. abundance around dusk ( Figure 5 and Table 3). The transcript for Expansin B (EB) exhibited no significant daily changes in both species.

Daily Patterns of Carbohydrate Accumulation in Cymodocea nodosa and Zostera marina
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).

General Overview of Daily Timing in Cymodocea nodosa and Zostera marina
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 FIGURE 5 | Daily transcript levels (as normalized -CT data) of GOIs with a function in cellular respiration and growth in C. nodosa (upper panels) and Z. marina (lower panels). Significant daily changes (P < 0.05) for each GOI based on SNK tests following one-way ANOVAs (Tables 3, 4), are indicated with different letters. Letters are coloured as the respective GOI. Values are mean ± SE for n = 3. 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 lightto-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 dawnrelated 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 FIGURE 6 | Daily patterns of TNC (soluble sugars and starch) accumulation in C. nodosa and Z. marina. Significant daily changes (P < 0.05) based on SNK tests following one-way ANOVAs (Table 5), are indicated with different letters. Values are mean ± SE for n = 3.
allocation of resources at night, coordinated by the circadian clock (Wiese et al., 2007;Graf and Smith, 2011).

Specific Daily Regulation of Individual Genes of Interests in Cymodocea nodosa and Zostera marina
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 triosephosphates 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  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 photosynthesisrelated 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).

CONCLUSION
Here, we provided a first overview of daily transcriptional timing of key metabolic pathways in two marine plants exposed to the same photoperiod. A strength of our study is that plants were exposed to natural irradiance conditions, with irregular fluctuations and gradual LD transitions at dusk and dawn. Experiments with terrestrial species are typically conducted exposing plants in controlled environments with constant (artificial) daily irradiance and abrupt LD transitions. Such studies are often not representative of natural settings, and diurnal variations of carbon metabolism can differ significantly between sunlight and artificial light conditions (Vialet-Chabrand et al., 2016Annunziata et al., 2017;Panter et al., 2019).
Overall, our results suggested a potential interplay between circadian clock and sugar levels in determining the periodicity of functions in both species, as demonstrated in terrestrial plants, although other factors could be at play. Large differences were observed between species in the daily regulation of transcripts related to photosynthesis and sugar metabolism, which together with the slightly different daily patterns of carbohydrate accumulation, highlight the importance of intrinsic biological and likely ecological attributes of the species. The two species could adopt e.g., different growth timing based on a differential strategy of resource allocation and mobilization, coordinated by the circadian clock. Such behavior could also derive from divergent evolutionary adaptations of the species to different photoperiodic conditions across their bio-geographical range of distributions. This could also partially explains the differential sensitivity of the two species during light-to-dark and dark-to-light transitions. Species-specific effects of twilight duration on circadian phases of animal behavior (Boulos et al., 1996) and physiology (Sosniyenko et al., 2009) have been well documented, as well as the occurrence of circadian phenotypes correlated with latitude and geographical distribution clines either in wild plant populations or within selectively bred crops (Greenham et al., 2017;Rees et al., 2021).
We are aware that our study has limitations due to a single day analysis. In addition, we can only infer, e.g., the circadian regulation of certain patterns, as this should be confirmed with ad hoc experiments under constant light or dark conditions. However, we believe our analysis addressing the diurnal response of two species along a single day is informative and broadens current knowledge of biological rhythms in seagrasses. Studies of circadian rhythms in marine plants are currently missing, although these could have wider implications, as the modulation of such rhythms could be used for improving marine plant performances and productivity (Kim et al., 2017;Steed et al., 2021) and may ultimately help in buffering the effects of climate changes on seagrass distribution and persistence (Resco et al., 2009).

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
ED and MR conceived the ideas and designed the study. ED, MR, JS, IB, JH, and KP participated in sample collection and mesocosm maintenance. MR and ED performed all geneexpression and biochemical analyses. MR led the writing of the manuscript. All authors contributed in drafting the work, revised it critically, and gave final approval for publication.

FUNDING
The research leading to these results received funding from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement no 730984, ASSEMBLE Plus project. This study also received Portuguese national funds from FCT -Foundation for Science and Technology through the projects GRASSMET (PTDC/MAR-EST/4257/2014) and UIDB/04326/2020. The work has been partially supported by the project Marine Hazard, PON03PE_00203_1, Italian Ministry of Education, University and Research (MIUR).