Eelgrass shoots were sampled in Norway (Røvik, 67°16′06.2” N, 15°15′38.4″ E) and France (Sète, Thau Lagoon, 43°25′08.0” N, 3°40′03.9″ E). The youngest two shoots of each of six plants were collected at each site at noon and the following midnight around summer solstice (June 21) and winter solstice (December 21) 2017. Intact shoots were placed in RNAlater (ThermoFisher) and stored after 1–2 days at −20°C.

Defining Zeitgeber Times ZTs (time after an environmental time cue that entrains circadian rhythms, such as sunlight) that correspond between the French and Norwegian samples is complicated by the fact that the two latitudes do not only differ in diurnal light rhythms but also in light quantity and quality (Taulavuori et al., 2010) (i.e., light conditions at dawn at the French site do not correspond to dawn light conditions at the Norwegian site and, thus may mark different ZTs in the circadian cycle). We chose midday and midnight as time points to compare samples from the two locations as these represent times at which the plants from both locations had the potential to receive most/least light of the day in both seasons, and the time that is at equal temporal distance from sunset and sunrise. This implies, however, that differences between the Norwegian and French samples at corresponding time points may be partly explained by them being at different ZTs within their circadian rhythm. However, differences in ZTs are less affecting seasonal differences between the Norwegian and French, because we compared daily averages that involved measurements from both midday and midnight samplings.

Snapshots of temperature (°C) and light (Lux; Table 1) were taken at each sampling event in ca. 20 cm water depth with Pendant data loggers (Onset HOBO UA-002). Annual variation of temperature (°C), salinity, and Lux conditions at the Norwegian sampling site were recorded from spring 2016 to spring 2017 with data loggers (HOBO U24-002-C for salinity/temperature with 30 min intervals, and HOBO UA-002 for light/temperature with 10 min intervals) anchored ca. 30 cm above sandy substrate during low tide. At low tide, the light/temperature loggers were falling dry. Temperature (°C) and salinity conditions in the Thau Lagoon (43° 22′ 44.4” N, 3° 34′ 17.4″ E) were downloaded for year 2016 from the marine environmental station database of Thau Lagoon (Messiaen et al., 2018) on July 31, 2020, which were recorded every 30 min with data loggers (WTW) at 1–2m water depth. Light (PAR μmol m² s⁻¹) conditions were recorded with dataloggers (SKUW 215 PAR Sensor, Skye Instruments) every 30 min for December 2015–July 2016 at the Thau Lagoon Observatory (43°24′53” N, 3°41′16″ E; Trombetta et al., 2019) and were downloaded from Sea scientific open data publication (www.seanoe.org, accessed on July 31, 2020). To make the light conditions comparable between the Norwegian and French sampling sites, Lux values were converted to approximate PAR (photosynthetically active radiation) values by applying a multiplication factor of 0.019, an approximation for the spectral power distribution of sunlight (6,000 K; Creswell, 2010).

Seasonal variation in light levels was greater in Norway, where PAR levels remained close to zero from mid-November to mid-January, but reached ca. 1,000 μmolm⁻² s⁻¹ in summer, similar to France (Figures 1C,D,E; Supplementary Table 1). In summer, light levels showed a diurnal rhythm in both Norway and France, but – as expected – reached slightly higher nocturnal levels in Norway.

Summer temperatures were ca. 10°C higher in France (20–25°C) compared to Norway (10–15°C; Supplementary Figure 1; Supplementary Table 1). Stronger Norwegian diurnal temperature variations likely reflected the differences between air and water temperatures around spring tides when the dataloggers may have been shortly exposed to air during low-tides but submerged during high-tides. In contrast, the eelgrass meadow in the Thau lagoon was constantly submerged. Accordingly, the data loggers revealed relatively constant salinity throughout the year at ca. 37 PSU (Supplementary Figure 2; Supplementary Table 6), while the high fluctuation (10–30 PSU) in Norway was due to seasonal rainfall, freshwater runoff, ice formation, and melting.

