# Epigenomics and bolting tolerance in sugar beet genotypes

Claire Hebrard, Daniel G. Peterson, Glenda Willems, Alain Delaunay, Béline

Jesson, Marc Lefebvre, Steve Barnes, Stéphane Maury

► To cite this version:

Claire Hebrard, Daniel G. Peterson, Glenda Willems, Alain Delaunay, Béline Jesson, et al. Epigenomics and bolting tolerance in sugar beet genotypes. Journal of Experimental Botany, 2016, 67 (1), pp.207 - 225. 10.1093/jxb/erv449. hal-02634949

HAL Id: hal-02634949

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

Submitted on 27 May 2020

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
RESEARCH PAPER

# Epigenomics and bolting tolerance in sugar beet genotypes

Claire Hébrard<fnref n="1" />,<fnref n="2" />,<fnref n="3" />*, Daniel G. Peterson<fnref n="4" />, Glenda Willems<fnref n="3" />, Alain Delaunay<fnref n="1" />,<fnref n="2" />, Béline Jesson<fnref n="5" />, Marc Lefèbvre<fnref n="3" />, Steve Barnes<fnref n="3" /> and Stéphane Maury<fnref n="1" />,<fnref n="2" />,†

Received 27 February 2015; Revised 17 September 2015; Accepted 21 September 2015

Editor: Christine Raines

## Abstract

In sugar beet (Beta vulgaris *altissima*), bolting tolerance is an essential agronomic trait reflecting the bolting response of genotypes after vernalization. Genes involved in induction of sugar beet bolting have now been identified, and evidence suggests that epigenetic factors are involved in their control. Indeed, the time course and amplitude of DNA methylation variations in the shoot apical meristem have been shown to be critical in inducing sugar beet bolting, and a few functional targets of DNA methylation during vernalization have been identified. However, molecular mechanisms controlling bolting tolerance levels among genotypes are still poorly understood. Here, gene expression and DNA methylation profiles were compared in shoot apical meristems of three bolting-resistant and three bolting-sensitive genotypes after vernalization. Using Cot fractionation followed by 454 sequencing of the isolated low-copy DNA, 6231 contigs were obtained that were used along with public sugar beet DNA sequences to design custom Agilent microarrays for expression (56k) and methylation (244k) analyses. A total of 169 differentially expressed genes and 111 differentially methylated regions were identified between resistant and sensitive vernalized genotypes. Fourteen sequences were both differentially expressed and differentially methylated, with a negative correlation between their methylation and expression levels. Genes involved in cold perception, phytohormone signalling, and flowering induction were over-represented and collectively represent an integrative gene network from environmental perception to bolting induction. Altogether, the data suggest that the genotype-dependent control of DNA methylation and expression of an integrative gene network participate in bolting tolerance in sugar beet, opening up perspectives for crop improvement.

**Key words:** Bolting tolerance, differentially expressed gene, differentially methylated region, epigenomics, microarray, sugar beet, vernalization.

Abbreviations: BD, bolting delay; BI, bolting index; CH3, methylation; Cy, cyanine; DEG, differentially expressed gene; DEMS, differentially expressed and differentially methylated sequence; DMR, differentially methylated region; EST, expressed sequence tag; FL1, FLOWERING LOCUS C-like; FT, FLOWERING LOCUS T; HAP, hydroxyapatite; ORF, open reading frame; QTL, quantitative trait locus; R, bolting resistant; RLGS, restriction landmark genome scanning; S, bolting sensitive.

<footnote n="1">¹ Université d'Orléans, Faculté des Sciences, Laboratoire de Biologie des Ligneux et des Grandes Cultures (LBLGC), UPRES EA 1207, 45067 Orléans, France</footnote>

<footnote n="2">² INRA, USC1328 Arbres et Réponses aux Contraintes Hydriques et Environnementales (ARCHE), 45067 Orléans, France</footnote>

<footnote n="3">³ SESVanderHave N.V./S.A., Soldatenplein Z2 nr15, Industriepark, B-3300 Tienen, Belgium</footnote>

<footnote n="4">⁴ Institute for Genomics, Biocomputing & Biotechnology, Mississippi State University, 2 Research Blvd., Box 9627, Mississippi State, MS 39762, USA</footnote>

<footnote n="5">⁵ IMAXIO/HELIXIO, Biopôle Clermont-Limagne, Saint-Beauzire, F-63360, France</footnote>

<footnote>* Present address: Université de Reims Champagne Ardenne, UFR Sciences Exactes et Naturelles, Unité de Recherche Vignes et Vins de Champagne, UPRES EA 4707, Laboratoire Stress Défenses et Reproduction des Plantes, 51687 Reims cedex 2, France.</footnote>

<footnote>† To whom correspondence should be addressed. E-mail: stephane.maury@univ-orleans.fr</footnote>
Introduction

Sugar beet (Beta vulgaris altissima) is a biennial root crop that provides ~15% of the world's sugar (FAO, 2012). This species is generally grown as a spring-sown crop and harvested the following autumn, the whole initial growing season corresponding to vegetative growth and sucrose production and storage. To induce bolting (i.e. rapid stem elongation followed by flowering), sugar beet requires a prolonged exposure to low temperatures between 2 °C and 10 °C (Lexander, 1980), a phenomenon called vernalization. This process requires mitotic activity of somatic cells in the shoot apical meristem for the subsequent formation of floral organs (Wellensiek, 1964). Acquisition of flowering ability is mitotically stable and constitutes a form of cell memory that can be maintained for weeks until optimal light conditions induce bolting (Lang, 1965; Metzger, 1988). Thus, under long-day conditions, bolting, associated with the use of stored sucrose, occurs in vernalized sugar beet and is usually followed by the development of an indeterminate inflorescence. In contrast to sugar beet, wild beets have a facultative requirement for vernalization; both annual and biennial forms have been described (Van Dijk, 2009). The annual habit in beet is under the genetic control of a dominant Mendelian factor termed *B* or 'bolting gene', which over-rides the need for vernalization by promoting flowering (Abegg, 1936; Boudry *et al.*, 1994; El-Mezawy *et al.*, 2002). Annual beets (*BB* or *Bb*) are non-responsive to vernalization, while biennial beets (*bb*) are responsive to vernalization. Recently, two additional loci affecting bolting, named *B2* and *B4*, were identified. The *B2* locus is unlinked to *B* and appears to act epistatically to *B*, while the *B4* locus is linked to the *B* locus and promotes annual bolting independently of *B* (Büttner *et al.*, 2010; Abou-Elwafa *et al.*, 2012).

Depending on the climate, sugar beet can be grown as a spring crop (temperate climates) or as a winter crop (warmer climates) sown in autumn (Mutasa-Göttgens *et al.*, 2012). The development of winter crop varieties in cooler temperate regions could result in at least a 26% increase in sugar yield (Jaggard and Werker, 1999). However, cool temperatures during the growing season may induce premature vernalization and subsequent early bolting in some genotypes, and reduce sugar production. There is a 12% loss in root sugar yield for one premature bolting plant per square metre in the field (Longden, 1989). Bolting tolerance is a quantitative trait reflecting the bolting response of biennial genotypes after a cold exposure. Several parameters are used to measure bolting tolerance, including bolting index (BI), corresponding to the percentage of bolting plants, and bolting delay (BD), which is the average number of days required for a visible initiation of bolting (Trap-Gentil *et al.*, 2011). At present, the selection of late bolting cultivars is based on phenotypic rather than genotypic criteria, and mechanisms controlling the variation in bolting tolerance between beet genotypes are still unclear.

The genetic control of flowering was first elucidated in
the annual model plant *Arabidopsis thaliana* (for a review,
see Andrés and Coupland, 2012). In this plant, vernaliza-
tion involves the epigenetic silencing of the *FLOWERING
LOCUS C* (*FLC*) gene. This MADS-box transcription factor

acts as a floral repressor inhibiting the expression of genes
required to switch the meristem to a floral fate, such as the
*FLOWERING LOCUS T (FT)* gene which encodes the
mobile floral-inductive signal (Michaels and Amasino, 1999;
Abe et al., 2005; Searle et al., 2006). The epigenetic silenc-
ing of *FLC* is stable and can be transmitted to daughter cells
by mitosis, but its expression is restored in each generation
(Michaels and Amasino, 2001; Henderson and Jacobsen,
2007). In the last few years, several vernalization genes were
identified in sugar beet by homology with counterparts of
*Arabidopsis*. Both partial conservation and divergence in the
genetic basis of bolting and flowering were found between
these two plants (Reeves et al., 2007; Pin et al., 2010, 2012).
Indeed, a reverse genetics approach allowed the identification
of a *FLOWERING LOCUS C*-like gene (*BvFL1*). The expres-
sion of *BvFL1* in sugar beet was shown to be down-regulated
during vernalization, but, unlike *AtFLC*, this repression
is not maintained after vernalization (Reeves et al., 2007).
In addition, two paralogues of the *Arabidopsis* *FT* gene,
named *BvFT1* and *BvFT2*, were isolated and shown to have
evolved antagonistic functions (Pin et al., 2010). In biennial
sugar beet, *BvFT1* acts as a floral repressor and is gradually
down-regulated during vernalization, enabling the induction
of *BvFT2*, which is the functional *FT* orthologue and pro-
motes flowering. Recently, the bolting locus *B* was map-based
cloned and sequenced, revealing multiple rearrangements
between the annual and the biennial alleles and the presence
of several linked genes (Pin et al., 2012). One is a pseudo-
response regulator gene called *BOLTING TIME CONTROL
1 (BvBTC1)* acting as an upstream regulator of *BvFT1* and
*BvFT2* and being up-regulated by long days in annuals. In
biennials, the rare recessive allele *Bvbtc1* retained a role as a
promoter of bolting, but imparts reduced sensitivity to pho-
toperiod which is restored by vernalization. In the same way,
the *B2* locus was map-based cloned and sequenced, reveal-
ing a zinc finger transcription factor gene termed *BvBBX19*
that is diurnally regulated and acts like *BvBTC1* upstream of
*BvFT1* and *BvFT2* (Dally et al., 2014). Taken together, these
data illustrate how evolutionary plasticity at a key regulatory
point can enable new life strategies.

Experimental observations, such as post-vernalization
bolting of some *Bvbtc1* and *BvFT2* gene knockdown (RNAi)
plants, indicate that additional mechanisms are involved in
initiation of bolting in sugar beet (Pin *et al.*, 2010, 2012).
In addition, the mechanisms controlling the expression
of these sugar beet orthologues are still poorly known.
As in *Arabidopsis* (for reviews, see Song *et al.*, 2012, 2013;
Baulcombe and Dean, 2014), an epigenetic control of ver-
alization has been shown in sugar beet (Trap-Gentil *et al.*,
2011). Specifically, the time course and amplitude of DNA
methylation variations in the shoot apical meristem have been
shown to correspond to critical points for the induction of
sugar beet bolting and represent an epigenetic component
of the genotypic bolting tolerance. Genotypes with differ-
ent bolting tolerance levels could be distinguished by their
global DNA methylation levels in the shoot apical meristem,
with sensitive genotypes being hypermethylated compared
with resistant genotypes. Recently, a few functional targets
of DNA methylation during vernalization such as RNA METHYLCYTOSINE TRANSFERASES (BvRNMT) were identified in the shoot apical meristem of sugar beet genotypes with distinct bolting tolerance (Hébrard et al., 2013). RNA methylation was shown to vary during cold exposure and between genotypes, to affect BvFL1 mRNA, and to delay bolting in an *Atrnmt* *Arabidopsis* mutant.

Here, high-throughput gene expression/methylation analysis of shoot apical meristems from six sugar beet genotypes that were submitted to 9 weeks of vernalization treatment is reported. A total of 169 differentially expressed genes (DEGs) and 111 differentially methylated regions (DMRs) were identified between bolting-resistant and bolting-sensitive genotypes. Two DEGs and 27 DMRs (all with a higher relative methylation level in bolting-resistant genotypes) mapping to the mitochondrial genome were also found. In addition, 14 sequences were found to be both differentially methylated and differentially expressed. It was shown that these DEGs, DMRs, and joint DEG-DMR sequences (DEMSs) belong to an integrated gene network involved in environment perception, phytohormone signalling, and flowering induction. Altogether, the data suggest that a genotype-dependent remodelling of DNA methylation and expression of an integrative gene network participate in bolting tolerance. These results also up open perspectives for breeding and provide the first candidate molecular markers for bolting tolerance in sugar beet.

## Materials and methods

### Plant material and growth conditions

