# The Seagrass Methylome Is Associated With Variation in Photosynthetic Performance Among Clonal Shoots

Alexander Jueterbock, Christoffer Bostrom, James A. Coyer, Jeanine L.

Olsen, Martina Kopp, Anusha K. S. Dhanasiri, Irina Smolina, Sophie

Arnaud-Haond, Yves van De Peer, Galice Hoarau

► To cite this version:

Alexander Jueterbock, Christoffer Bostrom, James A. Coyer, Jeanine L. Olsen, Martina Kopp, et al. The Seagrass Methylome Is Associated With Variation in Photosynthetic Performance Among Clonal Shoots. Frontiers in Plant Science, 2020, 11, pp.571646. 10.3389/fpls.2020.571646. hal-03410896

HAL Id: hal-03410896

https://hal.umontpellier.fr/hal-03410896

Submitted on 17 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for
updates

# The Seagrass Methylome Is Associated With Variation in Photosynthetic Performance Among Clonal Shoots

Alexander Jueterbock¹,²*, Christoffer Boström³, James A. Coyer²,⁴, Jeanine L. Olsen⁵,
Martina Kopp², Anusha K. S. Dhanasiri², Irina Smolina², Sophie Arnaud-Haond⁶,
Yves Van de Peer⁷,⁸,⁹ and Galice Hoarau²

**OPEN ACCESS**

**Edited by:**
Gidon Winters,
Dead Sea and Arava Science
Center, Israel

**Reviewed by:**
Jeroen T. F. Gillard,
California State University,
Bakersfield,
United States
Shin-Ichiro Ozawa,
Okayama University, Japan
Lazaro Marin-Guirao,
Spanish Institute of
Oceanography, Spain

*Correspondence:
Alexander Jueterbock
Alexander-Jueterbock@web.de

**Specialty section:**
This article was submitted to Marine and Freshwater Plants, a section of the journal Frontiers in Plant Science

**Received:** 11 June 2020
**Accepted:** 21 August 2020
**Published:** 04 September 2020

**Citation:** Jueterbock A, Boström C, Coyer JA, Olsen JL, Kopp M, Dhanasiri AKS, Smolina I, Arnaud-Haond S, Van de Peer Y and Hoarau G (2020) The Seagrass Methylome Is Associated With Variation in Photosynthetic Performance Among Clonal Shoots. *Front. Plant Sci.* 11:571646. doi: 10.3389/fpls.2020.571646

¹Algal and Microbial Biotechnology Division, Faculty of Biosciences and Aquaculture, Nord University, Bodo, Norway,
²Marine Molecular Ecology Group, Faculty of Biosciences and Aquaculture, Nord University, Bodo, Norway, ³Environmental and Marine Biology, Åbo Akademi University, Åbo, Finland, ⁴Shoals Marine Laboratory, University of New Hampshire, Durham, NH, United States, ⁵Ecological Genetics-Genomics Group, Groningen Institute of Evolutionary Life Sciences, University of Groningen, Groningen, Netherlands, ⁶UMR MARBEC, Université de Montpellier, Ifremer, IRD, CNRS, Sète, France, ⁷Department of Plant Biotechnology and Bioinformatics, Ghent University, Ghent, Belgium, ⁸Bioinformatics and Systems Biology, VIB Center for Plant Systems Biology, Ghent, Belgium, ⁹Center for Microbial Ecology and Genomics, Department of Biochemistry, Genetics and Microbiology, University of Pretoria, Pretoria, South Africa

Evolutionary theory predicts that clonal organisms are more susceptible to extinction than sexually reproducing organisms, due to low genetic variation and slow rates of evolution. In agreement, conservation management considers genetic variation as the ultimate measure of a population's ability to survive over time. However, clonal plants are among the oldest living organisms on our planet. Here, we test the hypothesis that clonal seagrass meadows display epigenetic variation that complements genetic variation as a source of phenotypic variation. In a clonal meadow of the seagrass *Zostera marina*, we characterized DNA methylation among 42 shoots. We also sequenced the whole genome of 10 shoots to correlate methylation patterns with photosynthetic performance under exposure to and recovery from 27°C, while controlling for somatic mutations. Here, we show for the first time that clonal seagrass shoots display DNA methylation variation that is independent from underlying genetic variation, and associated with variation in photosynthetic performance under experimental conditions. It remains unknown to what degree this association could be influenced by epigenetic responses to transplantation-related stress, given that the methylomes showed a strong shift under acclimation to laboratory conditions. The lack of untreated control samples in the heat stress experiment did not allow us to distinguish methylome shifts induced by acclimation from such induced by heat stress. Notwithstanding, the co-variation in DNA methylation and photosynthetic performance may be linked via gene expression because methylation patterns varied in functionally relevant genes involved in photosynthesis, and in the repair and prevention of heat-induced protein damage. While genotypic diversity has been shown to enhance stress resilience in seagrass meadows, we suggest that epigenetic variation plays a similar role in meadows dominated by a single genotype. Consequently,
conservation management of clonal plants should consider epigenetic variation as
indicator of resilience and stability.

Keywords: DNA methylation, ecological epigenetics, clonality, heat stress, seagrass, *Zostera marina* (eelgrass)

INTRODUCTION

Genetic variation is considered key to long-term survival of populations (Bijlsma and Loeschcke, 2012), and is recognized as a key for the conservation and restoration of biological diversity (Lande, 1988; Spielman et al., 2004; Laikre et al., 2009). In contrast, lack of genetic variation is regarded as an evolutionary dead-end (Lynch and Lande, 1998) but clonal growth challenges the expected relationship between genetic diversity and long-term survival. For example, conservation management follows the rough guideline that within a population >1,000 genetically different individuals must mate randomly to avoid inbreeding depression, and retain evolutionary potential (Frankham et al., 2014). Yet, roughly 40% of all plants can reproduce asexually (Tiffney and Niklas, 1985), mostly by clonal growth where parental genotypes (genets) grow vegetative modules (ramets) that often remain connected *via* underground stolons or rhizomes.

Although clones benefit from resource sharing, niche specialization, and rapid vegetative growth (Liu et al., 2016), they are predicted to survive only for short periods, and in stable environments (Silvertown, 2008). Asexual reproduction is assumed to lead to slow rates of genetic evolution, and the lack of DNA repair mechanisms afforded by meiosis (Muller's ratchet, mutational meltdown) (Muller, 1964; Gabriel et al., 1993; Lynch et al., 1993). Despite asexual reproduction, many of our most important crops are clones (McKey et al., 2010), including banana, garlic, hops, potatoes, and turmeric, as well as many of the earth's most invasive and oldest plants. For example, genets of Palmer's oak (*Quercus palmeri*) and seagrass (*Posidonia oceanica*) are estimated to be older than 10,000 years (May et al., 2009; Arnaud-Haond et al., 2012). Thus, long-term survival appears not to rely solely on sexual reproduction.

Somatic mutations can create a certain level of genetic diversity, and may explain some evolutionary potential of clonal organisms (Whitham and Slobodchikoff, 1981; Loxdale et al., 2003; Lushai et al., 2003; Reusch and Boström, 2011). For example, ~7,000 single nucleotide polymorphisms (SNPs), 597 in coding regions, and 432 non-synonymous, distinguish ramets of a large Finnish clone of the seagrass *Zostera marina* (Yu et al., 2020). To set this in perspective, 139,321 biallelic SNPs were reported in coding regions among four populations of the same species (Jueterbock et al., 2016). The degree to which epigenetic variation can contribute to phenotypic heterogeneity in ecologically relevant traits, independently from the underlying genetic variation, is a key question in assessing its contribution to stress tolerance, and long-term survival of clonal organisms.

The definition of epigenetics is currently heavily debated (Ptashne, 2007; Greally, 2018). Here, “epigenetics” implies molecular variations that do not alter the DNA sequence but have the potential to change gene expression, and include

non-coding RNAs (ncRNAs), histone modifications, and DNA methylation (Bossdorf et al., 2008). DNA methylation is, from an evolutionary perspective, the most relevant epigenetic mechanism because it can be independent from genetic variation (Bossdorf et al., 2008; Schmitz et al., 2013; Kilvitis et al., 2014), and transgenerationally stable (Boyko et al., 2010; Verhoeven et al., 2010; Ou et al., 2012; Bilichak et al., 2015; Williams and Gehring, 2017). DNA methylation involves the addition of a methyl-group to the C5 position of a cytosine in DNA sequence motifs (CG, CHG, and CHH in plants, where H stands for A, C, or T) (Kilvitis et al., 2014). Depending on sequence context, methylation can be associated with gene activation or silencing (Bossdorf et al., 2008; Niederhuth and Schmitz, 2017). While CG methylation in gene bodies often correlates with increased gene expression, methylation in promoters and repeat regions, such as transposable elements (TEs), silences expression (Feng et al., 2010; Seymour et al., 2014; Dubin et al., 2015; Bewick and Schmitz, 2017; Zhang et al., 2018).

The methylome, or set of DNA methylation modifications in an organism's genome, can change spontaneously at a rate of 2.5 × 10⁻⁴ to 6.3 × 10⁻⁴ methylation polymorphisms per CG site per generation, which is about 7 × 10⁴ higher than the genetic mutation rate of base substitutions per site per generation (Schmitz et al., 2011; van der Graaf et al., 2015). That methylome variation can enhance productivity, and pathogen resistance, has been shown in *Arabidopsis thaliana* plant populations (Latzel et al., 2013). Moreover, methylation variation explained the rapid invasive success of the Japanese knotweed (*Fallopia japonica*) by facilitating differentiation in response to new habitats despite decreased genetic variation (Richards et al., 2012). This suggests that methylation variation complements genetic variation as a source of phenotypic variation in plant populations deprived of genotypic diversity (Gao et al., 2010; Richards et al., 2012; Verhoeven and van Gurp, 2012; Latzel et al., 2013; Zhang et al., 2013; Vanden Broeck et al., 2018).

Unlike genetic variants, methylation variants can also switch state directly in response to the environment (Dowen et al., 2012) and, if stable enough, establish a molecular memory that can be involved in stress priming. Under clonal growth, DNA methylation patterns are expected to be more stably inherited than under sexual reproduction (Verhoeven and Preite, 2014), because clonal growth circumvents epigenetic reprogramming during gameto- and embryogenesis. Although stable transmission across asexual generations has been shown for environment-specific phenotypes and DNA methylation patterns in clonal plants (Verhoeven et al., 2010; Verhoeven and van Gurp, 2012; Vanden Broeck et al., 2018; Verhoeven et al., 2018), a link between methylation variation and fitness-related traits has yet only been demonstrated in non-clonal plants (Latzel et al., 2013). Thus, while epigenetic mechanisms have been suggested to
contribute to clonal plant success (Douhovnikoff and Dodd, 2014; Verhoeven and Preite, 2014; Dodd and Douhovnikoff, 2016; Latzel et al., 2016), empirical evidence is virtually lacking.

Clonal propagation is especially well-developed in aquatic plants (Barrett, 2015). Seagrasses, the only plants to inhabit the marine world, form the foundational basis of some of the most productive and highly diverse coastal marine ecosystems on the planet, and are essential for the health and abundance of economically exploited marine species (Costanza et al., 1997; Larkum et al., 2006; Orth et al., 2006). Ecosystem services are worth more than € 25,000 ha⁻¹ year⁻¹ (Costanza et al., 2014), including nursery grounds, habitat and food for fish and invertebrates, protection of the coastline from erosion, carbon sequestration of up to 186 g C m² yr⁻¹ (Fourqurean et al., 2012; Duarte et al., 2013), and nutrient fixation (Orth et al., 2006; Procaccini et al., 2007).

Over the last decades, losses of seagrass ecosystems have been documented worldwide due to increasing anthropogenic stressors such as invasive species, sediment and nutrient runoff, dredging, aquaculture, rising sea levels, and global warming (Orth et al., 2006; Chefaoui et al., 2018). Losses are expected to accelerate under projected global temperature increase, as ocean warming is considered the most severe threat among climate change factors (Repolho et al., 2017; Duarte et al., 2018) and seagrass meadows have tracked temperature changes in the past (Olsen et al., 2004). How accurate these shifts and losses can be predicted depends on our knowledge of drivers of adaptive potential, including epigenetic diversity (Duarte et al., 2018).

In this study, we characterize the functional relevance of epigenetic variation for heat stress resilience in the seagrass *Zostera marina*. *Z. marina* is the most widely distributed seagrass in the northern hemisphere, inhabiting highly contrasting habitats from sub-arctic to sub-tropical waters (Green and Short, 2004; Olsen et al., 2004; Boström et al., 2014). Few plants display such dramatic range in clonal diversity as *Z. marina* (Reusch et al., 2000; Olsen et al., 2004). Its clonal architecture varies from genetically diverse meadows with high levels of sexual reproduction, to meadows composed of a single large clone due to exclusive vegetative reproduction (Reusch et al., 1999; Olsen et al., 2004; Becheler et al., 2010; Reusch and Boström, 2011). Clonality of *Z. marina* peaks in the Baltic Sea, where large clones were estimated > 1,000 years old (Reusch et al., 1999). These meadows display remarkable phenotypic plasticity and persistence in time under drastic changes in ice cover, salinity, and temperature (Reusch et al., 1999), and under perturbations that represent environmental stress predicted elsewhere for the future (e.g. high temperatures) (Reusch et al., 2018). Moreover, extension of a single seagrass clone in space over ca. 160 × 40 m, and a water depth ranging from 1.5 to 4.5 m with an environmental gradient in light, sedimentation, wind exposure, and ice scour, strongly suggests niche differentiation among ramets of a single clone (Reusch et al., 1999). Thus, clonal meadows (such predominated by a single genet) confound experimental results showing a positive effect of genotypic (i. e. clonal) richness on the productivity and stress resilience of *Z. marina* (Hughes and Stachowicz, 2004; Reusch et al., 2005; Reusch and Hughes, 2006; Ehlers et al., 2008; Hughes et al., 2008). In other words, if stress resilience and tolerance would rely

strongly on genotypic diversity, clones would not be able to disperse widely in space or survive for such long time periods.

The resilience, longevity, and adaptive potential of clonal seagrass meadows remain unknown without a fundamental analysis of their epigenetic variation and its ecological relevance. Therefore, this study tests the hypothesis that variation in DNA methylation can promote functional phenotypic diversity and, thus, may explain how clonal seagrass meadows can persist over millennia.

Specific objectives were to: 1) Characterize DNA methylation variation in an ancient clonal meadow of *Z. marina* (>1,000-years-old); and 2) Identify the functional role of this variation in photosynthetic performance—specifically how photosynthetic performance under benign and stressful temperatures is linked to DNA methylation variation independently from underlying somatic mutations.

