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

HAL
open science

# DNA methylation dynamics during stress response in woodland strawberry (Fragaria vesca)

María-Estefanía López, David Roquis, Claude Becker, Béatrice Denoyes, Etienne Bucher

► To cite this version:

María-Estefanía López, David Roquis, Claude Becker, Béatrice Denoyes, Etienne Bucher. DNA methylation dynamics during stress response in woodland strawberry (Fragaria vesca). Horticulture research, 2022, online first, 10.1093/hr/uhac174. hal-03791445

HAL Id: hal-03791445

https://hal.inrae.fr/hal-03791445

Submitted on 29 Sep 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
DNA methylation dynamics during stress response in woodland strawberry (Fragaria vesca)

María-Estefanía López¹,², David Roquis¹, Claude Becker³, Béatrice Denoyes⁴ and Etienne Bucher¹*

¹Crop Genome Dynamics Group, Agroscope, 1260 Nyon, Switzerland, ²Department of Botany and Plant Biology, Faculty of Sciences, University of Geneva, Geneva, Switzerland, ³LMU BioCenter, Faculty of Biology, Ludwig-Maximilians-University Munich, D-82152 Martinsried, Germany and ⁴Univ. Bordeaux, INRAE, Biologie du Fruit et Pathologie, F-33140 Villenave d'Ornon, France

Running title: Stress-induced epigenetic changes in strawberry

Correspondence:

Etienne Bucher

E-mail: etienne.bucher@agroscope.admin.ch

Tel.: +41 58 483 97 53

Mobile: +41 79 158 54 65

María-Estefanía López (ORCID ID: 0000-0001-9068-4909; maria.lopez.1@etu.unige.ch)

David Roquis (ORCID ID: 0000-0001-5265-3132; david.roquis@agroscope.admin.ch)

Claude Becker (ORCID ID: 0000-0003-3406-4670; claude.becker@biologie.uni-muenchen.de)

Béatrice Denoyes (ORCID ID: 0000-0002-0369-9609; beatrice.denoyes@inrae.fr)

Etienne Bucher (ORCID ID: 0000-0002-3114-3763; etienne.bucher@agroscope.admin.ch)

© The Author(s) 2022. Published by Oxford University Press. All rights reserved. This is an Open Access article distributed under the terms of the Creative Commons Attribution License https://creativecommons.org/licenses/by/4.0/, which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
<table><tr><td>*Total word count (excluding summary, acknowledgements, data availability, author contribution, references, and legends)</td><td>6949</td><td>No. of figures:</td><td>6 (all in color)</td></tr><tr><td>Summary:</td><td>203</td><td>No. of Tables:</td><td>1</td></tr><tr><td>*Introduction:</td><td>865</td><td>No. of Supporting Information files:</td><td>(Fig. S1-S8; Table S1-S9)</td></tr><tr><td>*Materials and Methods:</td><td>1700</td><td></td><td></td></tr><tr><td>*Results:</td><td>2515</td><td></td><td></td></tr><tr><td>*Discussion:</td><td>1692</td><td></td><td></td></tr><tr><td>*Conclusions</td><td>177</td><td></td><td></td></tr><tr><td>Acknowledgements:</td><td>99</td><td></td><td></td></tr><tr><td>Data availability</td><td>59</td><td></td><td></td></tr><tr><td>Author Contribution:</td><td>56</td><td></td><td></td></tr></table>

## Summary

Environmental stresses can result in a wide range of physiological and molecular responses in plants. These responses can also impact epigenetic information in genomes, especially at the level of DNA methylation (5-methylcytosine). DNA methylation is the hallmark heritable epigenetic modification and plays a key role in silencing transposable elements (TEs). Although DNA methylation is an essential epigenetic mechanism, fundamental aspects of its contribution to stress responses and adaptation remain obscure.

We investigated epigenome dynamics of wild strawberry (*Fragaria vesca*) in response to variable ecologically relevant environmental conditions at the DNA methylation level. *F. vesca* methylome responded with great plasticity to ecologically relevant abiotic and hormonal stresses. Thermal stress resulted in substantial genome-wide loss of DNA methylation. Notably, all tested stress conditions resulted in marked hot spots of differential DNA methylation near centromeric or pericentromeric regions, particularly in the non-symmetrical DNA methylation context. Additionally, we identified differentially methylated regions (DMRs) within promoter regions of transcription factor (TF) superfamilies involved in plant stress-response and assessed the effects of these changes on gene expression.

These findings improve our understanding on stress-response at the epigenome level by highlighting the correlation between DNA methylation, TEs and gene expression regulation in plants subjected to a broad range of environmental stresses.
**Keywords:**

epigenetics, stress response, transposable elements, transcription factors, centromeres

**Introduction**

Plants in natural environments are exposed to multiple stimuli, including numerous biotic and abiotic stresses that make it necessary for plants to develop strategies to rapidly adapt. According to the Global Climate Report 2020, the past 10 years were the warmest recorded around the globe in our era. The greater temperature variability has resulted in both droughts and extreme precipitations, affecting not only natural plant populations but also crop production ¹. In order to face these challenges, we need to better understand the mechanisms which allow plants to rapidly adapt and evolve to better cope with increasing climate change-related stresses. Recent advances in genome sequencing have revealed how dynamic plant genomes can be under stressful scenarios ²,³. This dynamism can be attributed to both genetic and epigenetic mechanisms which can contribute to specific traits ⁴,⁵. However, how epigenetic information is influenced by stresses ⁵,⁶ and whether these can contribute to adaptation requires a better understanding. One key epigenetic mark is DNA methylation (5-methylcytosine) which exists in three sequence contexts in plants: CG, CHG, and CHH (H = A, C, or T). Each of them is regulated by distinct, but also interconnected silencing mechanisms ⁷. Symmetric methylation in the CG sequence context (mCG) has been found to be enriched in gene bodies but the biological function of gene body methylation (gbM) remains unclear ⁸. mCG is highly heritable and able to persist over many generations. Conversely, DNA methylation in the CHG (mCHG) and CHH (mCHH) sequence contexts show a lower stability ⁹. DNA methylation is a hallmark epigenetic modification, contributing to the regulation of many biological processes such as genome stability, definition of euchromatin and heterochromatin, control of gene expression, and, most importantly, silencing of transposable elements (TEs) ³,⁷. Studies of DNA methylation variability in natural *Arabidopsis* accessions have shown a clear correlation between epigenomic changes in coding and non-coding genomic regions and environmental stimuli, suggesting a role for DNA methylation in adaptation ¹⁰. More generally, it has been found that DNA methylation changes may be implicated in morphological changes in response to different climates in plants ¹¹. It has been suggested that DNA methylation might play a role in stress tolerance acquisition in plants. For instance, chilling stress in tea plant (*Camellia sinensis*) was found to cause
hypomethylation in cold-responsive genes greatly enhancing their transcription thereby contributing to low temperature tolerance¹². Similarly, low DNA methylation over gene bodies in seagrasses (*Posidonia oceanica* and *Cymodocea nodosa*) when sea-surface temperature increased showed transcriptional flexibility as an adaptive mechanism under changing environments¹³. Extreme water deficiency can result in DNA hypomethylation and upregulation of genes in the resurrection plant *Boea hygrometric* suggesting a rapid regulation of the stress-response machinery to contribute to dehydration tolerance¹⁴. In rice, it was show that a high affinity potassium transporter gene is regulated by DNA methylation in its promoter in order to maintain potassium and sodium homeostasis during salinity stress¹⁵. A notable example concerning environmental impacts on DNA methylation is the documented loss of methylation at TEs at high altitudes and how it may contribute to local adaptation in natural populations of wild strawberry¹⁶,¹⁷. However, whether stress-induced methylome alterations at individual loci or across the entire genome contribute to heritable changes in DNA methylation patterns and adaptation remains uncertain.

In the case of (octoploid) strawberry, it has previously been well described how abiotic stresses influence the normal development of the plant. To illustrate, in strawberry cultivars, drought stress was shown to affect chlorophyll, proline, and soluble carbohydrates levels. After rewatering, only proline and soluble carbohydrates went back to initial levels ¹⁸. Low and high temperatures are known to disturb strawberry yield and fruit quality. For example, at high temperatures, strawberries initiate the production of chemical compounds which stunt plant growth¹⁹. Similarly, cold negatively affected anthocyanin content and the accumulation of soluble sugars in strawberry fruits²⁰. Furthermore, hormonal levels (indole-acetic acid and abscisic acid) in seeds were also affected by low temperatures²⁰. Higher tolerance of salt stress in strawberry cultivars have been associated with low stomatal density. However, salt stress caused accumulation of toxic ions in the leaves of non-salt tolerant strawberry plants²¹,²². Plants in low light had low photosynthetic rate with an increase of photosynthetic pigments, reduced stomatal conductance and altered sugar accumulation²³. Conversely, at high light conditions, leaf height and photosynthetic rates augmented and as consequence, affected plant growth²⁴. Finally, it was previously shown that SA increased vegetative growth and early flowering²²,²⁵. However, little was so far known about how these stresses can influence the epigenome of strawberry.
Here, we wanted to assess how a broad spectrum of well-established, ecologically relevant and climate change-related stresses can influence the epigenome of *F. vesca* (diploid strawberry). From an epigenetic point of view, *F. vesca* is an interesting model plant as it can multiply clonally (via stolons) or sexually (via seeds)²⁶. Furthermore, in this species DNA methylation regulates key developmental traits such as seed dormancy²⁷ and fruit ripening²⁶. Overall, we identified genomic regions which may act as stress-responsive epigenetic rheostats as a strategy to maintain homeostasis during unfavorable growth conditions in *F. vesca*. Our findings can be used to target epigenetic regulated genes for improvement of adaptation or stress tolerance in crops through epigenome editing.

