Adaptation to local conditions is known to occur in seagrasses, however, knowledge of the genetic basis underlying this phenomenon remains scarce. Here, we analyzed Posidonia oceanica from six sites within and around the Stagnone di Marsala, a semi-enclosed coastal lagoon where salinity and temperature exceed the generally described tolerance thresholds of the species. Sea surface temperatures (SSTs) were measured and plant samples were collected for the assessment of morphology, flowering rate and for screening genome-wide polymorphisms using double digest restriction-site-associated DNA sequencing. Results demonstrated more extreme SSTs and salinity levels inside the lagoon than the outer lagoon regions. Morphological results showed significantly fewer and shorter leaves and reduced rhizome growth of P. oceanica from the inner lagoon and past flowering events were recorded only for a meadow farthest away from the lagoon. Using an array of 51,329 SNPs, we revealed a clear genetic structure among the study sites and confirmed the genetic isolation and high clonality of the innermost site. Fourteen outlier loci were identified and annotated with several proteins including those relate to plant stress response, protein transport, and regulators of plant-specific developmental events. Especially, five outlier loci showed maximum allele frequency at the innermost site, likely reflecting adaptation to the extreme temperature and salinity regimes, possibly due to the selection of more resistant genotypes and the progressive restriction of gene flow. Overall, this study helps us to disentangle the genetic basis of seagrass adaptation to local environmental conditions and may support future works on assisted evolution in seagrasses.

Ninety-five ddRAD-seq library construction and sequencing were conducted at IGATech (Udine, Italy) using an IGATech custom protocol, with minor modifications with respect to Peterson’s double digest restriction-site associated DNA preparation (Peterson et al., 2012). To ensure the quality of sequencing outcomes, for each site, one sample was randomly selected and sequenced twice. The final number of biological replicates for each site was n = 14 for OpenSea-C and n = 15 for the other sites (i.e. North-basin, South-basin, OpenSea-A, OpenSea-B, and OpenSea-D), respectively (i.e. 89 unique samples + 6 technical replicates). In short, gDNA was double digested with both SphI and MboI endonucleases (New England BioLabs). Fragmented DNA was purified with AMPureXP beads (Agencourt) and subsequently ligated with T4 DNA ligase (New England BioLabs). Samples were pooled on multiplexing batches and bead purified as before and then they were size-selected and underwent several purification steps. ddRAD-seq libraries were sequenced with 150 cycles in paired-end mode on NovaSeq 6000 instrument following the manufacturer’s instructions (Illumina, San Diego, CA).

Single nucleotide polymorphisms (SNPs) calling was performed de novo using Stacks software package v2.53 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013). First, raw Illumina reads were demultiplexed using the process_radtags utility (Catchen et al., 2013). The short reads of each sample were assembled into exactly matching stacks using the ustacks utility (Catchen et al., 2013). The creation of the loci catalog (i.e. a set of consensus loci from all the analyzed samples) was done using cstacks and matching each sample against the catalog using sstacks and tsv2bam utilities (Catchen et al., 2013). gstacks ultility (Catchen et al., 2013) was used to pull in paired-end reads, assemble the paired-end contigs and merge them with the single-end locus, align reads to the locus and ultimately call SNPs. Finally, detected loci were filtered using the populations program included in Stacks v2.53 (Catchen et al., 2013), with option –R=0.75 to retain only loci that were represented in at least the 75% of the whole metapopulation and with cutoff --max-obs-het=0.8, to process a nucleotide site at a locus with observed heterozygosity at a maximum of 80%.