# MATERIALS AND METHODS

## Sampling and Cultivation

In June 2015, we sampled seagrass shoots from a clonal meadow (Reusch et al., 1999) in the Baltic Sea (Åland Islands, 60°09'50.4''N, 19°31'48.1''E) (Figure 1) by collecting two to three shoots attached to the same rhizome (ramets) every 3 meters along a 250 m transect (Supplementary Table S1). Along the transect running perpendicular to the shore, water depth increased from 0.5 to 3 m. Shoots were transported in seawater-filled cooling boxes to the field station at Nord University, Norway. The most mature leaf of 42 shoots (field transect samples), randomly chosen from the two to three connected ramets, was flash frozen in liquid nitrogen for subsequent DNA extraction (Supplementary Table S1). Thus, recorded DNA methylation patterns are expected to reflect the average state of methylation of a mature leaf, along which global DNA methylation has been shown to vary with tissue age in the seagrass Posidonia oceanica (Ruocco et al., 2019b).

Ten ramets of the same genet (Supplementary Table S1), sampled at the same time and from the same transect as the field samples, were individually replanted in plastic containers (10 × 15 cm) filled with sand from the sampling site, placed in a 1,280-L aquarium at 15°C (corresponding to field temperature), and illuminated with a 16:8 h light–dark cycle, under light-saturating conditions (200–220 µmol m⁻² s⁻¹OSRAM Fluora, 150 W) (Dennison and Alberte, 1982; Alcoverro et al., 1999). Seawater at 5.5 PSU, corresponding to the salinity at the collecting site, was obtained by mixing freshwater and natural filtered seawater at 32 PSU. The water, filled to a 40-cm level, was kept in constant motion with airstones. Once a week, 50% of the water was renewed after removing epiphytic algae.

## Heat Stress Experiments

After two weeks of acclimation, the 10 clonal ramets were exposed to heat stress (Supplementary Table S1) in a climate chamber (Fitotron, weisstechnik), where they were distributed in three aquaria (60 L, 3–4 ramets in each aquarium) filled with
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Baltic Sea sampling site of a > 1,000 years old clonal meadow where 42 *Zostera marina* shoots were collected along a 250 m transect.</figcaption></figure>

aerated brackish water (5.5 PSU) to 30 cm, of which 50% was weekly renewed. Light was kept at 240 µmol m⁻² s⁻¹ and a 16:8-h light:dark cycle. The temperature in the climate chambers was increased with a daily increment of 3°C from 15°C to 27°C, which can be lethal for *Z. marina* (Greve et al., 2003) and exceeds maxima of sea surface temperatures recorded in the Baltic Sea by ca. 2°C (Reusch et al., 2005), but is likely to be reached in the Baltic Sea during future summer heat waves. After three weeks at 27°C, the temperature was decreased by a daily increment of 3°C from 27°C to 15°C, after which all shoots were returned to the 1,280 L aquarium at 15°C for a recovery period of 5.5 weeks where they were exposed to the same conditions as under acclimation.

At three time points (control, stress, recovery), one entire mature (outer-most) leaf of each shoot was excised and flash-frozen in liquid nitrogen. Thus, the methylation pattern of each sample represents the average methylation state across young and mature tissues, which have been shown to vary in their epigenetic heat-stress response in the seagrass *Posidonia oceanica* (Ruocco et al., 2019a). Control samples were collected at 15°C on the day before the temperature was increased. Stressed samples were collected at 27°C on the day before the temperature was decreased. Recovery samples were collected after the 5.5-week recovery period at 15°C.

## Photosynthetic Performance

At control, stress, and recovery time points, we measured for each shoot (two measurements for each of the two inner leaves) the increase in chlorophyll a fluorescence upon illumination after a 10 min dark period (OJIP curve) (Bussotti et al., 2010), with a PAM-Fluorometer (FluorPen FP100, Photon Systems Instruments) using a saturating pulse of 75% light intensity at

455 nm. From the measurements, we extracted the photosynthetic performance index PiABS (Strasser et al., 2000), reflecting the functionality of photosystem II (PSII) and photosynthetic performance in general (Živčák et al., 2008; Bussotti et al., 2010; Stefanov et al., 2011). PiABS combines three parameters: 1) the density of reaction centers, 2) the efficiency of electron transport beyond Qₐ at the onset of illumination, and 3) the probability for an absorbed photon to reach the reaction center in PS II. PiABS is calculated as: $PiABS = ((1 - F₀/Fₘ)/(M₀/Vⱼ)) × ((Fₘ - F₀)/F₀) × ((1 - Vⱼ)/Vⱼ)$. Here, F₀ is the minimum fluorescence intensity in a dark adapted leaf when all reaction centers are open (all quinone acceptors are oxidized), Fⱼ is the fluorescence intensity at 2 ms illumination, Fₘ is the maximum fluorescence intensity when all reaction centers are closed (all quinone acceptors are reduced), Vⱼ is the relative variable fluorescence at 2 ms calculated as Vⱼ = (FⱼF₀)/(FₘF₀), and M₀ reflects the initial slope of fluorescence kinetics, and is calculated as M₀ = 4 × (F₃₀₀ µsF₀)/(FₘF₀) (Živčák et al., 2008). We tested for correlation of PiABS values (mean values per shoot) between time points using a two-sided Pearson's product moment correlation test in R v3.4.4 (R Core Team, 2019).

## DNA Extraction

All flash-frozen tissue was freeze-dried for one day, then stored at -80°C until DNA extraction. DNA was extracted with the HP Plant DNA mini kit (Omega Bio-Tek, protocol version May 2013) from ≥ 5 mg of lyophilized leaf tissue after grinding with a mixer-mill (Retsch MM 400) in 2 mL Eppendorf tubes with tungsten beads supplied (60 s at 30 Hz). We added 10 µl beta-mercaptoethanol at step 2, and equilibrated the columns at step 8 of the standard protocol. The extracted DNA was eluted in 2 ×
100 µl EB buffer (Qiagen), cleaned, and concentrated with the Clean and Concentrator-5 kit (Zymo Research, protocol v1.2.1) using 15,000 g for all centrifugation steps, and finally eluted in 2 × 30 µl EB Buffer (Qiagen) at 60°C.

## Clonality

In order to identify ramets belonging to the same genet, we genotyped each shoot for seven microsatellite loci (ZosmarGA2, -GA3, -GA6, -CT3, -CT12, -CT19, -CT20) (Reusch and Boström, 2011). PCR was performed using a Veriti 96-Well Thermal Cycler (Applied Biosystems, Life Technologies) in two 10-µl multiplex reactions containing 2 µl cleaned genomic DNA, and 1× AccuStart II PCR ToughMix (Quanta bio). One multiplex reaction contained forward and reverse primers GA2, GA3, and GA6 at 0.5 µM each, and ran at 94°C for 4 min, followed by 30 cycles of 94°C for 60 s, 55°C for 90 s, and 72°C for 90 s, and a final extension at 72°C for 10 min. The other multiplex reaction contained forward and reverse primers CT3 (0.5 µM), CT12 (0.3 µM), CT19 (0.3 µM), and CT20 (0.3 µM), and ran at 94°C for 3 min, followed by 35 cycles of 94°C for 60 s, 57°C for 60 s, and 72°C for 60 s, and a final extension at 72°C for 10 min.

DNA fragment lengths were determined on an ABI 3500xl Genetic Analyzer from 1 µl of 1:99 diluted PCR products mixed with 8.9 µl of HiDi Formamide (Life Technologies), and 0.1 µl of Gene Scan 500 LIZ Size Standard (Life Technologies) after 5 min denaturation at 95°C. Alleles were called with the GeneMapper v4.1 Software (Applied Biosystems, Thermo Fisher Scientific). The shoots were assigned to multi-locus genotypes using the R package "RClone" (Bailleul et al., 2016).

## Whole Genome Sequencing, SNP Detection, and Genetic Distance in Experimental Samples

In order to detect SNPs resulting from somatic mutations in the 10 heat-stressed ramets, genomic DNA libraries were prepared according to the TruSeq DNA PCR-Free (Illumina) protocol, and sequenced on one Illumina HiSeq 3/4000 lane (2 × 150 bp) at the Norwegian Sequencing Centre (University of Oslo, Norway). Raw reads (25.7 million to 44.3 million per library, **Supplementary Table S2**, NCBI BioProject PRJNA575339) were quality-checked with FastQC v0.11.8¹ to control for aberrant read base content, length distribution, duplication, and over-representation. We used TrimGalore! v0.6.0² to remove adapter sequences with a stringency of 3 bp overlap, and low-quality bases with a Phred score Q < 20 (99% base call accuracy). The high-quality reads (25.5 to 44.0 million per library, **Supplementary Table S2**) were mapped to the *Z. marina* genome v2.1 (Olsen et al., 2016) with BWA v0.7.17 (Li and Durbin, 2009). Read duplicates were removed with MarkDuplicatesSpark within GATK v4.1.4.1 (Auwera et al., 2013). SNPs were called with HaplotypeCaller, followed by CombineGVCFs and GenotypeGVCFs within GATK v4.1.4.1 (Auwera et al., 2013).

Before filtering SNPs, we excluded indels, non-variant sites, and alternate alleles not present in any genotypes from the vcf file with SelectVariants within GATK v4.1.4.1 (Auwera et al., 2013). This set of 759,407 raw SNPs was reduced to 105,443 SNPs after hard-filtering with vcffilter from vcflib (Garrison, 2020) with thresholds that were based on density plots drawn with ggplot2 (Wickham, 2016): QualByDepth (QD < 15.0), FisherStrand (FS >12.0), RMSMappingQuality (MQ < 38), MappingQualityRankSumTest (MQRankSum < -1.5), ReadPositionRankSumTest (ReadPosRankSumand < -4.0), and Depth (DP > 4000.0, in order to remove SNPs potentially caused by genome duplication). Subsequently, we used VCFtools v0.1.15 to remove genotypes with genotype quality < 30 (-minGQ 30) or depth < 20 (-minDP 20), and to remove SNPs with more than 2 alleles (-min-alleles 2 -max-alleles 2), with a minor allele frequency of 0.01 (-maf 0.01), and with any missing genotype (-max-missing-count 0). From the remaining 15,508 high-quality SNPs (**Supplementary File S1**) we excluded all that shared the same genotypes among all 10 heat-stressed shoots, as these reflected genetic differences only to the reference genome. The remaining 1,079 SNPs (**Supplementary File S2**) were used to estimate Euclidean genetic distances among the 10 shoots using the R package vcfR v1.9.0 (Knaus and Grünwald, 2017) and the *dist* function of the R package "stats" v3.6.9 (R Core Team, 2019). We tested for correlation between genetic and physical distance among the 10 heat-stressed shoots with Mantel tests in the R package "vegan" v1.4-2 (Oksanen et al., 2019), using 1,000 permutations, and the Pearson's product moment correlation method.

## Methylome Characterization

