Global DNA methylation variations after short-term heat shock treatment in cultured microspores of Brassica napus cv. Topas

Heat stress can induce the cultured microspores into embryogenesis. In this study, whole genome bisulphite sequencing was employed to study global DNA methylation variations after short-term heat shock (STHS) treatments in cultured microspores of Brassica napus cv. Topas. Our results indicated that treatment on cultured Topas microspores at 32 °C for 6 h triggered DNA hypomethylation, particularly in the CG and CHG contexts. And the total number of T32 (Topas 32 °C for 6 h) vs. T0 (Topas 0 h) differentially methylated region-related genes (DRGs) was approximately two-fold higher than that of T18 (Topas 18 °C for 6 h) vs. T0 DRGs, which suggested that 32 °C might be a more intense external stimulus than 18 °C resulting in more changes in the DNA methylation status of cultured microspores. Additionally, 32 °C treatment for 6 h led to increased CHG differential methylations of transposons (DMTs), which were mainly constituted by overlaps between the hypomethylated differentially methylated regions (hypo-DMRs) and transposon elements (TEs). Further analysis demonstrated that the DRGs and their paralogs exhibited differential methylated/demethylated patterns. To summarize, the present study is the first methylome analysis of cultured microspores in response to STHS and may provide valuable information on the roles of DNA methylation in heat response.

involved in the control of cell growth during heat stress in tobacco BY-2 cells 2 . Although heat pre-treatment is one of the significant stresses for ME induction 7,15 , few studies have examined whether DNA methylation changes in cultured microspores after STHS treatment. Here, we tried to employ genome-wide bisulphite sequencing (GWBS) to decipher global DNA methylation variations after 32 °C and 18 °C treatments for 6 h in cultured microspores of B. napus cv. Topas at single-base resolution. And our results revealed that 32 °C heat treatment for 6 h was sufficient to induce global DNA hypomethylation in cultured Topas microspores. And 32 °C might be a more intense external stimulus than 18 °C generating more changes in the DNA methylation status of cultured microspores. To summarize, the present study is the first methylome analysis of cultured microspores in response to STHS and may provide valuable information on the roles of DNA methylation in heat response.

Results
Microspore collection and culture. A highly embryogenic cultivar Topas of B. napus was chosen for the analysis based on previous studies 8, 16 . We attempted to collect only late uninucleate microspores as the initial materials for in vitro heat treatment and culture by bud selection and mesh screening (Fig. 1A). The mean diameter of these isolated microspores was 19.58 ± 1.09 μ m (Fig. 1D). Then 32 °C and 18 °C treatments on the isolated microspores for 6 h were adopted in our experiments. Results showed that many enlarged microspores were formed after 6 h under 32 °C treatment rather than 18 °C treatment (Fig. 1B,C). The mean diameter and the frequency of the swollen microspores following the 32 °C treatment for 6 h were 26.22 ± 1.90 μ m and 57.50 ± 7.50%, respectively (Fig. 1D,E). However, the mean diameter of the unswollen microspores under this same heat treatment condition was 20.02 ± 1.62 μ m, which was almost identical to the size of the initial microspores (Fig. 1D). In order to decipher the global DNA methylation variations after treatments at 32 °C and 18 °C for 6 h in cultured microspores of B. napus cv. Topas, sample without treatment (Fig. 1A) was chosen as a control for DNA methylation comparisions.
The DNA methylation landscapes of baseline and STHS-treated microspores. Genomic DNA extracted from the three collected samples (T0, Topas 0 h; T18, Topas 18 °C for 6 h; T32, Topas 32 °C for 6 h) was  Table S1). The paired-end sequence files were subjected to multiple filtration steps and then aligned and deduplicated with Bismark (Supplementary Methods). The generated alignment reports indicated that the unique mapping efficiencies varied from 47.8% to 49.4% and the average proportions of methylation in three contexts (CG, CHG, and CHH) were 53.0%, 16.4%, and 3.5%, respectively (Supplementary Table S2). The level of symmetrical methylation was much higher compared with non-symmetrical methylation. Subsequent descriptive statistics by methylKit were shown in Supplementary Methods. Additionally, we depicted the identified methylated cytosines of the three samples on each chromosome ( Fig. 2A). And the amount of identified methylated cytosines of the three samples distributed on the C genome was obviously higher than that distributed on the A genome ( Fig. 2A).