Six biennial sugar beet genotypes (*Beta vulgaris* *altissima*, provided by SES-VanderHave, Tienen, Belgium) were germinated and grown for 8 weeks at 22 °C under a 16h photoperiod (700 µM m⁻² s⁻¹). Vernalization at 4 °C was then applied for 0, 3, 9, or 18 weeks, as described by Trap-Gentil *et al.* (2011). In this study, the genotypes vernalized for 9 weeks were used for only molecular biology analyses.

### DNA preparation, Cot fractionation, and 454 sequencing

Nuclear DNA was extracted from 500 g of sugar beet leaves from the R1 genotype according to Peterson *et al.* (1997). A Misonix Sonicator 3000 (Farmingdale, NY, USA) was used to shear sugar beet nuclear DNA into fragments with a mean length of 500 bp. The DNA fragments were dissolved in 0.03 M sodium phosphate buffer (SPB) containing 1 mM EDTA and loaded onto a hydroxyapatite (HAP) column equilibrated with 0.03 M SPB. SPB was prepared as detailed in Liu *et al.* (2011). The absorbance (260 nm) of the column eluant was monitored using an Agilent 8530 spectrophotometer (Massy, France). The HAP column was washed with 10 column volumes of 0.03 M SPB and 10 column volumes of 0.12 M SPB. These washes remove carbohydrate contaminants and secondary compounds that may affect DNA reassociation. The DNA was eluted from the column with the addition of 0.5 M SPB. The HAP-purified DNA was concentrated using a Centriplus YM-30 column (EMD Millipore Corporation, Billerica, MA, USA) and sealed in a glass microcapillary tube.

Cot filtration was used to fractionate the HAP-purified sugar beet DNA fragments into single/low-copy and repetitive DNA components. As detailed in Peterson (2005), the genome size of a eukaryote can be used to estimate the Cot1/2 value of its single-copy component using the following formula:

$$y = (\text{Cot1} / 2_{\text{Cot1}} \times \text{G}_{\text{org}}) \div \text{G}_{\text{coli}} \quad (1)$$

where y is the estimated Cot1/2 of single-copy DNA for the organism of interest (i.e. the Cot value at which half the DNA in the single-copy component has reassociated), G_coli is the genome size in base pairs of *Escherichia coli*, G_org is the 1C DNA content of the organism of interest in base pairs, and Cot1/2_coli is the Cot1/2 of *E. coli* DNA (Peterson, 2005). Inserting *E. coli*'s known genome size (4.639 Mbp; Blattner *et al.* 1997), its Cot1/2 value (4.545455 M s⁻¹; Zimmerman and Goldberg, 1977), and the 1C DNA content of sugar beet (1.223 Gbp; Bennett, 1972) into Equation 1 and solving for y yields an estimated Cot1/2 value for sugar beet single-copy DNA of 1200 M s⁻¹. Based on the two Cot decade principle (Peterson, 2005), 80% of the single-copy component of the sugar beet genome is expected to reassociate between 0.1y and 10y (i.e. Cot 120 and Cot 12 000 M s⁻¹). The concentrated sheared DNA sample in 0.5 M SPB (sealed in a glass microcapillary tube as described above) was heated at 95 °C for 10 min and allowed to renature to Cot 940 M s⁻¹. The slowest reassociating component of a Cot curve represents the reassociation of single-copy DNA (Supplementary Fig. S1 available at JXB online). It was chosen to reassociate the sample to the fairly stringent value of Cot 940 M s⁻¹ to ensure that most sequences recovered as unreassociated DNA were truly single copy in nature. The partially reassociated sample was fractionated into single-stranded and double-stranded DNA fractions by HAP chromatography as described by Peterson *et al.* (2002). Those fragments remaining single-stranded at Cot >940 M s⁻¹ were eluted in 0.5 M SPB, concentrated using a Millipore YM-30 column, and transferred into 10 mM aqueous TRIS (pH 8.0). Because recovery of Cot >940 M s⁻¹ was relatively low (i.e. 1–2 µg), the Cot >940 single-stranded DNA was amplified (and concomitantly made double-stranded) using the Qiagen Repli-g Midi Kit (Qiagen, Courtaboeuf, France). The Repli-g amplification resulted in a yield of 185 µg of double-stranded DNA with a concentration of 1840 µg ml⁻¹. The A₂⁶⁰/A₂⁸⁰ ratio for the Repli-g amplified DNA was 1.83. Gel electrophoresis indicates that the Repli-g amplification worked well, although it appears to have preferentially amplified the longer fragments in the reactions (> 2 kb; data not shown).

The Cot >940 DNA was processed for 454 sequencing at GATC-Biotech (Koblenz, Germany) using the GS-FLX technology (Roche, Boulogne Billancourt, France) following established manufacturer's protocols (Margulies *et al.*, 2005). Half an FLX run was performed. Sequence assembly was achieved using the Roche Diagnostic GS de novo Newbler Assembler, and repeated sequences were eliminated using Repeat Masker software (http://www.repeatmasker.org/).

## Microarray probe design