Sequencing libraries were prepared according to the MethylRAD protocol (Wang et al., 2015) with few adjustments. MethylRAD is a genome-reduction method based on the methylation-dependent restriction enzyme FspEI that targets fully methylated CCGG and CCWGG motifs, thus capturing methylation in CG and CHG sequence contexts. MethylRAD has the potential to reveal genome-wide DNA methylation patterns that are consistent with those generated from Whole Genome Bisulfite Sequencing (Wang et al., 2015). First, sense and anti-sense oligos of adapters A1 and A2 (**Supplementary File S3**) were annealed in 10 µl containing 10 µM of each oligo (Eurofins), 10 mM Tris HCl (Thermo Fisher), 50 mM NaCl (Thermo Fisher), and 1 mM EDTA (Thermo Fisher). Library preparation began with digestion of 100 ng cleaned genomic DNA at 37°C for 4 h in 15 µl containing 4 U FspEI (NEB), 1× CutSmart Buffer (NEB), and 30× Enzyme Activator Solution (NEB). Digestion was verified on a TapeStation 2200 with a D1000 ScreenTape. Second, adapters were ligated to the digested fragments over night at 4°C in 26 µl containing 13 µl digestion solution, 0.1 µM each of two annealed adapters, 1× T4 ligase buffer (NEB), 1.5 µM ATP (NEB) and 1040 U of T4 DNA ligase (NEB). Ligation products were amplified in 20-µl reactions containing 7 µl ligated DNA, 0.05 µM of each primer (P1 and P2, **Supplementary File S3**), 0.3 mM dNTP, 1× Phusion HF buffer (NEB) and 0.4 U Phusion high-fidelity DNA polymerase (NEB). PCR was conducted using a Veriti 96-Well Thermal

<footnote n="1">¹http://www.bioinformatics.babraham.ac.uk/projects/fastqc/</footnote>

<footnote n="2">²https://www.bioinformatics.babraham.ac.uk/projects/trim\_galore/</footnote>
Cycler (Applied Biosystems, Life Technologies) with 16 cycles of 98°C for 5 s, 60°C for 20 s, 72°C for 10 s, and a final extension of 5 min at 72°C. The target band (approx. 100 bp) was extracted from a 2% E-Gel (Thermo Fisher). For multiplex sequencing, shoot barcodes were introduced by means of PCR. Each 20 µl PCR reaction contained 12 µl of gel-extracted PCR product, 0.2 µM of each primer (P3 and index primer, **Supplementary File S3**), 0.3 mM dNTP, 1x Phusion HF buffer (NEB) and 0.4 U Phusion high-fidelity DNA polymerase (NEB). PCR was conducted with the same PCR cycling program outlined above. PCR products were purified using AMPURE XP beads (Beckman Coulter) using a 1.8:1 volume ratio of beads to product, and a final elution in 22 µl EB buffer (Qiagen). The purified fragments were sequenced on an Illumina NextSeq 500 (1 × 75 bp) using three high-output flow-cells: The 10 experimental samples were sequenced on a single flow-cell, while the field transect samples were split over two other flow-cells (detailed in **Supplementary Table S3**).

The sequenced reads were quality-trimmed with TrimGalore! v0.4.1¹ by removing the adapter sequences with a stringency of 3 bp overlap, low-quality bases with a Phred score Q < 20, and the terminal 2 bp from both ends in order to eliminate artifacts that might have arisen at the ligation position. Quality was checked with FastQC v0.11.8⁴ to control for aberrant read base content, length distribution, duplication and over-representation.

The high-quality reads were mapped with SOAP v1.11 (Li et al., 2008) to 628,255 *in silico* predicted MethylRAD tags that were extracted from the *Z. marina* genome v2.1 from ORCAE (Sterck et al., 2012) with the custom python script *InSilicoTypeIIbDigestion.py*<fnref n="5" />. For mapping, we allowed for two mismatches, filtered reads with >1 N, and used a seed size of 8 bp. Based on the uniquely mapped reads, we counted the coverage of each methylated site for each shoot using htseq-count (v0.7.2). Methylation calls were retained only for sites with ≥ 2× coverage, which reduced the false-positive rate from 1.10% to 0.23% for CG sites and from 2.50% to 0.89% for CHG sites. False-positive rates were estimated as the percentage of methylation sites supported by at least two reads in the generally unmethylated chloroplast genome. For each shoot, raw counts were normalized to reads-per-million by dividing reads per site through the total number of reads per shoot library, times one million.

The methylation sites were annotated with the v2.1 gff3 file from ORCAE (Sterck et al., 2012), and separated into genes, intergenic regions, and TEs. TEs were located in both genes and intergenic regions, and contained repeats classified by RepeatModeler<fnref n="6" /> as rnd-1/2/3/4/5 families, referring to the series of processing rounds of the *de novo* TE family identification program. The methylation sites were further separated into such containing CG and CHG recognition sites. Sequence contexts for all 628,255 *in silico* predicted MethylRAD tags are listed in **Supplementary Table S4**.

## Methylation Variation Between the Field Transect Shoots

In order to describe the level of intra- and inter-clonal epigenetic distance among the 42 transect shoots (of which shoot 27 belonged to a different genet than the 41 other ramets, **Supplementary Table S1**), we calculated the Euclidean distance between shoots based on their coordinates in the 2-dimensional PCA (Principal Component Analysis) plot using the *dist* function of the R package “stats” v3.6.0 (R Core Team, 2019). PCA was done on reads-per-million for all shoots with the PCA function of the R package “FactoMineR” (Lê et al., 2008). We tested for correlations between epigenetic and physical distance among the 41 ramets with Mantel tests using the R package “vegan” v1.4-2 (Oksanen et al., 2019) with 1,000 permutations, and the Pearson’s product moment correlation method. *P*-values were adjusted according the Benjamini-Hochberg method (Benjamini and Hochberg, 1995).

## Correlation Between Methylome Variation and Photosynthetic Performance

For the 10 experimental samples, we tested for correlations between epigenetic distance and photosynthetic performance difference (PiABS values) among shoots, while controlling for genetic distance with Partial Mantel tests using the function *partial.mantel.test* of the R package “ncf” v1.2.9 (Lê et al., 2008). Only significant correlations (p<0.05, adjusted for multiple comparisons according to (Bjornstad, 2020) in R (R Core Team, 2019)) with coefficients *R* > 0.65, when controlled for genetic distance, were considered strong enough to be reckoned as biologically linked. The same analysis ran Mantel tests for correlation between genetic and epigenetic distance, and between genetic distance and photosynthetic performance difference.

## Methylome Shift Under Experimental Conditions

The methylome of the 10 heat-stressed shoots was characterized under control, stress, and recovery conditions (3 shoots died from the stress), and compared with the methylome of field transect samples with PCA on reads-per-million data using the PCA function of the R package “FactoMineR” (R Core Team, 2019). Epigenetic distances between the samples were estimated for each sequence context (gene, intergene, TE, each in CG and CHG regions, respectively) as the Euclidean distances in the 2-dimensional PCA plot using the *dist* function of the R package “stats” v3.6.9 (Lê et al., 2008).

## Differential Methylation Analyses

To estimate the number of sites that changed in methylation state from control to stress and recovery conditions, we applied differential methylation analyses using the R package “edgeR” 3.20.9 (R Core Team, 2019) within the “SARTools” pipeline v1.6.6 (Robinson et al., 2009). Read counts were normalized using a trimmed mean of M-values (TMM) between each pair of samples (Varet et al., 2016). All samples taken under control, stress, and recovery served as replicates for the three different conditions. Methylation levels were considered significantly

<footnote n="3">³https://www.bioinformatics.babraham.ac.uk/projects/trim\_galore/</footnote>

<footnote n="4">⁴http://www.bioinformatics.babraham.ac.uk/projects/fastqc/</footnote>

<footnote n="5">⁵http://marinetics.org/2017/04/11/REdigestions.html</footnote>

<footnote n="6">⁶http://www.repeatmasker.org/RepeatModeler/</footnote>
increased (hyper-methylated) or decreased (hypo-methylated) when Benjamini-Hochberg adjusted *p*-values (Robinson and Oshlack, 2010) fell below 0.05.

To identify how the methylome differed between samples of high and low photosynthetic performance, we used differential methylation analyses using the R package “edgeR” 3.20.9 (Benjamini and Hochberg, 1995) within the “SARTools” pipeline v1.6.6 (Robinson et al., 2009). Read counts were normalized using a trimmed mean of M-values (TMM) between each pair of samples (Varet et al., 2016). Differential methylation analysis was done only for conditions and sequence contexts where a positive correlation was found between photosynthetic performance differences and epigenetic distances. For control conditions, we compared methylation levels in genes (CG regions) between the two samples of highest photosynthetic performance (79.1, and 13.2) and the two samples of lowest photosynthetic performance (63.1 and 59.2, **Figure 2**). For recovery conditions, we compared methylation levels in intergenic TEs (CHG regions) between the two samples of highest photosynthetic performance (79.1 and 17.1) and the two samples of lowest photosynthetic performance (57.1 and 59.2). Methylation levels were considered significantly different between shoots of high and low photosynthetic performance, when Benjamini-Hochberg adjusted *p*-values (Robinson and Oshlack, 2010) fell below 0.05.

For differentially methylated sites within gene bodies, we tested with Fisher's exact tests for enrichment of gene ontology terms of biological processes with the R package “topgo” (Benjamini and Hochberg, 1995) using Fisher's exact tests. GO terms were obtained from the v2.1 *Zostera* genome annotation from the ORCAE database (Alexa and Rahnenführer, 2010). To reduce redundancy in the significantly enriched GO terms (*p*-values < 0.05), we calculated “sim rel” scores (Sterck et al., 2012) (Allowed similarity=0.5), based on the *A. thaliana* GO-term database, using the REVIGO web server (Schlicker et al., 2006).

# RESULTS

## Methylome Characterization and Variation Among the Clonal Transect Shoots

On average, 74 million high-quality reads were obtained per sequencing library (*Methylome Characterization* section, ranging from 0.7 to 151 million **Supplementary Table S3**). DNA raw reads are accessible from NCBI under BioProject number PRJNA575339. On average, 35% of the high-quality reads mapped to the *in silico* digested *Z. marina* genome, and 11% mapped uniquely (annotated reads-per-million in **Supplementary File S4**). In total, 144,420 sites were methylated (covered by at least two reads) across all transect shoots. Across all transect shoots, 84,640 methylated CG sites and 59,780 methylated CHG sites were detected, which represents 59% and 41% of all methylated sites, respectively (**Supplementary Table S5**). About 23% of all methylated sites fell in gene bodies (of which 41% in TEs), and 77% in intergenic regions (of which 42% in TEs) (**Figure 3**). In gene bodies, 67% of the methylated sites fell in CG regions, in intergenic regions only 56%.

Based on their multi-locus microsatellite genotype, 41 of 42 shoots sampled along the transect (except sample 27 in **Supplementary Table S1**) belonged to the same genet (*Clonality* section; **Supplementary Table S6**). Methylome variation between these 41 clonal ramets (*Methylation Variation between the Field Transect Shoots*) exceeded the methylome shift over the course of the heat-stress experiment (*Methylome Shift Under Experimental Conditions* section; **Figure 4**). Epigenetic distances between the 41 ramets of one genet and the single ramet of the other genet was not higher than among the 41 clonal ramets (*Methylation Variation between the Field Transect Shoots* section; **Figures 5A, B, Supplementary Table S5**). Epigenetic distances were generally lower in CHG than in CG sequence contexts, and lowest in TEs in gene regions (**Figures 5A, B**). Epigenetic distances were not significantly correlated with physical distances in any sequence

<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Photosynthetic performance (PiAbs, absolute values) for all ten heat-stressed *Zostera marina* shoots at control, stress, and after recovery. Shoots 13.2, 15.2, and 27.2 did not recover from the stress.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | Proportion of methylated sites detected in genes, intergenic regions, and transposable elements (TEs), separated by CG and CHG sequence contexts.</figcaption></figure>

context (Methylation Variation between the Field Transect Shoots section; **Supplementary Table S7**).

## Genetic Variation Among Heat-Stressed Shoots

Based on whole genome sequences of the 10 heat-stressed shoots, we identified 15,508 high-quality SNPs (Whole Genome Sequencing, SNP Detection, and Genetic Distance in Experimental Samples section; **Supplementary File S1**). That all 10 shoots shared the same heterozygous state in 14,429 (93%) of all SNPs, and shared the same multi-locus microsatellite genotype (Clonality section; **Supplementary Table S6**), suggests that they were clones having originated from one ancestral zygote (**Supplementary Figure S1**) (Supek et al., 2011).