Identification of differentially methylated regions (DMRs).
To identify regions of the genome subjected to differential methylation, we used methylKit to calculate DMRs in a pairwise fashion ( Fig. 2B and  Supplementary Data S1). Then we randomly selected two DMRs for pyrosequencing validation, and the results confirmed the methylation status detected by GWBS ( Supplementary Fig. S3). Deeper analysis indicated that more DMRs were identified in T32 vs. T0 than those in T18 vs. T0, which was mainly due to more CG and CHG  DMRs were induced by 32 °C culture (Fig. 3A). Besides, more total DMRs were observed in the C genome than in the A genome in all pairwise samples ( Fig. 3B and Supplementary Data S1), particularly in the T32 vs. T0 comparison (Fig. 3B). In T32 vs. T0 and T32 vs. T18 comparisons, there were more hypomethylated DMRs (hypo-DMRs) than hypermethylated DMRs (hyper-DMRs) ( Fig. 3C and Supplementary Data S1). And these T32 vs. T0 and T32 vs. T18 hypo-DMRs were mainly occupied by CG and CHG DMRs ( Fig. 3C and Supplementary Data S1). Moreover, 222 and 116 T32 vs. T0 hypo-DMRs were distributed on the C and A genomes, respectively (Fig. 3D). Analysis of the percentages of hyper-and hypo-methylated regions per chromosome indicated that the 32 °C treatment on the cultured Topas microspores for 6 h resulted in a significant proportion of hypomethylation in the symmetrical CG and CHG contexts (Fig. 4).

Differential methylation in transposons.
Cytosine methylation is chiefly targeted towards transposon element (TE) silencing 17 . Stress-induced transposon activation has been confirmed by molecular data in many different hosts 14,18,19 . Therefore, it was necessary to investigate the overlapping information between the identified DMRs and TEs. Although Chalhoub et al. 20 previously analysed the TEs identified in the released B. napus genome, their position information is not available in the public database. We firstly employed RepeatScout and RepeatMasker to de novo identify TEs in the whole B. napus genome. A total of 146,998 TEs (43,083 in the A genome and 103,915 in the C genome) were identified. Of them, 102,624 were retrotransposons and 44,374 were DNA transposons. Then, the differential methylations of transposons (DMTs) were searched based on the position information of TEs and DMRs (Supplementary Data S2). Further analysis showed that the amount of CHG DMTs was higher than the CG and CHH DMTs in T32 vs. T0 and T32 vs. T18 comparisons; the greater number of CHG DMTs was attributed to an increased overlap between the TEs and hypo-DMRs rather than hyper-DMRs (Fig. 5A, Supplementary Data S2, and Supplementary Table S4). Additionally, the differential methylations of retrotransposons (DMRTs) in all contexts in the three pairwise samples were more than the differential methylations of DNA transposons (DMDTs); LINE, LTR/Copia, and LTR/Gypsy were the most abundant retrotransposable elements ( Fig. 5B and Supplementary Table S5). Subsequently, the distributions of DMTs on each chromosome were analysed. Results demonstrated that the numbers of DMTs located on A and C genomes were unequal (Supplementary Table S6).

