Front. Genet. Frontiers in Genetics Front. Genet. 1664-8021 Frontiers Media S.A. 866758 10.3389/fgene.2022.866758 Genetics Original Research 2b-RAD Genotyping of the Seagrass Cymodocea nodosa Along a Latitudinal Cline Identifies Candidate Genes for Environmental Adaptation Ruocco et al. 2b-RAD Genotyping in Cymodocea nodosa Ruocco Miriam 1 * Jahnke Marlene 2 * Silva João 3 Procaccini Gabriele 1 Dattolo Emanuela 1 1 Stazione Zoologica Anton Dohrn, Naples, Italy 2 Department of Marine Sciences, Tjärnö Marine Laboratory, University of Gothenburg, Gothenburg, Sweden 3 Centre of Marine Sciences, University of Algarve, Faro, Portugal

Edited by: Serena Aceto, University of Naples Federico II, Italy

Reviewed by: Lauric Reynes, UMI3614 Biologie évolutive et écologie des algues (EBEA), France

Aline Finger, Royal Botanic Garden Edinburgh, United Kingdom

 *Correspondence: Miriam Ruocco, miriam.ruocco@szn.it; Marlene Jahnke, marlene.jahnke@gu.se

This article was submitted to Evolutionary and Population Genetics, a section of the journal Frontiers in Genetics

 16 05 2022 2022 13 866758 31 01 2022 19 04 2022 Copyright © 2022 Ruocco, Jahnke, Silva, Procaccini and Dattolo. 2022 Ruocco, Jahnke, Silva, Procaccini and Dattolo

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Plant populations distributed along broad latitudinal gradients often show patterns of clinal variation in genotype and phenotype. Differences in photoperiod and temperature cues across latitudes influence major phenological events, such as timing of flowering or seed dormancy. Here, we used an array of 4,941 SNPs derived from 2b-RAD genotyping to characterize population differentiation and levels of genetic and genotypic diversity of three populations of the seagrass Cymodocea nodosa along a latitudinal gradient extending across the Atlantic-Mediterranean boundary (i.e., Gran Canaria—Canary Islands, Faro—Portugal, and Ebro Delta—Spain). Our main goal was to search for potential outlier loci that could underlie adaptive differentiation of populations across the latitudinal distribution of the species. We hypothesized that such polymorphisms could be related to variation in photoperiod-temperature regime occurring across latitudes. The three populations were clearly differentiated and exhibited diverse levels of clonality and genetic diversity. Cymodocea nodosa from the Mediterranean displayed the highest genotypic richness, while the Portuguese population had the highest clonality values. Gran Canaria exhibited the lowest genetic diversity (as observed heterozygosity). Nine SNPs were reliably identified as outliers across the three sites by two different methods (i.e., BayeScan and pcadapt), and three SNPs could be associated to specific protein-coding genes by screening available C. nodosa transcriptomes. Two SNPs-carrying contigs encoded for transcription factors, while the other one encoded for an enzyme specifically involved in the regulation of flowering time, namely Lysine-specific histone demethylase 1 homolog 2. When analyzing biological processes enriched within the whole dataset of outlier SNPs identified by at least one method, “regulation of transcription” and “signalling” were among the most represented. Our results highlight the fundamental importance signal integration and gene-regulatory networks, as well as epigenetic regulation via DNA (de)methylation, could have for enabling adaptation of seagrass populations along environmental gradients.

 seagrass SNPs 2b-RAD outlier loci latitude flowering Introduction

Environmental heterogeneity in space and time can impose a strong selective pressure driving adaptive divergence of populations (Linhart and Grant, 1996; Bergland et al., 2016; Urban et al., 2020). Quantifying the extent of this differentiation and identifying loci underlying such divergence is a major aim of evolutionary and ecological genetics (Feder and Mitchell-Olds, 2003; Savolainen et al., 2013) and seascape genomics (Selkoe et al., 2016). Environmental gradients, where environmental factors change along a geographic scale, offer a great opportunity for understanding patterns and processes responsible for phenotypic changes of populations (Endler, 2020). There are examples from a variety of species, where phenotypes tend to change in predictable ways across large-scale gradients such as latitude, altitude, or water depth (Endler, 2020). These geographical patterns can reflect genetic variation or phenotypic plasticity and eventually represent adaptive variation in response to selection gradients (Conover et al., 2009; Sanford and Kelly, 2011; Schneider and Meyer, 2017).

Variation in specific environmental variables (e.g., temperature or photoperiod) across geographical clines is directly related to aspects of phenotypic and genetic divergence among populations in land plants (Linhart and Grant, 1996; Hut and Beersma, 2011; Hut et al., 2013). For instance, climate is one of the most important drivers of adaptive phenotypic traits in forest trees (Richardson et al., 2009; Eckert et al., 2010). Fitness-related traits such as survival, growth and biomass partitioning have been shown to vary along temperature gradients associated with altitudinal and latitudinal clines or according to different precipitation and aridity regimes (Aitken and Hannerz, 2001; St Clair et al., 2005). Similarly, other traits associated with major life history events, e.g., seed dormancy and flowering time, were found to follow a clear latitudinal pattern in Arabidopsis thaliana accessions collected throughout the European range of the species (Debieu et al., 2013). In particular, timing of flowering across accessions was linked to molecular polymorphisms in key regulatory genes that control this trait, such as FRIGIDA and FLOWERING LOCUS C (Le Corre et al., 2002; Stinchcombe et al., 2004; Shindo et al., 2005).

In marine plants, these studies are still hampered due to the absence of genomic resources for most species (with few exceptions e.g., Zostera marina) and the general unfeasibility to perform molecular genetic studies across multiple generations. The reproduction of seagrasses under controlled conditions is indeed still challenging, and life cycles are often too long, which prevent identifying heritable variations (Hu et al., 2020; Pazzaglia et al., 2021). Seagrasses are the only group of flowering plants that has returned to the sea, (re)adapting their physiology and morphology to a completely submerged lifestyle (Wissler et al., 2011; Golicz et al., 2015; Lee et al., 2016; Olsen et al., 2016; Lee et al., 2018). Despite their low taxonomical diversity, they successfully colonized most coastal shores worldwide (Spalding, 2003; Larkum et al., 2006), providing fundamental ecosystem services (Ruiz-Frau et al., 2017) and contributing to climate change mitigation (Marbà et al., 2015; Stankovic et al., 2021). The current distribution of seagrass species along geographical and depth gradients is influenced by their tolerance window for environmental drivers, such as light, temperature and salinity (Short et al., 2007).

To date, population genetic studies in seagrasses have mostly relied on a limited number of nuclear DNA markers (e.g., simple sequence repeats, SSRs). Such studies have aided in resolving main geographic differentiation and structure of seagrass populations (e.g., Alberto et al., 2006; Arnaud-Haond et al., 2007; Alberto et al., 2008; Serra et al., 2010; Hernawan et al., 2017; Bijak et al., 2018; Jackson et al., 2021). Indices of genetic diversity based on SSRs have also been correlated with local environmental disturbances (Jahnke et al., 2015b), fundamental plant traits such as flowering synchronization (Jahnke et al., 2015a) and ecosystem functioning and resilience (Reusch et al., 2005; Hughes and Stachowicz, 2011). In a few studies, molecular polymorphisms have been functionally characterized through genome scan approaches and related to contrasting habitat types, depth and latitudinal gradients (Oetjen and Reusch, 2007; Oetjen et al., 2010; Jahnke et al., 2019). Yet, gene expression studies employing common-garden or reciprocal transplantation approaches have collectively suggested local adaptation of seagrass populations to thermal or light gradients (Franssen et al., 2011; Winters et al., 2011; Franssen et al., 2014; Jueterbock et al., 2016; Dattolo et al., 2017), even if the heritability of the observed phenotypic differences remain elusive.

Cymodocea nodosa (Ucria) Ascherson is a marine dioecious angiosperm present throughout the Mediterranean Sea and the adjoining Atlantic coasts and the dominant species structuring subtidal seagrass ecosystems along North West Africa and South West Europe (Green and Short, 2003; Alberto et al., 2008). Recent species distribution models predict the most relevant environmental variables defining its distribution to be sea surface temperature (SST) and salinity (Chefaoui et al., 2016). Large inter-regional and local-scale variations in abundance and structure (i.e., morphology and biomass allocation), as well as the extant of sexual reproduction (as seed production) in C. nodosa have been described and correlated with different environmental conditions [i.e., seasonal patterns of Photosynthetically Active Radiation (PAR) and SST], and evolutionary contexts (i.e., genetic diversity) (Máñez-Crespo et al., 2020). Populations inhabiting distinct biogeographical regions have also been shown to possess a differential resilience, performance and recovery capacity under local perturbations (e.g., shading) (Tuya et al., 2019; Tuya et al., 2021). The higher sensitivity of certain meadows to disturbance, as those at the peripheral distribution of the species in the Canary Islands, has been attributed to their genetic isolation and low genetic diversity (Tuya et al., 2019; Tuya et al., 2021). An important geographic heterogeneity was also observed in C. nodosa plants from contrasting thermal environments (warm-adapted vs. cold-adapted plants) in response to heatwaves, in terms of diversion of carbon reserves and biomass allocation (Marín-Guirao et al., 2018). However, molecular studies addressing how environmental gradients could affect local adaptation of natural populations through adaptive genetic variation at specific loci are currently missing in this species.

Restriction-site-associated DNA sequencing (RAD-Seq) techniques represent a family of cost-effective techniques compared to e.g., whole-genome sequencing that can be employed in non-model species without a reference genome and guarantee high-resolution population genomics data for demographic analyses (Andrews et al., 2016; Shafer et al., 2017). In addition, these approaches offer better opportunities in respect to those based on other markers (e.g., SSRs) to identify loci with a putative signal of selection from the background of neutral variation, and to test their functional importance by associating nucleotide variation in these genes with phenotypic variation in adaptive traits in natural populations (Storz, 2005; Biswas and Akey, 2006). RAD-seq approaches have been applied in marine engineering species (e.g., kelp, Guzinski et al., 2020; Fucales, Reynes et al., 2021a or Mediterranean corals, Pratlong et al., 2021) shedding first light on the detection of candidate SNPs for local adaptation of populations.