Despite being clonal, the 10 shoots differed in 1,079 SNPs resulting from somatic mutations (Whole Genome Sequencing, SNP Detection, and Genetic Distance in Experimental Samples section). Based on these SNPs, Euclidean genetic distances ranged from 9 to 35 (frequency distribution in **Supplementary Figure S2**), and were not significantly correlated (Whole Genome Sequencing, SNP Detection, and Genetic Distance in Experimental Samples section) with physical distances between shoots (R=0.18, p=0.12).

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | Methylation variation of field transect samples and of experimental samples of *Zostera marina* along the first two principal components. The samples (small symbols) are plotted along the first two principal components (Dim) based on methylation profiles across all sequence contexts. Circles represent 95% confidence intervals around group means (large symbols). Bracketed numbers represent the percentage of explained variation.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

**FIGURE 5** | Boxplot of relative epigenetic distances (distances divided by the maximum distance for each sequence context) showing the median, upper and lower quartiles (25th and 75th percentiles, box margins), minimum and maximum values falling within 1.5 X the interquartile range above and below the box (ends of the vertical lines), and all outlying points individually. Epigenetic distances in (A) CG and (B) CHG regions were comparable within ramets of the same genet (intraclone) and between ramets of different genets (interclone) of *Zostera marina* across all sequence contexts.

## Link Between Epigenetic Distance and Photosynthetic Performance Changes Under Heat Stress

Photosynthetic performance (PiABS, Photosynthetic Performance section) declined in all 10 shoots under heat stress (Figure 2, Supplementary Table S8, raw data in Supplementary Table S9). Seven of the 10 shoots recovered from the heat stress, yet photosynthetic performances did not reach pre-stress levels. Photosynthetic performances under control and recovery conditions were positively correlated (adjusted *p*<0.05, *R*=0.93, **Figure 6, Supplementary Table S10**).

Genetic distances correlated moderately (0.4<R<0.6, adjusted *p*<0.05) with photosynthetic performance differences in recovered samples, and with epigenetic distances (Correlation between Methylome Variation and Photosynthetic Performance section) in some sequence contexts of stressed and recovered samples (Figure 6, Supplementary Table S10). Nevertheless, epigenetic distances correlated strongly (R>0.65, adjusted *p*<0.05) with photosynthetic performance differences even after controlling for genetic distance based on 1,079 SNPs: 1) epigenetic distances among control samples in CG gene body regions correlated with photosynthetic performance differences prior to stress, and after recovery (stress resilience, Figure 6, Supplementary Table S10). However, some of the samples, including one that performed well before the stress (sample 13.2 in Figure 6), did not recover from the heat stress

(Photosynthetic Performance section). Epigenetic distances among recovered samples in CHG regions of intergenic TEs correlated with photosynthetic performance differences after recovery (Figure 6, Supplementary Table S10).

## Methylome Changes in the Course of the Heat-Stress Experiment

Over the course of the heat-stress experiment, methylation patterns in all sequence contexts changed and did not return to but instead diverged further from control (pre-stress) patterns during the recovery period (Methylome Shift Under Experimental Conditions section; Figure 7 for all sequence contexts combined, Supplementary Figure S3 for the different sequence contexts, Supplementary File S5 listing annotated reads-per-million for each sample). More sites became hyper- (increased in methylation) than hypo-methylated (decreased in methylation) in the course of the experiment (Differential Methylation Analyses section 2.11; Figures 8A, D). Methylation levels differed significantly at 437 sites between control and stress conditions (257 hyper- and 180 hypo-methylated), at 1788 sites between control and recovery conditions (1141 hyper-, and 647 hypo-methylated), and only at 39 sites between stress and recovery conditions (18 hyper-, and 21 hypo-methylated, Supplementary Table S11). After recovery, CG methylation had changed more strongly in gene body regions, and CHG methylation in intergenic regions (Figures 8B, E).
<figure><img src="image_8.png" /><figcaption>**FIGURE 6** | Correlation matrix of photosynthetic performance (PiABS) differences (first three columns) or of genetic distance (last column) among clonal *Zostera marina* shoots with epigenetic distances (first 24 rows) or with photosynthetic performance differences at control, stress, or recovery conditions (last three rows). All correlations between photosynthetic performance differences and epigenetic distances were controlled for genetic distances. Black squares represent untested correlations that were considered biologically not meaningful. Pearson product-moment correlation coefficients R are encoded by the color gradient explained in the bar on the right end, and by the shape of the ellipses. Narrow ellipses represent stronger correlations as compared with wide ones. Asterisks highlight strong (R>0.65) and significant (adjusted p<0.05) correlations. CG and CHG: sequence contexts of the methylated cytosine; TE: Transposable element.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**FIGURE 7** | Methylation patterns of all sequence contexts in shoots of a *Zostera marina* over the course of the stress experiment. See Supplementary Figure S3 for methylation shifts separated by sequence context. The samples (small symbols) are plotted along the first two principal components (Dim) based on methylation profiles across all sequence contexts. Circles represent 95% confidence intervals around group means (large symbols). Bracketed numbers represent the percentage of explained variation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 8** | Hyper- and hypo-methylated sites in CG (A–C) and CHG (D–F) regions in stressed vs. control samples (A, D), recovered vs. control samples (B, E), and in recovered vs. stressed samples (C, F). Hyper-methylation refers to higher methylation, hypo-methylation to lower methylation in samples taken during stress as compared with samples taken before stress (A, D) in samples taken after recovery as compared with samples taken before stress (B, E), and in samples taken after recovery as compared with samples taken during stress (C, F). TE, transposable element.</figcaption></figure>
After recovery, methylation increased (*Differential Methylation Analyses* section) in gene bodies with functions including DNA transcription and replication (GO:0006353, GO:0006261, GO:0006270); catabolism of misfolded proteins (GO:0006515), gamma-aminobutyric acid (GABA) (GO:0009450, GO:0009448), and neurotransmitter (GO:0042133, GO:0042135, GO:0001505); as well as amino acid metabolism (GO:0009072, GO:0009073, GO:0009063) (**Supplementary Table S12**). DNA methylation had decreased in gene bodies with functions including transmembrane transport of cations, ions, and protons (GO:0006811, GO:0055085, GO:0006810, GO:0034220, GO:0006812, GO:0098655, GO:0098660, GO:0015672, GO:0099132, GO:0015991, GO:0098662, GO:0099131, GO:0015988, GO:1902600); transport of ammonium and phospholipids (GO:0072488, GO:0015696, GO:0015914, GO:0006869); localization of lipids and organelles (GO:0051179, GO:0051234, GO:0010876, GO:0051640); exocytosis (GO:0006887, GO:0006904); and secretion (GO:0032940) (**Supplementary Table S12**).

## Differential Methylation Between Shoots of High and Low Photosynthetic Performance

Recovered shoots of high and low photosynthetic performance presented too similar methylation patterns to allow for detection of many regions showing significant differentiation in methylation levels (*Differential Methylation Analyses* section). Only four CHG sites in intergenic TEs were hyper-methylated in recovery samples of high photosynthetic performance (**Supplementary Table S13**). In contrast, 90 CG gene body sites were hyper-methylated (1 hypomethylated) in control samples of high photosynthetic performance (**Supplementary Table S13**). Enriched biological processes in the 90 hyper-methylated CG sites included “light harvesting in photosystem I” (GO:0009768, GO:0009765) and “protein folding” (GO:0006457, GO:0042026, GO:0006515) (**Supplementary Table S14**).

## DISCUSSION

Plant genes persisting >1,000 years challenge the positive correlation between genetic variation retained through recombination, and long-term survival. Our study shows for the first time that ramets of the same seagrass genet display DNA methylation variation that is independent from underlying genetic variation, and associated with phenotypic variation in the fitness-related trait of photosynthetic performance under experimental conditions. To what degree this association could have been affected by transplantation-induced stress responses in the methylome remains unknown, given the strong methylome shift under acclimation to laboratory conditions. The lack of untreated control samples in the heat stress experiment did not allow us to distinguish methylome shifts induced by acclimation from such induced by heat stress. Nevertheless, our findings support the hypothesis (Yu et al., 2020) that methylation variation, via variation in gene regulation, compensates potential costs of clonal reproduction (Dodd and Douhovnikoff, 2016), and contributes to

the long-term survival of clonal seagrass meadows by increasing variation in ecologically relevant traits that cannot be simply explained by the underlying somatic genetic variation.

## Methylome Variation of Functional Relevance

Methylome variation among ramets of the same genet resulted either from random epimutations or from microscale variations in the environment because its correlation with physical distance between shoots was weak. By sampling entire leaflets composed of young and mature tissue, we aimed to standardize variation in tissue maturity among shoots, since tissue maturity can affect methylation patterns (Lynch et al., 1993; Lynch and Lande, 1998). None of the shoots had been reproductive at the time of sampling, so that variation in reproductive status is unlikely to explain the recorded methylome variation among shoots (Ruocco et al., 2019b). The change in depth of 3 meters and, thus, gradual changes in environmental factors along the sampled transect may not have been extreme enough to impose a result. The disagreement between methylome variation and transect position may also result from recent uprooting and re-settling of some shoots. This was suggested to explain disagreement between genetic similarity of clonal shoots and their transect position in another clonal meadow of *Z. marina* (Marín-Guirao et al., 2019), and would further explain the absence of correlation between genetic and physical distance in our heat-stressed ramets. In such case, the re-settled shoots would display a methylome shaped by distant environmental conditions.

Epigenetic distance was comparable at the intra- and inter-clonal level at all sequence contexts (Figures 5A, B). It is important to note that the inter-clonal comparisons relied on 41 shoots of one genet versus a single shoot of another genet. That we could not identify more genets in our samples supports the finding that this meadow is dominated by a single genotype (Reusch et al., 1999). It is unlikely that this genet could belong to the same genet as the other 41 shoots with a somatic mutation explaining its different microsatellite genotype at locus *ZosmarCT19* (heterozygous with a 147bp fragment, and not homozygous for the 150bp fragment as the other shoots). This microsatellite locus had been developed as part of a set of seven less variable loci to discriminate genets, in contrast to other hypervariable microsatellite loci that allow to identify somatic mutations (Yu et al., 2020). However, the inter-clonal comparisons have to be taken with caution as it remains unknown how different the methylation profile of one or more additional genets from the same meadow could be.

The seagrass methylome is flexible and responds directly to environmental change, given its strong change from field transect samples to acclimated control samples within two weeks (Figure 4). Three factors may explain why the shift in methylation patterns appears stronger in response to acclimation to laboratory conditions than in response to heat stress (Figure 4). First, several environmental factors had likely changed from field to lab conditions, including light intensity, and nutrient composition of the water. The wider variation in methylation profiles among field samples than among lab samples may mirror microscale variations in the environment and appears
to have converged to more similar methylation patterns under controlled and uniform laboratory conditions. Although the applied heat stress (27°C) was strong, only temperature had changed over the course of the experiment, potentially requiring less changes in methylation state as compared with lab acclimation. Second, field and acclimation epigenotypes differed already before the acclimation. Thus, the methylation changes from field to lab conditions can not only be ascribed to changes in environmental conditions but also to intrinsic differences between the epigenotypes. This applies to a lower extent to the experimental samples, because the same epigenotypes had been sampled throughout the experiment. Third, since we had not run technical replicates on the three sequencing flow-cells, part of the methylome shift under acclimation may have been due to flow-cell batch effects. All experimental samples had been run on a different flow cell than the field transect samples, which had been run on two other flow-cells (Supplementary Table S3). Although the 95% confidence intervals of methylation patterns of the field transect samples from different flow-cells did not overlap, the patterns were still more similar between all field samples than between field and experimental samples (Supplementary Figure S4). This does not alter the conclusion that there has been a shift in methylation patterns from field to lab conditions, because it is supported by two independent sequencing runs. However, this shift may be weaker than shown in Figure 4, due to sequencing batch effects.

Methylome variation in CG gene regions predicted photosynthetic performance (Figure 6) and, thus, is likely functionally relevant for this fitness-relevant trait. Methylation variation did not predict photosynthetic performance under acute heat-stress, given that the photosynthetic performance of the different epigenotypes converged to low and highly similar photosynthetic performance values under stress (Figure 2). Although difficult to prove (Reusch and Boström, 2011), the correlations between methylome and photosynthetic performance differences before and after the heat-stress may be causal, given that control samples of high photosynthetic performance showed hyper-methylated CG sites in gene bodies with relevant functions: “light harvesting in PSI” (GO:0009768), “protein folding” (GO:0006457), “protein refolding” (GO:0042026), and “misfolded or incompletely synthesized protein catabolic process” (GO:0006515) (Supplementary Table S14). While increased light harvesting can enhance photosynthetic performance, the functions related to protein folding can prime against heat stress when increased methylation is associated with expression of the underlying genes and, thus, with accumulation of protective molecules like heat-shock proteins. Heat shock proteins are involved in repair and prevention of heat-induced protein damage (Crisp et al., 2016), and that can facilitate a fast stress response. Indeed, while genes are generally repressed by methylated promoter regions (Feder and Hofmann, 2002), genes are often activated by CG gene body methylation (Lisch, 2009) that prevents aberrant expression from intragenic promoters (Zhang et al., 2006; Schmitz et al., 2013; Yang et al., 2014; Dubin et al., 2015; Bewick and Schmitz, 2017; Niederhuth and Schmitz, 2017). Accordingly, in the seagrass