## Materials and Methods

### Plant growth and material

All strawberry plants used in this study were a homozygous cultivated near-isogenic line (NIL), Fb2:39-47, *F. vesca* cv. Reine des Vallées (RV), possessing the “r” locus on chromosome 2 which causes this accession to propagate vegetatively through stolon development²⁸. Seeds from a single founder plant were germinated in water over Whatman filter paper for two weeks and transferred to 50% Murashige & Skoog (MS) medium (Duschefa cat# M0222), 30% sucrose, and 2% phytagel (Sigma-Aldrich cat# P8169) and grown for 4 weeks prior to stress.

### Stress assays

One-month-old seedlings on agarose plates were exposed to different stresses under long-day conditions (16 h light 24°C/8 h dark 21°C) in plant growth chambers (Panasonic, phcbi: MLR-352/MLR-352H). The seedling age was optimized to assure stress tolerance. For salt and drought stress, one-month-old plantlets were transferred to MS media supplemented with 100 mM sodium chloride (NaCl) (Sigma-Aldrich, cat# S9888) and 5% polyethylene glycol (PEG-6000) (-0.05 MPa) (Sigma-Aldrich, cat# P7181), respectively. For cold and heat stresses, plants were initially grown as described above. The plates were then transferred to either 6°C or 37°C chambers. High light was induced by 20,000 lx of illuminance (460 µmol s⁻¹ m⁻²) and low light with 80% sunblock black net leading to 4,000 lx of illuminance (92 µmol s⁻¹ m⁻²). To simulate a hormone stress, MS medium was supplemented with 0.5 mM salicylic acid (SA) (Sigma-Aldrich, cat# 247588). All stress assays were carried out for 2 weeks with 2 recovery days after
one week. We sampled the following arial parts of plants: meristem, 2 cotyledon leaves, the first true leaf and the first trifoliate leaf. For molecular analysis, we chose 5 plants to reduce the variability among individual plants and three biological replicates for statistical analysis. The requirements for methylome sequencing has been well defined and standardized in different previous studies<sup>29,30</sup>. The experiments were designed to meet the requirements to run the necessary bioinformatic pipelines<sup>31,32</sup>. The samples were harvested in 1.5 mL tubes between 9:00 a.m. and 11:00 a.m. and immediately flash-frozen in liquid nitrogen and stored at -80°C.

## Genome sequencing and assembly NIL Fb2

Genomic DNA from strawberry plants was extracted by a Hexadecyltrimethylammonium bromide (Cetrimonium bromide, CTAB) modified protocol<sup>28</sup> and purified with Agencourt AMPure XP beads (cat# A63880). Long-read sequencing was performed for genome assembly; Genomic DNA by Ligation (Oxford Nanopore, cat# SQK-LSK109) library was prepared as described by the manufacturer and sequenced on a MinION for 72 h (Oxford Nanopore).

## Reference genome polishing

Reads obtained from nanopore were filtered with Filtlong v0.2.1 (https://github.com/rrwick/Filtlong) using --min_mean_q 80 and --min_length 200. Cleaned reads were then aligned to the most recent version of the *F. vesca* genome v4.0.a1<sup>33</sup>, with the annotation of *F. vesca* genome v4.0.a2 downloaded from the Genome Database for Rosaceae (GDR) (https://www.rosaceae.org/species/fragaria_vesca/genome_v4.0.a2), using minimap2 v2.21<sup>34</sup> with parameters -aLx map-ont --MD -Y. The generated BAM file was then sorted and indexed with samtools v1.11. We used mosdepth v0.3.1 to verify that coverage on chromosomal scaffolds was over 50 X. Sniffles v1.0.12a with parameters: -s 10 -r 1000 -q 20 --genotype -l 30 -d 1000 was used to detect structural variations larger than 30 bp. We observed that larger structural variants (SV) were most likely falsely identified due to misalignments in regions with gaps or Ns, therefore the VCF files obtained from Sniffles were sorted and filtered with BCFtools v1.14 to keep only SV with less than 200 kb, supported by 10 or more reads and with allelic frequencies above 0.8 to isolate homozygous changes. The complete filtering command used was "bcftools view -q 0.8 -Oz -i '(SVTYPE="DUP" || SVTYPE="INS" || SVTYPE="DEL" || SVTYPE="TRA" || SVTYPE="INV" || SVTYPE="INVDUP") && %FILTER="PASS" && FMT/DV>9 && SVLEN>29 && SVLEN<200000"". From the VCF listing all the structural variants that we detected in our *F. vesca* accession, we generated a substituted genome version
based on the reference *F. vesca* genome v.4.0.a2. The reference genome was first indexed with samtools faidx v1.11 and a sequence dictionary were generated with Picard CreateSequenceDictionary v2.25.6 (https://broadinstitute.github.io/picard). The VCF containing the SV produced from our Nanopore sequencing was also indexed with IndexFeatureFile v4.2.0.0 (https://gatk.broadinstitute.org/hc/en-us/articles/360037262651-IndexFeatureFile). FastaAlternateReferenceMaker v4.2.0.0 (https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker) was then run with the reference genome and the VCF file to generate a substituted genome representative of our *Fragaria* accession (Fb2). As substituting our genome with the detected structural variants changes genomic coordinates, we also corrected the public GFF genome annotation of *F. vesca* ³⁵ using liftoff v1.6.1. Liftoff also detects and annotates duplications within the substituted genome. Transposable elements annotation was carried out using the EDTA transposable element annotation pipeline v. 1.9.6 ³⁶ on the substituted genome using default parameters.

## Whole-genome bisulfite sequencing (WGBS)

A modified CTAB DNA extraction protocol was performed using frozen above-ground tissues ²⁸. DNA libraries were generated using the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs, cat# E7103S) according to the manufacturer's instructions with the following modification for bisulfite treatment. DNA was sheared to 350 bp using a Covaris S2 instrument. The bisulfite treatment step using the EZ DNA Methylation-Gold kit (Zymo Research, cat# D5007) was inserted after the adaptor ligation. After clean-up of the bisulfite conversion reaction, library enrichment was done using Kapa Hifi Uracil+ DNA polymerase (Kapa Biosystems, cat# KK1512) for 12 PCR cycles, using the 96 single-index NEBNext Multiplex Oligos for Illumina (New England Biolabs, cat# E7335S). Paired-end reads were obtained on an Illumina (150 bps) NovaSeq6000 instrument at Novogene (Hongkong, China).

## Processing and alignment of bisulfite-converted reads

Sequencing data was analyzed by data collection software from read alignment to DNA methylation analysis: Epidiverse/wgbs pipeline³¹. The pipeline included quality control using FastQC v.0.11.9 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Cutadapt v.3.5 (https://github.com/marcelm/cutadapt/). Genome mapping was performed using erne-bs5 v.2.1.1 with default parameters to generate the BAM files. Methylation calling and methylation bias correction was performed with Methyldackel v.0.6.1
(https://github.com/dpryan79/MethylDackel) with only uniquely-mapping reads. The pipeline used Nextflow v20.07.1 to run multitask in parallel. Because plant chloroplast DNA are not methylated³⁷, reads originating from those sequenced were used to evaluate the bisulfite conversion rate. The pipeline is available at https://github.com/EpiDiverse/wgbs. An average of 82,771,701 reads (~50X coverage) were produced per sample, of which 81% mapped properly to the *F. vesca* genome. The average non-bisulfite conversion rate among the samples was 0.10 (See Table S1 for more details). To calculate global methylation ratios, output files from wgbs pipeline were pre-filtered for a minimum coverage of 5 reads using awk command and only the common cytosine positions were kept among all the samples using bedtools 2.28.0. The data were tested for statistical significance with an unpaired Student's t-test, *p* < 0.05 was selected as the point of minimal statistical significance in all the analyses. R-packages ggplot2 v.3.3.5 and gplots v.3.1.1 were used for the visualization of the results. Global DNA methylation levels were computed by combing all bedGraph files of the individual samples into a unionbedg file and filtered for cytosine positions without sequencing coverage and average DNA methylation levels calculated.

## Identification of differentially methylated regions (DMRs)

First, bedGraph files from wgbs pipeline were pre-filtered for a minimum coverage of 5 reads using awk command. These output files were then used as input for the EpiDiverse/dmr bioinformatics analysis pipeline for non-model plant species to define DMRs ³² with default parameters (minimum coverage threshold 5; maximum q-value 0.05; minimum differential methylation level 10%; 10 as minimum number of Cs; Minimum distance (bp) between Cs that are not to be considered as part of the same DMR is 146 bp). The pipeline uses metilene v.0.2.6.1 (https://www.bioinf.uni-leipzig.de/Software/metilene/) for pairwise comparison between groups and R-packages ggplot2 v.3.3.5 and gplots v.3.1.1, for visualization results (Fig. S1). Based on our *F. vesca* genome transcript annotation and methylation data (overlapped regions with DNA methylation cytosines and DMRs), we detected the methylated genes, promoters, 3' UTRs, 5' UTR and transposable elements in strawberry. Global DNA methylation and DMR plots were performed with R-package ggplot2. Gene analyses by methylation patterns and analysis of per-family TE DNA methylation profiles were performed with deepTools v.3.5.0. DMRs comparison between treatments were done by the Venn diagram v.1.7.0 R-package. We
produced several genome browsers tracks with DMRs that we integrated in our local instance of
JBrowse available at the following url:

<https://jbrowse.agroscope.info/jbrowse/?data=fragaria_sub>

Heat stress and salt stress RNA-seq analysis

One-month-old seedlings on agarose plates were exposed to heat stress (37°C). For salt stress, plants were transferred to MS media supplemented with 100 mM sodium chloride (NaCl). The stresses were carried out for 2 weeks with 2 recovery days after one week. Following the recovery days, 5 plants were pooled to reduce the variability among individuals (three biological replicates per condition) and flash frozen in liquid nitrogen for transcriptome sequencing. RNA extractions were performed using NucleoSpin RNA Plus, Mini kit for RNA purification with DNA removal column (cat# 740984.50). Samples were sent in heat, salt, and control condition (n = 9) for Illumina paired-end reads sequencing (150 bps) to Novogene (Hongkong, China).

Differential gene expression analysis

RNA-seq analyses were performed as previously described in ³. Briefly FastQC v.0.11.9 and
trimmomatic v.0.39 packages were used for quality control and trimming. Salmon v1.4.0 was
used to map sequence reads to the *Fragaria vesca* reference transcript dataset. DESeq2 package
was use for the quantitative differential gene transcription analysis at the European Galaxy
platform with default parameters (Fig. S2).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of
differentially methylated and expressed genes

All methylated and differentially expressed genes were annotated based on GO and KEGG annotation downloaded from the Genome Database for Rosaceae (GDR) (https://www.rosaceae.org/species/fragaria_vesca/genome_v4.0.a2).

To better understand the potential function of the differentially methylated and expressed genes, GO functional classification of these genes was performed by AgriGOv1.2 ³⁸. Genes and promoters were classified by genes contained hypo and hypermethylated DMRs. Similarly, all DEGs were extracted with p-adj< 0.05. The GO slim library was used as reference GO reference type. Fisher's exact test p-values were calculated for over-representation of the differential methylated genes in all GO categories and Hochberg (FDR) as multi-test adjust method. GO terms with p < 0.05 were considered as significantly enriched. KEGG analysis was performed
with the R-package: `clusterProfiler`³⁹; `pvalueCutoff` = 0.05; `pAdjustMethod` = "fdr". Plots were
performed with R-package `ggplot2`.

## Results

### Stress-induced DNA methylation dynamics in *F. vesca*

To explore the ecological relevance and the dynamism of DNA methylation following stress exposure, we focused on multiple stimuli which plants may encounter in natural environments (Fig. 1a, see Materials and Methods for details). Here, *F. vesca* plants were submitted to two successive stress applications *in vitro*. To assess potentially stable DNA methylation changes we sampled the plants two days after the second stress. Compared to controls, stressed plants showed: reduce size after cold- and salt stress; tissue necrosis following heat stress; petiole elongation under low light; and reddish tissue coloration following SA exposure suggesting the accumulation of anthocyanins (Fig. 1a, Fig. S3). To assess genome-wide DNA methylation levels, DNA was extracted from these plants and submitted to whole genome bisulfite sequencing (WGBS, 20x genome coverage) (Table S1). We carried out a global quantification of DNA methylation in the three sequence contexts (CG, CHG and CHH). Overall, the global DNA methylation levels in all stress conditions were similar (Fig. 1b). *F. vesca* seedlings in control conditions had 40.11% mCG, 14.93% mCHG and 2.38% mCHH. We found significant decreases of 0.5%, 3.2% and 1.1% for global mCG after cold stress, heat stress and salt stress, respectively (Fig. 1b). In addition, we observed a significant global mCHG decrease of 1.8% in salt and an increase of 3.1% in the presence of SA. We did not detect significant changes in global mCHH level. To assess methylation variation in genic and non-genic sequences, we screened the methylome data in all three sequence contexts separately in three regions: 2 kb upstream of genes, along the gene body, and 2 kb downstream of genes using 100-bp sliding windows (Fig. 1c). In contrast to the global analysis, here we observed a higher local DNA methylation variability in the CHH context at the transcription start and end sites (TSS and TES, correspondingly) compared to the other two contexts. Notably for the CHH context, cold stress, heat stress and salt stress resulted in hypomethylation at the TSS, TES and over the gene body (Fig. 1c). We extracted genes with body methylation (gbM) similarly to the parameters defined previously ⁸ filtering for at least 20 CGs and a methylation level above the median value. Genes containing gbM showed low variability in all the conditions (Fig. S4).
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>(b)</figcaption></figure>

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

<figure><img src="image_6.png" /><figcaption>**Fig. 1 Effect of abiotic and hormone stresses on genome-wide DNA methylation levels in F. vesca.** (a) Top: Scheme of stress treatment time course. Bottom: Photographs of plates with one-month-old plants grown under different stress conditions (cold, drought, heat, high and low intensity of light, salt, and salicylic acid). (b) Average DNA methylation levels for each cytosine context (CG, CHG, CHH) between normal and stress conditions (only common cytosines positions among all samples were considered that had a minimum coverage of 5 reads). Asterisks indicate levels of significance between treated and control plants: *, p-value < 0.05; **, p-value < 0.01; ***, p-value < 0.001 (unpaired two-tailed Student's t-test). (c) Heat maps showing distribution of DNA methylation (top: CG, middle: CHG, and bottom: CHH context) around genes with and without stress (Control). Mean of the average methylation percentage (within a sliding 100-bp window) was plotted 2 kb upstream of TSS, over the gene body and 2 kb downstream of TES. Black arrows highlight the samples in which a reduction of DNA methylation can be observed in the vicinity of TSS and TES sites.</figcaption></figure>
Stress particularly affects DNA methylation variability in the non-CG contexts in F. vesca

To explore the dynamics of DNA methylation at specific loci in detail, we assessed differentially methylated regions (DMRs) for each sequence context. DMRs were defined using metilene v.0.2.6.1 which uses an algorithm to identify a base-pair window through sequence segmentation with significant methylation differences ³² (see Materials and Methods for parameters). We compared the methylomes of plants submitted to each stress condition with the methylomes from control plants. The majority of DMRs identified were in the CHH context for all conditions, with a maximum of 12,414 DMRs under heat stress. Fewer DMRs were detected in both CG sequence context, ranging from 104 (cold stress) to 3,016 (heat stress), and in the CHG context, ranging from 15 (cold stress) to 236 (heat stress) (Fig. 2a, Table S2). In line with our global DNA methylation analysis, most of the heat stress, cold stress, and salt stress DMRs were hypomethylated (hypoDMRs) relative to the control condition in all sequence contexts. For drought stress, high light stress, low light stress, and SA stress, most of the DMRs showed hypermethylation (hyperDMRs) in the CHH context (Fig. 2a). Next, to test whether genic and non-genic regions were rich in DMRs, we qualified DMRs based on their intersection with promoters (empirically defined as 2 kb to 50 bp upstream TSS), gene bodies, and intergenic regions. The minimum overlap required was 1bp. Many of the CG and CHG DMRs (30-44%, 45-60% respectively) were in genes, while most of the CHH DMRs (between 32-44%) were in promoters and intergenic regions (Fig. 2b). In summary, abiotic and hormone stresses led to DNA methylation changes primarily in the CHH context within promoters and intergenic regions. Overall, heat stress caused the most numerous DNA methylation changes.
<figure><img src="image_7.png" /><figcaption>Fig. 2 Differentially methylated regions in *F. vesca* seedlings grown at control and stress conditions.</figcaption></figure>

(a) Number of stress-induced hyper- (hyperDMRs) and hypomethylated DMRs (hypoDMRs) separated by sequence context. (b) Distribution of DMRs in genomic features: promoter (2 kb to 50 bp upstream TSS), gene body, and intergenic regions. Minimum overlap required: 1bp. (c) Venn diagrams of common promoter and genic regions containing hypo- and hyperDMRs per context among all the stress conditions. (d) Percentages of genic loci with hypo- and hyperDMRs which are shared among the stress conditions. Reading from left to the right (arrow). e. g. (Left block) 8.3% genic locations with cold-stressed hypoDMRs overlap with 38.7% genic locations with drought-stress induced hypoDMRs. (Right block) 24% genic loci with hyperDMRs from salt-stressed seedlings overlap with 11% genic loci with hyperDMRs
resulting from SA treatment. Ascending intensity of block colors correlate with overlap percentages.

<table><caption>**Table 1. Gene ontology (GO) enrichment analysis of genic regions (gene and promoter) with DMRs.** GO enrichment of genes with hypo- and hyperDMRs caused by thermal stress. The genes are arranged according to their DMRs patterns.</caption><thead><tr><th>Stress</th><th>Methylation</th><th>Description</th><th>FDR</th><th>Number in input list</th><th>Number in BG/Ref</th><th>Ontology</th></tr></thead><tbody><tr><td rowspan="2">Cold</td><td>Hypermethylation</td><td>transporter activity</td><td>0.024</td><td>23</td><td>838</td><td>F</td></tr><tr><td>Heat</td><td>Hypermethylation</td><td>generation of precursor metabolites and energy</td><td>0.021</td><td>8</td><td>153</td><td>P</td></tr><tr><td rowspan="6"></td><td rowspan="6"></td><td>Hypomethylation</td><td>transcription regulator activity</td><td>0.024</td><td>108</td><td>420</td><td>F</td></tr><tr><td></td><td>transcription factor activity</td><td>0.024</td><td>99</td><td>380</td><td>F</td></tr><tr><td></td><td>intracellular cell part</td><td>0.033</td><td>504</td><td>2276</td><td>C</td></tr><tr><td></td><td>cell</td><td>0.033</td><td>849</td><td>3930</td><td>C</td></tr><tr><td></td><td>cytoplasm</td><td>0.033</td><td>849</td><td>3930</td><td>C</td></tr><tr><td></td><td></td><td></td><td>0.033</td><td>239</td><td>1029</td><td>C</td></tr></tbody></table>

GO categories are listed with the false discovery rate-adjusted P-value <0.05. BG: background, Ref: reference, C: Cellular Component; F: Molecular Function; P: Biological process. Only well identified and annotated genes were included in the analysis.

## Identification of stress-induced DNA methylation changes in genic regions

To better understand the potential functional roles of the DMRs and the commonalities between the different stresses, we focused our analysis on DMRs located within promoters and gene bodies (Fig. 2c). We only identified one locus with a CG DMR (within *FvH4_6g40845* an unknown protein) and 44 loci with CHH DMRs that were in common to all stress conditions (Fig. 2c, Fig. S5). Comparing each stress data set of promoter and gene locations with hypo- and hyperDMRs we noted that heat stress shared more loci with hypoDMRs with the other conditions than all other comparisons (Fig. 2d, left). For example, 37.3% of loci with heat-stress hypoDMRs overlapped with 54.6% of all hypoDMRs found under cold-stress. Conversely, drought stress resulted in the highest number of genic loci with hyperDMRs shared with the other conditions. To illustrate, 33.4% of salt stress hyperDMRs overlapped with 14.6% drought stress hyperDMRs (Fig. 2d, right). In order to identify potential functional roles for these DMRs, we performed a singular enrichment analysis (SEA) using the AgriGO tool 38 (see Materials and Methods). *F. vesca* has around 34,000 genes but only 54% of which have been assigned a GO
number ³⁵. For this reason, “unknown” annotated genes were omitted for this analysis (Table S3). The analysis was based on the identified genic regions (gene and promoter) with DMRs for each stress condition according to their methylation change (hypo- or hypermethylated). Plants submitted to heat stress showed the largest variation in DNA methylation over genic regions. These were enriched with hyperDMR-associated genes involved in the generation of precursors of metabolites and energy. Genes with heat-stress induced hypoDMRs were enriched for transcription factor (TF) activity, transcriptional regulators, and genes related to cellular components (Table 1). We also found that cold stressed plants had hyperDMR-associated genes enriched for transporter activity (Table 1). Next, we wanted to correlate DNA methylation changes with transcriptional changes. For this analysis, we focused on heat stress and salt stress because of the significant loss of genome-wide DNA methylation that we observed for these conditions (Fig. 1b). To evaluate transcriptional changes, we performed RNA-sequencing on polyadenylated RNA extracted from the plants growing in control, heat stress, and salt stress conditions. We found 3966 and 5033 significantly differentially expressed genes (DEGs) in heat stress and salt stress, respectively (Table S4, Table S5). To identify the functions of the DEGs we performed an enrichment analysis using the aforementioned parameters in the AgriGO tool (Fig. S6a). Interestingly, for heat stress we found genes enriched in reproductive process, pollination, and protein modification process. On the other hand, salt stress resulted in a specific response in pyrophosphatase activity, transport, motor activity, homeostatic and hydrolase activity. Even though distinct stress-response genes were responding to the two stresses, heat stress and salt stress shared a substantial number of DEGs with similar expression patterns (Fig. 3a). To illustrate, 1860 down-regulated DEGs and 1049 up-regulated DEGs were in common for both conditions (Fig. 3b). This group of genes in heat stress and salt stress were related with molecular functions such as catalytic activity, transferase activity, nucleotide binding and transcription factors. Biological process genes were related with metabolic processes, and response to stress. Finally, cellular components included cell wall and membrane (Fig. 3c). In order to identify possible pathways playing a role in heat stress and salt stress response we performed a KEGG enrichment analysis³⁹. Genes related with enzymatic and transcription factor activity were enriched for both stresses. For example, peroxidases, reductases and MYB TFs following heat stress. Similarly, under salt stress, MYB TFs were the most enriched group of genes (Fig. S6b). Next, we wanted to identify if there was a relationship between the presence of
DMRs in genic regions (gene and promoter) and transcriptional changes (Fig. 3d). We identified 172 heat stress DEGs and 124 salt stress DEGs with DMRs located in the gene body; 637 and 366 DEGs with promoter DMRs; 58 and 14 DEGs with DMRs in both gene body and promoter in each stress respectively. However, the presence of hypo- or hyperDMRs did not correlate with the transcription patterns (up- or down-regulation) of these genes (Fig. 3d). For instance, hypoDMRs were equally distributed between up- and down-regulated genes in both conditions. Additionally, we identified DEG with heat-stress DMRs which were enriched in transcription factor activity, reproductive process, membrane, and cell wall (Fig. S6c). In the case of salt stress, the enrichment representation of DEG with salt-stress DMRs showed genes involved in catalytic activity (Fig. S6c). Taken together, heat stress caused more DNA methylation changes while salt stress caused more transcriptional changes.

<figure><img src="image_8.png" /><figcaption>**Fig. 3** Functional analysis of differentially expressed genes (DEGs) after heat stress and salt stress. (a) Heatmap of DEGs of plants grown under heat stress and salt stress. Heatmaps were generated with log₂ expression of counts normalized to transcript size and million mapped reads (DESeq2). DEGs are listed along the Y-axis in the order they clustered in as indicated by the colored line along the Y-axis. Each column contains expression values for individual genes, with groups indicated along the X-axis. Deeper red colors indicate higher expression while deeper blue indicates lower expression compared to the control. (b) Bar plots indicating the total number</figcaption></figure>
of DEGs in heat stress and salt stress. Intersected bars show the number of shared DEGs in heat stress and salt stress. (c) Singular enrichment analysis (SEA) of the total number of DEGs shared between heat stress and salt stress (results from AgriGOv2, p-adj:<0.05). The x-axis indicates the number of genes in a category. The y-axis indicates the most enriched GO terms in three categories: biological processes (BP, grey), cellular component (CC, orange), and molecular function (MF, dark blue). The x-axis indicates the number of genes in a category. (d) Scatter plot of DEGs related with DMRs (promoter and gene body) showing the relationship between transcript levels (fold change: log₂) and DNA methylation (meth ratio). Top: heat stress and bottom: salt stress. Blue: hypoDMRs; Red: hyperDMRs.

## Stress induces hypomethylation and its effect on the expression of transcription factor coding genes

Transcription factors (TFs) play key roles in plant growth, development, and stress responses. Remarkably, we observed that heat stress resulted in an enrichment of hypoDMRs at promoters of genes related to TFs. Furthermore, heat-stress DEGs associated with DMRs also showed an enrichment in TFs (Fig. S6c). We found that under heat stress, 99 genes from the GO category transcription factor activity (Table1) were associated with DMRs. Of these, 31% were *AP2/EREBP*, 19% *WRKY* genes and 4% heat shock transcription factors (*HSFs*). *AP2/EREBP* members (119 genes in total) have been characterized in the latest version of the *F. vesca* genome and are known to be involved in stress tolerance<sup>40,41</sup>. To evaluate the detailed DNA methylation changes at *AP2/EREBP* genes under different stress conditions, we looked at the distribution of DNA methylation at those loci (Fig. S7a). We observed a noted reduction in DNA methylation at the TSS in the CHH context for heat-, cold- and salt-stressed plants. Combining all stress conditions, a total of 74 DMRs were identified within the promoter regions of 44 *AP2/EREBP* genes (Table S6). In addition, 14 *HSFs* have been identified in *F. vesca* genome<sup>42</sup>; however, in the last genome annotation version<sup>35</sup>, we were able to identify 19 *HSFs* (Table S7). The distribution of methylation over *HSFs* showed high variability in TSS and TES in CHH context compared gene bodies (Fig. S7b). We identified 13 *HSFs* with DMRs within promoter regions mostly attributed to heat stress (Table S7). The marked presence of DMRs within TF promoters suggested a relationship between methylation and transcription. For this reason, we selected different members of TF families to evaluate changes in methylation and transcription (Fig.4). Four examples with clear changes in DNA methylation in the promoters of the *AP2/EREBP*, *HSF*, *WRKY*, and *MYB* TF gene families are shown in Fig. 4a. We therefore extracted the list of members of the TF superfamilies to verify their transcript levels in our RNA-seq data for heat stress and salt stress (Fig. 4b). Overall, we identified 13 *AP2/EREBP* genes that
were significantly up regulated in heat stress and 17 *AP2/EREBP* genes in salt stress (Fig. 4b).

We also obtained similar results by analyzing the relationship between methylation and expression of heat shock factors (*HSFs*). Only *FvHSFB1a* and *FvHSFB3a* showed significant changes in transcription in heat stress (Fig. 4b). Hypomethylation in promoter regions of *FvHSFB1a* and *FvHSB3a* showed a clear relationship with significantly increased transcript levels after heat stress. However, in salt stress *FvHSFB1a* and *FvHSB3* were up regulated without being associated with a DMR. Similarly, in *F. vesca*, 61 *WRKY* and 217 *MYB* genes have been identified<sup>35,43,44</sup> and described to be differentially expressed under abiotic stress conditions.

Here, as part of the DEGs analyzed, we found 16 *WRKY* genes to be up regulated by heat stress and 9 of which contained DMRs (Fig. 4b). On the contrary, 14 *WRKY* genes were differentially expressed after salt stress with almost no associated DMR. For the MYB superfamily, 48 genes were differentially expressed in heat stress and 64 in salt stress (TableS8). However, few of those were related with a DMR: 19 and 9 MYB TFs were related with a DMR in heat stress and salt stress respectively (Fig. 4c, TableS8). Collectively, our data show that heat stress and salt stress trigger the up- and down- regulation of TFs. Nonetheless, only heat-stress induces loss of DNA methylation mostly at promoter regions of genes related to specific TF families possibly influencing their expression. Overall, transcriptional changes seemed not to correlate with changes in DNA methylation.
<figure><img src="image_9.png" /><figcaption>**Fig. 4** Stress induced hypoDMRs at transcription factor coding genes. Genome browser views of DMRs located in promoter regions of (a) APETALA2/ethylene-responsive element binding protein (AP2/EREBP) superfamily genes: *ERF36* (ethylene-responsive binding factor 36); heat shock transcription factors (*HSF*) genes: *FvHSFB3a* (heat shock transcription factor B); *WRKY* superfamily: *WRKY 75.2*; *MYB* superfamily: *FvMYB85*. Depicted are genes structures (top panels, UTRs in light blue, exons in blue), TEs (red and dark blue) and DNA methylation levels (histograms). Boxes above the histograms indicate identified DMRs with methylation difference ratios (color codes for DNA methylation: red for CG, blue for CHG and yellow for CHH contexts). (b) Differentially expressed transcription factor (TF) families after heat stress. Transcription profiles of *AP2/EREBP*, *HSF*, *WRKY*, *MYB* genes. (c) selected differentially transcribed *MYB* genes related with a DMR (dashed box). Note that for space reasons all differentially expressed *MYBs* without DMRs have been omitted here. DMR location is indicated</figcaption></figure>
in the upper colored boxes: TES -2kb (grey), promoter (blue), none (green). The x-axis indicates the TF member per family. The y-axis indicates fold change (log₂) with p-adj. < 0.05.

## Stress leads to distinct methylation changes at transposable elements (TEs)

Since one of the most important functions of DNA methylation is to repress TE transcription and mobility⁷,⁴⁵, we investigated the effect of stress on DNA methylation at TEs. To describe the variation in DNA methylation in TE bodies and their flanking regions we plotted DNA methylation profiles in all three contexts. We used 50-bp sliding window: 2kb upstream, over the body and 2 kb downstream (Fig. 5a). In general, all TE families showed high methylation levels in CG context; however, DNA transposons such as the *Mariner* (DTT) and *Helitrons* superfamilies were characterized by low DNA methylation levels in the CHG context. Notably, *Miniature Inverted-Repeat Transposons* (MITEs) showed the highest DNA methylation levels in CHH context. For MITEs, mCHH was distinctly reduced under cold, heat and salt stress (Fig. 5a). To assess significant changes in DNA methylation at specific loci over TEs, we counted the number of DMRs within TE annotations for each stress condition (Fig. 5b). The minimum overlap DMRs/TEs was 1bp. As for genes, we identified the greatest number of DMRs within TEs in heat-stressed plants in all cytosine contexts (Fig. 5b). Cold stress, heat stress, and salt stress displayed more hypoDMRs in TEs compared to drought stress, low light stress, high light stress, and SA stress, which induced more hyperDMRs in TEs (Fig. 5c). Although there was no specific TE family significantly enriched in DMRs, heat stress resulted in at least 12% of all MITES acquiring hypoDMRs (Fig. 5c, Table S9) and most of them were close to genic regions (< 2kb upstream genes) (see Fig. S8 for examples). Overall, these results suggest that DNA methylation dynamics among all TE members changed proportionally in all superfamilies; nonetheless, the change regarding gain or loss of DNA methylation was clearly defined by the applied stresses.

## DMRs are enriched in distinct regions of the *F. vesca* in the genome

Even though genome-wide DNA methylation variation levels were low (Fig. 1b), we observed regions in the genome that were enriched in DMRs in all stress conditions (Fig. 6a). Indeed, we observed DMR hotspots in the *F. vesca* genome. Notably, the distribution showed a similar pattern on all chromosomes independently of the stress conditions. The DMR density was not clearly related to gene and TE density (Fig. 6b); nevertheless, when we analyzed individual TE families, *Helitron* density hotspots showed a clear correlation with DMR density (Fig. 6). The
*Helitron-rich regions may highlight centromeric regions, however little is currently known about the exact localization of centromeric and pericentromeric regions in strawberry* ⁴⁶. However, it has been reported previously that centromeric regions can be enriched in *Helitrons* in certain plant genomes*⁴⁷*. It is important to note here, that the *Helitrons* themselves were not enriched in DMRs.

<figure><img src="image_10.png" /><figcaption>**Fig. 5 Association of stress-induced differentially methylated regions with transposable elements in *F. vesca*.** (a) Heatmaps showing DNA methylation profiles for all the TE families separated by sequence context mCG (left), mCHG (center) and mCHH (right). The mean of the average DNA methylation percentage (within 50 bp sliding windows) was plotted for the TE</figcaption></figure>
bodies and 2 kb around the TSS and TES regions. (b) Number of stress induced DMRs in TEs per sequence context. (c) Number of hypoDMRs (left) and hyperDMRs (right) within different TE families. Class I elements (retrotransposons): *LTR-Copia*, *LTR-Gypsy*. Class II elements (DNA transposons): TIR: *Tcl-Mariner* (DTT), *hAT* (DTA), *Mutator* (DTM), *PIF-Harbinger* (DTH), *CACTA* (DTC); *Helitron*; *Miniature Inverted-Repeat Transposons* (MITEs). Upper char box indicates the colors which represent each treatment.

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

**Fig. 6 Genome-Wide distribution of DMR densities.** (a) Circa plot showing DMR densities on all 7 strawberry chromosomes (gray boxes) for each stress condition (2. Cold, 3. Drought, 4. Heat, 5. High light, 6. Low light, 7. SA, and 8. Salt). *Helitron* density (1) in the genome is depicted on the outer most circle in red with turquoise background. (b) DMR density depicted on chromosome 5 (Fvb5) displaying a common enrichment among all the stress conditions (left of the 20 Mb tick mark). Below, the chromosome is indicated in grey, green shows gene density, yellow the TE density and red Helitrons along the chromosome.

Discussion

Genome-wide DNA methylation patterns are altered under stressful environmental conditions

One genomic response of plants to abiotic and biotic stresses is to change their epigenome ⁵.
DNA methylation variation might trigger modifications in plant development and physiology
contributing to phenotypic variability, thereby presumably contributing to plant acclimation ⁴⁸.
Here, we wanted to investigate how these previously well-described stresses can affect DNA
methylation at the genome-wide level. Looking at global DNA methylation levels in the three
different sequence contexts, we observed only slight variations among the different stress
conditions. However, when performing the DMR analyses we detected extensive differences at
specific loci especially in the CHH context. This is in agreement with previous studies which reported that mCHH was most dynamic in response to different climates⁹. More specifically, in the case of *F. vesca*, altitude variations of natural strawberry populations was found to be correlated with a high variability in DNA methylation in the CHH context¹⁶,⁴⁹. In our study, cold stress, heat stress, and salt stress resulted in substantial local loss of DNA methylation in all sequence contexts, particularly in regions close to TSS and TES sites (Fig. 1c). Such DNA methylation changes have been associated with different degrees of cold tolerance in other plant species. For instance, chickpea presented an increased number of hypomethylated regions near abiotic stress response genes and transcription factors upon cold stress⁵⁰. Similarly, higher temperatures resulted in hypomethylation in rice seeds, soybean roots, and tobacco leaves, affecting plant growth by altering gene expression patterns in genes that control biosynthesis and catabolism of phytohormones⁵¹⁻⁵³. In addition, salt stress-induced epigenetic variation in Arabidopsis has been shown to be partially transmitted to offspring, primarily via the female germ line⁵⁴. On the other hand, we noticed some gain of global DNA methylation caused by different intensities of light or drought stress and SA treatments resulting a considerable number of hyperDMRs in the CHH context. Different results were found in tomato plants exposed to variable intensities of light which resulted in DNA hypomethylation and transcriptional changes causing male-sterility⁵⁵. Furthermore, drought stress induced mCG and mCHG hypermethylation and a slight decrease in mCHH methylation in mulberry⁵⁶. There is also evidence that drought, nitrogen-deficiency and heavy metal stresses can result in heritable changes in DNA methylation levels across generations in rice⁵⁷⁻⁵⁹. Another study showed that DNA methylation changes induced by hormone stresses via jasmonic and salicylic acid can be faithfully transmitted to offspring by asexual reproduction in dandelions⁶⁰. Taken together, a common pattern in the responses of DNA methylation changes between plant species does not seem to exist. Our results raise an important question that remains to be tested: Are stress-induced DNA methylation changes in *F. vesca* maintained over generations through asexual (stolon) and sexual (seeds) reproduction? This is an intriguing aspect we are currently investigating.

**Different TE superfamilies show contrasted responses to stresses**
DNA methylation plays a key role in limiting transcriptional activation and mobilization of TEs in order to ensure genome integrity⁴⁵,⁶¹,⁶². Indeed, TEs can be an important source of genetic and epigenetic variation that can influence stress-responses⁶³. Here, we wanted to better understand how different TE superfamilies respond to the stresses we applied. Using this approach, we observed distinct DNA methylation profiles at TEs depending on their superfamily and the applied stress conditions. Overall, we found that *F. vesca* TEs are highly methylated in both, CG and CHG sequence contexts. Similar results were observed in maize where LTR and TIR elements are highly methylated under normal conditions⁶⁴. On the other hand, we noticed that overall DNA methylation levels in the CHH context were lower and more variable among different TE superfamilies and stress conditions, implying contrasted responses of TEs to stresses and that these TEs are silenced by different transcriptional gene silencing pathways. These DNA methylation changes could have direct physiological impacts as non-CG methylation seems to create a boundary between genes and TEs⁹. For example, it has been found that TEs located close to stress-induced genes in Arabidopsis and rice are silenced by hypermethylation after phosphate starvation in order to prevent collateral activation of TE transcription during stress⁶⁵. In our study we found that heat stress affected mCHH in the TE body of all TE superfamilies. MITES had the highest percentage of members targeted by DMRs (Table S8 and Fig. S8). Previous studies have highlighted the importance of MITES in genome evolution and how MITE insertions in promoter regions can regulate the expression of genes in a wide variety plant species such as mulberry⁶⁶. Taken together, these observations suggest that specific TE superfamily members with dynamic DNA methylation levels may contribute to stress response strategies in plants.

## DMR location preferences for centromeric regions

It has long been established that DNA methylation is enriched in peri-/centromeric regions which follow the distribution of TEs over the chromosomes of genomes of the Brassicaceae family, as recently also confirmed for *Thlaspi arvense*⁶⁷,⁶⁸. In our study, while investigating the global distribution of DMRs over the *F. vesca* chromosomes, we found that regions with high DMR density correlated with regions enriched in *Helitrons*. Interestingly, we made this observation for all stress conditions (Fig. 6). Currently, *F. vesca* centromeres are not well defined; however, a genome-wide scan of the *F. vesca* genome for tandem repeats suggested the presence of
*Helitrons* near centromeres ⁴⁷. This supports our hypothesis that *F. vesca* centromeres or pericentromeric regions respond to stresses with DNA methylation changes. To further confirm the exact localization of *F. vesca* centromeres, immunoprecipitations using a CENH3 antibody followed by sequencing will be required ⁶⁹. Why centromere-associated regions are more prone to DNA methylation changes and the physiological relevance of this observation still needs to be determined.

## Heat stress and salt stress trigger similar transcriptional responses but different magnitudes of DNA methylation changes

In our study, the global analysis of the DMRs resulting from each tested stress condition highlighted numerous shared genomic patterns among the stresses but also interesting stress-specific genomic features. In the case of heat stress, we observed hypoDMRs predominantly in gene, promoter, and TE regions. Extreme temperatures are one of the main stresses affecting plants that particularly alter their development and potentially cause yield loss ⁷⁰. The functional analysis of DEGs after heat stress in *F. vesca* showed an enrichment in genes related to stress response, reproduction process, photosynthesis, and metabolic process (Fig. S6a) which is similar to what has been previously described as an effect of the high temperatures ⁷⁰. Similarly, DEGs resulting from salt stress were related to external stimuli response, transporter activity, transcription factors, catabolic and metabolic processes (Fig. S6a). It has previously been described that heat stress and salt stress trigger the activation of common genes related with osmotic stress, oxidative stress and leaf physiology, and particularly, photosynthetic activity in plants such as olive, rice, and coffee ⁷¹,⁷². In this study we found that DMRs were enriched in genes related to transcription factor regulation and activity as well as generators of metabolites and energy for heat stress. This result is consistent with the idea that transcription factors are required to reprogram stress-related genes ⁷³. There is evidence demonstrating that heat stress responses can also have an effect on epigenetic regulators and small RNAs to rapidly activate genes ⁷⁴. However, little is known about how these responses are directly influenced by changes in DNA methylation or *vice versa*. Transcription factor families such as *AP2/EREBP*, *WRKY*, *MYB*, *HSFs* play important roles in response to abiotic stresses in plants ⁴¹,⁴⁴,⁷⁵. To illustrate, in *Populus thricocarpa*, 1156 TFs (including MYB, AP2, WRKY, NAC, and bHLH) showed loss of methylation and changes in expression after drought stress ⁷⁶. Moreover, some *AP2/EREBP*
genes are known to be highly induced under heat stress conditions by HSFs through an
interconnected stress regulatory network<sup>77</sup>. Here, we provide epigenetic evidence suggesting
that members of the *AP2/EREBP* TF class might be regulated by changes in DNA methylation in
*F. vesca* after heat stress. On the other hand, salt stress resulted in the activation of similar TFs
without accompanying DMRs in the gene body or promoter regions. Among them, the promoter
region of the ethylene response factor (*ERF*) genes showed loss of methylation only after heat
stress but significant up-regulation in heat stress and salt stress. Recent studies have shown that
ERFs enhance basal thermotolerance by regulating heat-responsive genes and interacting with
HSF in Arabidopsis and tomato<sup>78,79</sup>. Similarly, hypomethylation in the TSS of genes involved in
the control of cell growth in tobacco and stress-tolerance genes in maize after heat stress
exposure is consistent with the increase in their transcription levels<sup>51,80</sup>. Correspondingly, we
identified 58% of *HSFs* genes with hypoDMRs in their promoter regions after heat stress. *HSFs*
are crucial for thermotolerance capacity and regulate the expression of several heat-stress
response genes such as heat shock proteins (HSPs)<sup>77</sup>. Here, we showed up-regulation of class B
*HSFs* in *F. vesca* after heat and salt stress. Comparable results were obtained in a transcriptome
analysis of the octoploid strawberry where *HSF* expression was induced by a heat shock
treatment<sup>81</sup>. Notably, we found 26% of all *WRKY* genes to be differentially expressed in heat
stress which is similar what was described by Wei *et al.*, (2016) who found 8 *WRKY* genes to be
up-regulated at 42°C in *F. vesca* by semi quantitative PCR<sup>43</sup>. At the same time, in soybean
*WRKY* genes were highly induced under salt stress conditions<sup>82</sup>. Concerning MYB TFs, it has
previously been reported that they can be regulated by gibberellic acid, abscisic acid, cold, and
heat treatments in *F. vesca*<sup>44</sup>. Markedly, differentially expressed MYB TFs were affected in
another way at the epigenetic level, while after heat stress 40% of those MYB TFs were
associated with a DMR, that was only the case for 14% of them after salt stress (TableS8). Taken
together, these findings provide insights into stress induced DNA methylation as being partially
independent from the transcriptional regulation. It will be of great interest to investigate the
direct role of DNA hypomethylation in promoter regions of genes in regulating or priming
transcription following heat stress and salt stress in *F. vesca*.

Conclusions

In summary, our data revealed how DNA methylation profiles at genes and transposable
elements can vary in response to stresses in wild strawberry. In addition, we observed changes in
DNA methylation and gene expression changes do not necessarily coincide following stress exposure. These results show how versatile and plastic the epigenome can be under unfavorable conditions. Furthermore, we provide insights into how specific chromosomal regions can vary at DNA methylation levels under stress conditions. These observations suggest that the epigenetic flexibility of centromeres may play an important role during plant stress response. These results obtained with *F. vesca* as a model plant will help to better understand the stress response of more complex genomes in the Rosacea family. Overall, this study with high-resolution methylome mapping of the *F. vesca* genome will contribute to a better comprehension of epigenetic responses under variable growth conditions. It remains to be tested if such epigenetic changes can be inherited during sexual or clonal propagation (which is common in *F. vesca*) and if such changes could contribute to adaptation to changing environments.

## DATA AVAILABILITY

The datasets generated and/or analyzed in this study are available in the Zenodo repository (DOI: 10.5281/zenodo.6141713). The sequencing data from this study have been submitted to European Nucleotide Archive (ENA, www.ebi.ac.uk/ena/, accessed on ERP135585) under the project PRJEB50996. The Bisulfite-sequencing raw read fastq accessions under ERR8684931-ERR8684954. The raw reads of the RNA-sequencing under the accessions: ERS11547432-ERS11547437 and ERS12247535-ERS12247537.

## Acknowledgements

The European Training Network “EpiDiverse” received funding from the EU Horizon 2020 program under Marie Skłodowska-Curie grant agreement No 764965; The European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program; No 725701 “BUNGEE” to E.B.. Funding for open access charge: Agroscope institutional funding. This study was supported by INRAE, Angers-Nantes, France and Agroscope, Nyon-Switzerland. We would like to thank all the members of the EpiDiverse consortium (www.epidiverse.eu) and the Crop Genome Dynamics research group for invaluable support, Katharina Jandrasits for preparing WGBS libraries and Dr. Marta Robertson for the careful reading of the manuscript.
**Author Contribution**

ME.L and E.B conceived the study. ME.L performed the experiments, analyzed sequencing data and wrote the manuscript. D.R. assembled the *F. vesca* genome and wrote the manuscript. C.B. performed experiments and wrote the manuscript. B.D provided plant material and wrote the manuscript. E.B. designed experiments, analyzed data, set up the genome browser and wrote the manuscript.

**Competing interests**

The authors declare they have no conflicts of interest.

**References**

1 World Meteorological Organization (WMO). State of the Global Climate 2021. 2021.

2 Kersey PJ. Plant genome sequences: past, present, future. *Curr Opin Plant Biol* 2019; **48**: 1–8.

3 Roquis D, Robertson M, Yu L, Thieme M, Julkowska M, Bucher E. Genomic impact of stress-induced transposable element mobility in Arabidopsis. *Nucleic Acids Res* 2021; **49**: 10431–10447.

4 Wang P, Zhao FJ, Kopittke PM. Engineering Crops without Genome Integration Using Nanotechnology. *Trends Plant Sci* 2019; **24**: 574–577.

5 Varotto S, Tani E, Abraham E et al. Epigenetics: Possible applications in climate-smart crop breeding. *J Exp Bot* 2020; **71**: 5223–5236.

6 MacKelprang R, Lemaux PG. Genetic Engineering and Editing of Plants: An Analysis of New and Persisting Questions. *Annu Rev Plant Biol* 2020; **71**: 659–687.

7 Zhang H, Lang Z, Zhu J-K. Dynamics and function of DNA methylation in plants. *Nat Rev Mol Cell Biol* 2018; **19**: 489–506.

8 Bewick AJ, Zhang Y, Wendte JM, Zhang X, Schmitz RJ. Evolutionary and experimental loss of gene body methylation and its consequence to gene expression. *G3 Genes, Genomes, Genet* 2019; **9**: 2441–2445.

9 Kenchanmane Raju SK, Ritter EJ, Niederhuth CE. Establishment, maintenance, and biological roles of non-CG methylation in plants. *Essays Biochem* 2019; **63**: 743–755.

10 Kawakatsu T, Huang S shan C, Jupe F et al. Epigenomic Diversity in a Global Collection of Arabidopsis thaliana Accessions. *Cell* 2016; **166**: 492–505.

11 Rendina González AP, Preite V, Verhoeven KJF, Latzel V. Transgenerational effects and epigenetic memory in the clonal plant trifolium repens. *Front Plant Sci* 2018; **871**: 1–11.

12 Tong W, Li R, Huang J et al. Divergent DNA methylation contributes to duplicated gene evolution and chilling response in tea plants. *Plant J* 2021; **106**: 1312–1327.

13 Entrambasaguas L, Ruocco M, Verhoeven KJF, Procaccini G, Marín-Guirao L. Gene body DNA methylation in seagrasses: inter- and intraspecific differences and interaction with transcriptome plasticity under heat stress. *Sci Rep* 2021; **11**: 1–15.

14 Sun RZ, Liu J, Wang YY, Deng X. DNA methylation-mediated modulation of rapid desiccation tolerance acquisition and dehydration stress memory in the resurrection plant
Boea hygrometrica. *PLoS Genet* 2021; **17**: 1–29.

15 Wang J, Nan N, Li N *et al*. A DNA methylation reader-chaperone regulator-transcription factor complex activates OsHKT1;5 expression during salinity stress. *Plant Cell* 2020; **32**: 3535–3558.

16 De Kort H, Panis B, Deforce D, Van Nieuwerburgh F, Honnay O. Ecological divergence of wild strawberry DNA methylation patterns at distinct spatial scales. *Mol Ecol* 2020; **29**: 4871–4881.

17 Sammarco I, Münzbergová Z, Latzel V. DNA Methylation Can Mediate Local Adaptation and Response to Climate Change in the Clonal Plant Fragaria vesca: Evidence From a European-Scale Reciprocal Transplant Experiment. 2022; **13**. doi:10.3389/fpls.2022.827166.

18 Ghaderi N, Siosemardeh A. Esponse to drought stress of two strawberry cultivars (cv. Kurdistan and Selva). *Hortic Environ Biotechnol* 2011; **52**: 6–12.

19 Turechek WW, Peres NA. Heat treatment effects on strawberry plant survival and angular leaf spot, caused by Xanthomonas fragariae, in nursery production. *Plant Dis* 2009; **93**: 299–308.

20 Lv J, Zheng T, Song Z *et al*. Strawberry Proteome Responses to Controlled Hot and Cold Stress Partly Mimic Post-harvest Storage Temperature Effects on Fruit Quality. *Front Nutr* 2022; **8**: 1–19.

21 Orsini F, Alnayef M, Bona S, Maggio A, Gianquinto G. Low stomatal density and reduced transpiration facilitate strawberry adaptation to salinity. *Environ Exp Bot* 2012; **81**: 1–10.

22 Jamali B, Eshghi S, Shahidi-Rad K. Growth and Fruit Characteristics of Strawberry cv. Selva as Affected by Different Application Timing of Salicylic Acid under Saline Conditions. *Int J Fruit Sci* 2015; **15**: 339–352.

23 Choi HG, Moon BY, Kang NJ. Correlation between strawberry (Fragaria ananassa Duch.) productivity and photosynthesis-related parameters under various growth conditions. *Front Plant Sci* 2016; 7. doi:10.3389/fpls.2016.01607.

24 Hidaka K, Dan K, Imamura H. Effect of supplementaHidaka, K., Dan, K., & Imamura, H. (2013). Effect of supplemental lighting from different light sources on growth and yield of strawberry. *Environmental Control in Biology*, 51(1), 41–47.
http://doi.org/10.2525/ecb.51.411 lighting from. *Environ Control Biol* 2013; **51**: 41–47.

25 Mohamed RA, Abdelbaset A-K, Abd-Elkader DY. Salicylic acid effects on growth, yield, and fruit quality of strawberry cultivars. *J Med Act Plants* 2018; **6**: 1–11.

26 Cheng J, Niu Q, Zhang B *et al*. Downregulation of RdDM during strawberry fruit ripening. *Genome Biol* 2018; **19**: 1–14.

27 Zhang L, Wang Y, Zhang X *et al*. Dynamics of phytohormone and DNA methylation patterns changes during dormancy induction in strawberry (Fragaria × ananassa Duch.). *Plant Cell Rep* 2012; **31**: 155–165.

28 Urrutia M, Bonet J, Arús P, Monfort A. A near-isogenic line (NIL) collection in diploid strawberry and its use in the genetic analysis of morphologic, phenotypic and nutritional characters. *Theor Appl Genet* 2015; **128**: 1261–1275.

29 Schmitz RJ, Schultz MD, Lewsey MG *et al*. Transgenerational epigenetic instability is a source of novel methylation variants. *Science (80- )* 2011; **334**: 369–373.

30 Becker C, Hagmann J, Müller J *et al*. Spontaneous epigenetic variation in the Arabidopsis thaliana methylome. *Nature* 2011; **480**: 245–249.
31 Nunn A, Otto C, Stadler PF, Langenberger D. Comprehensive benchmarking of software for mapping whole genome bisulfite data: from read alignment to DNA methylation analysis. *Brief Bioinform* 2021; **00**: 1–9.

32 Nunn A, Can SN, Otto C, Fasold M, Stadler PF, Langenberger D. EpiDiverse Toolkit: a pipeline suite for the analysis of bisulfite sequencing data in ecological plant epigenetics. 2021; **3**: 1–7.

33 Edger PP, VanBuren R, Colle M *et al*. Single-molecule sequencing and optical mapping yields an improved genome of woodland strawberry (Fragaria vesca) with chromosome-scale contiguity. *Gigascience* 2018; **7**: 1–7.

34 Li H. Minimap2: Pairwise alignment for nucleotide sequences. *Bioinformatics* 2018; **34**: 3094–3100.

35 Li Y, Pi M, Gao Q, Liu Z, Kang C. Updated annotation of the wild strawberry Fragaria vesca V4 genome. *Hortic Res* 2019; **6**: doi:10.1038/s41438-019-0142-6.

36 Ou S, Su W, Liao Y *et al*. Benchmarking transposable element annotation methods for creation of a streamlined, comprehensive pipeline. *Genome Biol* 2019; **20**: 1–18.

37 Fojtová M, Kovařík A, Matyášek R. Cytosine methylation of plastid genome in higher plants. *Fact or artefact? Plant Sci* 2001; **160**: 585–593.

38 Tian T, Liu Y, Yan H *et al*. AgriGO v2.0: A GO analysis toolkit for the agricultural community, 2017 update. *Nucleic Acids Res* 2017; **45**: W122–W129.

39 Wu T, Hu E, Xu S *et al*. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. *Innov* 2021; **2**: 100141.

40 Dong C, Xi Y, Chen X, Cheng ZM. Genome-wide identification of AP2/EREBP in Fragaria vesca and expression pattern analysis of the FvDREB subfamily under drought stress. *BMC Plant Biol* 2021; **21**: 295.

41 Xie Z, Nolan TM, Jiang H, Yin Y. AP2/ERF transcription factor regulatory networks in hormone and abiotic stress responses in Arabidopsis. *Front Plant Sci* 2019; **10**: 1–17.

42 Hu Y, Han YT, Wei W *et al*. Identification, isolation, and expression analysis of heat shock transcription factors in the diploid woodland strawberry Fragaria Vesca. *Front Plant Sci* 2015; **6**: 1–16.

43 Wei W, Hu Y, Han YT, Zhang K, Zhao FL, Feng JY. The WRKY transcription factors in the diploid woodland strawberry Fragaria vesca: Identification and expression analysis under biotic and abiotic stresses. *Plant Physiol Biochem* 2016; **105**: 129–144.

44 Li H, Zhou Y, Ma Z, Lu X, Li Y, Chen H. Genome-wide identification and expression analysis of myb transcription factors and their responses to abiotic stresses in woodland strawberry (Fragaria vesca). *Horticulturae* 2021; **7**: doi:10.3390/horticulturae7050097.

45 Bucher E, Reinders J, Mirouze M. Epigenetic control of transposon transcription and mobility in Arabidopsis. *Curr Opin Plant Biol* 2012; **15**: 503–510.

46 Li A, Chen L, Liu Z *et al*. Characterization of strawberry (Fragaria vesca) sequence genome. *bioRxiv* 2018. doi:10.1101/244418.

47 Xiong W, Dooner HK, Du C. Rolling-circle amplification of centromeric Helitrons in plant genomes. *Plant J* 2016; **88**: 1038–1045.

48 Zhang H, Zhu J, Gong Z, Zhu JK. Abiotic stress responses in plants. *Nat Rev Genet* 2021; **0123456789**: doi:10.1038/s41576-021-00413-0.

49 Kort H De, Toivainen T, Nieuwerburgh F Van, Panis B, Hytönen TP, Honnay O. Standing covariation between genomic and epigenomic patterns as source for natural selection in wild strawberry plants. *bioRxiv* 2021; : 2021.03.31.437859.
50 Rakei A, Maali-Amiri R, Zeinali H, Ranjbar M. DNA methylation and physio-biochemical analysis of chickpea in response to cold stress. *Protoplasma* 2016; **253**: 61–76.

51 Qian Y, Hu W, Liao J, Zhang J, Ren Q. The Dynamics of DNA methylation in the maize (Zea mays L.) inbred line B73 response to heat stress at the seedling stage. *Biochem Biophys Res Commun* 2019; **512**: 742–749.

52 Hossain MS, Kawakatsu T, Kim K Do et al. Divergent cytosine DNA methylation patterns in single-cell, soybean root hairs. *New Phytol* 2017; **214**: 808–819.

53 Suriyasak C, Hatanaka K, Tanaka H et al. Alterations of DNA Methylation Caused by Cold Plasma Treatment Restore Delayed Germination of Heat-Stressed Rice (Oryza sativa L.) Seeds. *ACS Agric Sci Technol* 2021; **1**: 5–10.

54 Wibowo A, Becker C, Marconi G et al. Elife-13546-V2. 2016; : 1–27.

55 Omidvar V, Fellner M. DNA methylation and transcriptomic changes in response to different lights and stresses in 7B-1 male-sterile tomato. *PLoS One* 2015; **10**: 1–23.

56 Li R, Hu F, Li B et al. Whole genome bisulfite sequencing methylome analysis of mulberry (Morus alba) reveals epigenome modifications in response to drought stress. *Sci Rep* 2020; **10**: 1–17.

57 Zheng X, Chen L, Li M et al. Transgenerational variations in DNA methylation induced by drought stress in two rice varieties with distinguished difference to drought resistance. *PLoS One* 2013; **8**: 1–13.

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

59 Kou HP, Li Y, Song XX et al. Heritable alteration in DNA methylation induced by nitrogen-deficiency stress accompanies enhanced tolerance by progenies to the stress in rice (Oryza sativa L.). *J Plant Physiol* 2011; **168**: 1685–1693.

60 Verhoeven KJF, Jansen JJ, van Dijk PJ, Biere A. Stress-induced DNA methylation changes and their heritability in asexual dandelions. *New Phytol* 2010; **185**: 1108–1118.

61 Slotkin RK, Martienssen R. Transposable elements and the epigenetic regulation of the genome. *Nat Rev Genet* 2007; **8**: 272–285.

62 Deniz Ö, Frost JM, Branco MR. Regulation of transposable elements by DNA modifications. *Nat Rev Genet* 2019. doi:10.1038/s41576-019-0106-6.

63 Naito K, Zhang F, Tsukiyama T et al. Unexpected consequences of a sudden and massive transposon amplification on rice gene expression. *Nature* 2009; **461**: 1130–1134.

64 Noshay JM, Anderson SN, Zhou P et al. Monitoring the interplay between transposable element families and DNA methylation in maize. *PLoS Genet* 2019; **15**: 1–25.

65 Secco D, Wang C, Shou H et al. Stress induced gene expression drives transient DNA methylation changes at adjacent repetitive elements. *Elife* 2015; **4**: 1–26.

66 Xin Y, Ma B, Xiang Z, He N. Amplification of miniature inverted-repeat transposable elements and the associated impact on gene regulation and alternative splicing in mulberry (Morus notabilis). *Mob DNA* 2019; **10**: 1–13.

67 Nunn A, Rodríguez-Arévalo I, Tandukar Z et al. Chromosome-level Thlaspi arvense genome provides new tools for translational research and for a newly domesticated cash cover crop of the cooler climates. *Plant Biotechnol J* 2022; : 1–20.

68 Naish M, Alonge M, Wlodzimierz P et al. The genetic and epigenetic landscape of the Arabidopsis centromeres. *Science (80–) 2021; **374**. doi:10.1126/science.abi7489.
69 Comai L, Maheshwari S, Marimuthu MPA. Plant centromeres. *Curr Opin Plant Biol* 2017; **36**: 158-167.

70 Janni M, Gullì M, Maestri E *et al*. Molecular and genetic bases of heat stress responses in crop plants and breeding for increased resilience and productivity. *J Exp Bot* 2020; **71**: 3780-3802.

71 dos Santos TB, Ribas AF, de Souza SGH, Budzinski IGF, Domingues DS. Physiological Responses to Drought, Salinity, and Heat Stress in Plants: A Review. *Stresses* 2022; **2**: 113-135.

72 Nahar L, Aycan M, Hanamata S, Baslam M, Mitsui T. Impact of Single and Combined Salinity and High-Temperature Stresses on Agro-Physiological, Biochemical, and Transcriptional Responses in Rice and Stress-Release. *Plants* 2022; **11**. doi:10.3390/plants11040501.

73 Ohama N, Kusakabe K, Mizoi J *et al*. The transcriptional cascade in the heat stress response of arabidopsis is strictly regulated at the level of transcription factor expression. *Plant Cell* 2016; **28**: 181-201.

74 Ohama N, Sato H, Shinozaki K, Yamaguchi-Shinozaki K. Transcriptional Regulatory Network of Plant Heat Stress Response. *Trends Plant Sci* 2017; **22**: 53-65.

75 Zhang CL, Wang YX, Hu X *et al*. An apple AP2/EREBP-type transcription factor, MdWRI4, enhances plant resistance to abiotic stress by increasing cuticular wax load. *Environ Exp Bot* 2020; **180**: 104206.

76 Liang D, Zhang Z, Wu H *et al*. Single-base-resolution methylomes of populus trichocarpa reveal the association between DNA methylation and drought stress. *BMC Genet* 2014; **15**: 1-11.

77 Liu HC, Liao HT, Charng YY. The role of class A1 heat shock factors (HSFA1s) in response to heat and other stresses in Arabidopsis. *Plant, Cell Environ* 2011; **34**: 738-751.

78 Klay I, Gouia S, Liu M *et al*. Ethylene Response Factors (ERF) are differentially regulated by different abiotic stress types in tomato plants. *Plant Sci* 2018; **274**: 137-145.

79 Huang J, Zhao X, Bürger M, Wang Y, Chory J. Two interacting ethylene response factors regulate heat stress response. *Plant Cell* 2021; **33**: 338-357.

80 Centomani I, Sgobba A, D'Addabbo P *et al*. Involvement of DNA methylation in the control of cell growth during heat stress in tobacco BY-2 cells. *Protoplasma* 2015; **252**: 1451-1459.

81 Liao WY, Lin LF, Zheng JL, Wang CC, Yang JH, Chou ML. Identification of heat shock transcription factor genes involved in thermotolerance of octoploid cultivated strawberry. *Int J Mol Sci* 2016; **17**: 1-21.

82 Song H, Wang P, Hou L *et al*. Global analysis of WRKY genes and their response to dehydration and salt stress in soybean. *Front Plant Sci* 2016; **7**: 1-15.