eArray software (https://earray.chem.agilent.com/earray/; Agilent technologies, Massy, France) was used to design transcriptomics and CH3 microarray probes corresponding to 60-mer oligomers based on the 42 004 predicted open reading frames (ORFs) and 6231 contigs, respectively, and on public sequence information from a large collection of sugar beet expressed sequence tags (ESTs; ~30 000) and the mitochondrial genome from the NCBI (http://www.ncbi.nlm.nih.gov/nucest/?term=beta%20vulgaris). Sugar beet genes known to be involved in the floral transition were added, such as FLC (Reeves *et al.*, 2007), FT (Pin *et al.*, 2010), or CONSTANS-LIKE (Chia *et al.*, 2008). Previous studies on sugar beet revealed 39 vernalization markers (Hébrard *et al.*, 2013) and 21 morphogenesis markers (Maury *et al.*, 2012) which were also included.
At least one transcriptomic probe was designed for each sugar beet EST, ORF, or gene. A second probe was designed when probes could not differentiate between sequences. A total of 53 752 probes were generated with this design, and 1227 probes allowed the differentiation between similar transcripts. Moreover, internal Agilent control probes and reproducibility controls, corresponding to 50 probes that were duplicated 10 times, were also used. For hybridization, 8×60K custom microarray slides (Agilent Technologies) were used.

The design of CH3 microarray probes was different depending on defined groups. The probe coverage was maximal, but the space between probes varied for each group. In the first and second groups, corresponding to public sugar beet EST data and the mitochondrial genome, respectively, probes were designed every 100 bp. The contigs obtained by sequencing of Cot-filtered DNA, vernalization, and morphogenesis markers were covered by 50 bp spaced probes. Candidate genes and ORF probes were designed every 30 bp. Moreover, internal Agilent control probes and reproducibility controls, corresponding to 50 probes that were duplicated 10 times, were also used. Finally, 227 242 probes were designed and 1×244K custom microarray slides were used (Agilent Technologies) for hybridization.

### RNA preparation and microarray hybridization

Total RNAs were independently isolated from shoot apical meristems after 9 weeks of vernalizing treatment in three bolting-sensitive (S1, S2, and S3) and three bolting-resistant (R1, R2, and R3) sugar beet genotypes using a Nucleospin® RNA Plant kit (Macherey-Nagel, Hoerdt, France). One shoot apical meristem per genotype was used. To confirm the absence of genomic DNA in all the RNA preparations, a control PCR was performed on RNA. The DNA contamination and RNA integrity with RNA integrity number (RIN) measurement (RIN threshold of 8) were estimated by the use of a bioanalyzer 2100 (Agilent Technologies).

Sugar beet total RNA was used to generate fluorescent cRNA (complementary RNA), according to Agilent Technologies' recommendations. Briefly, to be used as a positive control for labelling and data analysis, RNA spike mix, which is a mix of 10 polyadenylated RNAs synthesized *in vitro* from the *E1A* adenovirus gene, was obtained with the One-Color RNA Spike-In Kit (Agilent Technologies). For each sample, 80 ng of total RNAs were reverse transcribed. Then, cRNA was synthesized by T7 RNA polymerase, and cyanine 3 (Cy3)-labelled CTP was incorporated using an Input Quick Amp Labeling Kit, One Color (Agilent Technologies). Fluorescent cRNA was purified with the RNeasy Mini Kit (Qiagen). Quality control of synthesized fluorescent cRNA was based on cRNA and Cy3 concentration measurement as well as cRNA fragment size. For each sample, three quality thresholds were chosen: a yield of 825 ng of cRNA, a specific activity of 6 pmol<sub>Cy3</sub> g<sub>cRNA</sub><sup>-1</sup>, and a 200–2000 nucleotides fragment size range.

Imaxio (Clermont-Ferrand, France) conducted RNA hybridization following Agilent Technologies' instructions.

Briefly, the labelled cRNAs were broken up in Fragmentation mix (Agilent Technologies) for 30 min at 60 °C. A 2× GEx Hybridization Buffer HI-RPM (Agilent Technologies) was then added and cRNAs were loaded onto the arrays. The arrays were incubated at 65 °C for 17 h in Agilent hybridization chambers and then washed according to the Agilent protocol at room temperature.

### DNA preparation and microarray hybridization

Shoot apical meristems from three bolting-sensitive (S1, S2, and S3) and three bolting-resistant (R1, R2, and R3) sugar beet genotypes exposed to 9 weeks of vernalization were used for DNA isolation. For each genotype, shoot apical meristems from three individuals were pooled, ground into a fine powder in liquid nitrogen, and used for DNA extraction. DNA was extracted using Plant DNAzol Reagent (InVitrogen, Cergy-Pontoise, France) according to the manufacturer's recommendations. To avoid RNA contamination, 100 µg of RNase A ml<sup>-1</sup> Plant DNAzol was added at the beginning of the isolation procedure. DNA was precipitated at -20 °C for at least 30 min using 96% ethanol (2.5 vols) and 3 M Na-acetate (pH 5.2, 0.1 vol.). DNA was then washed using 70% ethanol, vacuum dried, and dissolved in ultrapure water. Concentrations and purity of DNA samples were estimated using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA).

Sugar beet DNA fragments ranging from 0.2 kb to 0.8 kb were obtained by sonication with a VC 505 Vibra-Cell sonicator (Fisher Scientific, Illkirch, France). DNA was isolated by immunoprecipitation according to Weng *et al.* (2009). One immunoprecipitation per DNA sample was realized. Briefly, 11 µg of sonicated DNA were diluted in 250 µl of 1× IP buffer composed of 10 mM Na<sub>3</sub>PO<sub>4</sub> (pH 7.0), 140 mM NaCl, 0.05% (v/v) Triton X-100, and denaturated by boiling for 10 min. Then DNA was cooled on ice for 10 min and 4 µg of DNA from each sample was kept at -20 °C to be used as the input control for microarray analysis. The remaining 7 µg of DNA of each sample was mixed with 10 µg of mouse monoclonal anti-5mC antibodies (Eurogentec, Angers, France) and incubated overnight on a rotator in a 4 °C cold room. The DNA-antibodies mixture was combined with Dynabeads Protein G (InVitrogen) during 2 h at 4 °C. Beads were washed four times with 1× IP buffer using a magnetic rack (Invitrogen). Antibodies were digested overnight on a rotator with 100 µg of proteinase K (Sigma-Aldrich, Saint-Quentin Fallavier, France) in 250 µl of a digestion buffer composed of 50 mM TRIS (pH 8.0), 10 mM EDTA, and 0.5% (w/v) SDS at 50 °C to elute methylated DNA. Isolated DNA was obtained with one phenol/chloroform/isoamyl alcohol (25:24:1) (v/v/v) and two chloroform/isoamyl alcohol (24:1) (v/v) extractions in phase lock gel heavy tubes (Eppendorf, Le Pecq, France) followed by ethanol precipitation. A NanoDrop 1000 spectrophotometer (Thermo Scientific) was used to estimate concentrations of DNA pellets that were solubilized in ultrapure water.

For each sample, labelled DNA was generated according to the recommendations of Agilent Technologies. Briefly, labeling of input control and immunoprecipitated methylated
DNA was performed using a Genomic DNA Enzymatic Labeling Kit (Agilent Technologies). Random primers and the exo-Klenow fragment were used to label input control and immunoprecipitated methylated DNA differentially with Cy3-labelled dUTP and Cy5-labelled dUTP, respectively. Labelled DNA was purified using 30 kDa Amicon columns (Millipore, Molsheim, France). Quality control of synthesized labelled DNA was based on DNA and Cy5/Cy3 concentration measurement. For each sample, three quality minimum thresholds were chosen: a yield of 2.5 µg of labelled DNA, a specific activity of 7–20 pmol$_{Cy5}$ g$_{DNA}^{-1}$, and 18–25 pmol$_{Cy3}$ g$_{DNA}^{-1}$.

Imaxio (Clermont-Ferrand, France) conducted DNA hybridization following Agilent Technologies' instructions. Briefly, Cy5- and Cy3-labelled DNAs were pooled in ddH₂O, and Hybridization mix (Agilent Technologies) was added before incubation at 95 °C for 3 min and then at 37 °C for 30 min. Samples were loaded onto the arrays, and hybridization was performed at 67 °C for 40 h at 20 rpm in Agilent hybridization chambers. Arrays were then washed according to the protocol at room temperature.

## Data collection, normalization, and cluster analysis

The single nucleotide polymorphism (SNP) collection was obtained in a previous study (Adetunji et al., 2014) by genotyping elite sugar beet breeding lines (from SESVanderHave) and wild beet accessions.

A G2505C scanner (Agilent Technologies) was used to scan microarray slides, and extraction of raw signal data was carried out using Feature Extraction 10.7 software. The reproducibility and reliability of each single microarray were evaluated using Quality Control report data (Feature Extraction, Agilent Technologies).

One-colour transcriptomics data were normalized per chip to the 75th percentile and per probe to medians and were transformed to log values. Feature Extraction flag criteria allowed data filtering: non-outlier population, not saturated; uniform and significant probe signal, probe signal above background. A probe signal was considered as 'well above background' if gBGSubSignal was >2.6 gBG_SD. A probe needed to pass all flag criteria for all replicates of at least one experimental condition (bolting-resistant or bolting-sensitive) to pass the filtering. A t-test was performed with MeV software (Saeed et al., 2003) to determine if genes were differentially expressed between experimental conditions. Statistical tests were considered significant at *P*<0.05. Benjamini-Hochberg's correction was used for false-positive discovery. The second filtering criterion which was used is a 2-fold change in transcript levels among conditions. Average linkage hierarchical clustering analysis was applied using Pearson correlation distance.

For the CH3 microarray, data normalization was performed by Feature Extraction software and used to correct cyanine bias with the linear method: the mean signal of probes selected for quality criteria (i.e. signal above background, non-outlier, and whose ranks are similar between the two dyes) allowed the normalization of all probes values. The Cy5/Cy3 log ratio was then calculated for each probe. Data filtering was performed using Feature Extraction flag

criteria: non-outlier population, not saturated; uniform and significant Cy3 and Cy5 probe signal, Cy3 and Cy5 probe signal above background. A probe signal was considered as 'well above background' if g/rBGSubSignal was >2.6 g/rBG_SD. A probe needed to pass all flag criteria for all replicates of at least one experimental condition (bolting-resistant or bolting-sensitive) to pass the filtering. A t-test was performed with MeV software (Saeed et al., 2003) to determine if genes were differentially methylated between experimental conditions. Statistical tests were considered significant at *P*<0.05. A probe signal fold change of two between the two experimental conditions was also chosen as a threshold.

## Bioinformatic analyses

Reads were mapped on the sugar beet genome (Dohm et al., 2014; The Beta Vulgaris Resource: bvseq.molgen.mpg.de) using CLC Genomics Workbench software (Qiagen). Contigs and ORFs were compared with the sugar beet genome using BLASTN. Differentially expressed and differentially methylated genes were annotated (i) by TBLATSX analyses using the *Arabidopsis* database from the National Center for Biotechnology Information (http://blast.ncbi.nlm.nih.gov/Blast.cgi); and (ii) by BLASTX analyses using The Beta Vulgaris Resource (bvseq.molgen.mpg.de). The Gene Ontology (GO) class enrichment was produced with Classification SuperViewer Tool w/ Bootstrap (http://www.blast2go.com/b2ghome; Götz et al., 2008).

## Quantitative PCR

Approximately 500 ng of total RNA were reverse transcribed using SuperScript™ III First-Strand Synthesis SuperMix for qRT-PCR kit (Invitrogen). Primers were designed with QuantPrime software (http://www.quantprime.de/) and are listed in Supplementary Table S1 at JXB online. Constitutively expressed genes encoding, respectively, a *UBIQUITIN* and an *ATP SYNTHASE BETA SUBUNIT* (accession nos BI544002 and BQ490219) were used as internal standards. To confirm the absence of genomic DNA in all the total RNA preparations, a control without reverse transcriptase during the cDNA synthesis was performed. PCRs were prepared using an automated system epMotion 5070 (Eppendorf). Amplification and analysis were performed using a Platinum SYBR Green qPCR Supermix-UDG kit (Eurogentec) on a Mastercycler Realplex (v1.6000; Eppendorf). Two biological and two technical replicates at least were performed for each gene and for six genotypes (three resistant and three sensitive to bolting). Quantitative PCR and microarrays results were compared calculating the Pearson correlation coefficient between qPCR amplification ratios and fold changes.

## Bisulphite sequencing

To determine the cytosine methylation status in CG, CHG, or asymmetric CHH (H could be A, T, or C) contexts, bisulphite
treatment (Epitect bisulfite kit, Qiagen), which results in the conversion of unmethylated cytosines to uracils but does not affect methylated cytosines, was applied on genomic DNA (500 ng) according to the manufacturer's recommendations. Controls for bisulphite conversion efficiency and amplification as well as primer design, PCR, and sequencing procedures have been previously detailed (Trap-Gentil *et al.*, 2011; Hébrard *et al.*, 2013). Primers used for bisulphite sequencing are listed in Supplementary Table S1 at JXB online. Duplicates were performed for each sequence and each genotype.

## Results

### Isolation of sugar beet single-copy sequences by Cot filtration and next-generation sequencing (Cot-SEQ)

In order to obtain more gene sequence data for sugar beet, which has a 1223 Mbp genome (Bennett, 1972) containing >63% repetitive DNA (Flavell *et al.*, 1974), Cot filtration was used to isolate single/low-copy sequences which were subsequently sequenced using 454 technology (Supplementary Fig. S2 at JXB online). An average read length of 353 bp was obtained for the 466 386 reads, corresponding to 117 949 495 nucleotides (Supplementary Fig. S3, Supplementary Table S2). The 466 386 reads were assembled into 6231 contigs, including 2211 large contigs (>500 bases, 832 bases average size, 14 860 bases largest size), and covered 2 953 454 bases (Supplementary Table S2). Large contigs have been deposited at DDBJ/EMBL/GenBank under the accession no. JMBQ00000000. About 78% of the total reads remained as singletons. Bioinformatic analyses revealed that single-copy sequences represented 98.6% of the covered bases while repeated sequences covered only 1.4%, mostly with small RNAs and simple repeats (Supplementary Table S3). Only 1.0% of the contigs (Supplementary Table S4) were mapped on the recent version of the reference sugar beet genome. Analysis of contigs identified thousands of potential ORFs that were blasted against public databases. These ORFs revealed very low sequence redundancy (1.5%), and 3.9% showed significant homology to previously annotated genes (Supplementary Tables S3, S5). GO annotations indicate that those sequences with homology to known genes represent a wide range of biological processes (Supplementary Fig. S4A) and molecular functions (Supplementary Fig. S4B). ORFs were then mapped on the sugar beet genome when it was released (Dohm *et al.*, 2014), and 0.52% of ORFs (i.e. 218) were mapped on the sugar beet genome (Supplementary Table S6). Among them, 184 ORFs showed sequence homology to genes, but only 84 distinct genes were identified. Using a Cot-SEQ strategy, information on single/low-copy sequences that were previously undescribed was generated and thus additional information was provided that was leveraged in designing microarrays for sugar beet.

### Bolting tolerance and genetic clustering of sugar beet genotypes

Six sugar beet genotypes were kept at 4 °C for up to 18 weeks and then characterized for bolting using the percentage of bolting plants (the bolting index, BI) after 0, 3, 9, or 18

weeks of cold and the average number of days required for bolting initiation (the bolting delay, BD; Fig. 1A). The BI increased with the duration of cold, depending on the genotypes. Genotypes S1, S2, and S3 started to bolt after 9 weeks at 4 °C and reached 75–100% of the BI at the end of the treatment, with BDs ranging from 17 d to 55 d. Conversely, genotypes R1, R2, and R3 did not bolt, or showed weak bolting responses during the treatment (R3 with a 15% BI and BD of 85 d after 18 weeks of cold exposure). Therefore, S1–S3 genotypes were considered as bolting sensitive compared with R1–R3 which were classified as bolting resistant. The shoot apical meristems were collected after 9 weeks of cold exposure for further microarray analyses. Indeed, this

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

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Characterization of the six sugar beet genotypes during vernalization. (A) Bolting index of genotypes after 0, 3, 9, or 18 weeks of vernalization treatment (4 °C). S1, S2, and S3 (open circles, squares, and triangles) were considered as bolting sensitive, while R1, R2, and R3 (filled circles, squares, and triangles) were considered as bolting resistant. On the right is indicated the bolting delay after 18 weeks of vernalization for each genotype. (B) Genetic clustering based on 708 single nucleotide polymorphisms for the three bolting-resistant (R1–R3) genotypes, the three bolting-sensitive (S1–S3) genotypes, and nine other genotypes with distinct bolting tolerance levels after 9 weeks of vernalization treatment.</figcaption></figure>
cold duration was sufficient to induce bolting in the three bolting-sensitive genotypes when transferred to optimal temperature and photoperiod conditions, but not in the three resistant genotypes.

A hierarchical genetic clustering (Fig. 1B) based on 708 SNPs (explaining ~46% of the variance) showed that the six genotypes were classified independently from their bolting characteristics. These SNPs were obtained from a previous study by genotyping elite sugar beet breeding lines (from SESVanderHave) and wild beet accessions (Adetunji et al., 2014).

### Transcriptomic microarray analysis of vernalized sugar beet genotypes with distinct bolting tolerance

To explore transcriptomic variations between sugar beet genotypes with distinct bolting tolerance levels, a custom long oligonucleotide (60-mer) microarray (8 × 56k) was generated using the Agilent SurePrint ink-jet technology and hybridized with total RNA extracted from the shoot apical meristems of each genotype after 9 weeks of cold exposure (Supplementary Fig. S2 at JXB online). Transcriptomic data are available in the GEO database (GEO accession no. GSE55951).

Of the 53 752 probes (29 094 from public databases and 24 658 from Cot-SEQ; Table 1), 22 448 (41.8%) remained after the filtration and normalization steps (Supplementary Fig. S5 at JXB online), mainly from the public databases probes. Using a fold change threshold of two, 1580 genes were found to be differentially expressed between S and R genotypes. According to a t-test (*P*<0.05), 169 of these 1580 genes revealed significant differences between genotypes and were considered as 'differentially expressed genes' (DEGs; Table 1; Supplementary Table S7). Microarray transcript abundance data were independently validated by quantitative reverse transcription-PCR (qRT-PCR) for a sample of genes (Supplementary Fig. S6), revealing a strong positive correlation between fold changes obtained from microarray and qPCR data (Pearson correlation coefficient *r*=0.986 at *P*<0.005, *n*=6). Figure 2A corresponds to a heatmap showing the 82 sequences with a higher relative expression level in S genotypes and the 87 sequences with a higher relative expression level in R genotypes. Genes with a lower relative

expression level exhibited homogeneous expression levels between genotypes of the same group (R or S), while genes with a higher relative expression level exhibited variations among genotypes, allowing a clustering in agreement with their bolting tolerance (Fig. 2A).

TBLASTX analysis of the 169 DEGs revealed annotations for 125 of them using the *Arabidopsis* gene database (Supplementary Table S7 at JXB online). Based on the TAIR percentage analysis, GO classification of these 125 sequences according to the 'Biological Process' component placed the largest proportion of genes in undefined categories such as 'Other metabolic processes' (31.9%), 'Unknown biological process' (31.9%), 'Other cellular processes' (27.9%), but also in 'Response to abiotic and biotic stimulus' (24.5%), 'Response to stress' (24.5%), and 'Developmental process' (12.0%; data not shown). These three last categories corresponded to 34 distinct sequences, including 17 sequences that were up-regulated in R genotypes and 17 in the S genotypes (Fig. 2B; Supplementary Table S7). Among them, eight genes could be regrouped in a 'Response to phytohormones' class (seven sequences up-regulated in R and one in S genotypes). Of the four sequences that were grouped together in the 'Response to cold' class, three had a higher relative expression level in S and one had a higher relative expression level in R. To validate the experimental approach, sequences directly involved in flowering regulation were studied, and three of them were found to have a higher relative expression level in S [SHORT VEGETATIVE PHASE (BvSVP) and AGAMOUS-LIKE 24 (BvAGL24)] and one a higher relative expression level in R [FRUITFULL (BvFUL); Table 2]. Altogether, the bolting tolerance transcriptomic analysis revealed an integrated network of genes that are involved in environmental perception (cold), phytohormone signalling, and flowering induction, and have a distinct relative expression level between R and S genotypes.

### Methylation microarray analysis of vernalized sugar beet genotypes with differing bolting tolerance

To explore epigenetic variation between sugar beet genotypes with distinct bolting tolerance levels, a custom

<table><caption>**Table 1.** Number of transcriptomic and methylation probes after each step of the microarray design and analysis</caption><thead><tr><th></th><th>Origin of the sequences</th><th>Submitted sequences</th><th>Designed probes</th><th>Filtered probes</th><th>DEGs or DMRs</th></tr></thead><tbody><tr><td rowspan="4">Transcriptomic</td><td rowspan="3">Public database</td><td>154 mitochondrial ORFs</td><td>132</td><td>106</td><td>2</td></tr><tr><td>30 023 ORFs from ESTs</td><td>28 918</td><td>21 043</td><td>159</td></tr><tr><td>22 candidate genes</td><td>44</td><td>26</td><td>3</td></tr><tr><td>Cot fraction</td><td>42 004 ORFs</td><td>24 658</td><td>1273</td><td>5</td></tr><tr><td>Total</td><td>72 203</td><td>53 752</td><td>22 448</td><td>169</td></tr><tr><td rowspan="4">Methylation</td><td rowspan="3">Public database</td><td>1 mitochondrial genome</td><td>3688</td><td>3211</td><td>27</td></tr><tr><td>29 585 ESTs</td><td>185 425</td><td>14 069</td><td>110</td></tr><tr><td>22 candidate genes</td><td>1426</td><td>109</td><td>1</td></tr><tr><td>Cot fraction</td><td>6231 contigs</td><td>36</td><td>2185</td><td>0</td></tr><tr><td>Total</td><td>35 839</td><td>227 242</td><td>19 574</td><td>138</td></tr></tbody></table>

DEG, differentially expressed gene; DMR, differentially methylated region.

The Cot fraction corresponds to single-copy sequences obtained after Cot filtration and 454 sequencing of sugar beet genomic DNA.

DEGs and DMRs correspond to sequences showing both significantly different signals (t-test analysis with *P*<0.05 and *P*<0.01, respectively) and fold change >2 between resistant and sensitive genotypes.
A

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

* Response to abiotic or biotic stimulus (1.54 E-03)

* Response to abiotic or biotic stimulus (4.85 E-04)

* Response to stress (2.05 E-04)

B

<table><thead><tr><th>GO categories</th><th>Up-regulated genes in R genotypes</th><th>Up-regulated genes in S genotypes</th></tr></thead><tbody><tr><td><strong>Response to abiotic stimulus/stress</strong></td><td></td><td></td></tr><tr><td>Response to phytohormon</td><td>7</td><td>1</td></tr><tr><td>Response to salt stress</td><td>2</td><td>3</td></tr><tr><td>Response to cold</td><td>1</td><td>3</td></tr><tr><td>Response to oxidative stress</td><td>2</td><td>2</td></tr><tr><td>Detoxification of nitrogen compound</td><td>1</td><td>1</td></tr><tr><td>Photoprotection</td><td>0</td><td>1</td></tr><tr><td>Response to water</td><td>0</td><td>1</td></tr><tr><td>Cellular response to phosphate starvation</td><td>0</td><td>1</td></tr><tr><td>Response to wounding</td><td>1</td><td>0</td></tr><tr><td>Response to cadmium ion</td><td>1</td><td>0</td></tr><tr><td><strong>Response to biotic stimulus/stress</strong></td><td></td><td></td></tr><tr><td>Response to other organism</td><td>1</td><td>1</td></tr><tr><td>Defense response to fungus</td><td>0</td><td>1</td></tr><tr><td><strong>Developmental process</strong></td><td></td><td></td></tr><tr><td>Regulation of floral development</td><td>0</td><td>2</td></tr><tr><td>Regulation of vernalization response</td><td>1</td><td>0</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Transcriptomic characterization of the shoot apical meristem in sugar beet genotypes after vernalization. (A) Heatmap representation of the 169 differentially expressed genes (fold change R versus S >2 and t-test P-value <0.05) between bolting-resistant (R1–R3) and bolting-sensitive (S1–S3) genotypes after 9 weeks of vernalization treatment. On the left are indicated the main biological processes in which sequences could be involved and the corresponding P-value. (B) GO categories for genes involved in response to abiotic or biotic stimulus/stress and developmental process. (This figure is available in colour at JXB online.)</figcaption></figure>

long-oligonucleotide (60-mer) microarray (1×244k) was generated using the Agilent SurePrint ink-jet technology and hybridized with either total or methylated DNA (isolated by methylated DNA immunoprecipitation; see the Materials and methods for details) extracted from the shoot apical meristems of each genotype after 9 weeks of cold

exposure (Supplementary Fig. S2 at JXB online). Epigenetic data are available in the GEO database (GEO accession no. GSE56415).

Of the 227 242 probes (190 539 from public databases and 36 703 from Cot-SEQ; Table 1), 19 574 (8.6%) remained after the filtration and normalization steps (Supplementary Fig. S7
<table><caption>**Table 2.** Sugar beet genes involved in the flowering pathway showing significant normalized values (≥[0.5]) for their methylation and expression levels</caption><thead><tr><th rowspan="2">Accession</th><th rowspan="2">Gene</th><th rowspan="2">Reference</th><th colspan="4">Methylation</th><th colspan="3">Expression</th></tr><tr><th>Position of probes groups</th><th>No. of probes per group</th><th>Methylation level in S<sup>a</sup></th><th>Methylation level in R<sup>a</sup></th><th>Expression level in S<sup>b</sup></th><th>Expression level in R<sup>a</sup></th></tr></thead><tbody><tr><td>HR807147</td><td>RNA METHYLCYTOSINE TRANSFERASE<br/>b (BvRNMTb)</td><td>Hébrard et al. (2013)</td><td>0-260</td><td>2</td><td>1.64</td><td>0.58</td><td>-</td><td>-</td></tr><tr><td>BQ592158</td><td>FLOWERING LOCUS VE (BvFVE)</td><td>Abou-Elwafa et al. (2010)</td><td>60-120</td><td>1</td><td>1.60</td><td>0.63</td><td>0.58</td><td>0.00</td></tr><tr><td>HM448911</td><td>FLOWERING LOCUS T2 (BvFT2)</td><td>Pin et al. (2010)</td><td>0-240</td><td>5</td><td>2.25</td><td>0.94</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td></td><td>210-630</td><td>9</td><td>-2.25</td><td>-2.24</td><td></td><td></td></tr><tr><td></td><td></td><td></td><td>1080-2070</td><td>17</td><td>1.42</td><td>0.75</td><td></td><td></td></tr><tr><td></td><td></td><td></td><td>2820-4260</td><td>6</td><td>1.31</td><td>0.41</td><td></td><td></td></tr><tr><td>HM448909</td><td>FLOWERING LOCUS T1 (BvFT1)</td><td>Pin et al. (2010)</td><td>60-260</td><td>3</td><td>1.60</td><td>1.39</td><td>-</td><td>-</td></tr><tr><td></td><td></td><td></td><td>1100-1773</td><td>8</td><td>1.48</td><td>0.59</td><td></td><td></td></tr><tr><td></td><td></td><td></td><td>2763-4893</td><td>6</td><td>0.92</td><td>0.28</td><td></td><td></td></tr><tr><td>DQ189214</td><td>FLOWERING LOCUS C (BvFL1)</td><td>Reeves et al. (2007)</td><td>790-1090</td><td>2</td><td>1.46</td><td>1.00</td><td>-0.17</td><td>0.30</td></tr><tr><td></td><td></td><td></td><td>5950-6130</td><td>5</td><td>-2.58</td><td>-2.56</td><td></td><td></td></tr><tr><td>BQ584677</td><td>FRUITFULL (BvFUL)</td><td>This study (e-value 6E-83)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-0.45</td><td>0.55</td></tr><tr><td>BQ588867</td><td>SHORT VEGETATIVE PHASE (BvSV/P)</td><td>This study (e-value 7E-22)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.00</td><td>-0.65</td></tr><tr><td>EU437782</td><td>CONSTANS (BvCOL1)</td><td>Chia et al. (2008)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>0.65</td><td>0.36</td></tr><tr><td>BQ488144</td><td>AGAMOUS-LIKE 24 (BvAGL24)</td><td>This study (e-value 4E-30)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>1.04</td><td>-0.75</td></tr></tbody></table>

<sup>a</sup>Average of the methylation levels for the indicated number of probes and for the three bolting-sensitive (S) or bolting-resistant (R) genotypes.

<sup>b</sup>Average of expression levels for the three S or R genotypes.

—, The corresponding probes were eliminated after filtration.
at JXB online), with 11.1% from Cot-SEQ probes. Using a fold change threshold of two, 1526 sequences were found to differ for their methylation levels between S and R genotypes. A *t*-test (*P*<0.05) revealed 707 out of these 1526 sequences that showed significant differences between genotypes, and 111 of them (with *P*<0.01) were considered as 'differentially methylated regions' (DMRs; Table 1; Supplementary Table S8). The DNA methylation ratios detected in this microarray experiment were independently validated by bisulphite sequencing (Supplementary Fig. S8). A heatmap is shown in Fig. 3A with the 111 sequences which all had a higher relative methylation level in S genotypes, allowing a clustering in agreement with the bolting tolerance of the six genotypes.

From the 111 DMRs, TBLASTX analysis detected 50 with annotations in the *Arabidopsis* gene database (46 with known function and four uncharacterized; Supplementary Table S8 at JXB online). Based on the TAIR percentage analysis, GO classification of these 50 sequences according to the ‘Biological Process’ component placed the largest proportions of genes in undefined categories such as ‘Other metabolic processes’ (18.9%) and ‘Other cellular processes’ (17.6%), but also in ‘Developmental process’ (8.8%), ‘Response to abiotic and biotic stimulus’ (8.1%), and ‘Response to stress’ (7.0%; data not shown). These three last categories corresponded to 17 distinct sequences, all with a higher relative methylation level in S genotypes (Fig. 3B). Among them, nine sequences could be grouped in a ‘Response to abiotic stimulus/stress’ class and six sequences in ‘Developmental process’. To validate the experimental approach, sequences directly involved in flowering regulation were studied, and five of them were found to be significantly hypermethylated in S [such as *FLOWERING LOCUS VE* (*BvFVE*); Table 2]. The methylation data indicated that bolting tolerance is associated with a gene network involved in the response to environment (cold) or in the induction of flowering, which exhibits distinct methylation states between S and R genotypes, these genes having a higher relative methylation level in S genotypes.

Transcriptomic and methylation microarray analysis of
the sugar beet mitochondrial genome

The sugar beet mitochondrial genome (368 kbp; GenBank
accession no. BA000009) was used to design oligonucleotides
for the custom transcriptomic (132 probes) and methylation
(3688 probes) microarrays. After normalization, 106 tran-
scriptomic probes and 3211 methylation probes were kept
(Table 1) and mapped on the mitochondrial genomic map for
R and S groups (Fig. 4). After filtration (t-test at P<0.05), only
two DEGs were identified: ORF104a and CYTOCHROME
OXIDASE C subunit 2 (COX2). DNA methylation probes
showed a succession of peaks all along the genome in both
S and R genotypes, but R genotypes had a significantly
higher relative methylation level compared with the S lines
(at P<0.001). The non-coding regions also had a significantly
higher relative methylation level compared with the coding
regions (at P<0.01). After filtration (t-test at P<0.05), 27
DMRs (Table 1) were found, with 18 of them mapping on
two main loci containing three ORFs (ORF152, ORF102b,

and ORF192) that exhibited lower normalized expression
values in R genotypes than in S genotypes. These data dem-
onstrated that the mitochondrial genome is also affected by
transcriptomic and epigenetic remodelling between sugar
beet genotypes.

Candidate sequences related to the genotypic bolting
tolerance

In order to identify new candidate genes for bolting tol-
erance, the 14 sequences found to be both differentially
methylated and differentially expressed between R and S
genotypes were selected (Table 3). TBLASTX analyses pro-
posed significant annotations for nine of them: four are
involved in biosynthesis of various compounds such as
suberin and cutin polymers (FG345845, a HXXXD-TYPE
ACYL-TRANSFERASE-LIKE PROTEIN), triterpenoids
(FG343208, a BETA-AMYRIN SYNTHASE), carot-
enoids (BQ582907, a CAROTENOID EPSILON-RING
HYDROXYLASE LUT1/CYP97C1), or ethylene (FG344953,
a 1-AMINOCYCLOPROPANE-1-CARBOXYLATE
OXIDASE-LIKE PROTEIN; Table 3); BQ586261, BQ489455,
and BQ585496 encode a GERMIN-LIKE PROTEIN
5 (GLP5/PDGLP1), an ATPASE E1-E2, and an SPX
DOMAIN CONTAINING PROTEIN 4 (SPX4), respectively,
and are involved in early development, endoplasmic reticulum
to Golgi vesicle-mediated transport, and response to phos-
phate starvation, respectively. Interestingly, one sequence is
involved in response to environmental stress including cold
(FG345694, a CYSTATIN 6 CYSB), and another one in the
transcriptional regulation of specific signalling pathways
such as flowering time (BQ582428, RNA POLYMERASE II
TRANSCRIPTION MEDIATOR; Table 3). A correlation
between methylation and expression values has been detected
for these 14 sequences (Pearson r = -0.34 and Spearman ρ=
-0.41 at P<0.01; Fig. 5; Supplementary Fig. S9A, B at JXB
online). Twenty-two candidate genes known to be involved
in sugar beet bolting and flowering were also studied. Nine
of them, displaying the highest variations of normalized
values for expression and/or methylation between R and S
genotypes, are detailed in Table 2. Four genes had a higher
relative expression level in S (SVP and the flowering activa-
tors FVE, AGL24, and CO), while two genes had a higher
relative expression level in R genotypes (the flowering repres-
sor FLC and FUL). In addition, five genes (RNMTb, FVE,
FT1, FT2, and FLC) had a higher relative methylation level in
S genotypes. Using the microarray data collected from sugar
beet ESTs and bolting/flowering candidate genes (Table 2), a
bolting tolerance model integrating phytohormone signalling
genes, response to environment genes, and key bolting/flower-
ing genes that were differentially expressed and/or methylated
between bolting-resistant and bolting-sensitive genotypes
after 9 weeks of vernalization is proposed (Fig. 6). The model
suggests that upstream regulators of flowering induction have
a higher relative methylation level after vernalization in the
shoot apical meristem of S genotypes and that their expres-
sion levels are in agreement with a potential induction state of
flowering in S genotypes and repression state in R genotypes.
A

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

B

<table><thead><tr><th>GO category</th><th>Hypermethylated genes in S genotypes</th></tr></thead><tbody><tr><td><strong>Response to abiotic stimulus/stress</strong></td><td></td></tr><tr><td>Response to abscisic acid stimulus</td><td>2</td></tr><tr><td>Response to oxidative stress</td><td>2</td></tr><tr><td>Response to cold</td><td>1</td></tr><tr><td>Response to salt stress</td><td>1</td></tr><tr><td>Response to cadmium ion</td><td>1</td></tr><tr><td>Response to heat shock</td><td>1</td></tr><tr><td>Cellular response to phosphate starvation</td><td>1</td></tr><tr><td><strong>Response to biotic stimulus/stress</strong></td><td></td></tr><tr><td>Response to xenobiotic stimulus</td><td>2</td></tr><tr><td><strong>Developmental process</strong></td><td></td></tr><tr><td>Vegetative to reproductive phase transition of meristem</td><td>3</td></tr><tr><td>Regulation of flower development</td><td>2</td></tr><tr><td>Leaf morphogenesis</td><td>1</td></tr></tbody></table>

**Fig. 3.** Methylation characterization of the shoot apical meristem in sugar beet genotypes after vernalization. (A) Heatmap representation of the 111 differentially methylated regions (fold change R versus S >2 and t-test P-value < 0.01) between bolting-resistant (R1–R3) and bolting-sensitive (S1–S3) genotypes after 9 weeks of vernalization treatment. On the left are indicated the main biological processes in which sequences could be involved and the corresponding P-value. (B) GO categories for genes involved in response to abiotic or biotic stimulus/stress and developmental process. (This figure is available in colour at JXB online.)

Discussion

*Epigenomic approaches in a crop plant*

Genomic resources for crop plants such as sugar beet are crucial to elucidate complex physiological processes such

as vernalization and bolting. When this study was started, the only public resource for sugar beet genomic sequences was a collection of ~30 000 ESTs in the NCBI dbEST database. However, the sugar beet genome has been released recently, opening up new perspectives (Dohm *et al.*, 2014).
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Expression and DNA methylation level analysis in the mitochondrial genome of sugar beet genotypes after vernalization. The three bolting-resistant (R, red lines) and the three bolting-sensitive (S, blue lines) genotypes have the same mitochondrial haplotype ('O' type, i.e. normal type mitochondrial genome leading to fertility, in contrast to the cytoplasmic male-sterile type) and were vernalized during 9 weeks at 4 °C. Green squares represent ORFs, purple squares represent rRNAs and tRNAs, and pink squares represent CpG islands. Black frames correspond to differentially methylated region-rich loci. Significantly different values for transcriptomic probes are marked with an asterisk (*P*<0.05).</figcaption></figure>

EST sequencing only provides data on the coding regions expressed in a given tissue or set of tissues. Therefore, Cot filtration, another reduced-representation sequencing technique without the bias for expressed sequences, was employed to obtain low-copy sequence regions including expressed genes, promoters, introns, and non-expressed genes (Peterson, 2005). In this regard, Cot filtration proved to be an effective tool for sugar beet genome analysis (567 Mbp and 42.3% repetitive sequences; Dohm et al., 2014). The efficiency of Cot filtration as a tool for gene enrichment has been demonstrated in several crop plants (Peterson et al., 2002; Yuan et al., 2003; Lamoureux et al., 2005; Paterson, 2006). In the present study, the sugar beet Cot >940 fraction contained >98% single-copy sequences, confirming the efficacy of the Cot filtration, and led to the assembly of 6231 contigs after sequencing. Only 5% of transcriptomic probes and 11% of methylation probes corresponding to Cot sequences remained after the filtration and normalization steps. Several explanations are possible, such as the low quality of the sequencing and assembly into contigs (1.0% mapped on sugar beet genome), the genetic distance between the genotype and the reference genome, or the quality of the actual version of the reference genome. In addition, fragments from the same gene have probably been re-assembled into several contigs, artificially increasing the number of ORFs used for probe design: 42 000 ORFs were predicted using the data, while only 27 421 protein-coding genes are predicted in the published genome assembly (Dohm et al., 2014). Therefore, a high proportion of predicted ORFs were 'unreal' and only 218 ORFs mapped on the genome, representing a very low proportion of potentially 'efficient' probes. In addition, 184 ORFs showed sequence homology to genes, but only 84 distinct genes were identified, revealing high redundancy. However, the aim of this study was to expand the sugar beet public database to enable development of the first custom sugar beet expression and methylation microarrays. Thus, 1273 additional filtered probes for expression analysis and 2185 additional filtered probes for methylation analysis were obtained from the sequencing data (Table 1) to complement the public genomic resources. In recent years, microarray technology has demonstrated the power of the high-throughput study of gene expression in unravelling key processes in plant biology (Galbraith, 2006; Clarke and Zhu, 2006). A number of strategies for high-throughput detection of DNA methylation based on microarrays have led to important discoveries (Zhang et al., 2006). In addition, microarrays have become especially relevant for species where little genome information is available and where intensive laboratory work is necessary to gain insight into a particular biological process, as is the case in sugar beet, as well as to identify candidate target genes for future breeding programmes (Martinez-Godoy et al., 2008; Mascarell-Creus et al., 2009). Recently, MeDIP-chip was compared with whole-genome bisulphite sequencing for the detection of DMRs and was also shown to be a reasonable alternative in population epigenetic studies (Wardenaar et al., 2013).
<caption>**Table 3.** Differentially expressed and differentially methylated genes between bolting-resistant and bolting-sensitive sugar beet genotypes</caption>

Genes were considered to be differentially expressed between genotypes when the fold change was >2 and the *t*-test *P*-value <0.05, and differentially methylated between genotypes when the fold change was >2.

<table><thead><tr><th rowspan="2">Sugar beet<br>accession</th><th rowspan="2">Sugar beet<br>gene</th><th rowspan="2">Arabidopsis thaliana<br>description</th><th rowspan="2">Biological process and/or<br>molecular function</th><th rowspan="2">Arabidopsis AGI</th><th rowspan="2">E-value</th><th colspan="2">Methylation</th><th colspan="2">Expression</th></tr><tr><th>Methylation<br>level in S<sup>a</sup></th><th>Methylation<br>level in R<sup>a</sup></th><th>Expression<br>level in S<sup>a</sup></th><th>Expression<br>level in R<sup>a</sup></th></tr></thead><tbody><tr><td>FG345845</td><td>Bv_39660_duks.t1</td><td>HXXXD-type acyl-transferase-like<br>protein</td><td>Suberin and cutin polymer<br>biosynthesis</td><td>AT5G07850</td><td>5.00E-39</td><td>-0.67</td><td>1.16</td><td>1.49</td><td>-1.62</td></tr><tr><td>BQ589458</td><td>No significant hit</td><td>No significant hit</td><td>NA</td><td>NA</td><td>NA</td><td>0.53</td><td>-0.52</td><td>0.85</td><td>-0.66</td></tr><tr><td>BQ586261</td><td>Bv6_125340_epqq.t1</td><td>Germin-like protein 5 (GLP5/<br>PDGLP1)</td><td>Early development, defence, apoptosis</td><td>AT1G09560</td><td>6.00E-71</td><td>0.72</td><td>-0.42</td><td>0.30</td><td>-0.75</td></tr><tr><td>FG343208</td><td>Bv7u_181680_oyth.t1</td><td>Beta-amyrin synthase</td><td>Triterpenoid synthesis</td><td>AT1G78950</td><td>4.00E-74</td><td>1.15</td><td>-0.04</td><td>-0.61</td><td>0.93</td></tr><tr><td>BQ582907</td><td>Bv7_160760_mmcw.<br>t1</td><td>Carotenoid epsilon-ring<br>hydroxylase (LUT1/CYP97C1)</td><td>Carotenoid biosynthesis, leaf<br>morphogenesis</td><td>AT3G53130</td><td>5.00E-40</td><td>1.22</td><td>-0.07</td><td>1.01</td><td>-0.89</td></tr><tr><td>FG344953</td><td>Bv3_049870_xxoz.t1</td><td>1-Aminocyclopropane-1-<br>carboxylate oxidase-like protein 9</td><td>Ethylene biosynthesis</td><td>AT1G06620</td><td>6.00E-74</td><td>1.40</td><td>0.17</td><td>0.63</td><td>-0.41</td></tr><tr><td>BQ582428</td><td>Bv4u_091120_gpyrt.t1</td><td>RNA polymerase II transcription<br>mediator activity</td><td>Transcriptional regulation of specific<br>signalling pathways</td><td>AT3G01435</td><td>2.00E-52</td><td>1.43</td><td>0.21</td><td>-1.01</td><td>0.88</td></tr><tr><td>CF543348</td><td>No significant hit</td><td>No significant hit</td><td>NA</td><td>NA</td><td>NA</td><td>1.59</td><td>0.53</td><td>-2.36</td><td>2.26</td></tr><tr><td>CK136737</td><td>Bv4_074360_ejeh.t1</td><td>Cystatin 6 (CYSB)</td><td>Cysteine proteinase inhibition, response<br>to stress</td><td>AT3G12490</td><td>4.00E-35</td><td>1.69</td><td>0.57</td><td>0.99</td><td>-0.93</td></tr><tr><td>BQ589845</td><td>Bv4u_091720_azjm.t1</td><td>No significant hit</td><td>NA</td><td>NA</td><td>NA</td><td>1.89</td><td>0.68</td><td>-2.35</td><td>2.03</td></tr><tr><td>FG345537</td><td>Bv9_208140_cjcm.t1</td><td>No significant hit</td><td>NA</td><td>NA</td><td>NA</td><td>1.99</td><td>0.29</td><td>-1.01</td><td>1.60</td></tr><tr><td>BQ489455</td><td>Bv7_169260_ewnm.<br>t1</td><td>ATPase E1-E2 type family protein</td><td>ER to Golgi vesicle-mediated transport</td><td>AT1G17500</td><td>2.00E-37</td><td>2.13</td><td>0.63</td><td>-2.59</td><td>1.44</td></tr><tr><td>BQ585496</td><td>Bv2_042300_ners.t1</td><td>SPX domain-containing protein<br>4 (SPX4)</td><td>Response to phosphate starvation</td><td>AT5G15330</td><td>8.00E-23</td><td>2.48</td><td>1.16</td><td>0.35</td><td>-0.68</td></tr><tr><td>BQ488791</td><td>No significant hit</td><td>No significant hit</td><td>NA</td><td>NA</td><td>NA</td><td>2.52</td><td>1.43</td><td>-0.45</td><td>1.03</td></tr></tbody></table>

R, bolting resistant; S, bolting sensitive; NA, not available.

<sup>a</sup>Average of the normalized values when several probes corresponded to the same accession.
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Relationship between methylation and expression in sugar beet genotypes after 9 weeks of vernalization treatment for the 14 sequences that are both differentially methylated (fold change R versus S >2) and differentially expressed (fold change R versus S >2 and *P*<0.05) between resistant and sensitive genotypes. The x- and y-axis correspond to expression and methylation normalized signals, respectively. Open circles, squares, and triangles represent the three bolting-sensitive genotypes, filled circles, squares, and triangles correspond to the bolting-resistant genotypes.</figcaption></figure>

*Bolting tolerance in sugar beet genotypes is associated with the differential expression of an integrative gene network in the shoot apical meristem*

The transcriptomic analysis revealed 169 DEGs between R and S vernalized genotypes. These genes are involved in environment perception, hormone signalling, and flowering control. In *Lolium perenne*, 70 cold response genes exhibited variation in expression during vernalization (*Ciannamea et al.*, 2006). Phytohormones, particularly those involved in gibberellin metabolism, are known to be involved in the vernalization response in sugar beet (*Mutasa-Göttgens et al.*, 2009). The present results showed that the expression of genes playing a role in the transition from the vegetative to the floral stage is remodelled in the shoot apical meristem during cold exposure, and in different ways according to the bolting tolerance of genotypes. These gene associations and their expression levels could correspond to a floral determination state in S genotypes, for which a 9 week cold exposure is sufficient to induce bolting under optimal flowering conditions (photoperiod and temperature). The R genotypes could remain in a vegetative state via the control of other genes. Indeed, some flowering activators, such as *BvFVE*, *BvAGL24*, and *BvCOL1*, had a higher relative expression level in S genotypes, while the floral repressor *BvFL1* had a higher relative expression level in R genotypes, corresponding to a suitable context for the bolting initiation. In a previous study, R and S genotypes also exhibited different expression patterns for key regulators of the vernalization in sugar beet (*Trap-Gentil et al.*, 2011; *Hébrard et al.*, 2013). Recently, the *BvFL1* gene was functionally characterized using RNAi and overexpression in transgenic sugar beet plants (*Vogt et al.*, 2014). *BvFL1* RNAi neither eliminated the requirement for vernalization of biennial beets nor had a major effect on bolting time after vernalization, while overexpression of *BvFL1* resulted in a moderate late-bolting phenotype. Neither down-regulation of *BvFL1* by RNAi nor its overexpression affected the

expression of the central floral regulators *BvBTC1*, *BvFT1*, and *BvFT2* in a consistent way. These data demonstrated that *BvFL1* is not a major regulator of the vernalization response in sugar beet and its role has to be clarified. In the present study, *BvFT1* and *BvFT2* were not retained after the filtration of the transcriptomic probes, indicating no expression or a low expression level in the shoot apical meristem, in agreement with the findings of Pin *et al.* (2010). In the same way, *BvRNMTa* and *BvRNMTb* seemed not to be expressed after 9 weeks of vernalization, in agreement with previous data (Hébrard *et al.*, 2013). In 2012, Mutasa-Göttgens *et al.* generated a sugar beet transcriptome by RNA sequencing (RNA-seq) and compared the expression level of genes in different genotypes before and after 21 weeks of vernalization and with or without gibberellin treatment. Some genes highlighted in that study were also evidenced in the present study, such as *BQ584677* corresponding to the *locus 6819* in Mutasa-Göttgens *et al.* (2012). This EST has strong homology with the *FUL* gene of *Arabidopsis*, which is involved in floral meristem identity, and was expressed differentially between R and S sugar beet genotypes. This RNA-seq analysis also revealed 20% of genes that were differentially expressed between genotypes with distinct vernalization requirements. These results indicated that such sugar beet genotypes are both genetically and transcriptionally distinct, in agreement with the present data. Similar divergences have been observed between Nipponbare and indica rice varieties, which revealed 1525 DEGs (Chodavarapu *et al.*, 2012).

*Bolting tolerance in sugar beet is associated with the hypermethylation of an integrative gene network in the shoot apical meristem of sensitive genotypes*

The methylation analysis revealed 111 DMRs between R and S vernalized genotypes. Therefore, microarray is an efficient tool to identify DMRs at a genomic scale using available public genomic resources. In a previous study, 39 DMRs between
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Bolting/flowering model in sugar beet. *BvAGL24*, AGAMOUS-LIKE 24; *Bvbtc1*, bolting time control 1; *BvCOL1*, CONSTANS-LIKE 1; *BvFL1*, FLOWERING LOCUS C; *BvFT1* and *BvFT2*, FLOWERING LOCUS T 1 and 2; *BvFVE*, FLOWERING LOCUS VE; *BvFUL*, FRUITFUL; *BvRNMT*, RNA METHYLTRANSFERASE. Arrows and a dotted line indicated a hypothetical control based on the *Arabidopsis* flowering pathway. Genes coloured in blue are floral inducers overexpressed in bolting-sensitive genotypes. *BvFL1*, in red, is a floral inhibitor overexpressed in bolting-resistant genotypes. Asterisks (*) indicate hypermethylated genes in bolting-sensitive genotypes.</figcaption></figure>

genotypes were identified during vernalization by the use of
the restriction landmark genome scanning (RLGS) method
(Hébrard *et al.*, 2013). Here, the 111 DMRs all had a higher
relative methylation level in S genotypes compared with the
R genotypes. These results are in agreement with previous
global DNA methylation data revealing the hypermethyl-
ated status of bolting-sensitive sugar beet genotypes after 9
weeks of vernalization (Trap-Gentil *et al.*, 2011). These ele-
ments indicated that R and S genotypes have distinct DNA
methylation patterns during vernalization, allowing their
early discrimination. DNA methylation divergences in the
same species have already been reported in two *Arabidopsis*
ecotypes, Landsberg *erecta* and C24, and their reciprocal
hybrids (Shen *et al.*, 2012). In rice, the Nipponbare and indica

varieties revealed 7% of epimutations (Chodavarapu *et al.*, 2012). The opposite DNA methylation and expression patterns observed in sugar beet genotypes during vernalization may reflect two meristematic determination states, one corresponding to a vegetative form and the other to a floral form. Indeed, in *Azalea japonica*, floral buds have been shown to be hypermethylated compared with vegetative buds (Meijon *et al.*, 2009).

Fourteen sequences that were both differentially expressed
and differentially methylated (DEMSs) between R and S
genotypes were also identified. In the same way, 45 DEMSs
have already been reported between the Nipponbare and
indica rice varieties (Chodavarapu *et al.*, 2012). Expression
and DNA methylation of these 14 sugar beet DEMSs were
negatively correlated, in agreement with other studies in *Arabidopsis* (Zhang *et al.*, 2006; Lister *et al.*, 2008), poplar (Vining *et al.*, 2012), and rice (He *et al.*, 2010). The biological significance of gene-body methylation is still not clear, although it has been proposed that it may suppress aberrant transcription from cryptic promoters inside the genes or participate in splicing variants (Zilberman *et al.*, 2007; Chodavarapu *et al.*, 2010; Lauria and Rossi, 2011). Nine of these 14 DEMSs exhibited putative functional annotations. Two sequences, encoding *CYSTATIN 6* (CYSB, FG345694) and a potential *RNA POLYMERASE II TRANSCRIPTION MEDIATOR* subunit (BQ582428), have already been shown to be involved in plant development and response to environmental stress. Cystatins are a group of proteins that are widespread in eukaryotes and inhibit cysteine proteases by direct interaction with the active site (Turk and Bode, 1991; Margis *et al.*, 1998). In plants, it has been suggested that cystatins may play a role in storage proteins, programmed cell death, defence mechanisms, regulation of endogenous proteases, and catabolism, germination, and seed maturation (Massonneau *et al.*, 2005). In *Arabidopsis*, seven *CYS* genes (*AtCYS1*–*AtCYS7*) have been identified (Martinez *et al.* 2005). The expression of *AtCYSa* and *AtCYSb* (synonymous with *AtCYS3* and *AtCYS6*, respectively) is strongly induced in cells and seedlings by multiple abiotic stresses including high salinity, drought, oxidative stress, and cold; overexpression of these genes in transgenic yeast and *Arabidopsis* plants increases resistance to these stresses (Zhang *et al.*, 2008). Other studies have reported an up-regulation of cystatin mRNA transcripts in organs submitted to adverse growth conditions, including low temperatures in amaranth, barley, maize, or wheat (Valdés-Rodriguez *et al.*, 2007; Gaddour *et al.*, 2001; Massonneau *et al.*, 2005, Christova *et al.*, 2006, respectively).

The Mediator complex is a large multiprotein complex conserved in all eukaryotes which acts as a bridge between the RNA polymerase II complex and transcription factors (Kim *et al.*, 1994; Koleske and Young, 1994). Developmental aberrations have been reported for *Arabidopsis* mutants in a number of genes subsequently identified as encoding Mediator subunits. For example, the loss of *SENSITIVE-TO-FREEZING6* (SFR6/MED16) function renders mutants unable to cold-acclimate and disrupts transcriptional outputs without low-temperature gene regulation, affecting expression of the flowering time pathway and circadian clock genes (Knight *et al.*, 2008). The *struwwelpeter* (*swp/med14*) mutants exhibited reduced cell numbers in all aerial organs, changes in the window of cell proliferation, and perturbations in the shoot apical meristem (Autran *et al.*, 2002). Recently, Zheng *et al.* (2013) showed that *Arabidopsis* MEDIATOR SUBUNIT 18 (MED18) affects flowering time and floral organ formation through *FLC* and *AG*. In the *MED18* loss-of-function mutant, flowering time was delayed and floral organ number as well as *FLC* and *AG* mRNA levels were altered. FG345694 and BQ582428 constitute interesting candidate genes for bolting tolerance.

Regarding another seven annotated sequences (FG345845, FG343208, BQ582907, FG344953, BQ489455, BQ586261,

and BQ585496), no relationship to bolting has been established according to the literature, and further functional studies are needed to show their putative involvement in this physiological process. The possible role of these 14 DEMSs in bolting tolerance has yet to be defined. The study of corresponding *Arabidopsis* mutants, as in Hébrard *et al.* (2013), or RNAi sugar beet mutants (Pin *et al.*, 2010, 2012; Vogt *et al.*, 2014) will be relevant to elucidate their implication.

The bolting tolerance model (Fig. 6) is in agreement with phenotypic plasticity resulting from the control of plasticity genes modulated by environmental conditions (DNA methylation and gene expression) on the activity of organogenesis genes. The model suggests that sugar beet bolting tolerance corresponds to a regulation pathway whose genes are activated precociously in bolting-sensitive genotypes (through DNA methylation and gene expression variations) compared with the bolting-resistant genotypes. To test this model, methylome and transcriptomic analyses (using whole-genome bisulphite sequencing and RNA-seq, respectively) could be performed using the sugar beet genome that is now released (Dohm *et al.*, 2014). This study could be applied on a larger set of genotypes vernalized or not during short cold periods (from 1 to 3 weeks) in order to identify the early mechanisms controlling bolting tolerance.

*The sugar beet mitochondrial genome is affected by transcriptomic and methylation changes after vernalization related to bolting tolerance*

The sugar beet mitochondrial genome was analysed for gene expression and DNA methylation in both R and S genotypes after 9 weeks of vernalization. The six genotypes had the same haplotype (ON type). Only two probes corresponding to *ORF104* and *COX2* were differentially expressed between R and S genotypes. Surprisingly, 3211 methylation probes were retained after the filtration step, and a mitochondrial methylome was established for the two types of genotypes. The methylation normalized values obtained were negative and very low for both groups of genotypes, indicating a very low methylation level. Methylation has already been detected in the mitochondrial DNA of several organisms, such as human, mouse, hamster, and a number of plant species (Nass, 1973; Shmookler Reis and Goldstein, 1983; Pollack *et al.*, 1984; Simkova, 1998; Bellizzi *et al.*, 2013). Other studies suggested that DNA methyltransferases could have access to different sites on the mitochondrial DNA, depending on the levels of protein occupancy (Rebelo *et al.*, 2009) and could be translocated into the mitochondria (Shock *et al.*, 2011). The sugar beet mitochondrial genome had a higher relative methylation level in R genotypes compared with the S genotypes, in contrast to the nuclear DNA. In *Sequoia sempervirens*, mitochondrial methylation is higher in juvenile shoots compared with adult shoots (6.6% versus 7.8–8.2%, respectively) and vice versa for the nuclear methylation (8% versus 6.5–7.5%, respectively; Huang *et al.*, 2012). The DNA methylation machinery could be regulated in different ways in the mitochondrion and in the nucleus, and according to the developmental stages. In addition, 510 differentially methylated
probes have been identified between R and S genotypes, and some of them were clustered in two loci of 10–15 kb. In these loci, only three ORFs were expressed (*ORF152*, *ORF102b*, and *ORF192*), and their expression level was higher in S compared with the R genotypes. This result suggests that the increased methylation level in these loci could reduce the expression of the corresponding coding sequences in R genotypes. This relationship has been demonstrated for some *Arabidopsis* nuclear genes (Zhang *et al.*, 2006; Lister *et al.*, 2008). Further analyses are needed to understand the role of these two mitochondrial loci that could be good candidates to distinguish S and R genotypes by assisted marker selection.

## Conclusion

Flowering is regulated by seasonal cues to help ensure reproductive success, synchronous flowering in agriculture, and maximal seed yield. Vernalization is the best known example of an environmentally induced epigenetic switch (Song *et al.*, 2012). Recently, using epigenetic recombinant inbred lines (*epiRILs*), DMRs have been shown to act as epigenetic quantitative trait loci (*QTLepi*), accounting for 60% of the heritability for flowering time in *Arabidopsis* (Cortijo *et al.*, 2014). However, the genotypic tolerance to bolting and flowering, particularly important for breeders, is still poorly understood in many crops. In *Arabidopsis*, quantitative modulation of chromatin silencing mechanisms through *cis* polymorphisms explains natural variation in vernalization (Coustham *et al.*, 2012). The present results showed that bolting tolerance is associated with a genotype-dependent remodelling of DNA methylation and expression of an integrative gene network involved in environmental perception, hormone signalling, and flowering in the shoot apical meristem. One perspective will be to construct the clusters of each bolting-tolerant group by development and use of new SNPs discovered in the DEG/DMR genes. In addition, the shoot apical meristem is a particularly relevant tissue in which to study the role of epigenetics in plant phenotypic plasticity and adaptation, and seems to be the site of a specific epigenetic control (Baubec *et al.*, 2014). Taken together, the present data, in agreement with previously published studies on sugar beet (Pin *et al.*, 2010; Trap-Gentil *et al.*, 2011; Hébrard *et al.*, 2013), open up new perspectives for sugar beet breeding. Moreover, a recent study revealed for the first time a major QTL for bolting tolerance in sugar beet, representing a major contribution towards the development of a winter sugar beet with controlled bolting behaviour (Pfeiffer *et al.*, 2014). The complete genomic sequence of sugar beet has been recently published (Dohm *et al.*, 2014) and will support molecular breeding, especially by the development of genome-wide genetic and epigenetic markers for relevant agronomic traits.

## Supplementary data

Supplementary data are available at JXB online.

Figure S1. A Cot curve, a graphic representation of a genome.

Figure S2. Diagram representing the strategy used to identify candidate genes for bolting tolerance.

Figure S3. Histogram representing the read length distribution after 454 high throughput sequencing.

Figure S4. TAIR percentage of gene ontology (GO) terms for ‘Biological Process’ and ‘Molecular Function’ of annotated sugar beet ORFs from the Cot fraction.

Figure S5. Box plots showing the distribution of the 22 448 transcriptomic filtered probes after the normalization step for the six sugar beet genotypes.

Figure S6. Fold change expression using qRT-PCR validates transcriptomic microarray data.

Figure S7. Histogram showing the distribution of the 19 574 methylation-filtered probes after the normalization step for the six sugar beet genotypes.

Figure S8. Bisulphite sequencing confirmation of DNA methylation changes detected with the methylation microarray.

Figure S9. Relationship between methylation and expression in sugar beet genotypes after 9 weeks of vernalization treatment.

Table S1. Sequence of primers.

Table S2. Assembly data summary.

Table S3. Functional analysis of the Cot fraction.

Table S4. Contigs with a blast hit in the sugar beet genome.

Table S5. ORF annotations.

Table S6. ORFs with a blast hit in the sugar beet genome.

Table S7. 169 differentially expressed genes (DEGs).

Table S8. 111 differentially methylated regions (DMRs).

## Acknowledgements

A PhD grant (CHd) was provided by SES-VanderHave (Tienen, Belgium) and Association Nationale de la Recherche et de la Technologie (France), convention CIFRE no. 0379/2008). The authors thank Valérie Vidal and Véronique Vidal from IMAXIO/HELIXIO (Clermont-Ferrand, France) for technical assistance regarding the microarrays. We acknowledge Martine Decoville from Orléans University/CNRS (France) for the use of sonication material. We thank Clément Lafon-Placette, Claude Joseph, and Franck Brignolas from Orléans University (France) for their discussions and continuous help.

## References

Abe M, Kobayashi Y, Yamamoto S, Daimon Y, Yamaguchi A, Ikeda Y, Ichinoki H, Notaguchi M, Goto K, Araki T. 2005. FD, a bZIP protein mediating signals from the floral pathway integrator FT at the shoot apex. Science **309**, 1052–1056.

Abegg FA. 1936. A genetic factor for the annual habit in beets and linkage relationships. Journal of Agricultural Research **53**, 493–511.

Abou-Elwafa SF, Buttner B, Kopisch-Obuch FJ, Jung C, Muller AE. 2012. Genetic identification of a novel bolting locus in *Beta vulgaris* which promotes annuality independently of the bolting gene B. Molecular Breeding **29**, 9.

Adetunji I, Willems G, Tschoep H, Bürkholz A, Barnes S, Boer M, Malosetti M, Horemans S, van Eeuwijk F. 2014. Genetic diversity and linkage disequilibrium analysis in elite sugar beet breeding lines and wild beet accessions. Theoretical and Applied Genetics **127**, 559–571.

Andrés F, Coupland G. 2012. The genetic basis of flowering responses to seasonal cues. Nature Reviews Genetics **13**, 627–639.

Autran D, Jonak C, Belcram K, Beemster GT, Kronenberger J, Grandjean O, Inzé D, Traas J. 2002. Cell numbers and leaf development
in *Arabidopsis*: a functional analysis of the STRUWWELPETER gene.
EMBO Journal **21**, 6036–6049.

Baubec T, Finke A, Mittelsten Scheid O, Pecinka A. 2014. Meristem-specific expression of epigenetic regulators safeguards transposon silencing in Arabidopsis. EMBO Reports **15**, 446–452.

Baulcombe DC, Dean C. 2014. Epigenetic regulation in plant responses to the environment. Cold Spring Harbor Perspectives in Biology **6**, a019471.

Bellizzi D, D'Aquila P, Scafone T, Giordano M, Riso V, Riccio A, Passarino G. 2013. The control region of mitochondrial DNA shows an unusual CpG and non-CpG methylation pattern. DNA Research **20**, 537–547.

Bennett MD. 1972. Nuclear DNA content and minimum generation time in herbaceous plants. Proceedings of the Royal Society B: Biological Sciences **181**, 109–135.

Blattner FR, Plunkett G, Bloch CA, et al. 1997. The complete genome sequence of *Escherichia coli* K-12. Science **277**, 1453–1462.

Boudry P, Wieber R, Saumitou-Laprade P, Pillen K, Dijk H, Jung C. 1994. Identification of RFLP markers closely linked to the bolting gene B and their significance for the study of the annual habit in beets (*Beta vulgaris* L.). Theoretical and Applied Genetics **88**, 852–858.

Büttner B, Abou-Elwafa S, Zhang W, Jung C, Müller A. 2010. A survey of EMS-induced biennial *Beta vulgaris* mutants reveals a novel bolting locus which is unlinked to the bolting gene B. Theoretical and Applied Genetics **121**, 1117–1131.

Chia TYP, Müller A, Jung C, Mutasa-Göttgens ES. 2008. Sugar beet contains a large CONSTANS-LIKE gene family including a CO homologue that is independent of the early-bolting (B) gene locus. Journal of Experimental Botany **59**, 2735–2748.

Chodavarapu RK, Feng S, Bernatavichute YV, et al. 2010. Relationship between nucleosome positioning and DNA methylation. Nature **466**, 388–392.

Chodavarapu RK, Feng S, Ding B, Simon SA, Lopez D, Jia Y, Wang GL, Meyers BC, Jacobsen SE, Pellegrini M. 2012. Transcriptome and methylome interactions in rice hybrids. Proceedings of the National Academy of Sciences, USA **109**, 12040–12045.

Christova PK, Christov NK, Imai R. 2006. A cold inducible multidomain cystatin from winter wheat inhibits growth of the snow mold fungus, *Microdochium nivale*. Planta **223**, 1207–1218.

Ciannamea S, Busscher-Lange J, de Folter S, Angenent GC, Immink RG. 2006. Characterization of the vernalization response in *Lolium perenne* by a cDNA microarray approach. Plant and Cell Physiology **47**, 481–492.

Clarke JD, Zhu T. 2006. Microarray analysis of the transcriptome as a stepping stone towards understanding biological systems: practical considerations and perspectives. The Plant Journal **45**, 630–650.

Cortijo S, Wardenaar R, Colomé-Tatché M, et al. 2014. Mapping the epigenetic basis of complex traits. Science **343**, 1145–1148.

Coustham V, Li P, Strange A, Lister C, Song J, Dean C. 2012. Quantitative modulation of polycomb silencing underlies natural variation in vernalization. Science **337**, 584–587.

Dally N, Xiao K, Holtgräwe D, Jung C. 2014. The B2 flowering time locus of beet encodes a zinc finger transcription factor. Proceedings of the National Academy of Sciences, USA **111**, 10365–10370.

Dohm JC, Minoche AE, Holtgräwe D, et al. 2014. The genome of the recently domesticated crop plant sugar beet (*Beta vulgaris*). Nature **505**, 546–549.

El-Mezawy A, Dreyer F, Jacobs G, Jung C. 2002. High resolution mapping of the bolting gene B of sugar beet. Theoretical and Applied Genetics **105**, 100–105.

FAO. 2012. FAOSTAT. http://faostat.fao.org. Accessed 10 January 2014.

Flavell RB, Bennett MD, Smith JB, Smith DB. 1974. Genome size and the proportion of repeated nucleotide sequence DNA in plants. Biochemical Genetics **12**, 257–269.

Gaddour K, Vicente-Carbojosa J, Lara P, Isabel-Lamoneda I, Diaz I, Carbonero P. 2001. A constitutive cystatin-encoding gene from barley (Icy) responds differentially to abiotic stimuli. Plant Molecular Biology **45**, 599–608.

Galbraith DW. 2006. DNA microarray analyses in higher plants. OMICS **10**, 455–473.

Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, Robles M, Talón M, Dopazo J, Conesa. 2008. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Research **36**, 3420–3435.

He G, Zhu X, Elling AA, et al. 2010. Global epigenetic and transcriptional trends among two rice subspecies and their reciprocal hybrids. The Plant Cell **22**, 17–33.

Hébrard C, Trap-Gentil M-V, Lafon-Placette C, Delaunay A, Joseph C, Lefebvre M, Barnes S, Maury S. 2013. Identification of differentially methylated regions during vernalization revealed a role of RNA methyltransferase in bolting. Journal of Experimental Botany **64**, 651–663.

Henderson IR, Jacobsen SE. 2007. Epigenetic inheritance in plants. Nature **447**, 418–424.

Huang LC, Hsiao LJ, Pu SY, Kuo CI, Huang BL, Tseng TC, Huang HJ, Chen YT. 2012. DNA methylation and genome rearrangement characteristics of phase change in cultured shoots of *Sequoia sempervirens*. Physiologia Plantarum **145**, 360–368.

Jaggard KW, Werker AR. 1999. An evaluation of the potential benefits and costs of autumn-sown sugarbeet in NW Europe. Journal of Agricultural Science **132**, 91–102.

Kim YJ, Bjorklund S, Li Y, Sayre MH, Kornberg RD. 1994. A multiprotein mediator of transcriptional activation and its interaction with the C-terminal repeat domain of RNA polymerase II. Cell **77**, 599–608.

Koleske AJ, Young RA. 1994. An RNA polymerase II hoLOenzyme responsive to activators. Nature **368**, 466–469.

Knight H, Thomson AJ, Mc Watters HG. 2008. Sensitive to freezing6 integrates cellular and environmental inputs to the plant circadian clock. Plant Physiology **148**, 293–303.

Lamoureux D, Peterson DG, Li W, Fellers JP, Gill BS. 2005. The efficacy of Cot-based gene enrichment in wheat (*Triticum aestivum* L.). Genome **48**, 1120–1126.

Lang A. 1965. Physiology of flower initiation. In: Lang A, ed. Encyclopedia of plant physiology. Berlin: Springer Verlag, 1380–1536.

Lauria M, Rossi V. 2011. Epigenetic control of gene regulation in plants. Biochimica et Biophysica Acta **1809**, 369–378.

Lexander K. 1980. Present knowledge of sugar bolting mechanisms. 43rd Winter Congress. Bruxelles: Institut International de Recherches Betteravières, 245–258.

Lister R, O'Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, Millar AH, Ecker JR. 2008. Highly integrated single-base resolution maps of the epigenome in *Arabidopsis*. Cell **133**, 523–536.

Liu W, Thummasuwan S, Sehgal SK, Chouvarine P, Peterson DG. 2011. Characterization of the genome of bald cypress. BioMed Central Genomics **12**, 553.

Longden PC. 1989. Effects of increasing weed-beet density on sugar-beet yield and quality. Annals of Applied Biology **114**, 527–532.

Margis R, Reis EM, Villeret V. 1998. Structural and phylogenetic relationships among plant and animal cystatins. Archives of Biochemistry and Biophysics **359**, 24–30.

Margulies M, Egholm M, Altman WE, et al. 2005. Genome sequencing in microfabricated high-density picolitre reactors. Nature **437**, 376–380.

Martinez M, Abraham Z, Carbonero P, Diaz I. 2005. Comparative phylogenetic analysis of cystatin gene families from arabidopsis, rice and barley. Molecular Genetics and Genomics **273**, 423–432.

Martinez-Godoy MA, Mauri N, Juarez J, Marques MC, Santiago J, Forment J, Gadea J. 2008. A genome wide 20 K citrus microarray for gene expression analysis. BioMed Central Genomics **9**, 318.

Mascarell-Creus A, Canizares J, Vilarrasa-Blasi J, et al. 2009. An oligo-based microarray offers novel transcriptomic approaches for the analysis of pathogen resistance and fruit quality traits in melon (*Cucumis melo* L.). BioMed Central Genomics **10**, 467.

Massonneau A, Condamine P, Wisniewski JP, Zivy M, Rogowsky PM. 2005. Maize cystatins respond to developmental cues, cold stress and drought. Biochimica et Biophysica Acta **1729**, 186–199.

Maury S, Trap-Gentil MV, Hébrard C, Weyens G, Delaunay A, Barnes S, Lefebvre M, Joseph C. 2012. Genic DNA methylation changes during *in vitro* organogenesis: organ specificity and conservation between parental lines of epialleles. Physiologia Plantarum **146**, 321–335.
Meijon M, Valledor L, Santamaria E, Testillano PS, Risueno MC, Rodriguez R, Feito I, Canal MJ. 2009. Epigenetic characterization of the vegetative and floral stages of azalea buds: dynamics of DNA methylation and histone H4 acetylation. Journal of Plant Physiology **166**, 1624–1636.

Metzger JD. 1988. Localization of the site of perception of thermoinductive temperatures in *Thlaspi arvense* L. Plant Physiology **88**, 424–428.

Michaels S, Amasino R. 2001. Loss of FLOWERING LOCUS C activity eliminates the late flowering phenotype of FRIGIDA and autonomous pathway mutations but not responsiveness to vernalization. The Plant Cell **13**, 935–941.

Michaels SD, Amasino RM. 1999. *FLOWERING LOCUS C* encodes a novel MADS domain protein that acts as a repressor of flowering. The Plant Cell **11**, 949–956.

Mutasa-Göttgens E, Qi A, Mathews A, Thomas S, Phillips A, Hedden P. 2009. Modification of gibberellin signalling (metabolism & signal transduction) in sugar beet: analysis of potential targets for crop improvement. Transgenic Research **18**, 301–308.

Mutasa-Göttgens ES, Joshi A, Holmes HF, Hedden P, Gottgens B. 2012. A new RNAseq-based reference transcriptome for sugar beet and its application in transcriptome-scale analysis of vernalization and gibberellin responses. BioMed Central Genomics **13**, 99.

Nass MM. 1973. Differential methylation of mitochondrial and nuclear DNA in cultured mouse, hamster and virus-transformed hamster cells. *In vivo* and *in vitro* methylation. Journal of Molecular Biology **80**, 155–175.

Paterson AH. 2006. Leafing through the genomes of our major crop plants: strategies for capturing unique information. Nature Reviews Genetics **7**, 174–184.

Peterson DG. 2005. Reduced representation strategies and their application to plant genomes. In: Meksem K, Kahl G, eds. The handbook of plant genome mapping: genetic and physical mapping. Weinheim: WILEY-VCH Verlag GmbH & Co. KGaA, 307–335.

Peterson DG, Boehm KS, Stack SM. 1997. Isolation of milligram quantities of nuclear DNA from tomato (Lycopersicon esculentum), a plant containing high levels of polyphenolic compounds. Plant Molecular Biology Reporter **15**, 148–153.

Peterson DG, Schulze SR, Sciara EB, Lee SA, Bowers JE, Nagel A, Jiang N, Tibbitts DC, Wessler SR, Paterson AH. 2002. Integration of Cot analysis, DNA cloning, and high-throughput sequencing facilitates genome characterization and gene discovery. Genome Research **12**, 795–807.

Pfeiffer N, Tränkner C, Lemnian I, Grosse I, Müller AE, Jung C, Kopisch-Obuch FJ. 2014. Genetic analysis of bolting after winter in sugar beet (*Beta vulgaris* L.). Theoretical and Applied Genetics **127**, 2479–2489.

Pin PA, Benlloch R, Bonnet D, Wremerth-Weich E, Kraft T, Gielen JJL, Nilsson O. 2010. An antagonistic pair of FT homologs mediates the control of flowering time in sugar beet. Science **330**, 1397–1400.

Pin PA, Zhang W, Vogt SH, et al. 2012. The role of a pseudo-response regulator gene in life cycle adaptation and domestication of beet. Current Biology **22**, 1095–1101.

Pollack Y, Kasir J, Shemer R, Metzger S, Szyf M. 1984. Methylation pattern of mouse mitochondrial DNA. Nucleic Acids Research **12**, 4811–4824.

Rebelo AP, Williams SL, Moraes CT. 2009. *In vivo* methylation of mtDNA reveals the dynamics of protein mtDNA interactions. Nucleic Acids Research **37**, 6701–6715.

Reeves PA, He Y, Schmitz RJ, Amasino RM, Panella LW, Richards CM. 2007. Evolutionary conservation of the FLOWERING LOCUS C-mediated vernalization response: evidence from the sugar beet (*Beta vulgaris*). Genetics **176**, 295–307.

Saeed AI, Sharov V, White J, et al. 2003. TM4: a free, open-source system for microarray data management and analysis. Biotechniques **34**, 374–378.

Searle I, He Y, Turck F, Vincent C, Fornara F, Krober S, Amasino RA, Coupland G. 2006. The transcription factor FLC confers a flowering response to vernalization by repressing meristem competence and systemic signaling *in Arabidopsis*. *Genes and Development* **20**, 898–912.

Shen H, He H, Li J, Chen W, Wang X, Guo L, Peng Z, He G, Zhong S, Qi Y, Terzaghi W, Deng XW. 2012. Genome-wide analysis of DNA

methylation and gene expression changes in two *Arabidopsis* ecotypes and their reciprocal hybrids. The Plant Cell **24**, 875–892.

Shmookler Reis RJ, Goldstein S. 1983. Mitochondrial DNA in mortal and immortal human cells. Genome number, integrity, and methylation. Journal of Biological Chemistry **258**, 9078–9085.

Shock, LS, Thakkar PV, Peterson EJ, Moran RG, Taylor SM. 2011. DNA methyltransferase 1, cytosine methylation and cytosine hydroxymethylation in mammalian mitochondria. Proceedings of the National Academy of Sciences, USA **108**, 3630–3635.

Simkova H. 1998. Methylation of mitochondrial DNA in carrot (Daucus carota L.). Plant Cell Reports **17**, 220–224.

Smit A. 1983. Influence of external factors on growth and development of sugar beet (*Beta vulgaris* L.). Agricultural Research Reports, PUDOC.

Song J, Angel A, Howard M, Dean C. 2012. Vernalization—a cold-induced epigenetic switch. Journal of Cell Science **125**, 3723–3731.

Song J, Irwin J, Dean C. 2013. Remembering the prolonged cold of winter—minireview. Current Biology **23**, 807–811.

Trap-Gentil M-V, Hébrard C, Lafon-Placette C, Delaunay A, Hagège D, Joseph C, Brignolas F, Lefebvre M, Barnes S, Maury S. 2011. Time course and amplitude of DNA methylation in the shoot apical meristem are critical points for bolting induction in sugar beet and bolting tolerance between genotypes. Journal of Experimental Botany **62**, 2585–2597.

Turk V, Bode W. 1991. The cystatins: protein inhibitors of cysteine proteinases. FEBS Letters **285**, 213–219.

Valdés-Rodriguez S, Guerrero-Rangel A, Melgoza-Villagomez C, Chagolla-Lopez A, Delgado-Vargas F, Martinez-Gallardo N, Sanchez-Hernandez C, Délano-Frier H. 2007. Cloning of a cDNA encoding a cystatin from grain amaranth (*Amaranthus hypochondracus*) showing a tissue-specific expression that is modified by germination and abiotic stress. Plant Physiology and Biochemistry **45**, 790–798.

Van Dijk H. 2009. Evolutionary change in flowering phenology in the iteroparous herb *Beta vulgaris* ssp. *maritima*: a search for the underlying mechanisms. Journal of Experimental Botany **60**, 3143–3155.

Vining KJ, Pomraning KR, Wilhelm LJ, Priest HD, Pellegrini M, Mockler TC, Freitag M, Strauss S. 2012. Dynamic DNA cytosine methylation in the *Populus trichocarpa* genome: tissue-level variation and relationship to gene expression. BioMed Central Genomics **13**, 27.

Vogt SH, Weyens G, Lefèbvre M, Bork B, Schechert A, Müller AE. 2014. The FLC-like gene BvFL1 is not a major regulator of vernalization response in biennial beets. Frontiers in Plant Science **5**, 146.

Wardenaar R, Liu H, Colot V, Colomé-Tatché M, Johannes F. 2013. Evaluation of MeDIP-Chip in the context of whole-genome bisulfite sequencing (WGBS-Seq) in *Arabidopsis*. In: Lee TL, Chun Shui Luk A, eds. Tiling arrays: methods and protocols. New York: Springer Science+Business Media, 203–224.

Wellensiek SJ. 1964. Dividing cells as the prerequisite for vernalization. Plant Physiology **39**, 832–835.

Weng Yi, Huang THM, Yan PS. 2009. Methylated DNA immunoprecipitation and microarray-based analysis: detection of DNA methylation in breast cancer cell lines. Methods in Molecular Biology **590**, 165–176.

Yuan Y, San Miguel PJ, Bennetzen JL. 2003. High-Cot sequence analysis of the maize genome. The Plant Journal **34**, 249–255.

Zhang X, Liu S, Takano T. 2008. Two cysteine proteinase inhibitors from *Arabidopsis thaliana*, AtCYSa and AtCYSb, increasing the salt, drought, oxidation and cold tolerance. Plant Molecular Biology **68**, 131–143.

Zhang X, Yazaki J, Sundaresan A et al. 2006. Genome-wide high-resolution mapping and functional analysis of DNA methylation in *Arabidopsis*. Cell **126**, 1189–1201.

Zheng Z, Guan H, Leal F, Grey PH, Oppenheimer DG. 2013. Mediator subunit18 controls flowering time and floral organ identity in *Arabidopsis*. PLoS One **8**, e53924.

Zilberman D, Gehring M, Tran RK, Ballinger T, Henikoff S. 2007. Genome-wide analysis of *Arabidopsis thaliana* DNA methylation uncovers an interdependence between methylation and transcription. Nature Genetics **39**, 61–69.

Zimmerman JL, Goldberg RB. 1977. DNA sequence organization in the genome of *Nicotiana tabacum*. Chromosoma **59**, 227–252.