*Posidonia oceanica* hyper-methylation was associated with transcriptionally active leaf segments (Takuno and Gaut, 2012; Neri et al., 2017), and appears to be a common response to stress (Greco et al., 2012, 2013; Ruocco et al., 2019a). Thus, *via* association with gene expression, the observed methylome variation in control samples may putatively create ecologically relevant phenotype variation that predicts photosynthetic performance. However, it is important to note that the methylome variation in control samples does not reflect the methylome variation present in the field but may have already been shaped by transplantation-related stress, given the strong methylome shifts under acclimation to laboratory conditions. It is also important to note that the methylome of a high-performing shoot does not always protect from stress, given that one of the highest-performing samples (shoot 13.2 in **Figure 2**) did not recover from the heat stress.

Only four CHG sites in intergenic TEs were differentially methylated between the recovered samples of high and low photosynthetic performance (Supplementary Table S13). This appears to contradict the correlation between methylation patterns and photosynthetic performance. However, especially the samples of the high-performance group (79.1 and 17.1) showed a big difference in photosynthetic performance values between themselves (Figure 2) and, thus, were likely also too different in methylation patterns to allow for the detection of sites being differently methylated as compared with the samples of the low-performance group.

Although our study was based on clonal ramets that presented identical seven-locus microsatellite genotypes, the 10 heat-stressed ramets varied genetically at 1,079 SNPs generated by somatic mutations. Epigenetic differentiation between populations can be tightly linked to genetic differentiation (Ruocco et al., 2019b). Accordingly, trans-acting SNPs appeared to explain the correlation between gene body methylation and latitude in *Arabidopsis* (Herrera and Bazaga, 2010; Gáspár et al., 2019). In contrast, part of the correlation between DNA methylation with ecological factors in landscape epigenomic studies on non-model plant species could not be predicted from the observed underlying patterns of genetic relatedness (Dubin et al., 2015). Genetic variation could also not explain the association between epigenetic variation and phenotypic variation (leaf, petiole and functional traits) in natural populations of the perennial herb *Helleborus foetidus* (Schulz et al., 2014; Foust et al., 2016; Gugger et al., 2016), and in clonal populations of the introduced clonal herb *Hydrocotyle vulgaris* (Medrano et al., 2014). In agreement, the correlation between methylome variation and photosynthetic performance in our study could not be predicted from the underlying genetic variation. This suggests that differences in performance among ramets in clonal plants can be explained not exclusively by somatic mutations (Wang et al., 2020) but at least partly by independent differences in their methylome. Independent from genetic variation was also the methylome shift in the experimental samples that provides potential for temporally stable phenotypic change at time scales unattainable by somatic mutations. Particularly across mitotically grown generations, epigenetic patterns can be expected to be more faithfully
inherited than across sexual generations, because clonal growth circumvents epigenetic reprogramming during meiosis and embryogenesis (Santelices et al., 2018; Simberloff and Leppanen, 2019; Yu et al., 2020). Thus, our results provide a first indication that DNA methylation variation provides a layer of ecologically relevant phenotypic variation that is independent from genetic variation in clonal seagrass meadows.

## Methylome Shift Over the Course of the Heat-Stress Experiment

Methylation profiles changed over the course of the heat stress experiment with a clear shift from control to stress and recovery profiles (Figure 7). Interestingly, methylation profiles were highly similar between the stress and recovery phase as shown by the overlap of the 95%-confidence intervals around the group means (Figure 7, Supplementary Figure S3), and the low number of differentially expressed sites between stress and recovery samples (Figures 8C, F, Supplementary Table S11). Nevertheless, recovery profiles showed an apparently stronger divergence from the control profiles (Figures 8B, E) than the stress profiles did (Figures 8A, D). This shift in methylation profiles may be due to different processes that likely acted in parallel.

A first explanation is that part of the recorded methylome shift was likely induced by heat stress, given that the recovered shoots showed hyper-methylation and, thus, potentially constitutive upregulation (Zhang et al., 2006; Schmitz et al., 2013; Yang et al., 2014; Dubin et al., 2015; Niederhuth and Schmitz, 2017) of genes involved in the functions of catabolism of misfolded proteins (GO:0006515) (Supplementary Table S12). This suggests increased investment in the breakdown of heat-denatured proteins, a common response to heat-stress (Feder and Hofmann, 2002). This function was also hyper-methylated in control samples of highest photosynthetic performance (Supplementary Table S14), that also recovered better from the stress (although some samples had died, Figure 2). Heat-responsive DNA methylation changes in plants appear not to show a consistent trend across different species, and little is yet known about their functional role (Liu et al., 2015). In Brassica rapa, heat-responsive DNA methylation was shown to be associated with differential expression of genes involved in RNA metabolic processing, and in heat stress signal transduction (Liu et al., 2018). Stress signal transduction may also have been affected by the experimental methylation changes in Z. marina as enriched biological processes in recovered shoots included the regulation of transmembrane transport of ions and protons (GO:0006811, GO: 0055085, GO:0006810, GO:0034220, GO:0006812, GO:0098655, GO:0098660, GO:0015672, GO:0099132, GO:0015991, GO:0098662, GO:0099131, GO:0015988, GO:1902600), and neurotransmitter levels (GO:0042133, GO:00042135, GO:0001505) (Supplementary Table S12). Experimental removal of DNA methylation, e.g. using zebularine or 5-Azacytidine (Griffin et al., 2016), or the targeted change of methylation patterns, e.g. via CRISPR (Xu et al., 2016), will ultimately allow us to identify the relationship between methylation changes under stress and adaptive phenotypic changes.

However, an alternative explanation for the shift in methylome profiles over the course of the experiment is the continued

acclimation of seagrass shoots to laboratory conditions. In order to correlate methylome differences with photosynthetic performance differences between the epigenotypes, we had to choose a longitudinal sampling design in which we followed single shoots/epigenotypes over the course of the experiment. To reduce the effect of acclimation on methylome profiles, we allowed for a 2-week acclimation period before the start of the experiment, but can not rule out that the methylome shift from field to lab conditions (Figure 4) was completed after this time period. Moreover, although we had planted the shoots in soil from the sampling site and observed continuous growth of the shoots, we cannot exclude the possibility that the dilution of natural seawater with freshwater may have resulted in lower nutrient contents as compared with field conditions, which may have induced a stress factor in addition to the applied heat. Most likely, heat-stress, acclimation, and nutrient depletion have contributed in parallel to the recorded methylome shift and can not be disentangled in this study. Thus, future studies need to run a cross-sectional sampling design with unstressed samples as control treatment in parallel to a longitudinal sampling design that traces methylome changes in single epigenotypes, in order to distinguish methylome changes induced by transplantation to laboratory conditions from such induced exclusively by the applied heat-stress.