Identification of DMR-related genes (DRGs).
We employed the methylKit package and an in-house R script to identify DRGs (Supplementary Data S3). Then the distribution of DRGs on each chromosome of B. napus was investigated (excluding the DRGs located on unassembled scaffolds) (Supplementary Table S7). The total number of DRGs in T32 vs. T0 (96) was approximately two-fold higher than T18 vs. T0 (52), mainly because of the increase in CG and CHG DRGs in T32 vs. T0 ( Fig. 6 and Supplementary Table S8). T32 vs. T0 DRGs could also be divided into 69 hypomethylated DRGs (hypo-DRGs) and 27 hypermethylated DRGs (hyper-DRGs), respectively. However, these two types of DRGs in T18 vs. T0 were almost equal (Table 1). Further analysis showed that only four common CG DRGs were identified between T32 vs. T0 and T18 vs. T0 ( Fig. 6 and Table 2). They exhibited similar tendencies of methylation change in both comparisons (Table 2). Functional annotation analysis indicated that BnaA03g36810D was similar to AT3G22840, which encodes an early light-inducible protein (ELIP) ( Table 2) and transiently accumulates in response to environmental stress 21 . Moreover, the mortality rates of plants lacking ELIPs are sometimes higher 22 . The in vitro culture itself was an environmental stress independent  of the different temperatures. Did the culture rather than temperature effects induce the similar methylation changes of the common CG DRGs identified between T32 vs. T0 and T18 vs. T0? To elucidate the differences of cultured microspores under different temperatures, subsequent analyses were mainly focused on the T32 vs. T18 DRGs ( Table 2, Table 3 and Supplementary Data S3). Among the total 77 T32 vs. T18 DRGs, 47 were hypo-DRGs and 30 were hyper-DRGs (Table 1). In addition, five common CG and two common CHH DRGs were sought between T32 vs. T18 and T18 vs. T0, respectively ( Fig. 6 and Table 2). Amazingly, the methylation/demethylation tendencies of these DRGs were totally opposite in two comparisons,  and only BnaA03g24920D and BnaC05g07550D were hypo-DRGs in T32 vs. T18 (Table 2). BnaA03g24920D was similar to proton gradient regulation 5-like 1B (PGRL1B). And PGRL1 has been proved to be involved in cyclic electron flow (CEF), which only generated ATP and was driven by photosystems I (PSI) 23 . Arabidopsis PSI CEF is abolished following thermal-stress 24 . We next asked whether PGRL1 maintained normal CEF functions during thermal-stress and generated sufficient energy for the survival of cultured microspores. BnaC05g07550D might encode LIM protein that regulating transcription or organizing the cytoskeleton by triggering the formation of actin bundles 25 . Actually, heat shock has been shown to cause changes in microtubule and cytoskeleton in cultured Topas microspores 26 . All the seven common CG and five common CHG DRGs further identified between T32 vs. T18 and T32 vs. T0 had negative meth.diff values in both comparisons ( Fig. 6 and Table 2). Among these DRGs, BnaA03g39290D was similar to UBC29 (ubiquitin-conjugating enzyme 29). UBCs participate in protein degradation via proteasome and may be involved in various biological processes 27 . As for the STHS treatment on cultured microspore can disrupt pollen differentiation 11 , we asked whether BnaA03g39290D functioned in sweeping the proteins that required for pollen differentiation pathway. Another common DRG BnaA04g24700D resembled PDF1 (PROTODERMAL FACTOR 1), which was involved in the fate determination of epidermal cell 28 . The CHG DRG BnaA06g03030D resembled a heat stress-responsive gene FUT12 (fucosyltransferase 12) 29 . In addition, BnaC01g24140D was similar to AtDOK1 (A. thaliana dolichol kinase 1), the expression of which could complement the temperature-sensitive growth and glycosylation defects of the Saccharomyces cerevisiaesec59 mutant 30 . Moreover, AtDOK1 is involved in the synthesis of dolichol phosphate (Dol-P), which serves as a carrier of complex polysaccharides during protein glycosylation 30 . Glycoproteins are confirmed to be involved in adaptation to biotic and abiotic stresses 31 .
Except for the common DRGs, there were 58 specific DRGs (33 hypo-DRGs and 25 hyper-DRGs) existing in T32 vs. T18, including 29 CG, 21 CHG, and eight CHH DRGs (Table 3). These specific DRGs contained several genes that might be involved in various biological processes (Table 3). For example, BnaA08g05750D (thioredoxin superfamily protein) and BnaC02g29760D (glutathione S-transferase family protein) might function in maintaining cell redox homeostasis (Table 3). High temperature provokes the accumulation of reactive oxygen species (ROS) in plants 32 . Heat stress-incuded ME also generates an oxidative burst and ROS 33 . Therefore, sustaining cell redox homeostasis was crucial for the survival of cultured microspore during heat stress. Additionally, two putative glycosyltransferases (BnaC06g12760D and BnaC06g04480D) and a putative fasciclin-like arabinogalactanprotein 2 (FLA2) (BnaA06g19130D) were found as specific hypo-DRGs in T32 vs. T18 (Table 3). Previous study has indicated that the protein encoded by AT5G39990 (similar to BnaC06g12760D) was involved in the biosynthesis of type II arabinogalactan (AG) and cell elongation during seedling growth 34 . Most AGPs are O-glycosylated at hydroxyproline residues by type II AG group 35 . And their pivotal roles in cell wall signal transduction, plant development and stress tolerance have also been discussed 36 . Intriguingly, FLA2 was found to be significant in several stress-related AFGC microarray experiments 37 . It could be part of the auxin transporters ABCB19/PINFORMED1 (PIN1) nanodomain 38 . In addition, BnaC03g42260D resembled AtSGP2 (A. thaliana G-protein) ( Table 3). Arabidopsis monomeric G-proteins are implied to be markers of early and late events in cell differentiation 39 . Intriguingly, a specific CHG hypo-DRG BnaC01g19320D coincidently matched with napin gene SESA4 (seed storage albumin 4) ( Table 3), which was identified as an early molecular marker for ME in B. napus 16,40 . Another specific CHG hypo-DRG BnaC01g21110D identified from T32 vs. T18 might be a putative SCL (scarecrow-like) gene (Table 3). Joosen et al. 9 previously identified a B. napus gene, which was homologous to Arabidopsis SCL11, as a robust marker for ME. Whether these STHS-induced DRGs are really related to ME is still an open question. Besides, other gene groups related to energy metabolism, protein degradation, transcription and translation, and signal transduction (Table 3) were also included in specific T32 vs. T18 DRGs. Schranz et al. 41 have built the A to X conserved collinear blocks (CCBs) of the reduced karyotype (n = 5) of A. thaliana, and these blocks were defined by their position in a proposed ancestral karyotype (n = 8). The genomes of Brassica rapa (AA) and its sister species Brassica oleracea (CC) were almost complete triplications of the genome of A. thaliana 42,43 . B. napus (AACC) was formed by recent allopolyploidy between the ancestors of B. oleracea and B. rapa 20,44 . Therefore, the CCBs of B. napus were constructed based on the position information of the A. thaliana A to X segments 41 . Then we investigated the number of DRGs in each of the A to X CCBs of B. napus. Results indicated that the U and F blocks possessed the most DRGs in the T32 vs. T0 and T18 vs. T0 comparisons, whereas no DRGs were distributed on the G, K, and S blocks (Supplementary Table 9). Due to the polyploid properties, it was also necessary to evaluate whether differential methylated An-Cn paralog gene pairs existed in B. napus. The An-Cn paralog gene pairs of B. napus were searched according to the methods described by Liu et al. 43 . Although almost all of the T32 vs. T0 and T18 vs. T0 DRGs possessed one or more paralog genes, none of these paralog genes were included in the identified DRGs ( Fig. 7 and Supplementary Data 4). This result suggested that distinct DNA methylation regulatory pathways might exist even for paralogs. Chang and Liao 45 demonstrated that DNA methylation could "rebalance" the overall expression dosage of paralogs.