In the present study, we used 2b-RAD genotyping (Wang et al., 2012) to characterize population differentiation and levels of genetic diversity among three populations of C. nodosa distributed along the Atlantic-Mediterranean transition region (i.e., Gran Canaria and Faro in the Atlantic and Ebro Delta in the Mediterranean; Figure 1A). Our main aim was to search for outlier loci that could underlie adaptive differentiation of populations across the latitudinal distribution of the species. The underlying hypothesis was that such genetic polymorphisms could be associated to local adaptation of populations to varying photoperiod-temperature regimes across latitudes. We also used SNP markers for assessing the clonality level of populations and compared these data, together with genetic diversity and differentiation, with estimates based on 7 SSR markers previously developed for the species.

Sampling scheme and environmental characteristics of the selected sites. (A) Sampling locations of C. nodosa: GC—Las Palmas de Gran Canaria (Lat: 27°N, Canary Islands); FA—Faro (Lat: 36°N, Portugal); EB—Ebro Delta (Lat: 40°N, Spain). (B) Annual patterns of photoperiod for the different sampling sites. Daylight hours were derived from annual data and averaged for every 15 days (https://www.timeanddate.com). (C) Long term monthly means of Sea Surface Temperature (SST) registered on the Atlantic and Mediterranean waters from NOAA derived from data for years 1971–2000 (https://psl.noaa.gov/data/gridded/data.noaa.ersst.v4.html).

 Materials and Methods Study Sites and Sample Collection

Individual shoots of Cymodocea nodosa were collected at three locations across the Atlantic-Mediterranean boundary from shallow-water meadows (1–4 m depth), thus encompassing most of the latitudinal distribution range of the species (Alberto et al., 2006; Cunha and Araújo, 2009). Collection sites were selected based on their geographic position along a latitudinal gradient and for sampling feasibility. Two sampling stations were located in the north-eastern Atlantic Ocean, i.e., Las Palmas de Gran Canaria (El pajar) (GC)—Canary Islands, Spain (27°45′14″N, 15°40′18″W) and in the Ria Formosa lagoon (Ilha da Culatra) (FA)—Faro, Portugal (36°59′23″N, 7°50′01″W), while the latter was located in the western Mediterranean Sea, i.e., Ebro Delta (EB)—Spain (40°34′48″N, 0°35′45″E) (Figure 1A). The Canary Islands represent the western distributional limit of the species, and are close to its southern distributional limit in Banc d’Arguin in Mauritania (Green and Short, 2003). The Ebro population represents the northernmost sampling site in our study, although the species extends up to the north Adriatic Sea (Venice lagoon, 45.60°N) (Rismondo et al., 2003; Sfriso and Facca, 2007).

In Gran Canaria, daylight hours oscillate between 10 h d−1–14 h d−1 throughout the year (Figure 1B). In the other two sites, there are greater fluctuations of daylight during seasons: the minimum is ca. 9 h d−1 in winter, while the maximum is 15 h d−1 during the summer solstice (Figure 1B). In Gran Canaria, the average seawater temperature (SST) during the year is about 19°C (Figure 1C). Faro, the coldest of our sites, displays an average SST during the year of around 16.5°C while Ebro Delta in the Mediterranean, despite its northern latitude, has an average seawater temperature of 18°C (Figure 1C). Twenty C. nodosa shoots were collected by snorkeling at each site at a reciprocal distance of around 3–5 m to minimize the risk of sampling within the same clonal patch (Arnaud-Haond et al., 2007; Serra et al., 2010; Jahnke et al., 2017; Jahnke et al., 2019). This number of samples can be considered representative of the meadows and it is comparable to what generally used in 2b-RAD studies (e.g., De Wit et al., 2020; Jahnke et al., 2022). Subsequently, leaf material (about 5–7 cm) was carefully cleaned from epiphytes and dried in silica gel prior to DNA isolation.

 2b-RAD Genotyping

Genomic DNA was extracted using the NucleoSpin® 96 Plant II kit (Macherey-Nagel) following the manufacturer’s instructions. DNA quality and quantity were checked through 1% agarose gel electrophoresis, and the QubitTM dsDNA BR assay kit (Thermo Fisher Scientific). 2b-RAD libraries were prepared following a modification of the protocol described by Wang et al., 2012 and available at https://github.com/z0on/2bRAD_denovo. Briefly, genomic DNA (∼100 ng) was digested using the type 2b restriction enzyme BcgI to produce uniform fragments of 36-bp, to which adaptors were ligated on the cohesive ends. The fragments were then amplified with barcoded adaptors and purification of the target bands was carried out by 2% agarose gel-electrophoresis. Subsequently, gel fragments were cleaned using a MinElute Gel Extraction Cleaning Kit (Qiagen) and pooled equimolarly into a single pool. Four individuals from each locality were used as technical replicates (i.e., replicated library preparation and sequencing). These technical replicates were used in the analysis step to assess the overall accuracy of genotyping, set quality filtering criteria, and quantify error rates between samples. In total, 70 C. nodosa samples (20 individuals +4 technical replicates for EB and FA; 18 individuals +4 technical replicates for GC) were successfully sequenced on two lanes of an Illumina NovaSeq platform, generating 50-bp paired-end reads, at the Science for Life Laboratory (SciLifeLab) Genomics, SNP&SEQ Technology Platform in Uppsala University, Sweden. Only reads 1 (“forward reads”) were used for this analysis.

Bioinformatic analyses were performed using the computer cluster of the Bioinforma Service of Stazione Zoologica Anton Dohrn (SZN), Naples (Italy). As no genome sequence is currently available for C. nodosa, the analysis followed a modified de novo pipeline available at https://github.com/z0on/2bRAD_denovo. Reads were first ‘demultiplexed’ based on barcodes and then adaptors were removed, then a quality filtering was performed using the Fastx-toolkit (Gordon and Hannon, 2010). Only reads containing 100% bases with a PHRED quality score of at least 20 were retained for downstream analysis. After trimming and quality filtering, we obtained a total of 82,883,209 reads. Individually trimmed fastq files were then merged to collect tags found in at least two individuals with a minimum depth of five for genotyping. Reads that had more than seven observations without reverse-complement were removed. Subsequently, tags were clustered with CD-hit (Li and Godzik, 2006) allowing for up to three mismatches, followed by the creation of a “reads-derived reference genome” based on 30 fake chromosomes, on which individual trimmed fastq files were mapped back using Bowtie 2 aligner (Langmead et al., 2009).

SNP-calling was performed using The Genome Analysis Toolkit (GATK) version 3.8 (McKenna et al., 2010). A first round of putative variants was generated using GATK’s UnifiedGenotyper, followed by base quality score recalibration (BQSR/BaseRecalibrator and PrintReads) based on a high confidence (>75th quality percentile) SNPs dataset. The realigned and recalibrated reads were then used to perform a second round of UnifiedGenotyper. We used the variant quality score recalibration (VQSR) step to generate an adaptive error model using the SNPs that were consistently genotyped across the technical replicates. A further filtration step was performed with vcftools (Danecek et al., 2011) to remove poorly-genotyped samples and to select only biallelic loci genotyped in at least 90% of individuals, and with a maximum heterozygosity of 50%. Harsh genotyping rate cut-off is recommended for best quality and to avoid RAD loci affected by null alleles because of mutations in the restriction site.

 Population Genomic Analyses

Our final dataset was thinned in order to keep one SNP per RAD fragment with maximal allele frequency (script thinner. pl with criterion = maxAF; https://github.com/z0on/2bRAD_denovo). Technical replicates were discarded from further analyses after examination for concordance in genotype estimates. Individual genetic variation was explored by a Principal Component Analysis (PCA) using the R package adegenet v. 2.1.3 (Jombart, 2008) and by using ADMIXTURE 1.3.0 (Alexander and Lange, 2011). To choose the best estimate of number of clusters (K), we used the ADMIXTURE’s cross-validation procedure with default settings. The hypothetical number of K was set from 1 to 15. PCA and ADMIXTURE analyses were re-calculated after removal of outlier loci confirmed by both employed approaches in the dataset (see 2.4), in order to determine whether outliers had disproportionate or distortive effects on the genetic structure analysis. As the analysis with all loci and only neutral loci showed similar patterns (data not shown), we kept all loci in the analyses. Significance levels of genic differentiation for each population pair (exact G test) across all loci was calculated with GENEPOP 4.7.5 (Rousset, 2008) with default Markov chain parameters (Dememorisation = 10000, Batches = 100, Iterations per batch = 5000). Pair-wise Weir and Cockerham mean and weighted FST estimates between C. nodosa populations were calculated with vcftools (Danecek et al., 2011). Observed heterozygosity (Ho) as well as FIS across all loci for each population were calculated using the R package diveRsity 1.9.90 (Keenan et al., 2013). The number of distinct Multi Locus Lineages (MLLs) for each population was calculated using the R package poppr (Kamvar et al., 2014; Kamvar et al., 2015). The genetic distance limit for setting delimitation of clones (Hamming distance >0.073) was determined based on the maximum genetic distance detected between technical replicates (i.e., “confirmed clones”). The genetic distance tree of individual samples obtained with poppr is depicted in Supplementary Figure S1. After clone delimitation, FST, FIS, and heterozygosity (Ho) were recalculated keeping only distinct MLLs per population.

Genetic diversity and differentiation, as well as clonal diversity results were compared with those obtained for the same populations using 7 highly polymorphic SSR markers (microsatellites) developed for C. nodosa (Ruggiero et al., 2004; Supplementary Table S1) and used in previous population genetics studies (Ruggiero et al., 2005a and b; Tuya et al., 2019). Only a subset of samples could be genotyped with both techniques (12 individuals for GC; 18 for FA and 12 for EB). Multiplex PCR amplifications were conducted in 25 μL reaction volumes containing 12.5 μL QIAGEN Multiplex PCR Master Mix (QIAGEN), and 0.5 μL DNA (6–10 ng). Thermal cycling consisted of 95°C for 15′, 35 cycles of 94°C for 60″, 58°C for 90″, and 72°C for 90″, followed by 72°C for 30'. PCR products were analyzed on an Automated Capillary Electrophoresis Sequencer 3730 DNA Analyzer (Applied Biosystems). Linkage disequilibrium (LD) and deviations from Hardy-Weinberg equilibrium (HWE) at each locus and across all loci in each population were tested with Genepop 4.7.5 (Rousset, 2008), using 1000 dememorisations, 100 batches and 1000 iterations per batch. LD was not detected among loci, indicating they behaved independently. The number of MLGs was identified using the GIMLET software (Valière, 2002). Clonal diversity was calculated as the R ratio: R = G-1/N-1, where G is the number of genotypes and N is the number of samples (Dorken and Eckert, 2001). Observed heterozygosity and inbreeding coefficient, as well as pairwise FST between populations were obtained with GenAlEx 6.503 (Peakall and Smouse, 2012).

 Outlier Loci Identification and Functional Annotation

To identify putative loci under selection, we used two statistical tools employing two different approaches based on FST and PCA, respectively, i.e., BayeScan v. 2.1 (Foll and Gaggiotti, 2008) and pcadapt v. 4 (Luu et al., 2017; Privé et al., 2020). BayeScan estimates FST for each SNP locus to perform a genomic scan for outlier FST values through a Bayesian method. It was used with prior odds set to 100, using thresholds of q = 0.3 and posterior probability P > 0.7 (Foll and Gaggiotti, 2008). The pcadapt method tests how much each variant is associated with population structure, assuming that outlier variants are indicative of local adaptation. Based on the scree plot (Jackson, 1993), we set 2 as the best number of principal components, and used the Manhattan plot and Histogram of p-values to graphically examine the presence of outliers. Significant outliers were then determined using the Bonferroni correction method for multiple comparisons with the R function p. adjust and α = 0.001. The outputs from the two methods were compared for overlap and shared loci were considered as the “best outliers”. Allele frequency of shared outlier loci was determined with GENEPOP 4.7.5 (Rousset, 2008). Outlier identification was repeated after removal of clone samples from the dataset (i.e., keeping only distinct MLLs for each population).

To determine if shared and non-shared outlier loci could be included in a potential coding sequence, chromosome regions of a length of 80 bp, corresponding to two adjacent 2b-RAD tags around each SNP of interest, were mapped against three available C. nodosa transcriptomes [Ruocco et al., 2017, Cymodocea nodosa Cn_1 Transcriptome (Project ID: 1264710) by Joint Genome Institute (JGI) and Dattolo et al. (unpublished)] by using the blastn algorithm (Altschul et al., 1997). Unambiguous positive-matches were determined if a percentage of identity ≥90% was found in stretches of sequences of at least 25 bp around the SNP position. This criterion for mapping was used to filter out possible errors that could result from non-specific similarities among homologous sites. Subsequently, the SNPs-carrying contigs were functionally annotated using the online Mercator-MapMan4 annotation tool (Schwacke et al., 2019).

 Results Accuracy of 2b-RAD Genotyping

The 2b-RAD sequencing generated an average of 1,184,045.85 reads ± SD 425,645.65 reads per sample (Supplementary Table S2). Using a set of 2,100 SNPs that were reproducible across technical replicates as a training set for the VQSR, we estimated the (true) transition/transversion ratio as Ti/Tv = 1.61. When applying this recalibration we chose a tranche with 99% truth sensitivity as the cut-off to call all SNPs from the overall dataset. After filtering out seven poorly genotyped samples as well as highly heterozygous sites (possible lumped paralogous), we obtained a genotype dataset of 7,562 SNPs for 63 individuals (Supplementary Table S3). Overall, transversions were more frequent than transitions and composed 56% of the identified changes (Figure 2). The most frequent changes were G ↔ T (1857) and A ↔ G (1686) (Figure 2; Supplementary Table S4). Genotyping correspondences among technical replicated individuals was 89 ± 6% (Supplementary Table S5). For population differentiation analysis, our final dataset was further thinned to remove technical replicates, two samples that were clear outliers in an exploratory PCA and to only keep one SNP per RAD fragment with maximum allele frequency, resulting in 4,941 loci for 49 individuals.

Distribution of transitions and transversions in the SNPs dataset. Percentages (%) of transitions and transversions of the whole set of 7,562 SNPs identified for 63 individuals is depicted using an alluvial diagram.

 Population Differentiation, Genetic and Genotypic Diversity

Principal Component Analysis (PC1 = 35.2% and PC2 = 25.6% total variance; Figure 3A) showed a strong genetic differentiation among the three sampling sites along the latitudinal gradient and this was further confirmed by results from ADMIXTURE analysis (Figure 3B). In this analysis, K = 4 was identified as the “optimal K” with the lowest cross-validation error of 0.239 (Supplementary Table S6).

Population structure analysis and clonal diversity. (A) Principal Component Analysis (PCA) displaying the three C. nodosa populations collected along a latitudinal gradient of distribution (i.e., Gran Canaria, Faro and Ebro). (B) Admixture plot of the C. nodosa populations at K = 4. (C) Number of distinct Multi Locus Lineages (MLLs) identified at each sampling site. (D) Weir and Cockerham FST weighted estimate between populations.

The three sampled populations had different levels of clonality. The C. nodosa population from Ebro Delta exhibited the highest genotypic richness (i.e., 11 distinct MLLs over 15 individuals, R = 0.71), in respect to Faro and Gran Canaria populations [i.e., only 2 (R = 0.05) and 4 (R = 0.23) distinct MLLs, respectively] (Figure 3C ,Table 1). In particular, in the C. nodosa population from Ria Formosa lagoon (Faro), one MLL largely dominated over the other (Figure 3C). To support clonal determination with SNP markers, we compared these data with those obtained for a subset of samples of the same populations using previously developed polymorphic SSR markers (Supplementary Table S7). Based on microsatellites, the C. nodosa population from Ebro Delta exhibited the highest genotypic richness with 11 MLGs over 12 samples (R = 0.91), while in Gran Canaria 6 distinct MLGs were identified (R = 0.45). In Faro, a single MLG was identified with microsatellite markers (Supplementary Table S8).

Information on sampling sites, sample details and genetic indices obtained with 2b-RAD genotyping. For each population, the following information are shown: acronym, geographic coordinates (latitude and longitude), number of sequenced individuals N and number of additional technical replicates in parenthesis, number of individuals retained after filtering procedures and removal of technical replicates, total number of distinct Multi Locus Lineages (MLLs), R value (G-1/N-1), observed heterozygosity (Ho) and inbreeding coefficient (FIS).

Sampling site Acronym Coordinates Total N N after filtering MLLs R Ho FIS
Las Palmas de Gran Canaria—Canary Islands GC 27°45′14″N, 15°40′18″W 18 (4) 14 4 0.23 0.13 −0.20
Ria Formosa lagoon—Faro, Portugal FA 36°59′23″N, 7°50′01″W 20 (4) 20 2 0.05 0.24 −0.52
Ebro Delta, Spain EB 40°34′48″N, 0°35′45″E 20 (4) 15 11 0.71 0.21 −0.13
Total — — 58 (12) 49 17 — — —

Genetic differentiation (based on allele frequencies across all loci) was significant for each population pair comparison (see exact G test in Supplementary Table S9). Global weighted pairwise FST values ranged from 0.31 to 0.43, with Faro and Gran Canaria exhibiting the highest differentiation (Figure 3D). Overall, inbreeding coefficient values (FIS) were negative for all populations, with the lowest value observed in Faro (Table 1). Mean observed heterozygosity (Ho) was highest in Faro, intermediate in Ebro and lowest in Gran Canaria (Table 1). Values of genetic differentiation (exact G test), FST, FIS, and observed heterozygosity (Ho) indicated above have been obtained following clones’ removal from the dataset. Genetic diversity and differentiation analyses repeated using SSR markers (Supplementary Table S8) gave comparable results across populations.

 Candidate Loci for Environmental Adaptation

BayeScan identified 9 FST outlier SNPs across all samples at q<0.3 and P >0.7 (Supplementary Figure S2A; Supplementary Table S10), while pcadapt identified 188 outliers based on PCA with a Bonferroni corrected p-value ≤0.001 (Supplementary Figure S2B,C; Supplementary Table S10). All 9 outliers identified by BayeScan were included among those identified by pcadapt (9 shared outliers; Supplementary Figure S2D). Interestingly, the frequency of the alternative allele for the 9-shared outliers was maximum in individuals from the Ebro population (Supplementary Table S11). The functional annotation of the 9-shared outlier SNPs is provided in Table 2, while full annotation of all outliers identified with one or both methods is provided in Supplementary Table S12.

Functional annotations of the 9 significant outlier SNPs shared between BayeScan and pcadapt. SNPs ID; E-value; identity with %; bp of the query match; contig name (representing the best hit of the blast analysis against available transcriptomes b ); name, description and related GO of the best hit of Arabidopsis thaliana according to TAIR (https://www.arabidopsis.org/index.jsp); bin codes of MapMan classification; and mutation type of the SNP (NS, nonsynonymous; SYN, synonymous), are outlined. Full annotation of all outlier SNPs obtained with one or both methods can be retrieved from Supplementary Table S12.

SNP ID E-value Identity (%) Query match Best hit C. nodosa Best hit A. thaliana Description TAIR GO biological process TAIR MapMan BIN Mutation
chr1_673698 a — — — No hits found — — — — —
chr1_421195 a 3E-10 36/37 (97) (11:47) TRINITY_DN164033_c0_g1_i1 AT3G13682 Lysine-specific histone demethylase 1 homolog 2. Involved in H3K4 methylation of target genes including the flowering loci FLC and FWA Histone H3-K4 methylation 26.7 NS: Glu/Gly
chr1_555347 a 3E-10 36/37 (97) (31:67) TRINITY_DN7669_c0_g1_i6 AT5G64060 NAC domain containing protein 103 Regulation of transcription, DNA-templated 27.3.27 NS: Ala/Thr
chr1_549579 a — — — No hits found — — — — —
chr2_19511 a 8E-11 34/34 (100) (7:40) TRINITY_DN5447_c0_g1_i3 — — Not assigned, unknown — —
chr1_379730 a — — — No hits found — — - — —
chr2_1253 a 0.005 27/29 (93) (20:48) TRINITY_DN226_c2_g1_i1 — — Not assigned, unknown 35.2 —
chr1_645789 a 1E-9 35/36 (97) (21:56) TRINITY_DN7020_c1_g2_i3 AT5G04840 bZIP protein Positive regulation of transcription, DNA-templated, regulation of transcription, DNA-templated 27.3.35 NS: Ala/Val
chr1_656786 — — — No hits found — — — — —

Outlier SNPs, retained after clone removal by at least one method.

Contig names assigned based on the Cymodocea nodosa transcriptome by JGI.

Five of the 9 shared outlier SNPs could be associated to specific contigs using available C. nodosa transcriptomes (Table 2). Two of these SNPs-carrying contigs (chr1_555347 and chr1_645789) encoded for proteins with a function in transcription regulation (“Regulation of transcription, DNA-templated”). One SNP (chr1_421195) was specifically associated to a gene involved in chromatin remodelling (“Histone H3-K4 methylation”) and photoperiodism/regulation of flowering time, namely Lysine-specific histone demethylase 1 homolog 2 (Martignago et al., 2019; Table 2). All three annotated SNPs involved non-synonymous aminoacid substitutions (Table 2). The other two SNPs associated with specific contigs had unknown functions (Table 2). When the analysis of the outlier SNPs was repeated keeping only distinct MLLs for each population, 8 over the 9 shared-outliers were significantly retained by at least one method (pcadapt at P ≤ 0.01) (Table 2). BayeScan identified the same loci as “top outliers” but with a lower (ns) posterior probability P >0.3.

Overall, 83 from the total of 188 outlier SNPs identified via Bayescan and/or pcadapt could be included in 14 different MapMan BINs (Figure 4; Supplementary Table S12). Among those, 27% SNPs-carrying contigs were associated to signalling (14%) and RNA metabolism (13%), while 54% of SNPs were not assigned (Figure 4). Within the “RNA” BIN (n. 27), 7 contigs encoded for transcription factors of various families (e.g., GARP, NAC, WRKY, bZIP) (Supplementary Table S12). In the “signalling” BIN (n. 30), 9 contigs were associated to receptor kinases (e.g., DUF 26), which play a critical role in plant response to stimuli and defence (Dievart et al., 2020), while 2 were wall associated kinases, which are involved in sensory and signal transduction pathways between the inner and outer surroundings of cell walls (Kohorn, 2001; Supplementary Table S12).

Functional annotation of 83 outlier loci displaying a positive match with transcribed regions. Functional annotation of SNPs with a positive blast hit against C. nodosa transcriptomes is depicted using a sunburst diagram. Each sequence has been assigned to a MapMan BIN category and sub-category based on its biological role or enzymatic activity. The definition of the BINs is included in the figure legend. A detailed description of BIN sub-categories can be retrieved from Supplementary Table S12.

 Discussion

Populations distributed over broad latitudinal gradients often show patterns of clinal variation in phenotype and/or genotype (Hut and Beersma, 2011; De Frenne et al., 2013; Bergland et al., 2016; Machado et al., 2016). Through a genome-wide SNPs analysis based on 2b-RAD genotyping, we here demonstrated the presence of differentiated polymorphisms in three populations of the seagrass C. nodosa along a latitudinal gradient encompassing the Atlantic-Mediterranean transition region. Confirmed by two differentiation-based outlier tests, we identified 9 outlier SNPs, three of which could be reliably associated to genes involved in functions relevant to environmental adaptation. Notably, one SNP-carrying contig encoded for a protein with a key role in chromatin remodelling via histone methylation, and specifically linked to the regulation of photoperiodism and flowering time in land plants (He, 2009; Zhou et al., 2018; Martignago et al., 2019). This supports our initial hypothesis that specific genetic polymorphisms linked to important traits (e.g., flowering) could be present in C. nodosa inhabiting different latitudes possibly underlying local adaptation of populations to different photoperiod-temperature cues or other environmental parameters.

Our population genomic data, based on both SSRs and SNPs markers, showed a clear genetic differentiation among the three populations analysed. Besides being genetically differentiated, they also exhibited variable levels of genetic and genotypic diversity. Previous studies on the genetic structure of C. nodosa along its distributional range based on SSRs, revealed a strong genetic discontinuity between Atlantic and Mediterranean regions and a certain divergence between western and eastern Mediterranean (Alberto et al., 2008). Yet, a great heterogeneity was found across populations in terms of clonal and genetic diversity (Ruggiero, 2004; Alberto et al., 2006; Alberto et al., 2008; Máñez-Crespo et al., 2020). A trend for decreasing allelic richness from the eastern Mediterranean towards the Atlantic was also advocated to explain the low ability of Atlantic populations to resist disturbances (Tuya et al., 2021). Our results are in line with these findings, as the Gran Canaria population (GC), exhibited the lowest genetic diversity (as observed heterozygosity, Ho) both with SSRs and SNPs markers. Genotypic richness (R) has been found to be highly variable across C. nodosa populations, even at small geographical scales (Alberto et al., 2006, 2008; Tuya et al., 2019). For instance, in the Canary Islands, the average R value previously calculated with available microsatellite markers was around 0.6, while it could vary from 0.3 to 0.9 for different meadows (Alberto et al., 2006; Máñez-Crespo et al., 2020). Genotypic richness obtained here for the Gran Canaria population with 7 SSR markers was within this range (0.4), while genotypic richness assessed with the SNP-dataset was lower (0.2). The C. nodosa population of the Ria Formosa lagoon (FA) has been historically considered monoclonal. Alberto et al., 2008 using SSR markers found only 5 genets over 220 ramets. Here using SSRs, we identified one single MLG over 18 samples analyzed, while using SNPs, 2 distinct MLLs could be detected. In general, we found a non-perfect match in the clonal discrimination power of SSRs and SNPs. SNPs have generally slower average mutation rates than SSRs; hence, each SNP is typically less informative (Allendorf, 2017). However, this is normally compensated by the much higher number of SNP markers, such that RAD-Seq can perform as well as, or better than SSRs in detecting population structure/divergence (Lemopoulos et al., 2019; Camacho-Sanchez et al., 2020; Sunde et al., 2020). However, SNP-datasets based on reduced representation libraries with moderate sequencing depth are also characterized by a substantial degree of missing data, error in SNP calling due to sequencing errors, lack of read depth or other sources of spurious allele calls (Mastretta-Yanes et al., 2015), which makes the designations of clones challenging (Kamvar et al., 2015). The challenge lies in particular in deciding on a threshold between designating a clone vs. a unique genotype. Several possible thresholds have been discussed (Kamvar et al., 2015), and we decided to use the difference between technical replicates as threshold for setting the limit we believe to be sequencing errors. While details differ slightly between SSRs and SNPs, general patterns are similar, and are also confirmed by other genetic parameters, particularly the inbreeding index FIS that is influenced by clonality. All three populations exhibited strongly negative FIS values based on the SNP genotyping, but the most clonal population (Faro) exhibited the lowest FIS (-0.51). Negative FIS has been shown to be strongly linked to the degree of clonality in facultative sexually reproducing species driven by genetic drift, both in theoretical (Prugnolle and De Meeûs, 2008) and population genetic assessments (Arnaud-Haond et al., 2020; Reynes et al., 2021b). Moreover, FIS has been shown to be not as sampling-sensitive as genotypic richness (R) and has been suggested as the measure of choice for assessing the importance of clonal reproduction in seagrasses (Arnaud-Haond et al., 2020).

While SNPs may have their challenges in clone detection, one of the main advantages of reduced representation sequencing is the possibility to identify locus-specific polymorphisms potentially responsible for local adaptation of populations to specific environmental conditions (Zimmerman et al., 2020). Here, we identified nine SNPs potentially involved in environmental adaptation of C. nodosa populations along its latitudinal distribution (Table 2). These nine outlier loci could either be under selection themselves, or they could be linked to loci involved in local adaptation. The high levels of linkage disequilibrium generated by the lack of recombination in clonal populations could conserve such linkage between neutral and selected sites to a higher degree than in sexually reproducing species. However, the consequences of clonal reproduction on genomic architecture and the effects on the theory of local adaptation are so far little investigated. The three loci with a functional annotation encoded for proteins with a role in regulation of transcription. In all three loci, the mutation at the protein level was non-synonymous, hence providing a change in the aminoacid sequence. Interestingly, when considering the functional annotation of the whole outliers’ dataset (Supplementary Table S12), the biological processes “signalling” and “regulation of transcription” constituted a large part of those associated with our sequences. In total, nine transcription factors (TF) of different families (e.g., NAC domain-containing TFs, WRKY, bZIP, and JUMONJI), were found across the outlier-dataset. Similarly, a number of cysteine-rich receptor-like protein kinases, cell wall-associated ser/thr kinases, and inositol 1,3,4-trisphosphate 5/6-kinase family proteins, were identified. This highlights the fundamental role signal integration and gene-regulatory networks could play for enabling adaptation to environmental conditions, as they allow to coordinate extracellular signals and intracellular regulatory machinery (Lopez-Maury et al., 2008).

Notably, one of the 9-shared outlier SNPs was specifically associated with photoperiodism and regulation of flowering time (chr1_421195), which is in an ecologically important life history trait. The timing of floral transition is critical to reproductive success in plants and varies across latitudes in response to changes in photoperiod and temperature patterns (Bäurle and Dean, 2006; Wilczek et al., 2010; Brambilla and Fornara, 2013). In Arabidopsis, this trait is genetically controlled by a network of flowering genes, whose expression is regulated by various chromatin modifications, among which is a central regulator of flowering, FLOWERING LOCUS C (FLC) (He, 2009; Wilczek et al., 2010). FLC inhibits floral transition largely by reducing the expression of flowering-time integrators, such as SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1) and FLOWERING LOCUS T (FT) (He et al., 2003; He, 2009). The regulation of FLC expression is associated with various “active” chromatin modifications including histone H3 lysine-4 (H3K4) methylation and various “repressive” histone modifications, as H3K4 demethylation and H3 lysine-27 (H3K27) methylation (He, 2009). Chr1_421195 was included in a contig encoding for Lysine-specific histone demethylase 1 homolog 2 (LDL2). In Arabidopsis, this gene is involved in H3K4 (de)methylation of target genes including FLC. FLC is down-regulated by LDL2 and ldl2 mutants display increased H3K4me3 levels at FLC compared to wild type (Martignago et al., 2019). LDL2 is additionally involved in the control of H3K4 methylation state of FWA, a homeodomain-containing transcription factor which interferes with floral transition (Jiang et al., 2007) and controls primary seed dormancy related-genes (Zhao et al., 2015). A similar function in regulating the histone methylation pattern of flowering genes has been also suggested for LDL homologs in other plant species (Gu et al., 2016). The mutation we identified for this locus in C. nodosa populations was non-synonymous, leading to a substitution of a Glutamic acid (Glu) in Glycine (Gly). We cannot establish the real effect of this substitution at the protein level, but as it will modify the charge of the peptide, it could potentially result in a loss of interactions with other molecules or could induce differences in the regulation of other genes related to the flowering pattern (Caicedo et al., 2004). It is worth noticing that the frequency of the alternative alleles of the locus chr1_421195 was maximum in individuals from the Ebro population (Supplementary Table S11), whereas the reference allele was present in Faro and Gran Canaria. Although we expected this locus to diverge in populations exposed to large differences in photoperiod, this was not the case as e.g. Faro and Ebro had the strongest similarity in daylight and seasonality variations among sampling sites. However, studies in terrestrial plants have demonstrated that even small differences in the duration of the light period might have large effects on flowering time (Roden et al., 2002). In addition, the strong level of the differentiation at this locus across populations could also be related to differences in the temperature regimes rather than photoperiodic cues. In Arabidopsis, vernalization-induced flowering is mediated by epigenetic regulations of FLC (Sharma et al., 2020); while other studies indicated several genes belonging to the FLOWERING LOCUS family influence flowering time in responses to the ambient growth temperature (Posé et al., 2013; Capovilla et al., 2015). We cannot exclude that this polymorphism could also be linked to differences in other environmental settings across sites or to a differential regulation of specific biological processes among the analysed populations.

Across the whole dataset of outlier SNPs (Supplementary Table S12), other loci with a functional annotation, were involved in the regulation of flowering time and circadian rhythms. Chr1_659852 was associated to a contig encoding for PWWP domain protein 3, which is specifically involved in the regulation of FLC and flowering time. PWWP domain proteins (PDPs) function together with other components (e.g., FLOWERING LOCUS VE and MULTICOPY SUPPRESSOR OF IRA 5) to regulate the function of the PRC2 histone methyltransferase complex, thereby facilitating the maintenance of H3K27me3 on FLC (Zhou et al., 2018). In addition, other two SNPs-carrying contigs (chr1_546361 and chr1_474761) encoded for transcription factors involved in flowering regulation and response to temperature stimulus. The first encoded for the protein JMJD5/JMJ30, which contains a jumonji-C (jmjC) domain, the second for EARLY FLOWERING MYB PROTEIN (EFM). In Arabidopsis, JMJD5 has a histone demethylase activity (H3-K36 specific) and interacts with EFM to repress the floral pathway integrator FT, thus negatively regulating flowering (Yan et al., 2014). In addition, EFM participate in the flowering thermosensory pathway of Arabidopsis (Gan et al., 2014).

Taken together, our results highlight that polymorphisms at flowering-related loci and transcription factors could make an important contribution to genetic adaptation of C. nodosa across latitude, although this should be further confirmed by larger-scale studies and ad hoc common-garden or reciprocal transplantation experiments (Reusch et al., 2014; Pazzaglia et al., 2021). Flowering-regulating loci, integrating annual responses to light and temperature patterns through complex gene-regulatory networks, determine the timing of reproduction, which is crucial for fitness and survival of populations (Austen et al., 2017), besides the plasticity of single long-lasting genotypes (Pazzaglia et al., 2021). This could also contribute to differences in flowering phenology observed in this species along the depth gradient (Buia and Mazzella, 1991) or across regions (Máñez-Crespo et al., 2020). Yet, epigenetic modifications, such as histone modifications, seem to play an important role in the regulation of flowering-related gene networks, and potentially other phenology-related pathways, in C. nodosa. In support of this, a recent work has demonstrated different levels of global genomic cytosine methylation (5-mC) in C. nodosa populations from contrasting thermal environments, as well as a low level of gene body methylation (predicted by high CpGO/E ratio) in transcripts with a differential expression depending on plants’ origin (Entrambasaguas et al., 2021).

While we acknowledge that our study including only three C. nodosa populations with a high level of clonality in at least two of them, has limitations, it sheds first light on genetic polymorphisms and related biological processes that could contribute to environmental adaptation of C. nodosa and potentially other seagrass species across their biogeographic range of distribution. In addition, we are aware that, although genome scan methods have been often applied in partially clonal organisms, there is a general lack of predictive works about the power of such methods in case of very strong clonality rate, as in our case. Future seascape genomic studies based on upgraded whole-genome information and employing a larger number of populations (or species) along latitudinal clines would allow a higher power for detecting signature of selection associated to local environmental settings in seagrasses.

 Data Availability Statement

The raw Illumina sequences used in the study are available in the Sequence Read Archive (SRA) repository under NCBI Bioproject PRJNA814006. Vcf files are available as Supplementary Material (Data Sheet 2).

 Author Contributions

ED and MR conceived the ideas and designed the study. ED, MR, and JS participated in the sampling activity. ED, MR, and MJ performed all molecular and bioinformatics 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. The work has been partially supported by the project Marine Hazard, PON03PE_00203_1, Italian Ministry of Education, University and Research (MIUR). Sequencing was financially supported by Centre for Marine Evolutionary Biology (http://www.cemeb.science.gu.se). MJ was supported by the Swedish Research Council Formas (Grant No. 2020-0008). This study received Portuguese national funds from FCT - Foundation for Science and Technology through projects UIDB/04326/2020, UIDP/04326/2020 and LA/P/0101/2020.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

We deeply thank all the staff of Tjärnö Marine Laboratory (Strömstad, Sweden) for the hospitality, and for logistical and technical supports. Sequencing was performed at Science for Life Laboratory (SciLifeLab) - Genomics, SNP&SEQ Technology Platform in Uppsala University, Sweden. Part of sequence data used for outlier SNPs functional annotation were produced by the US Department of Energy Joint Genome Institute https://www.jgi.doe.gov/in collaboration with the user community. We also thank Fernando Tuya for the collection of C. nodosa samples in Las Palmas de Gran Canaria. We are grateful to the SZN-Bioinforma service and in particular to Marco Miralto for the support with computer cluster usage and software installation. We deeply thank Arturo Santaniello for his invaluable help with the 2b-RAD bioinformatic analyses.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.866758/full#supplementary-material

References Aitken S. N. Hannerz M. (2001). “Genecology and Gene Resource Management Strategies for Conifer Cold Hardiness,” in Conifer Cold Hardiness (Springer), 23–53. 10.1007/978-94-015-9650-3_2 Alberto F. Arnaud-Haond S. Duarte C. Serrão E. (2006). Genetic Diversity of a Clonal Angiosperm Near its Range Limit: the Case of Cymodocea Nodosa at the Canary Islands. Mar. Ecol. Prog. Ser. 309, 117–129. 10.3354/meps309117 Alberto F. Massa S. Manent P. Diaz-Almela E. Arnaud-Haond S. Duarte C. M. (2008). Genetic Differentiation and Secondary Contact Zone in the seagrassCymodocea Nodosaacross the Mediterranean-Atlantic Transition Region. J. Biogeogr. 35, 1279–1294. 10.1111/j.1365-2699.2007.01876.x Alexander D. H. Lange K. (2011). Enhancements to the ADMIXTURE Algorithm for Individual Ancestry Estimation. BMC Bioinforma. 12, 246. 10.1186/1471-2105-12-246 Allendorf F. W. (2017). Genetics and the Conservation of Natural Populations: Allozymes to Genomes. Wiley Online Library. Altschul S. Madden T. L. Schäffer A. A. Zhang J. Zhang Z. Miller W. (1997). Gapped BLAST and PSI-BLAST: a New Generation of Protein Database Search Programs. Nucleic acids Res. 25, 3389–3402. 10.1093/nar/25.17.3389 Andrews K. R. Good J. M. Miller M. R. Luikart G. Hohenlohe P. A. (2016). Harnessing the Power of RADseq for Ecological and Evolutionary Genomics. Nat. Rev. Genet. 17, 81–92. 10.1038/nrg.2015.28 Arnaud‐Haond S. Stoeckel S. Bailleul D. (2020). New Insights into the Population Genetics of Partially Clonal Organisms: when Seagrass Data Meet Theoretical Expectations. Mol. Ecol. 29, 3248–3260. 10.1111/mec.15532 Arnaud-Haond S. Migliaccio M. Diaz-Almela E. Teixeira S. Van De Vliet M. S. Alberto F. (2007). Vicariance Patterns in the Mediterranean Sea: East-West Cleavage and Low Dispersal in the Endemic Seagrass Posidonia Oceanica . J. Biogeogr. 34, 963–976. 10.1111/j.1365-2699.2006.01671.x Austen E. J. Rowe L. Stinchcombe J. R. Forrest J. R. K. (2017). Explaining the Apparent Paradox of Persistent Selection for Early Flowering. New Phytol. 215, 929–934. 10.1111/nph.14580 Bäurle I. Dean C. (2006). The Timing of Developmental Transitions in Plants. Cell 125, 655–664. 10.1016/j.cell.2006.05.005 Bergland A. O. Tobler R. González J. Schmidt P. Petrov D. (2016). Secondary Contact and Local Adaptation Contribute to Genome‐wide Patterns of Clinal Variation in Drosophila melanogaster . Mol. Ecol. 25, 1157–1174. 10.1111/mec.13455 Bijak A. L. van Dijk K.-j. Waycott M. (2018). Population Structure and Gene Flow of the Tropical Seagrass, Syringodium Filiforme, in the Florida Keys and Subtropical Atlantic Region. PloS one 13 (9), e0203644. 10.1371/journal.pone.0203644 Biswas S. Akey J. M. (2006). Genomic Insights into Positive Selection. TRENDS Genet. 22, 437–446. 10.1016/j.tig.2006.06.005 Brambilla V. Fornara F. (2013). Molecular Control of Flowering in Response to Day Length in Rice. J. Integr. Plant Biol. 55, 410–418. 10.1111/jipb.12033 Buia M. C. Mazzella L. (1991). Reproductive Phenology of the Mediterranean Seagrasses Posidonia Oceanica (L.) Delile, Cymodocea Nodosa (Ucria) Aschers., and Zostera Noltii Hornem. Aquat. Bot. 40 (4), 343–362. 10.1016/0304-3770(91)90080-o Caicedo A. L. Stinchcombe J. R. Olsen K. M. Schmitt J. Purugganan M. D. (2004). Epistatic Interaction between Arabidopsis FRI and FLC Flowering Time Genes Generates a Latitudinal Cline in a Life History Trait. Proc. Natl. Acad. Sci. U.S.A. 101, 15670–15675. 10.1073/pnas.0406232101 Camacho‐Sanchez M. Velo‐Antón G. Hanson J. O. Veríssimo A. Martínez‐Solano Í. Marques A. (2020). Comparative Assessment of Range‐wide Patterns of Genetic Diversity and Structure with SNPs and Microsatellites: A Case Study with Iberian Amphibians. Ecol. Evol. 10, 10353–10363. 10.1002/ece3.6670 Capovilla G. Schmid M. Posé D. (2015). Control of Flowering by Ambient Temperature. J. Exp. Bot. 66 (1), 59–69. 10.1093/jxb/eru416 Chefaoui R. M. Assis J. Duarte C. M. Serrão E. A. (2016). Large-Scale Prediction of Seagrass Distribution Integrating Landscape Metrics and Environmental Factors: The Case of Cymodocea Nodosa (Mediterranean-Atlantic). Estuaries Coasts 39, 123–137. 10.1007/s12237-015-9966-y Conover D. O. Duffy T. A. Hice L. A. (2009). The Covariance between Genetic and Environmental Influences across Ecological Gradientsfluences across Ecological Gradients. Ann. N. Y. Acad. Sci. 1168, 100–129. 10.1111/j.1749-6632.2009.04575.x Cunha A. H. Araújo A. (2009). New Distribution Limits of Seagrass Beds in West Africa. J. Biogeogr. 36, 1621–1622. 10.1111/j.1365-2699.2009.02135.x Danecek P. Auton A. Abecasis G. Albers C. A. Banks E. Depristo M. A. (2011). The Variant Call Format and VCFtools. Bioinformatics 27, 2156–2158. 10.1093/bioinformatics/btr330 Dattolo E. Marín-Guirao L. Ruiz J. M. Procaccini G. (2017). Long-term Acclimation to Reciprocal Light Conditions Suggests Depth-Related Selection in the Marine Foundation speciesPosidonia Oceanica. Ecol. Evol. 7, 1148–1164. 10.1002/ece3.2731 Debieu M. Tang C. Stich B. Sikosek T. Effgen S. Josephs E. (2013). Co-variation between Seed Dormancy, Growth Rate and Flowering Time Changes with Latitude in Arabidopsis thaliana . PloS one 8, e61075. 10.1371/journal.pone.0061075 De Wit P. Jonsson P. R. Pereyra R. T. Panova M. André C. Johannesson K. (2020). Spatial Genetic Structure in a Crustacean Herbivore Highlights the Need for Local Considerations in Baltic Sea Biodiversity Management. Evol. Appl. 13 (5), 974–990. 10.1111/eva.12914 Dievart A. Gottin C. Périn C. Ranwez V. Chantret N. (2020). Origin and Diversity of Plant Receptor-like Kinases. Annu. Rev. Plant Biol. 71, 131–156. 10.1146/annurev-arplant-073019-025927 Dorken M. E. Eckert C. G. (2001). Severely Reduced Sexual Reproduction in Northern Populations of a Clonal Plant, Decodon Verticillatus (Lythraceae). J. Ecol. 89, 339–350. 10.1046/j.1365-2745.2001.00558.x Eckert A. J. Bower A. D. González-martínez S. C. Wegrzyn J. L. Coop G. Neale D. B. (2010). Back to Nature: Ecological Genomics of Loblolly Pine (Pinus Taeda, Pinaceae). Mol. Ecol. 19, 3789–3805. 10.1111/j.1365-294x.2010.04698.x Endler J. A. (2020). Geographic Variation, Speciation and Clines.(MPB-10), Vol. 10. Princeton University Press. Entrambasaguas L. Ruocco M. Verhoeven K. J. F. Procaccini G. Marín-Guirao L. (2021). Gene Body DNA Methylation in Seagrasses: Inter- and Intraspecific Differences and Interaction with Transcriptome Plasticity under Heat Stress. Sci. Rep. 11, 14343. 10.1038/s41598-021-93606-w Feder M. E. Mitchell-Olds T. (2003). Evolutionary and Ecological Functional Genomics. Nat. Rev. Genet. 4, 649–655. 10.1038/nrg1128 Foll M. Gaggiotti O. E. (2008). A Genome-Scan Method to Identify Selected Loci Appropriate for Both Dominant and Codominant Markers: A Bayesian Perspective. Genetics 180, 977–993. 10.1534/genetics.108.092221 Franssen S. U. Gu J. Bergmann N. Winters G. Klostermeier U. C. Rosenstiel P. (2011). Transcriptomic Resilience to Global Warming in the Seagrass Zostera Marina , a Marine Foundation Species. Proc. Natl. Acad. Sci. U.S.A. 108, 19276–19281. 10.1073/pnas.1107680108 Franssen S. U. Gu J. Winters G. Huylmans A.-K. Wienpahl I. Sparwel M. (2014). Genome-wide Transcriptomic Responses of the Seagrasses Zostera Marina and Nanozostera Noltii under a Simulated Heatwave Confirm Functional Types. Mar. Genomics 15, 65–73. 10.1016/j.margen.2014.03.004 Frenne P. Graae B. J. Rodríguez‐Sánchez F. Kolb A. Chabrerie O. Decocq G. (2013). Latitudinal Gradients as Natural Laboratories to Infer Species' Responses to Temperature. J. Ecol. 101, 784–795. 10.1111/1365-2745.12074 Gan E. S. Xu Y. Wong J. Y. Goh J. G. Sun B. Wee W. Y. (2014). Jumonji Demethylases Moderate Precocious Flowering at Elevated Temperature via Regulation of FLC in Arabidopsis . Nat. Commun. 5, 5098. 10.1038/ncomms6098 Golicz A. A. Schliep M. Lee H. T. Larkum A. W. D. Dolferus R. Batley J. (2015). Genome-wide Survey of the Seagrass Zostera Muelleri Suggests Modification of the Ethylene Signalling Network. J. Exp. Bot. 66, 1489–1498. 10.1093/jxb/eru510 Gordon A. Hannon G. (2010). Fastx-toolkit. FASTQ/A Short-Reads Preprocessing Tools. Available at: http://hannonlab. cshl. edu/fastx_toolkit Green E. P. Short F. T. (2003). World Atlas of Seagrasses. University of California Press. Gu T. Han Y. Huang R. Mcavoy R. J. Li Y. (2016). Identification and Characterization of Histone Lysine Methylation Modifiers in Fragaria Vesca . Sci. Rep. 6, 23581. 10.1038/srep23581 Guzinski J. Ruggeri P. Ballenghien M. Mauger S. Jacquemin B. Jollivet C. Coudret J. Jaugeon L. Destombe C. Valero M. (2020). Seascape Genomics of the Sugar Kelp Saccharina Latissima along the North Eastern Atlantic Latitudinal Gradient. Genes 11 (12), 1503. 10.3390/genes11121503 He Y. (2009). Control of the Transition to Flowering by Chromatin Modifications. Mol. plant 2, 554–564. 10.1093/mp/ssp005 He Y. Michaels S. D. Amasino R. M. (2003). Regulation of Flowering Time by Histone Acetylation in Arabidopsis . Science 302, 1751–1754. 10.1126/science.1091109 Hernawan U. E. van Dijk K.-j. Kendrick G. A. Feng M. Biffin E. Lavery P. S. (2017). Historical Processes and Contemporary Ocean Currents Drive Genetic Structure in the seagrassThalassia Hemprichiiin the Indo-Australian Archipelago. Mol. Ecol. 26 (4), 1008–1021. 10.1111/mec.13966 Hu Z.-M. Zhong K.-L. Weinberger F. Duan D.-L. Draisma S. G. A. Serrão E. A. (2020). Linking Ecology to Genetics to Better Understand Adaptation and Evolution: a Review in Marine Macrophytes. Front. Mar. Sci. 7, 545102. 10.3389/fmars.2020.545102 Hughes A. R. Stachowicz J. J. (2011). Seagrass Genotypic Diversity Increases Disturbance Response via Complementarity and Dominance. J. Ecol. 99, 445–453. 10.1111/j.1365-2745.2010.01767.x Hut R. A. Beersma D. G. M. (2011). Evolution of Time-Keeping Mechanisms: Early Emergence and Adaptation to Photoperiod. Phil. Trans. R. Soc. B 366, 2141–2154. 10.1098/rstb.2010.0409 Hut R. A. Paolucci S. Dor R. Kyriacou C. P. Daan S. (2013). Latitudinal Clines: an Evolutionary View on Biological Rhythms,. Proc. R. Soc. B 280, 20130433. 10.1098/rspb.2013.0433 Jackson D. A. (1993). Stopping Rules in Principal Components Analysis: a Comparison of Heuristical and Statistical Approaches. Ecology 74, 2204–2214. 10.2307/1939574 Jackson E. L. Smith T. M. York P. H. Nielsen J. Irving A. D. Sherman C. D. (2021). An Assessment of the Seascape Genetic Structure and Hydrodynamic Connectivity for Subtropical Seagrass Restoration. Restor. Ecol. 29 (1), e13269. 10.1111/rec.13269 Jahnke M. Casagrandi R. Melià P. Schiavina M. Schultz S. T. Zane L. (2017). Potential and Realized Connectivity of the seagrassPosidonia Oceanicaand Their Implication for Conservation. Divers Distrib. 23 (12), 1423–1434. 10.1111/ddi.12633 Jahnke M. D’Esposito D. Orrù L. Lamontanara A. Dattolo E. Badalamenti F. (2019). Adaptive Responses along a Depth and a Latitudinal Gradient in the Endemic Seagrass Posidonia Oceanica . Heredity 122, 233–243. 10.1038/s41437-018-0103-0 Jahnke M. Moknes P.-O. Le Moan A. Martens G. A. Jonsson P. R. (2022). Seascape Genomics Identify Adaptive Barriers Correlated to Tidal Amplitude in the Shore Crab Carcinus maenas. Mol. Ecol. 31, 1980–1994. 10.1111/mec.16371 Jahnke M. Olsen J. L. Procaccini G. (2015b). A Meta-Analysis Reveals a Positive Correlation between Genetic Diversity Metrics and Environmental Status in the Long-Lived seagrassPosidonia Oceanica. Mol. Ecol. 24, 2336–2348. 10.1111/mec.13174 Jahnke M. Pagès J. F. Alcoverro T. Lavery P. S. McMahon K. M. Procaccini G. (2015a). Should We Sync? Seascape-Level Genetic and Ecological Factors Determine Seagrass Flowering Patterns. J. Ecol. 103, 1464–1474. 10.1111/1365-2745.12470 Jiang D. Yang W. He Y. Amasino R. M. (2007). ArabidopsisRelatives of the Human Lysine-specific Demethylase1 Repress the Expression ofFWAandFLOWERING LOCUS Cand Thus Promote the Floral Transition. Plant Cell 19, 2975–2987. 10.1105/tpc.107.052373 Jombart T. (2008). Adegenet: a R Package for the Multivariate Analysis of Genetic Markers. Bioinformatics 24, 1403–1405. 10.1093/bioinformatics/btn129 Jueterbock A. Franssen S. U. Bergmann N. Gu J. Coyer J. A. Reusch T. B. H. (2016). Phylogeographic Differentiation versus Transcriptomic Adaptation to Warm Temperatures inZostera Marina, a Globally Important Seagrass. Mol. Ecol. 25, 5396–5411. 10.1111/mec.13829 Kamvar Z. N. Brooks J. C. Grünwald N. J. (2015). Novel R Tools for Analysis of Genome-wide Population Genetic Data with Emphasis on Clonality. Front. Genet. 6, 208. 10.3389/fgene.2015.00208 Kamvar Z. N. Tabima J. F. Grünwald N. J. (2014). Poppr: an R Package for Genetic Analysis of Populations with Clonal, Partially Clonal, And/or Sexual Reproduction. PeerJ 2, e281. 10.7717/peerj.281 Keenan K. Mcginnity P. Cross T. F. Crozier W. W. Prodöhl P. A. (2013). diveRsity: AnRpackage for the Estimation and Exploration of Population Genetics Parameters and Their Associated Errors. Methods Ecol. Evol. 4, 782–788. 10.1111/2041-210x.12067 Kohorn B. D. (2001). WAKs; Cell Wall Associated Kinases. Curr. Opin. cell Biol. 13, 529–533. 10.1016/s0955-0674(00)00247-7 Langmead B. Trapnell C. Pop M. Salzberg S. L. (2009). Ultrafast and Memory-Efficient Alignment of Short DNA Sequences to the Human Genome. Genome Biol. 10, R25. 10.1186/gb-2009-10-3-r25 Larkum A. W. D. Orth R. J. Duarte C. M. (2006). SEAGRASSES: BIOLOGY, ECOLOGYAND CONSERVATION. Springer Netherlands. Le Corre V. Roux F. Reboud X. (2002). DNA Polymorphism at the FRIGIDA Gene in Arabidopsis thaliana: Extensive Nonsynonymous Variation Is Consistent with Local Selection for Flowering Time. Mol. Biol. Evol. 19, 1261–1271. 10.1093/oxfordjournals.molbev.a004187 Lee H. Golicz A. A. Bayer P. E. Severn-Ellis A. A. Chan C.-K. K. Batley J. (2018). Genomic Comparison of Two Independent Seagrass Lineages Reveals Habitat-Driven Convergent Evolution. J. Exp. Bot. 69, 3689–3702. 10.1093/jxb/ery147 Lee H. Golicz A. A. Bayer P. Jiao Y. Tang H. Paterson A. H. (2016). The Genome of a Southern Hemisphere Seagrass Species (Zostera Muelleri). Plant physiol. 172, 00868.2016. 10.1104/pp.16.00868 Lemopoulos A. Prokkola J. M. Uusi‐Heikkilä S. Vasemägi A. Huusko A. Hyvärinen P. (2019). Comparing RADseq and Microsatellites for Estimating Genetic Diversity and Relatedness - Implications for Brown Trout Conservation. Ecol. Evol. 9, 2106–2120. 10.1002/ece3.4905 Li W. Godzik A. (2006). Cd-hit: a Fast Program for Clustering and Comparing Large Sets of Protein or Nucleotide Sequences. Bioinformatics 22, 1658–1659. 10.1093/bioinformatics/btl158 Linhart Y. B. Grant M. C. (1996). Evolutionary Significance of Local Genetic Differentiation in Plants. Annu. Rev. Ecol. Syst. 27, 237–277. 10.1146/annurev.ecolsys.27.1.237 López-Maury L. Marguerat S. Bähler J. (2008). Tuning Gene Expression to Changing Environments: from Rapid Responses to Evolutionary Adaptation. Nat. Rev. Genet. 9, 583–593. 10.1038/nrg2398 Luu K. Bazin E. Blum M. G. B. (2017). Pcadapt: anRpackage to Perform Genome Scans for Selection Based on Principal Component Analysis. Mol. Ecol. Resour. 17, 67–77. 10.1111/1755-0998.12592 Machado H. E. Bergland A. O. O’Brien K. R. Behrman E. L. Schmidt P. S. Petrov D. A. (2016). Comparative Population Genomics of Latitudinal Variation in Drosophila simulans and Drosophila melanogaster . Mol. Ecol. 25, 723–740. 10.1111/mec.13446 Máñez-Crespo J. Tuya F. Fernández-Torquemada Y. Royo L. Pilar-Ruso Y. d. Espino F. (2020). Seagrass Cymodocea Nodosa across Biogeographical Regions and Times: Differences in Abundance, Meadow Structure and Sexual Reproduction. Mar. Environ. Res. 162, 105159. 10.1016/j.marenvres.2020.105159 Marbà N. Arias-Ortiz A. Masqué P. Kendrick G. A. Mazarrasa I. Bastyan G. R. (2015). Impact of Seagrass Loss and Subsequent Revegetation on Carbon Sequestration and Stocks. J. Ecol. 103, 296–302. 10.1111/1365-2745.12370 Marín-Guirao L. Bernardeau-Esteller J. García-Muñoz R. Ramos A. Ontoria Y. Romero J. (2018). Carbon Economy of Mediterranean Seagrasses in Response to Thermal Stress. Mar. Pollut. Bull. 135, 617–629. 10.1016/j.marpolbul.2018.07.050 Martignago D. Bernardini B. Polticelli F. Salvi D. Cona A. Angelini R. (2019). The Four FAD-dependent Histone Demethylases of Arabidopsis Are Differently Involved in the Control of Flowering Time. Front. Plant Sci. 10, 669. 10.3389/fpls.2019.00669 Mastretta-Yanes A. Arrigo N. Alvarez N. Jorgensen T. H. Piñero D. Emerson B. C. (2015). Restriction Site-Associated DNA Sequencing, Genotyping Error Estimation and De Novo Assembly Optimization for Population Genetic Inference. Mol. Ecol. Resour. 15, 28–41. 10.1111/1755-0998.12291 Mckenna A. Hanna M. Banks E. Sivachenko A. Cibulskis K. Kernytsky A. (2010). The Genome Analysis Toolkit: a MapReduce Framework for Analyzing Next-Generation DNA Sequencing Data. Genome Res. 20, 1297–1303. 10.1101/gr.107524.110 Oetjen K. Ferber S. Dankert I. Reusch T. B. H. (2010). New Evidence for Habitat-specific Selection in Wadden Sea Zostera Marina Populations Revealed by Genome Scanning Using SNP and Microsatellite Markers. Mar. Biol. 157, 81–89. 10.1007/s00227-009-1297-8 Oetjen K. Reusch T. B. H. (2007). Genome Scans Detect Consistent Divergent Selection Among Subtidal vs. Intertidal Populations of the Marine angiospermZostera Marina. Mol. Ecol. 16, 5156–5157. 10.1111/j.1365-294x.2007.03577.x Olsen J. L. Rouzé P. Verhelst B. Lin Y.-C. Bayer T. Collen J. (2016). The Genome of the Seagrass Zostera Marina Reveals Angiosperm Adaptation to the Sea. Nature 530, 331–335. 10.1038/nature16548 Pazzaglia J. Reusch T. B. H. Terlizzi A. Marín‐Guirao L. Procaccini G. (2021). Phenotypic Plasticity under Rapid Global Changes: The Intrinsic Force for Future Seagrasses Survival. Evol. Appl. 14 (5), 1181–1201. 10.1111/eva.13212 Peakall R. Smouse P. E. (2012). GenAlEx 6.5: Genetic Analysis in Excel. Population Genetic Software for Teaching and Research-Aan Update. Bioinformatics 28, 2537–2539. 10.1093/bioinformatics/bts460 Posé D. Verhage L. Ott F. Yant L. Mathieu J. Angenent G. C. (2013). Temperature-dependent Regulation of Flowering by Antagonistic FLM Variants. Nature 503 (7476), 414–417. 10.1038/nature12633 Pratlong M. Haguenauer A. Brener K. Mitta G. Toulza E. Garrabou J. Aurelle D. (2021). Separate the Wheat from the Chaff: Genomic Scan for Local Adaptation in the Red Coral Corallium Rubrum . Peer Community J. 1, 1–None. 10.24072/pcjournal.12 Privé F. Luu K. Vilhjálmsson B. J. Blum M. G. B. (2020). Performing Highly Efficient Genome Scans for Local Adaptation with R Package Pcadapt Version 4. Mol. Biol. Evol. 37, 2153–2154. 10.1093/molbev/msaa053 Prugnolle F. De Meeûs T. (2008). The Impact of Clonality on Parasite Population Genetic Structure. Parasite 15, 455–457. 10.1051/parasite/2008153p455 Reusch T. B. (2014). Climate Change in the Oceans: Evolutionary Versus Phenotypically Plastic Responses of Marine Animals and Plants. Evol. appl. 7 (1), 104–122. 10.1111/eva.12109 Reusch T. B. H. Ehlers A. Hämmerli A. Worm B. (2005). Ecosystem Recovery after Climatic Extremes Enhanced by Genotypic Diversity. Proc. Natl. Acad. Sci. U.S.A. 102, 2826–2831. 10.1073/pnas.0500008102 Reynes L. Aurelle D. Chevalier C. Pinazo C. Valero M. Mauger S. (2021a). Population Genomics and Lagrangian Modeling Shed Light on Dispersal Events in the Mediterranean Endemic Ericaria Zosteroides (= Cystoseira Zosteroides)(Fucales). Front. Mar. Sci. 8. 10.3389/fmars.2021.683528 Reynes L. Thibaut T. Mauger S. Blanfuné A. Holon F. Cruaud C. (2021b). Genomic Signatures of Clonality in the Deep Water kelpLaminaria Rodriguezii. Mol. Ecol. 30, 1806–1822. 10.1111/mec.15860 Richardson B. A. Rehfeldt G. E. Kim M. S. (2009). Congruent Climate‐Related Genecological Responses from Molecular Markers and Quantitative Traits for Western White Pine (Pinus Monticola). Int. J. Plant Sci. 170, 1120–1131. 10.1086/605870 Rismondo A. Curiel D. Scarton F. Mion D. Caniglia G. (2003). “A New Seagrass Map for the Venice Lagoon,” in Proceedings of the Sixth International Conference on the Mediterranean Coastal Environment (Medcoast), 843–852. Roden L. C. Song H.-R. Jackson S. Morris K. Carre I. A. (2002). Floral Responses to Photoperiod Are Correlated with the Timing of Rhythmic Expression Relative to Dawn and Dusk in Arabidopsis . Proc. Natl. Acad. Sci. U.S.A. 99 (20), 13313–13318. 10.1073/pnas.192365599 Rousset F. (2008). genepop'007: a Complete Re-implementation of the Genepop Software for Windows and Linux. Mol. Ecol. Resour. 8, 103–106. 10.1111/j.1471-8286.2007.01931.x Ruggiero M. V. Capone S. Pirozzi P. Reusch T. B. H. Procaccini G. (2005a). Mating System and Clonal Architecture: a Comparative Study in Two Marine Angiosperms. Evol. Ecol. 19 (5), 487–499. 10.1007/s10682-005-0292-x Ruggiero M. V. Reusch T. B. H. Procaccini G. (2005b). Local Genetic Structure in a Clonal Dioecious Angiosperm. Mol. Ecol. 14 (4), 957–967. 10.1111/j.1365-294x.2005.02477.x Ruggiero M. V. Reusch T. B. H. Procaccini G. (2004). Polymorphic Microsatellite Loci for the Marine Angiosperm Cymodocea Nodosa . Mol. Ecol. Notes 4, 512–514. 10.1111/j.1471-8286.2004.00709.x Ruggiero M. V. (2004). Spatial Genetic Structure in the Clonal Marine Angiosperm Cymodocea Nodosa: The Influence of Dispersal Potential, Mating System and Species Interactions. UK: The Open University. Ruiz-Frau A. Gelcich S. Hendriks I. E. Duarte C. M. Marbà N. (2017). Current State of Seagrass Ecosystem Services: Research and Policy Integration. Ocean Coast. Manag. 149, 107–115. 10.1016/j.ocecoaman.2017.10.004 Ruocco M. Musacchia F. Olivé I. Costa M. M. Barrote I. Santos R. (2017). Genomewide Transcriptional Reprogramming in the seagrassCymodocea Nodosaunder Experimental Ocean Acidification. Mol. Ecol. 26, 4241–4259. 10.1111/mec.14204 Sanford E. Kelly M. W. (2011). Local Adaptation in Marine Invertebrates. Annu. Rev. Mar. Sci. 3, 509–535. 10.1146/annurev-marine-120709-142756 Savolainen O. Lascoux M. Merilä J. (2013). Ecological Genomics of Local Adaptation. Nat. Rev. Genet. 14, 807–820. 10.1038/nrg3522 Schneider R. F. Meyer A. (2017). How Plasticity, Genetic Assimilation and Cryptic Genetic Variation May Contribute to Adaptive Radiations. Mol. Ecol. 26, 330–350. 10.1111/mec.13880 Schwacke R. Ponce-Soto G. Y. Krause K. Bolger A. M. Arsova B. Hallab A. (2019). MapMan4: a Refined Protein Classification and Annotation Framework Applicable to Multi-Omics Data Analysis. Mol. plant 12, 879–892. 10.1016/j.molp.2019.01.003 Selkoe K. D’Aloia C. Crandall E. Iacchei M. Liggins L. Puritz J. von der Heyden S. Toonen R. (2016). A Decade of Seascape Genetics: Contributions to Basic and Applied Marine Connectivity. Mar. Ecol. Prog. Ser. 554, 1–19. 10.3354/meps11792 Serra I. A. Innocenti A. M. Di Maida G. Calvo S. Migliaccio M. Zambianchi E. (2010). Genetic Structure in the Mediterranean seagrassPosidonia Oceanica: Disentangling Past Vicariance Events from Contemporary Patterns of Gene Flow. Mol. Ecol. 19, 557–568. 10.1111/j.1365-294x.2009.04462.x Sfriso A. Facca C. (2007). Distribution and Production of Macrophytes and Phytoplankton in the Lagoon of Venice: Comparison of Actual and Past Situation. Hydrobiologia 577, 71–85. 10.1007/s10750-006-0418-3 Shafer A. B. A. Peart C. R. Tusso S. Maayan I. Brelsford A. Wheat C. W. (2017). Bioinformatic Processing of RAD‐seq Data Dramatically Impacts Downstream Population Genetic Inference. Methods Ecol. Evol. 8, 907–917. 10.1111/2041-210x.12700 Sharma N. Geuten K. Giri B. S. Varma A. (2020). The Molecular Mechanism of Vernalization in Arabidopsis and Cereals: Role of Flowering Locus C and its Homologs. Physiol. Plant. 170, 373–383. 10.1111/ppl.13163 Shindo C. Aranzana M. J. Lister C. Baxter C. Nicholls C. Nordborg M. (2005). Role of FRIGIDA and FLOWERING LOCUS C in Determining Variation in Flowering Time of Arabidopsis . Plant physiol. 138, 1163–1173. 10.1104/pp.105.061309 Short F. Carruthers T. Dennison W. Waycott M. (2007). Global Seagrass Distribution and Diversity: a Bioregional Model. J. Exp. Mar. Biol. Ecol. 350, 3–20. 10.1016/j.jembe.2007.06.012 Spalding M. (2003). The Distribution and Status of Seagrasses. University of California Press Berkeley, USA: World atlas of seagrasses, 5–27. St Clair J. B. Mandel N. L. Vance-Borland K. W. (2005). Genecology of Douglas Fir in Western Oregon and Washington. Ann. Bot. 96, 1199–1214. 10.1093/aob/mci278 Stankovic M. Ambo-Rappe R. Carly F. Dangan-Galon F. Fortes M. D. Hossain M. S. (2021). Quantification of Blue Carbon in Seagrass Ecosystems of Southeast Asia and Their Potential for Climate Change Mitigation. Sci. Total Environ. 783, 146858. 10.1016/j.scitotenv.2021.146858 Stinchcombe J. R. Weinig C. Ungerer M. Olsen K. M. Mays C. Halldorsdottir S. S. (2004). A Latitudinal Cline in Flowering Time in Arabidopsis thaliana Modulated by the Flowering Time Gene FRIGIDA. Proc. Natl. Acad. Sci. U.S.A. 101, 4712–4717. 10.1073/pnas.0306401101 Storz J. F. (2005). INVITED REVIEW: Using Genome Scans of DNA Polymorphism to Infer Adaptive Population Divergence. Mol. Ecol. 14, 671–688. 10.1111/j.1365-294x.2005.02437.x Sunde J. Yıldırım Y. Tibblin P. Forsman A. (2020). Comparing the Performance of Microsatellites and RADseq in Population Genetic Studies: Analysis of Data for Pike (Esox lucius) and a Synthesis of Previous Studies. Front. Genet. 11, 218. 10.3389/fgene.2020.00218 Tuya F. Fernández‐Torquemada Y. Zarcero J. del Pilar‐Ruso Y. Csenteri I. Espino F. (2019). Biogeographical Scenarios Modulate Seagrass Resistance to Small‐scale Perturbations. J. Ecol. 107, 1263–1275. 10.1111/1365-2745.13114 Tuya F. Fernández-Torquemada Y. Del Pilar-Ruso Y. Espino F. Manent P. Curbelo L. (2021). Partitioning Resilience of a Marine Foundation Species into Resistance and Recovery Trajectories. Oecologia 196, 515–527. 10.1007/s00442-021-04945-4 Urban M. C. Strauss S. Y. Pelletier F. Palkovacs E. P. Leibold M. A. Hendry A. P. (2020). Evolutionary Origins for Ecological Patterns in Space. Proc. Natl. Acad. Sci. U.S.A. 117, 17482–17490. 10.1073/pnas.1918960117 Valière N. (2002). Gimlet: a Computer Program for Analysing Genetic Individual Identification Data. Mol. Ecol. Notes 2, 377–379. 10.1046/j.1471-8286.2002.00228.x Wang S. Meyer E. Mckay J. K. Matz M. V. (2012). 2b-RAD: a Simple and Flexible Method for Genome-wide Genotyping. Nat. Methods 9, 808–810. 10.1038/nmeth.2023 Wilczek A. M. Burghardt L. T. Cobb A. R. Cooper M. D. Welch S. M. Schmitt J. (2010). Genetic and Physiological Bases for Phenological Responses to Current and Predicted Climates. Phil. Trans. R. Soc. B 365 (1555), 3129–3147. 10.1098/rstb.2010.0128 Winters G. Nelle P. Fricke B. Rauch G. Reusch T. (2011). Effects of a Simulated Heat Wave on Photophysiology and Gene Expression of High- and Low-Latitude Populations of Zostera Marina . Mar. Ecol. Prog. Ser. 435, 83–95. 10.3354/meps09213 Wissler L. Codoñer F. M. Gu J. Reusch T. B. Olsen J. L. Procaccini G. (2011). Back to the Sea Twice: Identifying Candidate Plant Genes for Molecular Evolution to Marine Life. BMC Evol. Biol. 11, 8. 10.1186/1471-2148-11-8 Yan Y. Shen L. Chen Y. Bao S. Thong Z. Yu H. (2014). A MYB-Domain Protein EFM Mediates Flowering Responses to Environmental Cues in Arabidopsis . Dev. cell 30, 437–448. 10.1016/j.devcel.2014.07.004 Zhao M. Yang S. Liu X. Wu K. (2015). Arabidopsis Histone Demethylases LDL1 and LDL2 Control Primary Seed Dormancy by Regulating DELAY of GERMINATION 1 and ABA Signaling-Related Genes. Front. Plant Sci. 6, 159. 10.3389/fpls.2015.00159 Zhou J.-X. Liu Z.-W. Li Y.-Q. Li L. Wang B. Chen S. (2018). Arabidopsis PWWP Domain Proteins Mediate H3K27 Trimethylation on FLC and Regulate Flowering Time. J. Integr. Plant Biol. 60, 362–368. 10.1111/jipb.12630 Zimmerman S. J. Aldridge C. L. Oyler-Mccance S. J. (2020). An Empirical Comparison of Population Genetic Analyses Using Microsatellite and SNP Data for a Species of Conservation Concern. BMC genomics 21, 382. 10.1186/s12864-020-06783-9