If the recorded methylation shift over the experiment has been primarily a response to the applied heat-stress, the similarity in methylation profiles between stress and recovery phases can be interpreted as the formation of an epigenetic stress memory (Crisp et al., 2016). Continued divergence of methylation profiles from stress to recovery conditions (Figures 8B, D) agrees with increasing divergence of gene transcription profiles from heat-stress to recovery conditions in a Danish *Z. marina* population (Franssen et al., 2011). We speculate that gene expression changes, and other molecular mechanisms involved in the heat-stress response, could have triggered additional methylation changes after the stress was removed (e.g. Secco et al., 2015). A stress-memory lasting longer than 5 weeks, would be long enough to potentially heat-harden the same generation of previously exposed shoots. In agriculture, hardening/priming of seeds is a long-standing practice to enhance crop resistance to environmental challenges, including hot, cold, dry, or saline conditions, or pathogen infections (Ibrahim, 2016; Wojtyla et al., 2016; Pawar and Laware, 2018). Priming refers to the plants' ability to acquire a stress memory that enhances performance under second stress exposure by responding faster, stronger, or in response to a lower threshold as compared with naïve plants (Balmer et al., 2015; Lämke and Bäurle, 2017). Molecular mechanisms involved in forming a stress memory include stalled RNA polymerase II, storage of chemical signaling factors, accumulation and phosphorylation of transcription factors, and epigenetic mechanisms such as microRNAs, histone modifications, and DNA methylation (Iwasaki and Paszkowski, 2014; Crisp et al., 2016; Hilker et al., 2016; Gallusci et al., 2017; Lämke and Bäurle, 2017). Heat-priming has only very recently been described in seagrasses (Nguyen et al., 2020). *Zostera muelleri* and *Posidonia australis* both performed better under a second heat-wave, in terms of photosynthetic capacity, leaf growth, and chlorophyll a content, when they had been previously exposed to a first heat-wave as compared with naïve
controls (Nguyen et al., 2020). This could explain why no mortality was reported for the seagrass *Posidoinia oceanica* after intense and long-lasting heat-waves in 2012, 2015, and 2017 (Darmaraki et al., 2019), although it had suffered high mortality rates after the 2006 heatwave ((Marbà and Duarte, 2010), as discussed in Nguyen et al. (2020). However, whether the experimental shift in methylation patterns in *Z. marina* may be involved in heat-priming can only be answered *via* exposure to a second heat-wave.

## Conclusions and Perspectives

Our study suggests that DNA methylation is functionally relevant for photosynthetic performance, independent from underlying somatic mutations. In seagrass meadows composed of several genotypes, stress resilience, growth, and associated invertebrate species diversity is enhanced by genotypic variation (Hughes and Stachowicz, 2004; Reusch et al., 2005; Ehlers et al., 2008). In clonal meadows, epigenetic variation may play a similar role in the potential to secure function and resilience not only of *Z. marina* plants, but also of the entire associated ecosystem.

Due to anthropogenic stressors, nearly one-third of global seagrass area has disappeared over the last 100 years, and the rate of loss accelerated from ca. 1% yr⁻¹ before 1940 to 7% yr⁻¹ since 1990 (Waycott et al., 2009). At the same time, rising temperatures open up new thermally suitable habitat in the Arctic (Krause-Jensen and Duarte, 2014). How fast and far warm-temperate and subarctic range edges will move polewards depends on the ability of seagrass to rapidly acclimate and adapt to rising temperatures and other environmental changes (Duarte et al., 2018). Thus, future studies are needed to assess the adaptive value and transgenerational stability of the epigenetic stress response, and to compare the ability to build up epigenetic variation between seagrass meadows composed of a single or multiple clones, as well as between range center *versus* edge populations.

The functional role of methylation variation in plant genets is not only of fundamental interest but also of applied interest for management programs of clonal organisms designed to assess evolutionary potential and population stability, and to minimize the loss of biodiversity. Our results can be relevant for restoration of seagrass ecosystems that largely depends on the success of replanted shoots to overcome natural variability/stress (Suykerbuyk et al., 2016; van Katwijk et al., 2016). Given that 40% of all plant species can reproduce clonally (Tiffney and Niklas, 1985), our findings are further important to other fields, such as invasion biology and crop breeding strategies (Bilchak and Kovalchuk, 2016) of clonal plants.

## REFERENCES

* Alcoverro, T., Zimmerman, R. C., Kohrs, D. G., and Alberte, R. S. (1999). Resource allocation and sucrose mobilization in light-limited elgrass Zostera marina. *Mar. Ecol. Prog. Ser.* 187, 121–131. doi: 10.3354/meps187121

* Alexa, A., and Rahnenführer, J. (2010). *Enrichment analysis for Gene Ontology. R Package.* "topGO" v2.16.0. Available at: http://www.bioconductor.org/packages/2.11/bioc/html/topGO.html.

* Arnaud-Haond, S., Duarte, C. M., Diaz-Almela, E., Marbà, N., Sintes, T., and Serrão, E. A. (2012). Implications of extreme life span in clonal organisms: Millenary clones in meadows of the threatened seagrass *Posidonia oceanica*. *PloS One* 7, e30454. doi: 10.1371/journal.pone.0030454

## DATA AVAILABILITY STATEMENT

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ **Supplementary Material**.

## AUTHOR CONTRIBUTIONS

GH (project leader) and AJ were planning the project and designing the experiments. AJ, CB, IS, and GH collected the shoots. AJ analyzed the data and wrote the manuscript. AJ, MK, AD, JC, and IS performed the DNA extraction, library preparation, and sequencing. SA-H, JO, and YP were involved in data interpretation. All authors contributed to the article and approved the submitted version.

## FUNDING

This work was supported by the Norwegian Research Council (Havkyst project 243916), the Åbo Akademi University Foundation sr to CB, and a personal research talents grant from Nord University to AJ. Open access publication fees were covered by Nord University's Opean Access fund.

## ACKNOWLEDGMENTS

We thank Marvin Choquet from Nord University for help with re-planting the shoots in the wet lab. We acknowledge Steinar Johnsen from Nord University for setting up the aquaria and water flow in the wet lab. We are grateful to the four reviewers, and the editor, Gidon Winter, whose suggestions and comments significantly improved the article. The content of this manuscript has previously appeared online as a preprint in bioRxiv (Jueterbock et al., 2019).

## SUPPLEMENTARY MATERIAL

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

* Auwera, G. A., Carneiro, M. O., Hartl, C., Poplin, R., del Angel, G., Levy-Moonshine, A., et al. (2013). From FastQ Data to High-Confidence Variant Calls: The Genome Analysis Toolkit Best Practices Pipeline. *Curr. Protoc. Bioinform.* 43, 11.10.1–11.10.33. doi: 10.1002/0471250953.bi1110s43

* Bailleul, D., Stoeckel, S., and Arnaud-Haond, S. (2016). RClone: a package to identify MultiLocus Clonal Lineages and handle clonal data sets in r. *Methods Ecol. Evol.* 7, 966–970. doi: 10.1111/2041-210X.12550

* Balmer, A., Pastor, V., Gamir, J., Flors, V., and Mauch-Mani, B. (2015). The "prime-ome": Towards a holistic approach to priming. *Trends Plant Sci.* 20, 443–452. doi: 10.1016/j.tplants.2015.04.002

* Barrett, S. C. H. (2015). Influences of clonality on plant sexual reproduction. *Proc. Natl. Acad. Sci.* 112, 8859–8866. doi: 10.1073/pnas.1501712112
Becheler, R., Diekmann, O., Hily, C., Moalic, Y., and Arnaud-Haond, S. (2010). The concept of population in clonal organisms: Mosaics of temporally colonized patches are forming highly diverse meadows of *Zostera marina* in Brittany. *Mol. Ecol.* **19**, 2394–2407. doi: **10.1111/j.1365-294X.2010.04649.x**

Benjamini, Y., and Hochberg, Y. (1995). Controlling the False Discovery Rate - a Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society Series B-Methodological 1995.pdf. *J. R. Stat. Soc. Ser. B* **57**, 289–300. doi: **10.2307/2346101**

Bewick, A. J., and Schmitz, R. J. (2017). Gene body DNA methylation in plants. *Curr. Opin. Plant Biol.* **36**, 103–110. doi: **10.1016/j.pbi.2016.12.007**

Bijlsma, R., and Loeschcke, V. (2012). Genetic erosion impedes adaptive responses to stressful environments. *Evol. Appl.* **5**, 117–129. doi: **10.1111/j.1752-4571.2011.00214.x**

Bilichak, A., and Kovalchuk, I. (2016). Transgenerational response to stress in plants and its application for breeding. *J. Exp. Bot.* **67**, 2081–2092. doi: **10.1093/jxb/erw066**

Bilichak, A., Ilnytskyy, Y., Wóycicki, R., Kepeshchuk, N., Fogen, D., and Kovalchuk, I. (2015). The elucidation of stress memory inheritance in *Brassica rapa* plants. *Front. Plant Sci.* **6**:5. doi: **10.3389/fpls.2015.00005**

Bjornstad, O. N. (2020). *ncf: Spatial Covariance Functions*. Available at: https://cran.r-project.org/package=ncf

Bossdorf, O., Richards, C. L., and Pigliucci, M. (2008). Epigenetics for ecologists. *Ecol. Lett.* **11**, 106–115. doi: **10.1111/j.1461-0248.2007.01130.x**

Boström, C., Baden, S., Bockelmann, A.-C., Dromph, K., Fredriksen, S., Gustafsson, C., et al. (2014). Distribution, structure and function of Nordic eelgrass (*Zostera marina*) ecosystems: Implications for coastal management and conservation. *Aquat. Conserv. Mar. Freshw. Ecosyst.* **24**, 410–434. doi: **10.1002/aqc.2424**

Boyko, A., Blevins, T., Yao, Y., Golubov, A., Bilichak, A., Ilnytskyy, Y., et al. (2010). Transgenerational adaptation of *Arabidopsis* to stress requires DNA methylation and the function of dicer-like proteins. *PloS One* **5**, e9514. doi: **10.1371/journal.pone.0009514**

Bussotti, F., Desotgiu, R., Pollastrini, M., and Cascio, C. (2010). The JIP test: a tool to screen the capacity of plant adaptation to climate change. *Scand. J. For. Res.* **25**, 43–50. doi: **10.1080/02827581.2010.485777**

Chefaoui, R. M., Duarte, C. M., and Serrão, E. A. (2018). Dramatic loss of seagrass habitat under projected climate change in the Mediterranean Sea. *Glob. Change Biol.* **24** (10), 4919–4928. doi: **10.1111/gcb.14401**

Costanza, R., Arge, R., Groot, R., De, Farberk, S., Grasso, M., Hannon, B., et al. (1997). The value of the world's ecosystem services and natural capital. *Nature* **387**, 253–260. doi: **10.1038/387253a0**

Costanza, R., de Groot, R., Sutton, P., van der Ploeg, S., Anderson, S. J., Kubiszewski, I., et al. (2014). Changes in the global value of ecosystem services. *Glob. Environ. Change* **26**, 152–158. doi: **10.1016/j.gloenvcha.2014.04.002**

Crisp, P. A., Ganguly, D., Eichten, S. R., Borevitz, J. O., and Fogson, B. J. (2016). Reconsidering plant memory: Intersections between stress recovery, RNA turnover, and epigenetics. *Sci. Adv.* **2**, e1501340–e1501340. doi: **10.1126/sciadv.1501340**

Darmaraki, S., Somot, S., Sevault, F., and Nabat, P. (2019). Past Variability of Mediterranean Sea Marine Heatwaves. *Geophys. Res. Lett.* **46**, 9813–9823. doi: **10.1029/2019GL082933**

Dennison, W. C., and Alberte, R. S. (1982). Photosynthetic Responses of *Zostera marina* L. (Eelgrass) to in situ Manipulations of Light Intensity. *Oecologia* **55**, 137–144. doi: **10.1007/BF00384478**

Dodd, R. S., and Douhovnikoff, V. (2016). Adjusting to Global Change through Clonal Growth and Epigenetic Variation. *Front. Ecol. Evol.* **4**:86. doi: **10.3389/fevo.2016.00086**

Douhovnikoff, V., and Dodd, R. S. (2014). Epigenetics: a potential mechanism for clonal plant success. *Plant Ecol.* **216**, 227–233. doi: **10.1007/s11258-014-0430-z**

Dowen, R. H., Pelizzola, M., Schmitz, R. J., Lister, R., Dowen, J. M., Nery, J. R., et al. (2012). Widespread dynamic DNA methylation in response to biotic stress. *Proc. Natl. Acad. Sci. U. S. A.* **109**, E2183–E2191. doi: **10.1073/pnas.1209329109**

Duarte, C. M., Kennedy, H., Marbà, N., and Hendriks, I. (2013). Assessing the capacity of seagrass meadows for carbon burial: Current limitations and future strategies. *Ocean Coast. Manage.* **83**, 32–38. doi: **10.1016/j.ocecoaman.2011.09.001**

Duarte, B., Martins, I., Rosa, R., Matos, A. R., Roleda, M. Y., Reusch, T. B. H., et al. (2018). Climate change impacts on seagrass meadows and macroalgal forests:

an integrative perspective on acclimation and adaptation potential. *Front. Mar. Sci.* **5**:190-190. doi: **10.3389/FMARS.2018.00190**

Dubin, M. J., Zhang, P., Meng, D., Remigereau, M. S., Osborne, E. J., Casale, F. P., et al. (2015). DNA methylation in *Arabidopsis* has a genetic basis and shows evidence of local adaptation. *Elife* **4**, 1–23. doi: **10.7554/eLife.05255**

Ehlers, A., Worm, B., and Reusch, T. B. H. (2008). Importance of genetic diversity in eelgrass *Zostera marina* for its resilience to global warming. *Mar. Ecol. Prog. Ser.* **355**, 1–7. doi: **10.3354/meps07369**

Feder, M. E., and Hofmann, G. E. (2002). Heat-shock proteins, molecular chaperones, and the stress response: Evolutionary and ecological physiology. *Annu. Rev. Physiol.* **61**, 243–282. doi: **10.1146/annurev.physiol.61.1.243**

Feng, S., Cokus, S. J., Zhang, X., Chen, P.-Y. P.-Y., Bostick, M., Goll, M. G., et al. (2010). Conservation and divergence of methylation patterning in plants and animals. *Proc. Natl. Acad. Sci.* **107**, 8689–8694. doi: **10.1073/pnas.1002720107**

Fourqurean, J. W., Duarte, C. M., Kennedy, H., Marbà, N., Holmer, M., Mateo, M. A., et al. (2012). Seagrass ecosystems as a globally significant carbon stock. *Nat. Geosci.* **5**, 505–509. doi: **10.1038/ngeo1477**

Foust, C. M., Preite, V., Schrey, A. W., Alvarez, M., Robertson, M. H., Verhoeven, K. J. F. F., et al. (2016). Genetic and epigenetic differences associated with environmental gradients in replicate populations of two salt marsh perennials. *Mol. Ecol.* **25**, 1639–1652. doi: **10.1111/mec.13522**

Frankham, R., Bradshaw, C. J. A., and Brook, B. W. (2014). Genetics in conservation management: Revised recommendations for the 50/500 rules, Red List criteria and population viability analyses. *Biol. Conserv.* **170**, 56–63. doi: **10.1016/j.biocon.2013.12.036**

Franssen, S. U., Gu, J., Bergmann, N., Winters, G., Klostermeier, U. C., Rosenstiel, P., et al. (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. doi: **10.1073/pnas.1107680108**

Gabriel, W., Lynch, M., and Bürger, R. (1993). Muller's Ratchet and Mutational Meltdowns. *Evol. (N. Y.)* **47**, 1744–1757. doi: **10.1111/j.1558-5646.1993.tb01266.x**

Gallusci, P., Dai, Z., Génard, M., Gauffretau, A., Leblanc-Fournier, N., Richard-Molard, C., et al. (2017). Epigenetics for Plant Improvement: Current Knowledge and Modeling Avenues. *Trends Plant Sci.* **22**, 610–623. doi: **10.1016/j.tplants.2017.04.009**

Gao, L., Geng, Y., Li, B., Chen, J., and Yang, J. (2010). Genome-wide DNA methylation alterations of *Alternanthera philoxeroides* in natural and manipulated habitats: implications for epigenetic regulation of rapid responses to environmental fluctuation and phenotypic variation. *Plant Cell Environ.* **33**, 1820–1827. doi: **10.1111/j.1365-3040.2010.02186.x**

Garrison, E. (2020). C++ library and cmdline tools for parsing and manipulating VCF files. Available at: https://github.com/vcflib/vcflib.

Gáspár, B., Bossdorf, O., and Durka, W. (2019). Structure, stability and ecological significance of natural epigenetic variation: a large-scale survey in *Plantago lanceolata*. *New Phytol.* **221**, 1585–1596. doi: **10.1111/nph.15487**

Greally, J. M. (2018). A user's guide to the ambiguous word "epigenetics." *Nat. Rev. Mol. Cell Biol.* **19**, 207–208. doi: **10.1038/nrm.2017.135**

Green, E. P., and Short, F. (2004). World atlas of seagrasses (Berkeley, USA: University of California Press). doi: **10.1515/BOT.2004.029**

Greco, M., Chiappetta, A., Bruno, L., and Bitonti, M. B. (2012). In Posidonia oceanica cadmium induces changes in DNA methylation and chromatin patterning. *J. Exp. Bot.* **63**, 695–709. doi: **10.1093/jxb/err313**

Greco, M., Chiappetta, A., Bruno, L., and Bitonti, M. B. (2013). Effects of light deficiency on genome methylation in Posidonia oceanica. *Mar. Ecol. Prog. Ser.* **473**, 103–114. doi: **10.3354/meps09955**

Greve, T. M., Borum, J., and Pedersen, O. (2003). Meristematic oxygen variability in eelgrass (*Zostera marina*). *Limnol. Oceanogr.* **48**, 210–216. doi: **10.4319/lo.2003.48.1.2010**

Griffin, P. T., Niederhuth, C. E., and Schmitz, R. J. (2016). A comparative analysis of 5-azacytidine-and zebularine-induced DNA demethylation. *G3 Genes Genomes Genet.* **6**, 2773–2780. doi: **10.1534/g3.116.030262**

Gugger, P. F., Fitz-Gibbon, S., PellEgrini, M., and Sork, V. L. (2016). Species-wide patterns of DNA methylation variation in *Quercus lobata* and their association with climate gradients. *Mol. Ecol.* **25**, 1665–1680. doi: **10.1111/mec.13563**

Herrera, C. M., and Bazaga, P. (2010). Epigenetic differentiation and relationship to adaptive genetic divergence in discrete populations of the violet *Viola cazorlensis*. *New Phytol.* **187**, 867–876. doi: **10.1111/j.1469-8137.2010.03298.x**
Hilker, M., Schwachtje, J., Baier, M., Balazadeh, S., Bäurle, I., Geiselhardt, S., et al. (2016). Priming and memory of stress responses in organisms lacking a nervous system. *Biol. Rev.* 91, 1118–1133. doi: 10.1111/brv.12215

Hughes, A. R., and Stachowicz, J. J. (2004). Genetic diversity enhances the resistance of a seagrass ecosystem to disturbance. *Proc. Natl. Acad. Sci.* 101, 8998–9002. doi: 10.1073/pnas.0402642101

Hughes, A. R., Inouye, B. D., Johnson, M. T. J., Underwood, N., and Vellend, M. (2008). Ecological consequences of genetic diversity. *Ecol. Lett.* 11, 609–623. doi: 10.1111/j.1461-0248.2008.01179.x

Ibrahim, E. A. (2016). Seed priming to alleviate salinity stress in germinating seeds. *J. Plant Physiol.* 192, 38–46. doi: 10.1016/j.jplph.2015.12.011

Iwasaki, M., and Paszkowski, J. (2014). Epigenetic memory in plants. *EMBO J.* 33, 1987–1998. doi: 10.15252/embj.201488883

Jueterbock, A., Franssen, S. U., Bergmann, N., Gu, J., Coyer, J. A., Reusch, T. B. H. H., et al. (2016). Phylogeographic differentiation versus transcriptomic adaptation to warm temperatures in *Zostera marina*, a globally important seagrass. *Mol. Ecol.* 25, 5396–5411. doi: 10.1111/mec.13829

Jueterbock, A., Boström, C., James, A. C., Olsen, J., Kopp, M., Dhanasiri, A., et al. (2019). Methylation variation promotes phenotypic diversity and evolutionary potential in a millenium-old clonal seagrass meadow. *bioRxiv* 787754. doi: 10.1101/787754

Kilvitis, H. J., Alvarez, M., Foust, C. M., Schrey, A. W., Robertson, M., and Richards, C. L. (2014). "Ecological Epigenetics," in *Ecological Genomics: Ecology and the Evolution of Genes and Genomes Advances in Experimental Medicine and Biology*. Eds. C. R. Landry and N. Aubin-Horth (Dordrecht: Springer Netherlands) 191–210. doi: 10.1007/978-94-007-7347-9

Knaus, B. J., and Grünwald, N. J. (2017). vcfr : a package to manipulate and visualize variant call format data in R. *Mol. Ecol. Resour.* 17, 44–53. doi: 10.1111/1755-0998.12549

Krause-Jensen, D., and Duarte, C. M. (2014). Expansion of vegetated coastal ecosystems in the future Arctic. *Front. Mar. Sci.* 1:77. doi: 10.3389/fmars. 2014.00077

Laikre, L., Allendorf, F. W., Aroner, L. C., Baker, C. S., Gregovich, D. P., Hansen, M. M., et al. (2009). Neglect of genetic diversity in implementation of the Convention on Biological Diversity. *Conserv. Biol.* 24, 86–88. doi: 10.1111/j.1523-1739.2009.01425.x

Lämke, J., and Bäurle, I. (2017). Epigenetic and chromatin-based mechanisms in environmental stress adaptation and stress memory in plants. *Genome Biol.* 18, 1–11. doi: 10.1186/s13059-017-1263-6

Lande, R. (1988). Genetics and demography in biological conservation. *Sci. (80.-).* 241, 1455–1460. doi: 10.1126/science.3420403

Larkum, A. W. D. D., Orth, R. J., and Duarte, C. M. (2006). *Seagrasses: Biology, Ecology, and Conservation*. Eds. A. W. D. D. Larkum, R. J. Orth and C. M. Duarte (Dordrecht: Springer Netherlands). doi: 10.1111/j.1439-0485.2006.00138.x

Latzel, V., Allan, E., Bortolini Silveira, A., Colot, V., Fischer, M., and Bossdorf, O. (2013). Epigenetic diversity increases the productivity and stability of plant populations. *Nat. Commun.* 4, 2875. doi: 10.1038/ncomms3875

Latzel, V., Rendina González, A. P., and Rosenthal, J. (2016). Epigenetic Memory as a Basis for Intelligent Behavior in Clonal Plants. *Front. Plant Sci.* 7:1354:1354. doi: 10.3389/fplms.2016.01354

Lé, S., Josse, J., and Husson, F. (2008). FactoMineR: An R Package for Multivariate Analysis. *J. Stat. Software* 25, 1–18. doi: 10.1016/j.envint.2008.06.007

Li, H., and Durbin, R. (2009). Fast and accurate short-read alignment with Burrows-Wheeler transform. *Bioinformatics* 25, 1754–1760. doi: 10.1093/bioinformatics/btp324

Li, R., Li, Y., Kristiansen, K., and Wang, J. (2008). SOAP: Short oligonucleotide alignment program. *Bioinformatics* 24, 713–714. doi: 10.1093/bioinformatics/btn025

Lisch, D. (2009). Epigenetic Regulation of Transposable Elements in Plants. *Annu. Rev. Plant Biol.* 60, 43–66. doi: 10.1146/annurev.arplant.59.032607.092744

Liu, J., Feng, L., Li, J., and He, Z. (2015). Genetic and epigenetic control of plant heat responses. *Front. Plant Sci.* 06, 267. doi: 10.3389/fplms.2015.00267

Liu, F., Liu, J., and Dong, M. (2016). Ecological Consequences of Clonal Integration in Plants. *Front. Plant Sci.* 7, 770. doi: 10.3389/fplms.2016.00770

Liu, G., Xia, Y., Liu, T., Dai, S., and Hou, X. (2018). The DNA methy lome and association of differentially methylated regions with differential gene expression during heat stress in *Brassica rapa*. *Int. J. Mol. Sci.* 19, 1414. doi: 10.3390/ijms19051414

Loxdale, H. D., Lushai, G., Aafc, A. C., Cereal, E., and Ecorc, C. (2003). Rapid changes in clonal lines: the death of a 'sacred cow'. *Biol. J. Linn. Soc.* 79 (1), 3–16. doi: 10.1046/j.1095-8312.2003.00177.x

Lushai, G., Loxdale, H. D., and Allen, J. A. (2003). The dynamic clonal genome and its adaptive potential. *Biol. J. Linn. Soc.* 79 (1), 193–208. doi: 10.1046/j.1095-8312.2003.00189.x

Lynch, M., and Lande, R. (1998). The critical effective size for a genetically secure population. *Anim. Conserv.* 1, 70–72. doi: 10.1017/S136794309822110X

Lynch, M., Butcher, D., Bürger, R., and Gabriel, W. (1993). The mutational meltdown in asexual populations. *J. Hered.* 84, 339–344. doi: 10.1093/oxfordjournals.jhered.a111354

Marbà, N., and Duarte, C. M. (2010). Mediterranean warming triggers seagrass (*Posidonia oceanica*) shoot mortality. *Glob. Change Biol.* 16, 2366–2375. doi: 10.1111/j.1365-2486.2009.02130.x

Marín-Guirao, L., Entrambasaguas, L., Ruiz, J. M., and Procaccini, G. (2019). Heat-stress induced flowering can be a potential adaptive response to ocean warming for the iconic seagrass Posidonia oceanica. *Mol. Ecol.* 28, 2486–2501. doi: 10.1111/mec.15089

May, M. R., Provance, M. C., Sanders, A. C., Ellstrand, N. C., and Ross-Ibarra, J. (2009). A pleistocene clone of palmer's oak persisting in Southern California. *PloS One* 4, 1–5. doi: 10.1371/journal.pone.0008346

McKey, D., Elias, M., Pujol, M. E., and Duputié, A. (2010). The evolutionary ecology of clonally propagated domesticated plants. *New Phytol.* 186, 318–332. doi: 10.1111/j.1469-8137.2010.03210.x

Medrano, M. M., Herrera, C. M., and Bazaga, P. (2014). Epigenetic variation predicts regional and local intraspecific functional diversity in a perennial herb. *Mol. Ecol.* 23, 4926–4938. doi: 10.1111/mec.12911

Muller, H. J. (1964). The relation of recombination to mutational advance. *Mutat. Res. - Fundam. Mol. Mech. Mutagen.* 1, 2–9. doi: 10.1016/0027-5107(64)90047-8

Neri, F., Rapelli, S., Krepelova, A., Incarnato, D., Parlato, C., Basile, G., et al. (2017). Intragenic DNA methylation prevents spurious transcription initiation. *Nature* 543, 72–77. doi: 10.1038/nature21373

Nguyen, H. M., Kim, M., Ralph, P. J., Marín-Guirao, L., Pernice, M., and Procaccini, G. (2020). Stress memory in seagrasses: first insight into the effects of thermal priming and the role of epigenetic modifications. *Front. Plant Sci.* 11:494:494. doi: 10.3389/FPLS.2020.00049

Niederhuth, C. E., and Schmitz, R. J. (2017). Putting DNA methylation in context: from genomes to gene expression in plants. *Biochim. Biophys. Acta - Gene Regul. Mech.* 1860, 149–156. doi: 10.1016/j.bbagrm.2016.08.009

Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). *vegan: Community Ecology Package*. Available at: https://cran.r-project.org/package=vegan.

Olsen, J. L., Stam, W. T., Coyer, J. A., Reusch, T. B. H. H., Billingham, M., Boström, C., et al. (2004). North Atlantic phylogeography and large-scale population differentiation of the seagrass *Zostera marina* L. *Mol. Ecol.* 13, 1923–1941. doi: 10.1111/j.1365-294X.2004.02205.x

Olsen, J. L., Rouzé, P., Verhelst, B., Lin, Y., Bayer, T., Collen, J., et al. (2016). The genome of the seagrass *Zostera marina* reveals angiosperm adaptation to the sea. *Nature* 530, 331–335. doi: 10.1038/nature16548

Orth, R. J., Carruthers, T. J. B., Dennison, W. C., Duarte, C. M., Fourqurean, J. W., Heck, K. L., et al. (2006). A Global Crisis for Seagrass Ecosystems. *Bioscience* 56, 987–996. doi: 10.1641/0006-3568(2006)56[987:agcfse]2.0.co;2

Ou, X., Zhang, Y., Xu, C., Lin, X., Zang, Q., Zhuang, T., et al. (2012). Transgenerational Inheritance of Modified DNA Methylation Patterns and Enhanced Tolerance Induced by Heavy Metal Stress in Rice (*Oryza sativa* L.). *PloS One* 7, e41143. doi: 10.1371/journal.pone.0041143

Pawar, V. A., and Laware, S. L. (2018). Seed Priming to A Critical Review. *Int. J. Sci. Res. Biol. Sci.* 5, 94–101. doi: 10.26438/ijsrbs/v5i5.94101

Procaccini, G., Olsen, J. L., and Reusch, T. B. H. (2007). Contribution of genetics and genomics to seagrass biology and conservation. *J. Exp. Mar. Bio. Ecol.* 350, 234–259. doi: 10.1016/j.jembe.2007.05.035

Ptashne, M. (2007). On the use of the word "epigenetic." *Curr. Biol.* 17, 233–236. doi: 10.1016/j.cub.2007.02.030

R Core Team (2019). *R: A Language and Environment for Statistical Computing*. Available at: https://www.r-project.org/.

Repolho, T., Duarte, B., Dionisio, G., Paula, J. R., Lopes, A. R., Rosa, I. C., et al. (2017). Seagrass ecophysiological performance under ocean warming and acidification. *Sci. Rep.* 7, 41443. doi: 10.1038/srep41443
Reusch, T. B. H., and Boström, C. (2011). Widespread genetic mosaicism in the marine angiosperm *Zostera marina* is correlated with clonal reproduction. *Evol. Ecol.* 25, 899–913. doi: 10.1007/s10682-010-9436-8

Reusch, T. B. H., and Hughes, A. R. (2006). The Emerging Role of Genetic Diversity for Ecosystem Functioning: Estuarine Macrophytes as Models. *Estuaries Coasts* 29, 159–164. doi: 10.1007/BF02784707

Reusch, T. B. H., Boström, C., Stam, W. T., and Olsen, J. L. (1999). An ancient eelgrass clone in the Baltic. *Mar. Ecol. Prog. Ser.* 183, 301–304. doi: 10.3354/meps183301

Reusch, T. B. H., Stam, W. T., and Olsen, J. L. (2000). A microsatellite-based estimation of clonal diversity and population subdivision in *Zostera marina*, a marine flowing plant. *Mol. Ecol.* 9, 127–140. doi: 10.1046/j.1365-294x.2000.00839.x

Reusch, T. B. H., Ehlers, A., Hammerli, A., and Worm, B. (2005). Ecosystem recovery after climatic extremes enhanced by genotypic diversity. *Proc. Natl. Acad. Sci.* 102, 2826–2831. doi: 10.1073/pnas.0500008102

Reusch, T. B. H., Dierking, J., Andersson, H. C., Bonsdorff, E., Carstensen, J., Casini, M., et al. (2018). The Baltic Sea as a time machine for the future coastal ocean. *Sci. Adv.* 4, eaar8195. doi: 10.1126/sciadv.aar8195

Richards, C. L., Schrey, A. W., and Pigliucci, M. (2012). Invasion of diverse habitats by few Japanese knotweed genotypes is correlated with epigenetic differentiation. *Ecol. Lett.* 15, 1016–1025. doi: 10.1111/j.1461-0248.2012.01824.x

Robinson, M. D., and Oshlack, A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data. *Genome Biol.* 11, R25. doi: 10.1186/gb-2010-11-3-r25

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2009). edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics* 26, 139–140. doi: 10.1093/bioinformatics/btp616