Real-time PCR to analyse DRG expression.
To examine the relationships between DNA methylation and gene expression levels, 16 DRGs (13 hypo-DRGs and three hyper-DRGs) were randomly selected for real-time PCR verification based on the UniGene information (Supplementary Data S5). Results indicated that seven genes were up regulated, two genes were down regulated, and the expression levels of seven genes remained unchanged (Fig. 8). DNA demethylation and methylation are considered to be associated with elevating and suppressing gene expressions, respectively. However, the expressions of four hypo-DRGs and two hyper-DRGs did not follow this principle ( Fig. 8 and Table 4). More dynamic and complex relationships between DNA methylation and expression have been illustrated in other studies 46,47 . Among these 16 selected genes, BnaA09g40710D was a T32 vs. T0 CG hypo-DRG and its expression was up-regulated by nearly two-fold in T32 compared with that of T0 (Fig. 8). Furthermore BnaA09g40710D was identical to Arabidopsis ELF3 (early flowering 3) (Table 4), which controlled elongation growth in response to temperature 48 . Similarly, another hypo-DRG BnaC05g29060D was up-regulated in the 32 °C treatment (Fig. 8). This gene resembled Arabidopsis SCAMP5 (secretory carrier membrane protein 5) (Table 4). Whether BnaC05g29060D functions in vesicles transportation and protein trafficking during heat treatment remains unknown.