At each sampling event, we assessed leaf photochemistry using a Chlorophyll a Pulse Amplitude Modulated (PAM) Fluorometer (FluorPen FP100, Photon Systems Instruments), by measuring the increase in chlorophyll a fluorescence upon illumination after a 10-min dark period induced by leaf-clips (OJIP curve; Bussotti et al., 2010; (Duarte et al., 2017). using a saturating pulse at 455 nm from 10 seagrass shoots.

When a dark-adapted leaf is illuminated with a saturating light intensity, it exhibits a so-called OJIP curve from a state where all reaction centers are opened (all quinone acceptors are oxidized and can accept electrons) to a state where all quinone acceptors are reduced and chlorophyll fluorescence reaches its maximum. The rise in fluorescence transient along this polyphasic OJIP curve depicts the rate of reduction kinetics of various components of Photosystem II (PSII). The OJIP curve provides information on the processes involved in photon absorption at the PSII antennae (ABS), transformation into an exciton and its free energy flux (RE), and transduction of its energy in the form of electronic potential vial electron transport (ET) from the donor side of the PSII to the electron acceptor of the PSI. The analysis of this OJIP transient (JIP test) allows to estimate the energy fluxes, involving absorption, entrapment in the PSII, transport to the quinone pool, and dissipation in the form of heat or fluorescence that helps to avoid photoinhibition damages. The formulae and definitions of terms derived from the JIP-test (chlorophyll a induction curves, or Kautsky curves measurements; Supplementary Table 2) are listed in Table 2.

Differences in these photosynthetic characteristics between populations (day/night/season) were tested for significance with Wilcoxon rank sum tests followed by correcting the values of p for multiple pairwise comparisons using the Benjamini-Hochberg method with the R package “rstatix” v0.6.0 (Kassambara, 2020). Multivariate differences in photochemistry were estimated from photochemical raw data using a canonical analysis of principal coordinates (CAP) with Euclidean distances with the Primer 6 software (Clarke and Gorley, 2005). This multivariate approach is insensitive to heterogeneous data and frequently used to compare different sample groups (Duarte et al., 2017; Cruz de Carvalho et al., 2020). Differences in the overall photochemical metabolism were estimated from normalized photochemical raw data with a canonical analysis of principal coordinates (Duarte et al., 2017). Multivariate differences were visualized as Euclidean distances along the first two principal coordinates.

RNA was extracted from small pieces of stored samples (RNAlater; ThermoFisher) that were homogenized (PowerLyzer 24; Qiagen) and processed according to the standard protocol of the RNeasy PowerPlant Kit (Qiagen). We consistently excised the mid portion of the youngest fully mature leaf in the shoot (usually the second-rank leaf), so that tissue-age dependent within- or among-leaf variations were unlikely to explain differences in gene expression (Ruocco et al., 2019a,b). RNA concentrations were measured with Qubit™ RNA HS Assays (ThermoFisher) using a Qubit 4 fluorometer. RNA integrity was verified on TapeStation 2200 with RNA ScreenTape. RIN values were >3.6.

RNA libraries were prepared from 1 μg of extracted RNA with NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina (New England Biolabs) following the polyA mRNA workflow of the kit protocol. Libraries were verified on TapeStation 2200 using High Sensitivity D1000 ScreenTape. Paired-end (2 × 150 bp) RNA sequencing (RNASeq) data were generated using a NextSeq 500/550 High Output Kit v2 (Illumina) on an Illumina NextSeq 500 platform at Nord University.

Raw reads were trimmed for low quality (Phred score Q<20, 99% base call accuracy), and TruSeq adapter sequences (overlap > 3bp) before reads <25bp were removed with TrimGalore! v 0.6.0.¹ Read quality was checked with FastQC v 0.11.8 to control for aberrant read base content, length distribution, duplication and over-representation.²

Trimmed reads were mapped to the genome of Z. marina (Olsen et al., 2016) v 2.2 from ORCAE (Sterck et al., 2012) and the organellar genomes (chloroplast and mitochondria) of Z. marina v2.1 from ORCAE (Sterck et al., 2012) with the splice-aware RNA-seq aligner STAR v 2.7.2a (sjdbOverhang 149; Dobin et al., 2013, 2016), guided by gff3 gene annotations v2.2 from ORCAE (Sterck et al., 2012), converted to gtf with the program gffread.³ Alignments that contained noncanonical splice junctions were filtered out. Duplication rates were checked with the R package “dupRadar” v1.14.0 (Sayols et al., 2016).

Read counts for 21,069 genes (Supplementary Table 3) were obtained with featureCounts in the R package “Rsubread” v1.34.6 (Liao et al., 2019) with the non-default settings minOverlap = 3, MinMQS = 20, strandSpecific = 2, isPairedEnd = TRUE. Genes of low expression (library average<5, n = 7,023) were excluded from subsequent analyses to avoid potential artefacts from sequencing errors.

Gene IDs of the Z. marina genome v2.2 (Olsen et al., 2016) were annotated with gene descriptions, Gene Ontology (GO) terms, KEGG Orthology groups (KO), and Enzyme Commission (EC) numbers (based on enzyme nomenclature) from the ORCAE database (Sterck et al., 2012). Gene ontology terms assign genes to one or more categories related to the molecular function, cellular component or biological process of the gene products. KEGG Orthology groups assign gene products to their function in a biological pathway of molecular interactions. The EC number classifies genes numerically based on the chemical reactions they catalyze.

To characterize overall transcriptomic differentiation, we clustered the samples hierarchically by the first five principal components of gene coverage (14,046 genes), which together explained 73.4% of the variation, using the R package “FactoMineR” v. 1.42 (Lê et al., 2008). For Principal Component Analysis (PCA), the scale.unit was set to FALSE not to scale the expression values to unit variance. Prior to PCA, we regularized log-transformed gene count values (Supplementary Table 4) with DESeq2 v 1.24.0 (Love et al., 2014) in order to assure homoscedasticity and to account for differences in sequencing depth.

We tested for differential gene expression between summer and winter samples (day and night samples combined) and between day and night samples in France and Norway, respectively, with DESeq2 v1.24.0 (Love et al., 2014). The values of p were adjusted with Independent Hypothesis Filtering based on the R package “IHW” v1.14.0 (Ignatiadis et al., 2016). We plotted overlaps of differentially expressed genes between seasons and locations with the R package “eulerr” v6.1.0 (Larsson, 2019) and “nVennR” v0.2.1 (Pérez-Silva et al., 2018).

For each population, we tested (Fisher’s exact test) for enriched gene ontology terms (biological process, cellular component, and molecular function) in genes that were differentially expressed between day and night samples or between summer and winter samples with the R package “topgo” (Alexa and Rahnenführer, 2010). Significance levels were adjusted with the Benjamini and Hochberg correction (Benjamini and Hochberg, 1995), using the p.adjust function in R (R Core Team, 2019), to control for the false discovery rate when conducting multiple comparisons.

The adjusted-log10 of the values of p of significantly enriched GO-terms were visualized for each population in heatmaps using the heatmap.2 function of the R package “gplots” v.3.0.3 (Warnes et al., 2020). GO terms were hierarchically clustered with the R package “goSTAG” v1.12.0 (Bennett and Bushel, 2017) based on euclidean distances and Ward’s minimum variance clustering method with Ward’s (Ward, 1963) clustering criterion implemented (ward.D2 option; Murtagh and Legendre, 2014). The sample dendrogram was based on Euclidean distances and a complete linkage clustering method.

Enriched/perturbed KEGG metabolic pathways in the fold-change matrices of differentially expressed genes were identified using the R package “gage” v2.37.0 (Luo et al., 2009). Nuclear genes were mapped to KO (Kegg Orthology)-terms by translating EC numbers from the genome annotation to corresponding KO-terms based on the KEGG API list http://rest.kegg.jp/link/ec/ko. For organellar genes, EC numbers and, thus KO-terms were not annotated. For pathways perturbed (p < 0.05) in at least one of the populations and seasons, we plotted the significant fold-changes of molecular functions (KO-terms) for all populations and seasons with the R package “pathview” v3.11 (Luo and Brouwer, 2013). When several genes referred to the same KO-term, we plotted their average fold-change in the respective pathway.

Overall photosynthetic performance, represented by the photochemical performance index (PI_ABS; Figure 2A) showed inverse seasonal shifts in Norway and France. While the French population performed significantly better during summer, the Norwegian population performed significantly better during winter. The reduced summer performance of the Norwegian population was associated with peculiar photosynthetic characteristics. In contrast to the French population and the Norwegian winter samples, light harvesting ability increased in the Norwegian summer samples under non-saturating light levels during night-time, as represented by the increase in absorbed (ABS/CS; Figure 2C) and trapped (TR₀/CS; Figure 2D) energy flux, and in the number of oxidized reaction centers per leaf cross section (RC₀/CS; Figure 2G). The contribution of dark reactions to the overall photochemistry [Ψ₀/(1-Ψ₀); Figure 2I] decreased significantly during summer nights. However, despite readiness of the electron transport chain (ETC) for electron transport, as indicated by the significant increase in the size of the oxidized quinone pool during Norwegian summer nights (Area; Figure 2B), the electron transport energy flux (ET₀/CS; Figure 2E) did not increase significantly. During winter, despite a significant increase in the contribution of light and dark reactions to the primary photochemistry [TR₀/DI₀, and Ψ₀/(1−Ψ₀); Figures 2H,I], the Norwegian population showed a significant reduction in all energy fluxes, particularly during night-time (ABS/CS, TR₀/CS, and ET₀/CS; Figures 2C–E), including dissipated energy loss (DI₀/CS; Figure 2F).

The samples clustered well by sampling location and daylength (78 of 80 samples correctly identified) along the first two principal components showing that the photochemical traits from the chlorophyll a induction curves could well classify the samples (Figure 3). The Norwegian summer samples grouped into one cluster, irrespective of whether they were collected at day or night. The Norwegian winter samples collected during nighttime formed another distinct cluster. Winter daytime samples from France and Norway grouped close together, indicating that winter had a prevalent effect over location.

On average, 76% of 1 million raw reads (average per library) mapped to exons (Supplementary Table 5). Based on the expression of 13,932 genes covered >5 times, the samples grouped into four clusters, separated by location (France/Norway) and season (Summer/Winter; Figure 4). The most distinct gene expression patterns were recorded in the winter samples from Norway (NW). Gene expression was most similar between winter samples from France (FS) and summer samples from Norway (NS; Figure 4). Strictly distinct clusters of gene expression between day and night were only observed in France during summer (FS).

Seasonal differences in gene expression were found in 8652 genes in the French population and in 8569 genes in the Norwegian population (Supplementary Table 6). In summer, 4,274 genes were upregulated in France, and 4,084 genes in Norway, of which 1,475 overlapped between the two populations (Figure 1A). In winter, 2,691 genes were upregulated in France, and 2,798 genes in Norway, of which 1,687 overlapped between the two populations (Figure 1A).

Diurnal differences in expression were strongest (>3,000 genes) in, and most similar (>500 genes shared) between the French summer samples and the Norwegian winter samples (Figure 1B). In contrast, the number of genes with diurnal expression differences was an order of magnitude lower (ca. 600 genes) in the French winter samples or the Norwegian summer samples (Figure 1B). All genes with diurnal differences in gene expression are characterized in Supplementary Table 7.

Genes involved in photosynthesis, carbon fixation (dark reactions of photosynthesis; Supplementary Figure 3), and sugar metabolism (fructose and mannose metabolism: Supplementary Figure 4 and glycolysis: Supplementary Figure 5) showed contrasting seasonal patterns between the Norwegian and French seagrass populations. Those are likely to explain the functional basis of the strong transcriptome similarity between Norwegian winter samples and French summer samples (Figure 4). Summer upregulation/winter downregulation of most of the organellar functions, both related to the chloroplast and the mitochondrion, characterizes the functional transcriptome specificities of the Norwegian samples (Figure 5). All pathways enriched for seasonal or diurnal differentially expressed genes are listed in Supplementary Table 8.

Upregulation of photosynthesis-related genes in the Norwegian population during summer (Figure 6) was also shown in the seagrass P. oceanica when exposed to high light levels (Dattolo et al., 2014) and is likely to enable the plants to process the increased amount of captured light energy (Walters, 2005). Upregulation of genes with functions in carbon capture during summer in the Norwegian (but not the French) samples (Supplementary Figure 3) is necessary for the rapid build-up of stored carbohydrates during the short summer period, and provides a molecular explanation for increasing autumn starch content with latitude (in Zostera noltii; Soissons et al., 2018).

Surprisingly, the light intensity at the Norwegian sampling site showed a clear day-night rhythm during summer. Even though the summer sun does not set below the horizon, shading by mountain ridges and the low sun angle during night-time appears to create a short, but clearly noticeable drop in available photosynthetically active radiation. This diurnal light rhythm was paralleled by diurnal changes in gene expression. For example, the daytime upregulation of arginine catabolism (Figure 7) may mobilize nitrogen storages and thus, allow for active growth under nitrogen limitation (Winter et al., 2015). In addition, the nighttime upregulation of trehalose biosynthesis (Figure 7), in which the intermediate compound trehalose-6-phosphate acts as a sensor for available sucrose, is likely involved in the regulation of sucrose-levels that have built up over the day (Lunn et al., 2014; Figueroa and Lunn, 2016) and plays a role in plant growth, development, and protection against various environmental stressors (Lunn et al., 2014; Joshi et al., 2020). Thus, the diurnal changes in gene expression during summer are likely related to active growth and processing of assembled sugars. Local geography, e.g., the presence of mountain ridges, may play an important role in allowing eelgrass to process a dark phase in high latitudes.

Despite diurnal fluctuations in daylight levels, the lack of complete darkness during summer significantly lowers the photosynthetic performance (PI_ABS) of the Norwegian samples (Figure 2A), in contrast to the seasonal change in photosynthetic performance in the French population (with significantly higher values during summer: Figure 2A). Particularly during summer nights, the Norwegian population showed a significant increase in absorbed (ABS/CS; Figure 2C), trapped (TR₀/CS; Figure 2D), and electron transport (ET₀/CS) energy fluxes, as well as the size of the oxidized quinone pool (ABS) and the number of oxidized reaction centers per leaf cross section (RC₀/CS; Figure 2G). Although these photosynthetic adjustments increase the light harvesting ability and efficiency to use low-light, the Norwegian samples also showed a significant increase in energy dissipation (DI/CS), suggesting an inability to use excessive amounts of light energy in summer.

Low light intensities during the polar summer nights appear to have allowed the Norwegian plants to process a dark phase to recycle intermediates of the Krebs cycle (such as malic acid) that are consumed during the light phase of photosynthesis. The dark phase is triggered by an inhibitor (2-carboxylarabinitol-1,5-bisphosphate, CAP) that is only produced under low light when the chloroplast pH shifts from acid to basic. It was not found in any of the genes showing diurnal differential expression (Supplementary Table 7). A crucial limit to the northward shift of seagrass may be set by the latitude at which the sunlight is too strong to allow the procession of a dark phase. The existence and identification of this latitude, however, relies on future studies involving transplant experiments to a range of increasing latitudes.

Our findings suggest that eelgrass can store sufficient energy reserves to survive dark periods lasting longer than those at 67°N (our Norwegian sampling site). In fact, if sugar levels were depleted in mid-winter, we would have expected to see an upregulation of genes involved in protein degradation, as autophagic protein degradation has been suggested as an alternative energy source under sugar depletion in eelgrass (Mochida et al., 2019). As the expected upregulation was not observed, a poleward shift of eelgrass should not be hindered by the potential to store sufficient photosynthates. In order to test this hypothesis, future studies must evaluate carbohydrate reserves in autumn, a good indicator of winter survival in the seagrass Zostera noltii (Govers et al., 2015), and reserve depletion over long dark periods directly in below-ground rhizomes.

While energy storing rhizomes allow eelgrass to survive dark periods, they also represent respiratory burdens that can increase light requirements, potentially preventing poleward migration (Olesen and Sand-Jensen, 1993; Lee et al., 2007). This is particularly true as respiration increases under rising temperatures (Lee et al., 2007; McMinn and Martin, 2013), faster than photosynthesis (higher Q10; Thom et al., 2008; Staehr and Borum, 2011). The respiratory demands of seagrass roots and rhizomes explain the relatively high amounts of light (10–30% surface irradiance; Ochieng et al., 2010) required for long-term survival, as compared with 1–5% for many species of phytoplankton and algae (Kenworthy and Fonseca, 1996). Below-ground rhizomes in eelgrass consume up to 15% of the carbon fixed, and up to 25% at the deepest range-margin (Kraemer and Alberte, 1993). Moreover, even if carbohydrate concentrations are not depleted, seagrass originating from 36°N cannot survive extreme experimental winter light limitation (2 h irradiance saturated photosynthesis), presumably through the effects of anoxia on increased translocation of energy reserves (Alcoverro et al., 1999). Thus, increasing carbohydrate storage does not guarantee long-term survival in the dark.

The need of eelgrass to elevate carbon storage in warming Arctic oceans may be lowered by thermal acclimatization that increases respiratory energy efficiency, as has been shown in terrestrial Arctic tundra species (Kornfeld et al., 2013), or by increasing CO₂ levels (Zimmerman et al., 1997). Eelgrass is generally carbon limited, because it cannot efficiently use the inorganic carbon from the water column (HCO3-; Lee et al., 2007). Thus, increasing CO₂ levels are expected to increase photosynthetic rates without affecting respiration (Zimmerman et al., 1997). If turbidity does not hinder long-term survival (Moore et al., 1997) in the future Arctic, increasing CO₂ levels will reduce the need for producing high below-ground biomass by shortening the seasonal period during which eelgrass experiences a negative carbon balance.

In this study, we show that distinct seasonal expression patterns of an eelgrass population from the Norwegian subarctic are driven by genes encoding functions of the organelles (chloroplast and mitochondrion). These are likely to provide key physiological adaptations that allow eelgrass to endure the polar dark season, and suggest metabolic strategies to enhance restoration success in turbid waters. The two most prominent features were:

This study alone, based on two meadows and an assessment of metabolic and transcriptomic responses, does not allow us to predict whether eelgrass will be able to migrate poleward into a warming, ice-free Arctic. Although this study reveals specificities that are functionally relevant at the species’ current northernmost range of distribution, its potential to migrate further poleward depends on the plastic and/or adaptive potential of the identified key functions. While genetic adaptations would rely on standing genetic diversity that allows for evolutionary change over several generations, acclimatization – which can build up within the lifetime of an organism (Kelly et al., 2012; Donelson et al., 2019) – would allow eelgrass to immediately occupy niches that open up in the Arctic. For example, in the seagrass Thalassia testudinum, phenotypic variation among subpopulations across an environmental gradient in the Florida Bay were mainly ascribed to phenotypic plasticity because gene flow was high within 50 km distance (Bricker et al., 2011). Further studies are needed to identify to what degree the north–south differences that we identified result from acclimatization or genetic adaptation (Box 1):

In order to test whether the reaction norm of eelgrass is wide enough to acclimatize to a warming Arctic, we need to expose eelgrass to a simulation of near-future polar conditions. For example, seagrass meadows in Australia (Halodule pinifolia and Halophila ovalis) were shown to die at >38 days at 0.1mol photons m⁻² d⁻¹ (0.35% SI; Longstaff and Dennison, 1999). To determine the latitude at which eelgrass would lack sufficient daylight throughout the year, we would need an equivalent study that imitates future temperature and CO₂ levels at different daylight lengths in a common-garden setup and, ideally, includes a variety of genotypes – since these can differ in their dark tolerance (Salo et al., 2015).