Ruocco, M., De Luca, P., Marin-Guirao, L., and Procaccini, G. (2019a). Differential Leaf Age-Dependent Thermal Plasticity in the Keystone Seagrass Posidonia oceanica. *Front. Plant Sci.* 10, 1556. doi: 10.3389/fpls.2019.01556

Ruocco, M., Marin-Guirao, L., and Procaccini, G. (2019b). Within- and among-leaf variations in photo-physiological functions, gene expression and DNA methylation patterns in the large-sized seagrass Posidonia oceanica. *Mar. Biol.* 166, 1–18. doi: 10.1007/s00227-019-3482-8

Santelices, B., Gallegos Sánchez, C., and González, A. V. (2018). Intraorganismal genetic heterogeneity as a source of genetic variation in modular macroalgae. *J. Phycol.* 54, 767–771. doi: 10.1111/jpy.12784

Schlicker, A., Domingues, F., Rahnenfuhrer, J., and Lengauer, T. (2006). A new measure for functional similarity of gene products based on Gene Ontology. *BMC Bioinf.* 7, 302. doi: 10.1186/1471-2105-7-302

Schmitz, R. J., Schultz, M. D., Lewsey, M. G., O'Malley, R. C., Urich, M., Libiger, O., et al. (2011). Transgenerational epigenetic instability is a source of novel methylation variants. *Science* 334, 369–373. doi: 10.1126/science.1212959