Discussion
Exposure of Arabidopsis plants to heat stress results in an increased global methylation 49 . Nevertheless, in cotton (Gossypium hirsutum) anthers, high temperature leads to the genome-wide hypomethylation at the tetrad stage and the tapetal degradation stage 50 . Our results showed that the T32 vs. T0 and T32 vs. T18 DMRs were mainly comprised of hypo-DMRs rather than hyper-DMRs ( Fig. 3C and Supplementary Data S1). The percentages of the hyper-and hypo-methylated regions per chromosome clearly indicated that the 32 °C treatment on cultured microspores generated a significant proportion of hypomethylation in the symmetrical CG and CHG contexts (Fig. 4). This result demonstrated that STHS initiated the demethylation process, giving rise to a decrease in global DNA methylation in cultured microspores. Actually, the percentage of methylated cytosine increased from 3 to about 11% from microspore to mature stage in pollen of B. napus 51 . Therefore, it seems that there is no consistent trend in the changes of DNA methylation under heat in different species or in different cell types. Whether plant male germ cell tended to be global hypomethylated rather than hypermethylated after heat stress? If so, why the DNA methylation tendencies were different? Additionally, we also found that the total number of DRGs in T32 vs. T0 was approximately two-fold higher than T18 vs. T0 suggesting 32 °C may be a more intense external stimulus than 18 °C and ultimately resulted in more changes in the DNA methylation status of cultured microspores (Table 1 and Fig. 6).
In plants, DNA methylation occurs frequently in all three sequence contexts: the symmetric CG and CHG contexts and the asymmetric CHH context. Each type of DNA methylation is vital for development and responses to environmental stresses 52 . Here, our data revealed that the 32 °C treatment brought about increased CHG DMTs in the T32 vs. T0 and T32 vs. T18 comparisons due to abundant overlaps between the TEs and hypo-DMRs (Fig. 5A). Moreover, DMRTs in all contexts in the three pairwise samples were higher than DMDTs. LINE, LTR/ Copia, and LTR/Gypsy were the most abundant retrotransposable elements (Fig. 5B and   Yang et al. 47 previously observed that retrotransposons were more tightly controlled by methylation than DNA transposons during the floral development of Arabidopsis. ONSEN, an LTR/copia type retrotransposon was found activated by heat stress in Arabidopsis 18,19 . Deeper studies are needed to answer whether and how these DMRTs function during STHS treatment in cultured microspores. Intriguingly, some molecular markers, such as napin, G-protein, SCL and AGP, for early ME in B. napus identified via transcriptome and proteome analysis were also existed among the 58 specific T32 vs. T18 DRGs (Table 3). Additional experiments are required to decipher how DNA methylation functions as an epigenetic response to heat stress in cultured microspores, whether DNA methylation status is associated with the expression of the overlapped markers for early ME in B. napus and whether the micropsores from different cultivars that with different embryogenic potentials possess different epigenetic response to heat stress. Analysis of the overlaps between TEs and DMRs. First, RepeatScout (http://bix.ucsd.edu/repeatscout/) was used to construct a repeat library of the B. napus genome using an ab initio approach. Then this library was used to screen DNA sequences for interspersed repeats and low complexity DNA sequences using RepeatMasker (http://www.repeatmasker.org/). Retrotransposon and DNA transposons position information was retrieved from the out file generated by RepeatMasker. An in-house Perl script was used to calculate the overlaps between TEs and DMRs.