Schmitz, R. J., Schultz, M. D., Urich, M. A., Nery, J. R., Pelizzola, M., Libiger, O., et al. (2013). Patterns of population epigenomic diversity. *Nature* 495, 193–198. doi: 10.1038/nature11968

Schulz, B., Eckstein, R. L., and Durka, W. (2014). Epigenetic variation reflects dynamic habitat conditions in a rare floodplain herb. *Mol. Ecol.* 23, 3523–3537. doi: 10.1111/mec.12835

Secco, D., Wang, C., Shou, H., Schultz, M. D., Chiarenza, S., Nussaume, L., et al. (2015). Stress induced gene expression drives transient DNA methylation changes at adjacent repetitive elements. *Elife* 4, e09343. doi: 10.7554/eLife.09343.001

Seymour, D. K., Koenig, D., Hagmann, J., Becker, C., and Weigel, D. (2014). Evolution of DNA Methylation Patterns in the Brassicaceae is Driven by Differences in Genome Organization. *PloS Genet.* 10, e1004785. doi: 10.1371/journal.pgen.1004785

Silvertown, J. (2008). The Evolutionary Maintenance of Sexual Reproduction: Evidence from the Ecological Distribution of Asexual Reproduction in Clonal Plants. *Int. J. Plant Sci.* 169, 157–168. doi: 10.1086/523357

Simberloff, D., and Leppanen, C. (2019). Plant somatic mutations in nature conferring insect and herbicide resistance. *Pest Manage. Sci.* 75, 14–17. doi: 10.1002/ps.5157

Spielman, D., Brook, B. W., and Frankham, R. (2004). Most species are not driven to extinction before genetic factors impact them. *Proc. Natl. Acad. Sci. U. S. A.* 101, 15261–15264. doi: 10.1073/pnas.0403809101

Stefanov, D., Petkova, V., and Denev, I. D. (2011). Screening for heat tolerance in common bean (*Phaseolus vulgaris* L.) lines and cultivars using JIP-test. *Sci. Hortic.* (*Amsterdam*) 128, 1–6. doi: 10.1016/j.scienta.2010.12.003

Sterck, L., Billiau, K., Abeel, T., Rouzé, P., and Van de Peer, Y. (2012). ORCAE: online resource for community annotation of eukaryotes. *Nat. Methods* 9, 1041. doi: 10.1038/nmeth.2242

Strasser, R. J., Srivastava, A., and Tsimilli-Michael, M. (2000). “The fluorescence transient as a tool to characterize and screen photosynthetic samples,” in *Probing Photosynthesis: Mechanism, Regulation & Adaptation*. Eds. M. Yunu, U. Pathre and P. Mohanty (London and New York: Taylor and Francis), 445–483. Available at: http://www.hansatech-instruments.com/docs/thefluorescencetransient.pdf.

Supek, F., Bošnjak, M., Škunca, N., and Šmuc, T. (2011). Revigo summarizes and visualizes long lists of gene ontology terms. *PloS One* 6, e21800. doi: 10.1371/journal.pone.0021800

Suykerbuyk, W., Govers, L. L., Bouma, T. J., Giesen, W. B. J. T., de Jong, D. J., van de Voort, R., et al. (2016). Unpredictability in seagrass restoration: analysing the role of positive feedback and environmental stress on *Zostera noltii* transplants. *J. Appl. Ecol.* 53, 774–784. doi: 10.1111/1365-2664.12614

Takuno, S., and Gaut, B. S. (2012). Body-methylated genes in *Arabidopsis thaliana* are functionally important and evolve slowly. *Mol. Biol. Evol.* 29, 219–227. doi: 10.1093/molbev/msr188

Tiffney, B., and Niklas, K. (1985). “Clonal growth in land plants - a paleobotanical perspective,” in *Population biology and evolution of clonal organisms*. Eds. R. E. Cook, J. B. C. Jackson and L. W. Buss (New Haven: Yale University Press).

van der Graaf, A., Wardenaar, R., Neumann, D. A., Taudt, A., Shaw, R. G., Jansen, R. C., et al. (2015). Rate, spectrum, and evolutionary dynamics of spontaneous epimutations. *Proc. Natl. Acad. Sci. U. S. A.* 112, 6676–6681. doi: 10.1073/pnas.1424254112

van Katwijk, M. M., Thorhaug, A., Marbà, N., Orth, R. J., Duarte, C. M., Kendrick, G. A., et al. (2016). Global analysis of seagrass restoration: The importance of large-scale planting. *J. Appl. Ecol.* 53, 567–578. doi: 10.1111/1365-2664.12562

Vanden Broeck, A., Cox, K., Brys, R., Castiglione, S., Cicatelli, A., Guarino, F., et al. (2018). Variability in DNA Methylation and Generational Plasticity in the Lombardy Poplar, a Single Genotype Worldwide Distributed Since the Eighteenth Century. *Front. Plant Sci.* 9:1635. doi: 10.3389/fpls.2018.01635

Varet, H., Brillet-Guéguen, L., Coppée, J. Y., and Dillies, M. A. (2016). SARTools: A DESeq2- and edgeR-based R pipeline for comprehensive differential analysis of RNA-Seq data. *PloS One* 11, 1–8. doi: 10.1371/journal.pone.0157022

Verhoeven, K. J. F., and Preite, V. (2014). Epigenetic variation in asexually reproducing organisms. *Evol. (N. Y.)* 68, 644–655. doi: 10.1111/evo.12320

Verhoeven, K. J. F., and van Gurp, T. P. (2012). Transgenerational effects of stress exposure on offspring phenotypes in apomictic dandelion. *PloS One* 7, e38605. doi: 10.1371/journal.pone.0038605

Verhoeven, K. J. F., Jansen, J. J., van Dijk, P. J., and Biere, A. (2010). Stress-induced DNA methylation changes and their heritability in asexual dandelions. *New Phytol.* 185, 1108–1118. doi: 10.1111/j.1469-8137.2009.03121.x

Verhoeven, K. J. F., Verbon, E. H., van Gurp, T. P., Oplaat, C., Ferreira de Carvalho, J., Morse, A. M., et al. (2018). Intergenerational environmental effects: functional signals in offspring transcriptomes and metabolomes after parental jasmonic acid treatment in apomictic dandelion. *New Phytol.* 217, 871–882. doi: 10.1111/nph.14835

Wang, S., Lv, J., Zhang, L., Dou, J., Sun, Y., Li, X., et al. (2015). MethylRAD: a simple and scalable method for genome-wide DNA methylation profiling using methylation-dependent restriction enzymes. *Open Biol.* 5:150130. doi: 10.1098/rsob.150130

Wang, M. Z., Li, H. L., Li, J. M., and Yu, F. H. (2020). Correlations between genetic, epigenetic and phenotypic variation of an introduced clonal herb. *Heredity* (*Edinb.*), 124, 146–155. doi: 10.1038/s41437-019-0261-8

Waycott, M., Duarte, C. M., Carruthers, T. J. B., Orth, R. J., Dennison, W. C., Olyarnik, S., et al. (2009). Accelerating loss of seagrasses across the globe threatens coastal ecosystems. *Proc. Natl. Acad. Sci. U. S. A.* 106, 12377–12381. doi: 10.1073/pnas.0905620106

Whitham, T. G., and Slobodchikoff, C. N. (1981). Evolution by individuals, plant-herbivore interactions, and mosaics of genetic variability: The adaptive significance of somatic mutations in plants. *Oecologia* 49, 287–292. doi: 10.1007/BF00347587
Wickham, H. (2016). *ggplot2* *Elegant Graphics for Data Analysis* (New York: Springer New York). doi: 10.1007/978-0-387-98141-3

Williams, B. P., and Gehring, M. (2017). Stable transgenerational epigenetic inheritance requires a DNA methylation-sensing circuit. *Nat. Commun.* 8, 2124. doi: 10.1038/s41467-017-02219-3

Wojtyla, Ł., Lechowska, K., Kubala, S., and Garnczarska, M. (2016). Molecular processes induced in primed seeds—increasing the potential to stabilize crop yields under drought conditions. *J. Plant Physiol.* 203, 116–126. doi: 10.1016/j.jplph.2016.04.008

Xu, X., Tao, Y., Gao, X., Zhang, L., Li, X., Zou, W., et al. (2016). A CRISPR-based approach for targeted DNA demethylation. *Cell Discov.* 2, 16009. doi: 10.1038/celldisc.2016.9

Yang, X., Han, H., DeCarvalho, D. D., Lay, F. D., Jones, P. A., and Liang, G. (2014). Gene body methylation can alter gene expression and is a therapeutic target in cancer. *Cancer Cell* 26, 577–590. doi: 10.1016/j.ccr.2014.07.028

Yu, L., Boström, C., Franzenburg, S., Bayer, T., Dagan, T., and Reusch, T. B. H. (2020). Somatic genetic drift and multilevel selection in a clonal seagrass. *Nat. Ecol. Evol.* 3, 833335. doi: 10.1038/s41559-020-1196-4

Zhang, X., Yazaki, J., Sundaresan, A., Cokus, S., Chan, S. W. L., Chen, H., et al. (2006). Genome-wide High-Resolution Mapping and Functional Analysis of DNA Methylation in *Arabidopsis*. *Cell* 126, 1189–1201. doi: 10.1016/j.cell.2006.08.003

Zhang, Y.-Y. Y., Fischer, M., Colot, V., and Bossdorf, O. (2013). Epigenetic variation creates potential for evolution of plant phenotypic plasticity. *New Phytol.* 197, 314–322. doi: 10.1111/nph.12010

Zhang, H., Lang, Z., and Zhu, J. K. (2018). Dynamics and function of DNA methylation in plants. *Nat. Rev. Mol. Cell Biol.* 19, 489–506. doi: 10.1038/ s41580-018-0016-z

Živčák, M., Brestič, M., Olšovská, K., and Slamka, P. (2008). Performance index as a sensitive indicator of water stress in *Triticum aestivum* L. *Plant Soil Environ.* 54, 133–139. doi: 10.17221/392-PSE

**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.

Copyright © 2020 Jueterbock, Boström, Coyer, Olsen, Kopp, Dhanasiri, Smolina, Arnaud-Haond, Van de Peer and Hoarau. 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.