Multigenerational Exposure to Heat Stress Induces Phenotypic Resilience, and Genetic and Epigenetic Variations in Arabidopsis thaliana Offspring

Plants are sedentary organisms that constantly sense changes in their environment and react to various environmental cues. On a short-time scale, plants respond through alterations in their physiology, and on a long-time scale, plants alter their development and pass on the memory of stress to the progeny. The latter is controlled genetically and epigenetically and allows the progeny to be primed for future stress encounters, thus increasing the likelihood of survival. The current study intended to explore the effects of multigenerational heat stress in Arabidopsis thaliana. Twenty-five generations of Arabidopsis thaliana were propagated in the presence of heat stress. The multigenerational stressed lineage F25H exhibited a higher tolerance to heat stress and elevated frequency of homologous recombination, as compared to the parallel control progeny F25C. A comparison of genomic sequences revealed that the F25H lineage had a three-fold higher number of mutations [single nucleotide polymorphisms (SNPs) and insertions and deletions (INDELs)] as compared control lineages, suggesting that heat stress induced genetic variations in the heat-stressed progeny. The F25H stressed progeny showed a 7-fold higher number of non-synonymous mutations than the F25C line. Methylome analysis revealed that the F25H stressed progeny showed a lower global methylation level in the CHH context than the control progeny. The F25H and F25C lineages were different from the parental control lineage F2C by 66,491 and 80,464 differentially methylated positions (DMPs), respectively. F25H stressed progeny displayed higher frequency of methylation changes in the gene body and lower in the body of transposable elements (TEs). Gene Ontology analysis revealed that CG-DMRs were enriched in processes such as response to abiotic and biotic stimulus, cell organizations and biogenesis, and DNA or RNA metabolism. Hierarchical clustering of these epimutations separated the heat stressed and control parental progenies into distinct groups which revealed the non-random nature of epimutations. We observed an overall higher number of epigenetic variations than genetic variations in all comparison groups, indicating that epigenetic variations are more prevalent than genetic variations. The largest difference in epigenetic and genetic variations was observed between control plants comparison (F25C vs. F2C), which clearly indicated that the spontaneous nature of epigenetic variations and heat-inducible nature of genetic variations. Overall, our study showed that progenies derived from multigenerational heat stress displayed a notable adaption in context of phenotypic, genotypic and epigenotypic resilience.

Introduction

Plants are continuously exposed to their surrounding environment which effects their growth and impacts agricultural yields (Wei et al., 2020). These interactions result in changes in gene expression and physiological and biomolecular responses that alter a plant’s phenotype (Zheng et al., 2013; Wei et al., 2020).

Plants respond to the environment through the adaptation and tolerance that ensure the survival of both, the plant itself and its progeny (Huang et al., 2012). This phenomenon is feasible due to the plasticity of the plant genome and epigenome, and it evolves through the environmentally induced alterations and the development of phenotypic resilience (Nicotra et al., 2010; Lind and Spagopoulou, 2018; Zhang et al., 2018; Miryeganeh and Saze, 2020). The phenotypic plasticity in response to environmental conditions is coordinated through the stress perception and activation of signaling pathways and responsive genes (Chinnusamy et al., 2004; Zheng et al., 2013). The adaptation to the environment is not restricted to changes in the physiology of the plant, but it also results in changes in the genome and epigenome (Gutzat and Scheid, 2012; Zhang et al., 2018; Miryeganeh and Saze, 2020).

Environmental stresses have an impact on the directly exposed generation as well as on their offspring via parental effects and/or transgenerational effects (Herman and Sultan, 2011). Parental effects, which are also referred as intergenerational effects, are described by phenotypic changes in the immediate progeny of stressed plants. In contrast, transgenerational effects refer exclusively to changes in the phenotype of the progeny that is separated from the stressed parent by at least one generation, thus representing true inheritance of acquired changes. However, many of the described intergenerational effects share their mechanisms with transgenerational effects, and many publications do not distinguish between the two, often referring to both as transgenerational (Perez and Lehner, 2019). The RNA interference-related mechanisms can regulate the transgenerational inheritance of a specific chromatin or DNA modifications (Duempelmann et al., 2020). Very recently, Yang et al. (2020) have reported that the segregation of the MSH1 RNAi transgene produces the heritable non-genetic memory in association with methylome reprogramming. The msh1 reprogramming is dependent on functional HDA6 and DNA-methyltransferase MET1, and the transition to memory requires the RdDM pathway. This system of phenotypic plasticity may play a vital role for plant adaptation to the changing environment (Yang et al., 2020).

While the exact mechanisms of stress adaption and transgenerational stress memory remain unclear; the understanding of these processes has far-reaching implications for plant breeding, genetic engineering, and the development of stress tolerant crops (Bilichak and Kovalchuk, 2016; Zhang et al., 2018). A substantial number of studies demonstrate the ability to maintain the memory of stress exposure throughout ontogenesis and transmit this memory to the following generation (Boyko and Kovalchuk, 2010; Suter and Widmer, 2013a; Bilichak and Kovalchuk, 2016; Ramírez-Carrasco et al., 2017; Zheng et al., 2017).

The mechanism of stress memory most likely includes genetic and epigenetic changes, with the latter ones being more prevalent. The genetic variations within species are identified as single nucleotide polymorphisms (SNPs) and insertions and deletions (INDELs). It has been reported that stress induced the accumulation of genetic mutations in stressed lineages of bacteria (Al Mamun et al., 2012), yeast (Shor et al., 2013), and Arabidopsis (Jiang et al., 2014; Belfield et al., 2020). SNPs may cause the phenotypic diversity such as flowering time adaptation, the color of plant fruit and tolerance to various biotic and abiotic factors (Vidal et al., 2012; Jang et al., 2015).

In the case of epigenetics, phenotype traits are transferred to the progeny without changes in gene sequence (Eichten et al., 2014; Bilichak and Kovalchuk, 2016; Zheng et al., 2017; Lind and Spagopoulou, 2018; Miryeganeh and Saze, 2020). Epigenetic regulation consists of three primary mechanisms – DNA methylation, histone modifications and small RNAs (sRNAs) activity. DNA methylation is the most established process that alters gene expression utilized by plants for changing their phenotypes in stressful environment, acting as a complementary mode of transferring heritable information (Mirouze and Paszkowski, 2011; Schmitz et al., 2011). Many recent reports have demonstrated that plant DNA methylation can be altered at individual loci or across the entire genome under environmental stress conditions (Zhang et al., 2018 and reference therein). In rice, the naturally occurring variations in DNA methylation are well studied epimutations at single gene loci that have been shown to result in heritable morphological variations without altering the DNA sequence of rice genes (Miura et al., 2009; Zheng et al., 2017). The chilling-induced loss of tomato flavor is associated with changes in DNA methylation, including the RIN gene (Zhang et al., 2016). In Arabidopsis, quantitative resistance to clubroot infection is mediated by transgenerational epigenetic variations (Liégard et al., 2019).

Changes in DNA methylation include differentially methylated positions (DMPs), where single cytosines are involved, and differentially methylated regions (DMRs), where multiple cytosines in a given region (arbitrarily defined) are considered. DMPs and DMRs are specific epimutations that can be associated with the epigenetic inheritance of stress response (Zhang et al., 2018). Exposure to environmental factors can induce changes in the specific DMRs, which suggest the existence of the exposure specific DMRs (Manikkam et al., 2012; Haque et al., 2016). Some examples include DMRs associated with the differentially expressed genes in the regulation of the Arabidopsis thaliana immune system against Pseudomonas syringae (Dowen et al., 2012) and the location-specific DMRs at specific cis-regulatory sites in the Hevea brasiliensis tree crop in response to cold and drought stresses (Uthup et al., 2011).

Previous research from our laboratory demonstrated that the progeny of plants exposed to various abiotic and biotic stressors exhibited changes in phenotypic traits such as leaf size and number, flowering time, and seed size (Boyko et al., 2010; Bilichak et al., 2012; Migicovsky et al., 2014). The progeny of stressed plants also displayed a certain degree of stress tolerance and cross-tolerance lasting for one or two generations without the maintenance of stress factor (Rahavi and Kovalchuk, 2013; Migicovsky et al., 2014). Our previous studies also demonstrated changes in the genome – rearrangements at the resistance genes loci in the progeny of plants infected with tobacco mosaic virus (TMV) (Boyko et al., 2007) and in the epigenome – changes in DNA methylation and histone modifications in the progeny of plants exposed to various stresses such as salt (Bilichak et al., 2012), heat (Migicovsky et al., 2014), cold (Migicovsky and Kovalchuk, 2013), UVC (Migicovsky and Kovalchuk, 2014), and viral infection (Boyko et al., 2007). Although these studies did not explicitly analyze the mechanism of transgenerational inheritance, they proposed the role of epigenetic factors and demonstrated genetic changes relevant to stress exposure (Boyko and Kovalchuk, 2010; Boyko et al., 2010; Bilichak et al., 2015).

Based on studies conducted in our lab and several other labs, we hypothesized that multigenerational exposure to heat stress would be better in conveying stress memory. We also hypothesized that the plant population propagated under multigenerational heat stress would show a broader genotypic and epigenotypic diversity compared with the population of plants propagated without stress.

In this study, we have shown the ability of Arabidopsis plants to form heat stress memory using the repeated heat stress exposure over 25 consecutive generations. The F25H generation of stressed progeny was compared with the parental F2C and parallel F25C control progenies for heat stress resilience and the presence of heat stress-induced genetic and epigenetic changes.

Materials and Methods

Parental Generation and Progeny Plants

Plant seeds used for the parental generation (P) were obtained from single homozygous Arabidopsis thaliana (Columbia ecotype, Col-0) plant transgenic for the luciferase (LUC) recombination reporter gene (15d8 line) carrying a copy of a direct repeat of the luciferase recombinant transgene. All parental plants were grown to flower, and seeds were collected and pooled from approximately 16–20 plants per each experimental line. The subsequent progeny plants were obtained from these seeds; specifically, ∼100 seeds were sown, and ∼20 randomly selected plants were either exposed to heat or grown under normal conditions, both allowed to set seeds, thus forming the next generation. The same process was repeated consecutively until 25 generations representing the progeny plant lines (Figure 1A). The seeds were kept in storage under dry conditions at room temperature.

FIGURE 1

Figure 1. The heat-stressed progeny exhibited an increase in heat stress resilience and higher frequency of recombination events. (A) Schematic diagram for propagations of 25 generations of heat stressed and parallel control lineages. Single Arabidopsis thaliana plant (Col-0, 15d8 line) was used for seed origin. Plants were grown for 25 generations either exposed to heat stress (42°C) or in a control ambient environment (22°C). The ‘H’ group included the heat-stressed progenies, while the ‘C’ group included plants grown under control conditions acting as a parallel control. Seeds of second-generation F2C were used as parental control because seeds of F0 and F1C were no longer viable. Lineages used in this study are shown in the red box. (B) Trends in the germination of the F25H progeny under heat stress. Representative photographs of germination under control (0 h) and heat-stress (6 h) conditions (the left panel). Seeds were either from the parallel control (F25C) or the heat-stressed progeny (F25H). Each plate represents one set, and in total, five sets were used for each treatment. The right panel: The germination index graph. 0 h: non-stress control (0 h at 45°C); 6 h: heat-stress (6 h at 45°C). Heat-stress was applied at 45°C for 6 h, just 4-days after stratification at 4°C. The germination index was measured on the 6th day. The data were analyzed by two-way ANOVA with the Tukey’s multiple comparisons test using GraphPad Prism version 8.4.2 for Windows. The data are shown as an average mean ± SD. Mean values that were significantly different from each other are indicated by * (p < 0.05), ** (p < 0.01) and those that were non-significant are indicated by ns (p > 0.05). (C) Left panel shows survival of mature plants of stressed progeny (F25H vs. F25C) under heat-stress. A survival rate was calculated after 10-days of recovery and shown as an average percentage (with ± SD) of plants surviving the stress from three independent experiments. Right panel shows the average fresh weight (FW) per plant (with ± SE) of survived plants. The data were analyzed by one-way ANOVA with Tukey’s multiple comparisons test using GraphPad Prism version 8.4.2 for Windows. The asterisk shows a significant difference between the stressed progeny and parallel control progeny (***p < 0.001). (D) Somatic HRF under non-stress conditions. Progeny of stressed plants exhibited a higher somatic homologous recombination (HRF) without exposure to heat-stress. Left panel shows the scheme of the HRF analysis. Right panel shows the recombination events in F25C and F25H plants. The upper panel is a superimposed photograph of HR events and plants, and the lower panel is an original florescence photograph. Left panel also shows the average no. of HR spots per plant (with ± SE). The data were analyzed by one-way ANOVA with Tukey’s multiple comparisons test using GraphPad Prism version 8.4.2 for Windows. The asterisk shows a significant difference between the stressed progeny and parallel control progeny (*** p < 0.001).

The progeny plants obtained from the progeny seeds through the process described above for every generation had two groups: the stressed (H) and non-stressed groups (C) which acted as the parallel control group within each generation. All progeny plants representing the ‘H’ group were subjected to 42°C heat stress for 4 h per day for four consecutive days, while the trays with plants were covered with plastic domes to prevent loss of moisture, starting at day 10 post germination, and they were grown to flower and set seeds. Seeds were then collected and propagated to obtain the next successive generation, thus creating 25 generations of the heat stressed ‘H’ group. On the other hand, the progeny plants representing the ‘C’ group were grown under normal conditions and set seeds; this was also done for 25 generations.

Plant Growth Conditions

The seeds were kept at 4°C for 7 days to initiate stratification, and then sown in all-purpose potting soil prepared with water containing a generic fertilizer (Scotts Miracle-Gro, Marysville, OH, United States). The moist soil containing seeds was further stratified for 48 h at 4°C. Plants were grown in a growth chamber (BioChambers, Winnipeg, MB, Canada) at 22°C under the extended day conditions of 16 h light (120 μmol photons m^–2 s^–1) and 8 h dark. Approximately 3–5 days post germination (dpg), plant seedlings were transplanted into individual pots containing a 9:1 ratio of soil to vermiculite composition prepared with a fertilizer. Pots containing plants were placed in trays and watered from below.

Lineages Used for Heat Stress Phenotyping

The seeds of F0 (parental) and F1 generations were no longer viable, hence F2C seeds were used as parental control. Overall, the lineages F2C (parental control), F25H (the 25th generation of stressed progeny) and F25C (the 25th generation of parallel control progeny) were used for this study.

Phenotypic Analysis Under Heat Stress Treatment

The F25H stressed progenies were studied for heat stress memory and response. The response to heat stress was determined at the germination or mature plant stage. For germination experiment, three replicates were used for the stressed progeny, and two replicates were used for the parallel control progeny in one MS media plate. Each plate represented one set, and in total, five sets were used for each treatment. Heat stress was applied at 45°C for 6 h, just 4 days after stratification at 4°C. The percentage of germination and the speed of germination were qualitatively and quantitatively evaluated between the groups. To analyze the heat stress tolerance at the seedling stage, we exposed F25C and F25H plants to 42°C for 4 h for four consecutive days (keeping plants under plastic domes – in the same way as they were propagated under stress for 25 generations). To test the heat stress phenotype at the mature plant stage, the plants were grown under normal condition, and heat stress was applied just before bolting for 96 h at 40°C. The qualitative and quantitative assessment of survived plants has been done after 10-days of recovery from stress. Phenotype changes such as the survival rate and fresh weight of the survived plants were measured.

The Analysis of Homologous Recombination Frequency

The Arabidopsis thaliana #15d8 line transgenic for the luciferase (LUC) recombination reporter gene enables the analysis of Homologous Recombination Frequency (HRF). Recombination events occur via the rearrangement between two homologous inactive copies of the luciferase transgene. The scheme of the HRF analysis is shown in Figure 1D (left panel). The F25H and F25C plants were analyzed for HRF under non-stress conditions. HR events are analyzed by scoring bright sectors on a dark background with a CCD camera after spraying with luciferin (Ilnytskyy et al., 2004). These sectors represent cells in which recombination events occurred. For the HRF analysis, we used at least 50 3-week-old plants per lineages. Each experiment was repeated at least twice.

Lineages Used for Genomic and Epigenomic Study

The parental F2C plants and the progeny of 25th generation (both stressed progeny F25H and parallel control progeny F25C) plants were grown to about 21 days post germination (dpg), rosette leaves tissues were harvested from individual plants and snap frozen in liquid nitrogen, then they were stored at –80°C for DNA extraction for sequence analysis. From each generation, five individual plants were sequenced, resulting in a total of 15 samples for whole genome sequencing (WGS) and whole genome bisulfite sequencing (WGBS) analysis.

Whole Genome Sequencing and Whole Genomic Bisulfite Sequencing

The total genomic DNA was extracted from approximately 100 mg of leaf tissue homogenized in liquid nitrogen using a CTAB protocol. The isolated genomic DNA was used for both WGS and WGBS to assist in identifying the genomic and epigenomic variations. The WGS and WGBS libraries construction and sequencing have been carried out at the Centre d’expertise et de service Génome, Montreal, QC, Canada.

Whole Genome Sequencing Libraries Construction and Sequencing

gDNA was quantified using the Quant-iT™ PicoGreen^® dsDNA Assay Kit (Life Technologies, Carlsbad, CA, United States). Libraries were generated using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA, United States) as per the manufacturer’s recommendations. Adapters and PCR primers were purchased from IDT. Size selection of libraries contained the desired insert size has been performed using SparQ beads (Qiagen, Hilden, Germany). Libraries were quantified using the Kapa Illumina GA with Revised Primers-SYBR Fast Universal kit (Kapa Biosystems, Woburn, MA, United States). Average size fragment was determined using a LabChip GX (PerkinElmer, Waltham, MA, United States) instrument.

Whole Genome Bisulfite Sequencing Libraries Construction

gDNA was quantified using the Quant-iT™ PicoGreen^® dsDNA Assay Kit (Life Technologies, Carlsbad, CA, United States). Libraries were generated with the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, MA, United States) using 250 ng of input genomic DNA. Adapters were purchased from NEB. Size selection of libraries containing the desired insert size have been performed using sparQ PureMag Beads (Quantabio, Beverly, MA, United States). Bisulfite conversion had been carried out with the EZ DNA Methylation-Lightning Kit (Zymo Research, Irvine, CA, United States). Libraries were quantified using the Kapa Illumina GA with Revised Primers-SYBR Fast Universal kit (Kapa Biosystems, Woburn, MA, United States) and average size fragment was determined using a LabChip GX (PerkinElmer, Waltham, MA, United States) instrument.

Whole Genome Sequencing and Whole Genomic Bisulfite Sequencing

The libraries were normalized and pooled and then denatured in 0.05N NaOH and neutralized using HT1 buffer. ExAMP was added to the mix following the manufacturer’s instructions. The pool was loaded at 200 pM on an Illumina cBot and the flowcell was ran on a HiSeq X for 2 × 151 cycles (paired-end mode). A phiX library was used as a control and mixed with libraries at 1% level. The Illumina control software was HCS HD 3.4.0.38, the real-time analysis program was RTA v. 2.7.7. Program bcl2fastq2 v2.20 was then used to demultiplex samples and generate fastq reads.

The Computation and Analysis of Genome Sequence Data

Raw sequencing reads were quality controlled and trimmed using Trim Galore software (version 0.4.4). The trimmed reads were aligned to the TAIR10 reference genome, and the duplicates were marked using Picard tool. Local realignments around SNPs and INDELs were performed using GATK (genome analysis toolkit) which accounts for genome aligners and mapping errors and identifies the consistent regions that contain SNPs and INDELs. The resulting reads were quality controlled with Haplotype scores and variant sample sites were called individually and jointly using the HaplotypeCaller with GATK. The sites marked as those that had a low-quality score by GATK were filtered out. The effects of variants in the genome sequences were classified using the SnpEff program (Cingolani et al., 2012). According to the SnpEff program (Cingolani et al., 2012) utilized in this study, upstream was defined as 5 kb region upstream of the distal transcription start site (TSS), and downstream was defined 5 kb region downstream of the most distal polyA addition site. Variants affecting the non-coding regions were expounded, and biotypes were identified with the available information after comparing with the TAIR 10 reference Arabidopsis genome.

The genomation tools were used to obtain a biological understanding of genomic intervals over the pre-defined functional regions such as promoters, exons, and introns and the Functional Classification SuperViewer to create gene association profiles and show the difference between samples. The genes nearest to the non-overlapping SNP and INDEL sites were annotated.

The Computation and Analysis of Whole Genome Bisulfite Sequencing Data

Whole genome bisulfite sequencing allows for the investigation of genome-wide patterns of DNA methylation at a single-base resolution. It involves the sodium bisulfite conversion of unmethylated cytosine into uracil, with the resulting cytosine residues in the sequence representing methylated cytosine in the genome which is then mapped to a reference genome (Susan et al., 1994). Binomial tests were applied and used to determine the observed methylation frequency against the bisulfite conversion reaction, and the percentage of methylation levels were calculated at each base (Schultz et al., 2012).

The WGBS raw sequencing data were analyzed using tools found in the methylKit package. Raw sequencing reads were quality controlled and trimmed using Trim Galore software (version 0.4.4). The trimmed reads were then aligned to the TAIR10 reference genome using the bisulfite mapping tool Bismark (Krueger and Andrews, 2011). The methylated cytosines (mCs) were extracted from the aligned reads with the Bismark methylation extractor on default parameters followed by the computation of methylation frequency using the R package software, methylKit. The percentage of methylation was calculated by counting the frequency ratio of Cs divided by reads with C or a T at each base and computed at bases with coverage ≥ 10 (Akalin et al., 2012).

$\begin{array}{c}\mathrm{The}\mathrm{percentage}(\%)\mathrm{of}\mathrm{methylation}\hfill \\ \mathrm{\hspace{1em}}=\{\mathrm{Frequency}\mathrm{of}\mathrm{C}÷\mathrm{readcoverage}\}\times 100\hfill \end{array}$

Common bases covered across all samples were extracted and compared, and the differential hyper- and hypo- methylated positions in each chromosome were extracted. The DMPs overlapping with genomic regions were assessed (in the preference of promotor > exon > intron), and the average percentage of methylation of DMPs around genes with the distances of DMPs to the nearest TSSs were also calculated.

The annotation analysis was performed with genomation package within methylKit to obtain a biological understanding of genomic intervals over the pre-defined functional regions such as promoters, exons, and introns (Akalin et al., 2014). The functional commentary of the generated gene methylation profiles was performed using the SuperViewer tool with Bootstrap to show the difference between samples (Provart and Zhu, 2003). Hierarchical clustering of samples was used to analyze similarities and detect sample outliers based on percentage methylation scores. Principal Component Analysis (PCA) was utilized for variations and any biological relevant clustering of samples. Scatterplots and bar plots showing the percentage of hyper-/hypo- methylated bases, the overall chromosome and heatmaps were used to visualize similarities and differences between DNA methylation profiles.

The Differentially Methylated Regions

A comparison of differential DNA methylation levels between samples reveals the locations of significant differential changes in the epigenome. The obtained information of DMRs was investigated over the predefined regions in all contexts: CG, CHG, and CHH for 100 and 1,000 bp tiles across the genome to identify both stochastic and treatment associated DMRs (Akalin et al., 2012).

The differential hyper-/hypo- methylated regions were extracted and compared across the samples. By default, DMRs were extracted with q-values < 0.01 and the percentage of methylation difference > 25%. The differential methylation patterns between the sample groups and the differences in methylation events of per chromosome were also extracted. The methylation profiles of the sample groups used were F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C. In summary, the sliding windows of 100 and 1,000 bp were considered for both DMRs and DMPs, and extractions were made based on at least 25 and 50% differences (q-values < 0.01) to assess the significant differences among samples.

Biological Enrichment Analysis

The functional classification of variants (DMRs and DMPs) unique to each test group was interpreted using SuperViewer to identify regions with the statistically significant number of over- or under-represented genes and genomic features. Biological processes that might be enriched or under-represented within and between generations were assessed. All values were normalized by bootstrap x100, and p-values < 0.05 were retrieved as significant.

Statistical Analysis and Quality Control Values

The mapped reads were obtained with a quality score of <30, the differential hyper- and hypomethylated bases were extracted with q-values < 0.01 and the percentage of methylation difference larger than 25% in methylKit. The Heatmaps of differentially methylated bases were quantified at q-values < 0.01, and the percentage of methylation difference was more significant than 50%. The distances of DMPs to the nearest TSSs obtained from genomation at both >25 and >50% of methylation change. The distance between TSSs and DMPs was extracted within ±1,000 bp and annotated at DMPs > 50% methylation difference. DNA methylation profiles obtained from the melthylKit used the obtained pairwise correlation coefficients of methylation levels (in %) and the 1-Pearson’s correlation coefficients for hierarchical clustering of samples. Logistic regression and Fisher’s exact test were used for the determination of differential methylation with calculations of q-values and the Benjamini–Hochberg procedure for the correction of p-values. The T-test for mean difference between groups was calculated with p-values < 0.05. The results of global genome methylation were analyzed by one-way ANOVA with the Tukey’s multiple comparisons test using GraphPad Prism version 8.4.2 for Windows. The data are shown as the average percentage (with ± SD) of methylated cytosines from five individual methylomes in each of the progenies. The asterisks show a significant difference between the stressed progeny and the parental control progeny (*p < 0.1, ^**p < 0.05). The phenotypic data were analyzed by the unpaired t-test with Welch’s correction using GraphPad Prism version 9.1.1 (255) for Windows. The data were shown as mean ± SD. A p-value less than 0.05 (p ≤ 0.05) was considered statistically significant.

Results

The Progeny of Stressed Plants Showed the Heat Stress Tolerance Phenotype

The F25H progeny was studied for heat stress memory and response. The response to heat stress was determined at the germination, seedling or mature-plant stage. The F25H stressed progeny did not display any significant difference compared to F25C in germination index represented by the speed of germination and the total germination percentage (Figure 1B). When we exposed 10 days old seedlings of F25C and F25H groups to the same stress that their ancestral progenies were exposed to, the F25H showed slightly higher fresh weight than F25C under heat stress, however, it was not significantly different (Supplementary Figure 1); it should be noted that fresh weight of F25H plants was significantly higher than F25C plants without heat stress. Under heat stress, at the mature-plant stage, the F25H stressed progeny survived better than their respective non-stressed parallel control progenies (p < 0.001) (Figure 1C). The F25H heat-stressed progeny plants also displayed a higher fresh weight when compared with parallel progeny plants.

The Progeny of Stressed Plants Exhibited a Higher Homologous Recombination Frequency

Our previous reports demonstrated that plants exposed to stress and their progeny exhibit higher HRF (Kathiria et al., 2010; Migicovsky and Kovalchuk, 2013). Here we also found that the F25H stressed progeny exhibited a higher HRF as compared with parallel control (Figure 1D, right panel).

The Multigenerational Heat-Stressed Progeny Exhibited Higher Single Nucleotide Polymorphisms and Insertions and Deletions

Genomic variants consisting of SNPs and INDELs were determined by mapping and comparing F25H, F2C and F25C sequencing reads to TAIR 10 reference genome. The total SNPs identified for F25H, F2C and F25C were 53,678, 15,599, and 15,526, respectively (Figure 2A). The total numbers of INDELs identified for F25H, F2C and F25C were 17,669, 7,986 and 8,143, respectively (Figure 2A). After combining SNPs and INDELs, the total genetic variations in F25H, F2C and F25C were 71,347, 23,585, and 23,669, respectively (Figure 2A). We have also determined the SNPs and INDELs by using stricter parameters, which displayed an overall significantly reduced number of SNPs (12,000 in F25H; 5,190 in F2C; 5,188 in F25C) and INDELs (6,177 in F25H; 4,643 in F2C; 4,781 in F25C), however, comparative patterns and conclusion remained almost similar (Supplementary Table 1). Relaxed parameters can result in a higher number of false positive genetic variations, and stricter parameters can miss true positives. However, both relaxed and strict parameters outputs resulted in almost similar genetic variations patterns between groups. The parallel F25C and the parental F2C controls exhibited a comparable number of mutations with a slight difference between the numbers of SNPs and INDELs (Figure 2A). Thus, the F25H heat-stressed progeny showed over three times more mutations than the parental and parallel controls (in case of stricter parameters around two times; Supplementary Table 1), suggesting that heat-stress induced genetic variations in the stressed progeny (Figure 2A). In addition, INDELs were of a similar size in F25C and F2C, whereas they were significantly larger in F25H (Figure 2B); this is likely caused by multigenerational heat-stress exposure.

FIGURE 2

Figure 2. The heat-stressed progeny showed a higher number of genetic variations (with relaxed parameters). (A) The total number of SNPs and INDELs in the genomes of stressed progeny and in the parallel and parental control progenies when all samples were jointly considered. Genomic variants consisting of SNPs and INDELs were determined by mapping and comparing F25H, F2C, and F25C sequencing reads to TAIR 10 reference genome. (B) The box plot for the length of insertions and deletions of variants within F2C, F25C, and F25H. The asterisk above bracket (*) shows a significant difference between F25H vs. F2C and F25H vs. F25C where *** indicates p < 0.001. Wilcoxon rank-sum text was used to determine a statistical difference (p < 0.05). (C) A Venn diagram of jointly called variants. The non-overlapping area of Venn indicates the number of unique SNPs and INDELs in each sample group. Sample variants were jointly called by the “GATK haplotype Caller” using five samples together. (D) The total number of variant effects by major genomic regions. The number of variant effects in the major genomic regions are shown above each bar, and the genomic regions are indicated on the x-axis.

The GATK Haplotype caller calls SNPs and INDELs via an assembly of haplotypes in the active region of the genome which enables the identification of variants that are unique to each sample genome by extracting the non-overlapping sites specific for F25H, F25C, and F2C. The analysis of non-overlapping site-specific variants revealed 39,712, 1,313, and 1,308 unique SNPs for the F25H, F2C, and F25C plants, respectively; 11,478, 1,576, and 1,801 unique INDELs were identified for F25H, F2C, and F25C, respectively (Figure 2C). In case of analysis with stricter parameters, we observed an overall significantly reduced number of unique SNPs (7,474 in F25H; 613 in F2C; 629 in F25C) and INDELs (2,467 in F25H; 794 in F2C; 953 in F25C), however, comparative patterns remained almost similar to relaxed parameters (Supplementary Table 1). Thus, the F25H plants showed a dramatically higher number of unique SNPs and INDELS than the control plants, suggesting that exposure to multigenerational heat predominantly led to the generation of SNPs and INDELs (Figure 2C). The parental F2C and parallel F25C controls showed a very similar number of unique SNPs (Figure 2C).

The Heat-Stressed Progeny Exhibited a Higher Number of the Effects of Genetic Variants

The total number of genetic variant effects calculated by the corresponding genomic locations of SNPs and INDELs such as the intronic, exonic, untranslated regions (5′UTR or 3′UTR), upstream and downstream of gene regions and intergenic regions were calculated. Among SNPs, the potential effect on coding sequences was analyzed by calculating the number of synonymous or non-synonymous mutations. Positional profiles of SNPs and INDELs revealed a higher number of variant effects in the F25H heat-stressed progeny compared with the non-stressed parallel F25C and parental F2C progenies (Figure 2D). The F25H plants showed the highest total number of variant effects (293,405) compared with the F25C (82,179) and F2C (81,657) controls (Figure 3). The classification of variant effects by genomic regions revealed the largest number of effects observed in upstream of the gene (F25H – 110,837; F25C – 28,762; F2C – 28,403) followed by downstream (F25H – 103,301; F25C – 28,035; F2C – 27,947) and intergenic regions (F25H – 40,941; F25C – 15,184; F2C – 15,058) (Figures 2D, 3). Variant effects were also found within the 5′ untranslated (F25H – 1,817; F25C – 438; F2C – 442) and the 3′ untranslated regions (F25H – 2,312; F25C – 461; F2C – 442) of the genome (Figure 3). It appears that the difference in variants was most drastic in 5′ and 3′ untranslated regions, with ∼4.2- and 5.0-fold difference between F25H and either of the control samples, followed by upstream (∼3.9), downstream (∼3.7) and intergenic regions (∼2.7).

FIGURE 3

Figure 3. The distribution of variants across the genome by their potential effects based on genomic positions (with relaxed parameters). (A) F25H, (B) F25C, and (C) F2C. The number of variant effects is shown above each bar, and genomic positions effects from SnpEff are listed on the x-axis.

The F25H plants showed 12x and 6x higher number of synonymous and non-synonymous mutations, 3,891 and 4,510, respectively, than the F25C plants (313 - synonymous, 752 – non-synonymous) and the F2C plants (316 - synonymous, 742 – non-synonymous) (Figure 3). In F25H, 18 SNP missense and 66 SNP nonsense mutations were identified (Figure 3C). These numbers were much lower in F25C and F2C where 5 and 5 stop codon losses and 17 and 10 stop codon gains were observed, respectively (Figures 3A,B). The non-synonymous mutations contribute to biological variations in the living organism, consequently, they are subjected to natural selection.

The Analysis of Epigenetic Variations

Multigenerational Exposure to Heat Stress Reduces the Percentage Methylation of Total Global DNA in the CHG and CHH Contexts

The bisulfite sequencing data revealed that the average percentage of global genome methylation in the CG context did not show any significant difference between the F25H plants and either of the control samples (Figure 4A). However, in the CHG context, the F25H plants showed a significantly lower global methylation level (p < 0.10) than the control plants (Figure 4A). Similar data were found for the CHH context; the F25H plants exhibited a significantly lower methylation level (p < 0.05) compared with the controls (Figure 4A). The reduction of methylation levels in the CHG and CHH contexts might be a part of heat stress adaptation strategy.

FIGURE 4

Figure 4. Global methylation changes, differentially methylated positions (DMPs) and differentially methylated regions (DMRs) in F25H, F25C, and F2C plants. (A) The average percentage of methylated cytosines in F2C, F25C, and F25H in the CG, CHG, and CHH sequence contexts (H = A, T, C and mean, n = 5). Methylation levels were determined from reads with the minimum coverage ≥ 10 mapped to the TAIR 10 reference genome by using Bismark software. The data were analyzed by one-way ANOVA with the Tukey’s multiple comparisons test using GraphPad Prism version 8.4.2 for Windows. The data are shown as the average percentage (with ± SD) of methylated cytosines from five individual methylomes in each of the progenies. The asterisk shows a significant difference between the stressed progeny and the parental control progeny (*p < 0.1, **p < 0.05). (B) The total number of DMPs in F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C. (C) The total number of DMPs in the CG, CHG, and CHH contexts in F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C. (D) The total number of DMRs (100 and 1000 base pair window) in F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C comparison groups. (E) The total number of DMRs (100 base pair window) in CG, CHG, and CHH contexts in F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C.

The Analysis of the Number of Differentially Methylated Positions and Differentially Methylated Regions Reveals a Higher Number of Differentially Methylated Cytosines in the CHG and CHH Contexts in the F25H Plants

The total number of DMPs for F25H vs. F2C, F25H vs. F25C and F25C vs. F2C were 66,491, 78,412, and 80,464, respectively, suggesting that there is no effect of heat stress on the total number of DMPs (Figure 4B). The analysis of DMPs in a specific sequence context revealed a different pattern. While DMPs in the CG context were higher in the F25C vs. F2C control comparison group, in the CHG and CHH contexts, DMPs were much higher in comparison groups involving the F25H stressed progeny (Figure 4C). The difference was especially obvious for hypomethylated cytosines at the CHG context; the F25H vs. F2C group showed 299 hypomethylated DMPs, while the F25C vs. F2C group – only 230. A similar pattern was observed in the CHH context, nearly a twofold larger number of hypomethylated DMPs was observed in the F25H vs. F2C group in comparison with the F25C vs. F2C group – 342 vs. 176, respectively; the difference for hypermethylated DMPs was also substantial, 133 vs. 87 (Figure 4C). These data suggest that epimutations induced by multigenerational heat stress are primarily associated with changes in the CHG and CHH contexts, and that hypomethylation is likely a prevalent mechanism.

The total DMRs in the 100 bp window for F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C were 10,557, 13,068, 12,810, respectively (Figure 4D). In the case of the 1,000-bp window, the total DMRs for F25H vs. F2C, F25H vs. F25C and F25C vs. F2C were 49, 36, 32, respectively. As to the CHG and CHH context, there was a larger number of DMRs in the F25H vs. F2C group compared with other groups (Figure 4E).

The Progeny Exposed to Multigenerational Heat-Stress Did Not Show Any Difference in the Methylation Levels

The percentage of methylation (the analysis of the percentage of methylated cytosines at specific sites) in F2C, F25C, and F25H showed that most of the cytosines had either a high (70–100%) or a low (0–10%) level of methylation in the methylated CG context, however, we did not see any difference in methylation levels between F2C, F25C, and F25H progeny (Supplementary Figure 2). Cytosine methylation in the CHG context was much lower compared with the CG context; most sites had 0–60% methylation levels including very high number of cytosines methylated at 0-10%, and sites with high levels of methylation (70–100%) were not observed (Supplementary Figure 3). Similarly, methylation in the CHH context was much lower than in the CG or CHG contexts; there was an even lower frequency of occurrence of sites with the 0–20% methylation level (Supplementary Figure 4). When each sequence context was considered, there was no significant difference in the methylation levels between the F2C, F25C, and F25H groups. In general, F2C, F25C, and F25H showed a similar percentage of methylation in the CG, CHG, and CHH contexts (Supplementary Figures 2–4).

The Heatmaps Based Hierarchical Clustering Analysis of Differentially Methylated Positions and Differentially Methylated Regions Revealed Directional Changes of Heat-Induced Epimutations in the Stressed Progeny

We analyzed the relatedness of the control and stressed progeny samples using the hierarchical clustering heatmap analysis (Figure 5). The heatmap analysis of changes in both hypo- and hyper-methylated DMPs revealed a clear separation of five samples of F25H from five samples of the F2C parental progeny in all three CG, CHG, and CHH cytosine methylation contexts (Figure 5A). Similarly, a comparison between F25H and F25C also showed separate clustering (Supplementary Figure 5A). A comparison between F25C and F2C samples showed more similarity, with two F25C samples clustering together with F2C samples (Supplementary Figure 5B). Heatmap clustering of CG-DMRs also revealed a clear separation between F25H and F2C samples (Figure 5B). Considering that CHH methylation happens in plants de novo, separate clustering of stressed samples from control progeny suggests the importance of CHH methylation in stress adaptation and indicates that multigenerational exposure to heat stress likely induces non-random changes in the methylome in the CHH context.

FIGURE 5

Figure 5. A hierarchical clustering heatmap analysis. (A) Heat maps of DMPs for hypermethylated cytosines (the upper panel) and hypomethylated cytosines (the lower panel) in CG, CHG, and CHH contexts in F25H vs. F2C. Differentially methylated cytosines in the genome with differences of > 50% in F25H vs. F2C. (B) Heat maps of DMRs of hypermethylated cytosines (the upper panel) and hypomethylated cytosines (the lower panel) in CG contexts in F25H vs. F2C. Differentially methylated positions in the genome with differences > 50% in F25H vs. F2C. In ‘the upper panel’, the red section indicates the larger percentage of methylation, and the yellow section indicates the lower percentage, and in “the lower panel,” the green section indicates the larger percentage of methylation and the yellow one indicates the lower percentage, q-value < 0.01.

Annotation of Genes With Differentially Methylated Positions and Differentially Methylated Regions

Differentially methylated positions and DMRs were then characterized to determine whether they were preferably located near genes. The location of hypo- or hyper- methylated DMRs was compared to the annotated Arabidopsis genes using genomation Bioconductor package (Akalin et al., 2014). Both hyper- and hypo-methylated CG-DMRs were predominantly located in the gene body, while were less frequent at TSS (Transcription Start Site) and TTS (Transcription Termination Site) regions (Figures 6A,C). In contrast, in the case of TEs, F25H progeny displayed less pronounced hypermethylation in CG-DMRs located in TEs body as compared to upstream and downstream regions of TEs (Figures 6B,D). Overall, methylation changes in the gene body were much more pronounced than in TEs in F25H group as compared to other groups (Figure 6).

FIGURE 6

Figure 6. The distribution of hyper-and hypo-methylated CG-DMRs in the gene body and TEs in F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C comparison groups. (A) Distribution of hypermethylated CG-DMRs in gene bodies. (B) Distribution of hypermethylated CG-DMRs in TEs. (C) Distribution of hypomethylated CG-DMRs in gene bodies. (D) Distribution of hypomethylated CG-DMRs in TEs.

In the CG context, the percentage of hypermethylated and hypomethylated DMPs in all groups (F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C) was the highest in the exon region, followed by the promoter, intron and intergenic regions (Figures 7A–F). In the case of the CHG context, all comparison groups (F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C) showed the highest percentage of hypermethylated DMPs in the promoter, followed by intergenic, exon, and intron regions (Figures 7A,C,E). Interestingly, F25H vs. F2C groups showed the highest percentage of hypomethylated DMPs in the intergenic regions (48%) followed by the promoter (44%), while in the F25C vs. F2C group, the highest percentage of DMPs was in the promoter region (52%) followed by the intergenic regions (44%) (Figure 7F). In the CHH context, in all groups, the percentage of both hyper- and hypo- DMPs was higher in the promoter, followed by the intergenic, exon, and intron regions (Figure 7).

FIGURE 7

Figure 7. The distribution of hyper-and hypo-methylated DMPs in the genome in the comparison groups F25H vs. F2C (A – hyper, B – hypo), F25H vs. F25C (C – hyper, D – hypo) and F25C vs. F2C (E – hyper, F – hypo). The data are shown as an average percentage of differentially methylated cytosines at the promoter, exon, intron and intergenic regions. DMPs were considered as regions of methylation with > 25% differences.

Overall, in the case of the CG context, the highest percentage of DMPs was observed in the gene body, whereas in the case of the CHG and CHH contexts, it was highest in the promoter and intergenic regions.

Comparative Analysis Revealed the Heat-Inducible Nature of Genetic Variants and the Spontaneous Nature of Epigenetic Variants

Strict genetic analysis showed that the ratio between the number of epigenetic and genetic variants was 6.27x, 7.3x, and 35.9x in the comparison groups of F25H vs. F2C, F25H vs. F25C, and F25C vs. F2C, respectively (Table 1), indicating that epigenetic variations were more prevalent than genetic variations. The largest difference in the ratio of epigenetic to genetic variations was observed between F25C and F2C control plants, which indicated the spontaneous nature of epigenetic variations. In the case of epigenetic variations, F25C parallel control progeny showed 1.2 times (80,464 vs. 66,491) more variation than F25H stressed progeny when both were compared with the parental control progeny F2C. In contrast, the F25H stressed progeny displayed five times (10,597 vs. 2,238) more genetic variations than the F25C parallel control progeny when both were compared with the parental control progeny F2C. These comparisons suggest that heat stress predominantly leads to the accumulation of genetic rather than epigenetic variations.

TABLE 1

Table 1. Comparison of epigenetic and genetic variations in F25H, F25C, and F2C plants.

There Was No Difference in the Enrichment of Genetic Variants for Any Biological Process in F25H

The most enriched biological process for SNPs was classed as “unknown biological processes” regardless of the tested group (Supplementary Figure 6). The only over-represented group was “transport” found in F2C INDELs sample. All other biological processes were underrepresented for all samples, and no significant difference between the F25H group and any other group was observed (Supplementary Figure 7 and Supplementary Table 2).

Epimutation-Associated Genes Participate in Stress Response Pathways

Gene Ontology (GO) analysis of epimutations (DMRs) in the hypermethylated CG context in F25H group showed a significant difference in the enrichment of biological processes such as “response to abiotic or biotic stimuli,” “cell organization and biogenesis” and “DNA or RNA metabolism” compared to F2C and F25C groups (Figure 8 and Supplementary Table 3). In the hypomethylated CG context, no significant difference was observed between the F25H vs. F2C group and other groups (Supplementary Figure 9 and Supplementary Table 4).

FIGURE 8

Figure 8. Enrichment analysis of hypermethylated DMRs on CG sites and their classification based on biological processes. Y-axis shows normalized Class Score with binomial coefficients as calculated by SuperViewer. To calculate enrichment, *p-value of <0.05, ±bootstrap StdDev were used.

In the hypomethylated CHG context, the F25H vs. F2C group exhibited the enrichment in “protein metabolism process” which was not observed in other groups (Supplementary Figure 9A and Supplementary Table 5). In the hypomethylated CHH context, the F25H vs. F2C group exhibited the enrichment in “transport,” “other metabolic processes,” and “developmental” processes (Supplementary Figure 9B and Supplementary Table 6).

Discussion

The current study aimed to understand the effects of multigenerational exposure to heat stress on plant phenotypes, genotypes, and epigenotypes. Specifically, we compared the progeny of plants exposed to heat stress over 25 consecutive generations with their counterparts (the parental and parallel control progenies). We studied the genetic variations (SNPs and INDELs) by using the whole genome sequencing analysis as well as the epigenetic variations (DMPs and DMRs) by using the whole genome bisulfite sequencing analysis of the heat-stressed and control progenies. Overall, this study has deciphered the phenotypic, genomic and epigenomic changes in the heat-stressed progeny in response to multigenerational exposure to heat stress.

Multigenerational Heat Stress Induced the Phenotypic Resilience in the Stressed Progeny

The phenotype of the progeny can be affected by the parental environment (Groot et al., 2016). Multigenerational stress exposure elicits different outcomes compared to single generation exposure (Groot et al., 2016). Our study showed that the progeny derived from 25 generations of heat stress displayed a notable stress tolerance at the mature stage (Figure 1C). Like our study, Wibowo et al. (2016) reported that the direct progeny of G2–G5 stressed plants exhibited the higher germination and survival rates under high-salinity stress. The authors concluded that stress priming requires repetitive exposure to stress. In some cases, it was evident that parental heat stress affected the response of offspring generation because heat stress accelerated leaf production to speed up the development of the progeny of stressed plants (Porter, 2005; Liu et al., 2015). Our previous studies showed that the heat-induced modifications of gene expression might have a role in phenotype changes such as leaf number and flowering time (Migicovsky et al., 2014). Zheng et al. (2017) also found that multigenerational exposure to drought improved the drought adaptability of rice offspring. The analysis of plant weight is useful for determining the active growth and development of a plant (Yadav et al., 2012). F25H heat-stressed progeny plants displayed a higher fresh weight when compared with parallel progeny plants.

The degree or severity of high temperature stress can impact cellular homeostasis (Kotak et al., 2007), and exposure to above the optimum growth temperature over time can result in the acquired thermotolerance (Bray et al., 2000; Sung et al., 2003). Moreover, the expression of transgenerational effects in the form of changes of phenotype responses strongly depends on the offspring environment (Boyko and Kovalchuk, 2010; Groot et al., 2016) and varies between the genotype and trait (Verhoeven and van Gurp, 2012; Suter and Widmer, 2013a). Some reports have mentioned the effects of epigenetic inheritance or transference of transgenerational stress memory to be restricted to a single generation without the maintenance of stress (Pecinka et al., 2009; Suter and Widmer, 2013b). Also, it was shown that multigenerational exposure often reduces the expression of parental effects compared to single-generation exposure (Groot et al., 2016).

Both intergenerational and transgenerational changes in response to stress include alterations in HRF. The progeny of plants that were stressed show the elevated HRF (Migicovsky and Kovalchuk, 2013). In the current report, under non-stress conditions, the stressed progeny F25H exhibited the elevated HRF compared with their parallel control progenies (Figure 1D). It has been reported that plants exposed to salt and UV stresses show a much higher HRF than control plants (Boyko et al., 2010). The progeny of pathogen-infected plant has a higher level of somatic HRF, and the progeny also show the decreased methylation at loci for resistance to TMV and they have global hypermethylation (Boyko et al., 2007). Likewise, the progeny of plants exposed to heavy metal stress showed high HRF (Rahavi et al., 2011). If stresses persist through multiple generations, it is possible that the epigenetic memory will be stronger and persist longer in the absence of stress, thus exhibiting changes in HRF for a longer period of time (Migicovsky and Kovalchuk, 2013).

The Analysis of Genomic Variations Induced by Multigenerational Exposure to Heat Stress

In our study, we have shown that the progeny of heat-stressed plants (F25H) had 3x higher frequency of genomic variations than the non-stressed control plants. It is well established that environmental stresses can be mutagenic and can cause genome instability (Boyko et al., 2010; Boyko and Kovalchuk, 2011; Gill et al., 2015). Previously, Jiang et al. (2014) reported around 2x more genetic variations in 10th generation of salinity stress lineages compared to control lineage. Environmental changes may also increase homologous recombination events or can facilitate the mobilization of transposable elements (Long et al., 2009; Boyko et al., 2010; Migicovsky and Kovalchuk, 2013). Eventually, environmental conditions may induce mutations and increase the chances of genomic diversity leading to the adaptation to cope up with the ever-changing environmental conditions (Filichkin et al., 2010).

The position of SNPs may influence the rate of transcription directly or indirectly depending on their corresponding genomic location; SNPs can also directly alter the coding sequence leading to a synonymous or non-synonymous amino acid replacement (Cingolani et al., 2012). SNPs and INDELs are vital genetic variations that can directly disrupt gene function and affect plant adaptability in the changing environment (Shastry, 2009). For instance, SNPs can affect light response and flowering time by changing amino acids in phytochromes A and B (Maloof et al., 2001; Filiault et al., 2008). Jain et al. (2014) indicated that INDELs occur in different genomic regions, and they further suggested that positions of INDELs affect the function and expression of genes in a different manner. In the current study, genetic variants were mostly found in the upstream and downstream region of the gene and were significantly higher (3x) in the F25H stressed progeny than in F25C parallel and F2C parental control progenies (Figures 2A,D). In contrast, the number of genetic variations was almost similar between the F25C parallel and F2C parental control progenies (Figure 2A). Very recently, Belfield et al. (2020) also showed higher mutation frequency in 11th generation of mild heat stress lineage (29°C) than control lineage. Positional profiles of genetic variations in our work revealed the higher number of variant effects in the F25H heat-stressed progeny compared with the non-stressed parallel F25C and parental F2C progenies (Figure 3). The F25H heat-stressed progeny exhibited higher rate of synonymous (12x) and non-synonymous (6x) mutations than the parallel and parental control progenies (Figure 3). SNPs can create new splice sites and alter gene functions (Guyon-Debast et al., 2010). Therefore, the higher number of genetic variations in the heat-stressed progeny may be an indication of the adaptive processes occurring in the progeny of stressed plants.

Interestingly, the INDELS were of a similar size in F25C and F2C, whereas they were significantly larger in F25H (Figure 2B), suggesting that multigenerational exposure to heat stress resulted in an increase in the size of INDELs in the heat-stressed progeny. Larger deletions may have a greater effect on gene expression and changes in the phenotype. It has been reported that small and large INDELs can cause pathogen sensitivity (Mindrinos et al., 1994; Kroymann et al., 2003).

The Analysis of Epigenetic Variations Induced by Multigenerational Exposure to Heat Stress

The mechanism behind intergenerational and transgenerational stress memory in plants likely includes changes in the DNA sequence, DNA methylation and chromatin structure (Crisp et al., 2016; Lind and Spagopoulou, 2018; Zhang et al., 2018; Perez and Lehner, 2019; Miryeganeh and Saze, 2020). The inheritance of changes in DNA methylation is one of the mechanisms of adaptation – epialleles often lead to the appearance of new alleles (Becker et al., 2011; Paszkowski and Grossniklaus, 2011; Zhang et al., 2018). It has been reported that abiotic stresses can cause hyper- or hypomethylation in the genome after either short-term or long-term exposure (Uthup et al., 2011; Zhang et al., 2018). Several studies suggested that multigenerational stress exposure could cause substantially higher heritable epigenetic variations compared with single-generation exposure (Remy, 2010; Groot et al., 2016; Zheng et al., 2017). This phenomenon was explained by the gradual acclimatization of epigenetic effects (Groot et al., 2016).

In our study, the bisulfite sequencing data revealed that the average percentage of global genome methylation in the CG context was similar among all three groups (Figure 4A). However, in the CHH contexts, the F25H stressed progeny showed a lower global methylation than the control progenies (Figure 4A). The data suggest that the reduction of methylation in the CHG and CHH contexts might be a part of adaptation strategies to heat stress. In plants, de novo DNA methylation at all three cytosine sequence contexts is established by RNA-directed DNA methylation (RdDM) pathway via DRM2 (DOMAINS REARRANGED METHYLTRANSFERASE 2) methyltransferase that is directed to the target loci by siRNAs (mainly small 24 nucleotide short interfering RNAs). Once the DNA methylation is established, it is maintained by various methyltransferases, including MET1 (METHYLTRANSFERASE 1; mainly CG context), CMT3 (CHROMOMETHYLASE 3; mainly CHG context), and DRM2 (mainly CHH context) (Law and Jacobsen, 2010; Stroud et al., 2013). Thus, CHH methylation is exclusively established and re-established by the de novo RdDM pathway (Erdmann and Picard, 2020; Wassenegger and Dalakouras, 2021). Therefore, upon meiosis/mitosis, and in the absence of RNA triggers, CG and CHG methylation can be maintained by MET1 and CMT3, respectively, while CHH methylation is lost. In the current study, the global CG methylation level was almost unaffected in F25H stressed progeny, but the global CHH methylation level was significantly reduced, which strongly suggests that under heat stress, the RdDM pathway was affected, however, the MET1 pathway was not impaired. It has been reported that abiotic stresses can cause hyper or hypomethylation in the genome after either short-term or long-term exposure (Uthup et al., 2011; Thiebaut et al., 2019). High temperature (HT) induced a global disruption of DNA methylation in cotton anther, mainly causing CHH hypomethylation in HT-sensitive cotton varieties (Ma et al., 2018). A decrease in global DNA methylation might result in the interruption of glucose- and reactive oxygen species (ROS)- producing metabolic pathways which may lead to microspore sterility (Ma et al., 2018). Min et al. (2014) also reported that in cotton anthers, HT decreased the expression of S-ADENOSYL-L-HOMOCYSTEINE HYDROLASE1 (SAHH1) and methyltransferases DRM1/DRM3, resulting in the genome-wide hypomethylation. Another report revealed that in soybean, heat stress induced the global DNA hypomethylation, especially in the CHH context, in both root hairs and stripped roots (Hossain et al., 2017). In cultured microspores of Brassica napus cv. Topas, heat-shock treatment (32°C for 6 h) triggers DNA hypomethylation, particularly in the CG and CHG contexts (Li et al., 2016). However, Gao et al. (2014) reported that under heat stress (37°C for 2 h, and then 45°C for 3 h), global genome methylation levels were significantly higher in the heat-sensitive genotype of rapeseed than in the heat-tolerant genotype. In Brassica rapa, the global DNA methylation levels are relatively stable under heat stress, while changes in CHH methylation at TE suggest that CHH methylation may be important for heat-stress response and tolerance (Liu et al., 2018). The methylome analysis revealed that water deficit was associated with a decrease in CHH methylation in apple cultivars, which might result in the hypomethylated status of TEs (Xu et al., 2018). Wibowo et al. (2016) also showed that methylation changes in CHG and CHH were well correlated with hyperosmotic stress treatment, whereas changes in CG methylation did not show any correlation, suggesting that CG methylation patterns occurred stochastically in the stressed and non-stressed samples. Further, they confirmed that hyperosmotic stress directed DNA methylation changes primarily at non-CG sites, and these epigenetic modifications were associated with an acquired transient adaptation to stress (Wibowo et al., 2016). Exposure to abiotic stressors such as salt also causes DNA hypomethylation in the progeny; Jeremias et al. (2018) reported the inheritance of DNA hypomethylation in response to salinity stress in Daphnia magna. Also, Jiang et al. (2014) showed that lineages exposed to soil salinity stress accumulated more methylation at the CG sites than the control progenies. Recently, Atighi et al. (2020) showed a massive and genome-wide hypomethylation as a crucial plant defense mechanism in response to nematode or bacterial pathogen infection in rice and tomato. Further, they demonstrated that DNA hypomethylation in the CHH context was associated with a reduced susceptibility to root-parasitic nematodes in rice.

Hierarchical clustering of epimutations and the heat map analysis of changes in both hypo- and hypermethylated DMPs revealed that the F25H progeny of stressed plants showed a clear separation from plants of the parental and parallel control progenies (Figure 5 and Supplementary Figure 4), indicating directional epimutations in the F25H stressed progeny due to multigenerational exposure to heat stress. Ganguly et al. (2017) also found similar clustering patterns in response to heat. They showed that in Arabidopsis, the progeny of stressed plants clustered separately from the progeny of non-stressed plants in response to the treatment with heat stress (Ganguly et al., 2017). Becker et al. (2011) compared global DNA methylation patterns among 10 Arabidopsis thaliana lines propagated for 30 generations from a common ancestor. Hierarchical clustering based on DMPs grouped siblings of the 3rd and 31st generation lines together in separate groups and suggested that epimutations did not distribute randomly throughout the genome (Becker et al., 2011). Similar to our study, Zheng et al. (2017) showed that the multigenerational drought stress induced directional epimutations in the methylome of rice offspring. Eichten and Springer (2015) reported that the separate hierarchical clustering of epimutations was associated with cold stress treatment.

We observed that the total number of DMPs was higher between the control groups F25C vs. F2C (80,464), compared with the stressed progeny groups F25H vs. F2C (66,491) and F25H vs. F25C (78,412) (Figures 4B,C). It was somewhat surprising because we expected an opposite result; we hypothesized that heat stress would induce epimutations and expected to see a higher number of epigenomic variations in the progeny of heat stressed plants compared with the controls. It is possible that heat stress does not result in an overall increase in epimutations, but it rather causes epimutations to occur at the specific, heat stress-related sites in the genome. The fact that the observed progenies of heat-stressed plants cluster together confirms this hypothesis. Becker et al. (2011) have compared the genome wide DNA methylation patterns between the 31st generation (that grew under control conditions) and the ancestor generation (the 3rd) and deciphered 30,000 differentially methylated cytosines (DMPs). More than one third of recurrent events were observed in the whole DMPs and DMRs, indicating the non-random nature of these epimutations. Further, they reported an average of 3,300 DMPs between immediate generations; this was around three times more than what would have been expected from 30,000 DMPs accumulated between the 3rd and the 31st generation lines (Becker et al., 2011). This data suggests that the number of epimutations does not increase linearly with time, indicating that only a fraction of new DMPs is maintained by the transgenerational inheritance over the longer term (Becker et al., 2011).

The analysis of DMPs in a specific sequence context revealed a different pattern. While DMPs in the CG context were higher in the F25C vs. F2C control group, in the CHG and CHH context, DMPs were higher in comparison groups involving the F25H stressed progeny (Figure 4C). The F25H vs. F2C comparison group showed the higher number of hypomethylated and hypermethylated DMPs in the CHG and CHH contexts compared with the control comparison groups (Figure 4C). These data suggest that the epimutations induced by multigenerational exposure to heat stress are primarily associated with changes in the CHG and CHH contexts, and that hypomethylation is a prevalent mechanism. A similar phenomenon has been narrated in response to heavy metal stresses where it has been reported that hypomethylation prevails at several loci in hemp and clover (Aina et al., 2004). In epigenetic mutants, methylation patterns of DMPs link DNA methylation to the response to the environmental signals like heat stress (Popova et al., 2013) and salt stress (Yao et al., 2012). Several reports also suggest that it is common for plants to methylate cytosines in the sequence contexts of CHG and CHH, but the reestablishment of epigenetic modifications is mostly guided by small RNAs or the heterochromatin-directed methylation pathways (Feng et al., 2010; Stroud et al., 2013). Further, we observed that in the case of the CG context, the highest percentage of DMPs was observed in gene body, whereas in the case of the CHG and CHH contexts, it was the highest in the promoter and intergenic regions. Changes in the CG sequences were usually found to cluster around the regulatory region of genes (Zheng et al., 2013), while changes in CHG and CHH in the gene body, likely contributing to changes in phenotypes (Bewick and Schmitz, 2017). It was reported before that the promoters of stress-responsive genes may be hypomethylated under stress conditions (Yao and Kovalchuk, 2011; Bilichak et al., 2012).

In our study, F25H stressed progeny displayed more pronounced methylation changes in gene body while less pronounced changes in the body of TEs (Figure 6). This is an interesting finding, as it is believed that TEs are quite unstable and are activated in response to stress. It remains to be shown whether differences in methylation of TEs in these groups would also result in lower transposon activity in F25H plants. Our previous work demonstrated that the immediate progeny of heat-stressed plants showed higher activity of TEs (Migicovsky et al., 2014). It is possible that multigenerational exposure to heat stabilizes transposons. Indeed, previously, Becker et al. (2011) also reported the transgenerational maintenance of CG methylation of TEs in response to propagation of Arabidopsis for 30 generations under normal conditions.

The GO analysis identified that epimutations in the F25H stressed progeny showed mainly enrichment in the processes such as response to abiotic and biotic stimulus, cell organizations and biogenesis, and DNA or RNA metabolism (Figure 8). Liu et al. (2018) found the enrichment of heat-induced DMRs in the RNA metabolism and modification process, which may suggest the importance of RNA metabolism for stress adaptation. Zheng et al. (2017) reported that in rice, the long-term drought stress-induced epimutations were directly involved in the stress-responsive pathways. Wibowo et al. (2016) showed that hyperosmotic stress-induced DMRs were enriched with functions related to metabolic responses and ion transport; it was also suggested that exposure to hyperosmotic-stress targeted the discrete and epigenetically labile regions of the genome. The DMR-associated genes in the drought-tolerant introgression line DK151 were mainly involved in stress response, programmed cell death, and the nutrient reservoir activity, which may contribute to the constitutive drought tolerance (Wang et al., 2016). Stress priming can be utilized to develop epi-populations with epigenetic variations in various crops (Lämke and Bäurle, 2017; Dalakouras and Vlachostergios, 2021). Some of these stress associated epigenetic variations (DMPs/DMRs) are stable and can faithfully transfer to the offspring (Wibowo et al., 2016; Zheng et al., 2017). These heritable stressed induced epigenetic variations epi-QTLs (DMPs and DMRs) can be used in epigenetic-assisted breeding and gene editing programs to develop superior next-generation stress tolerant crops (Dalakouras and Vlachostergios, 2021; Gogolev et al., 2021).

In conclusion, the stressed progeny derived from 25 generations of heat stress displayed a notable heat stress tolerance and higher genetic and epigenetic variations. Hence, our study suggests that the genetic and epigenetic plasticity induced by heat stress in the progeny might result in a phenotypic resilience to adverse environments and potentially trigger the process of microevolution under long-term stress conditions. Further study is required to discover the stable heritable genetic and epigenetic variations associated with heat stress resilience; a topic currently studied in our laboratory. The discovery of stress-associated stable heritable genetic and epigenetic markers could provide a potential source to develop a molecular breeding program for stress-tolerant crops.

Data Availability Statement

The sequencing data files from this study are deposited in the NCBI Sequence Read Archive (NCBI SRA) under BioProject – PRJNA802915 with accession number SRP357958. SRA identifier link: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA802915.

Author Contributions

NY designed and performed the experiments, analyzed the data, and wrote the manuscript. VT and IA performed the experiments and analyzed the data. BB and YI conducted the bioinformatics analysis. IK came up with the concept of the study, analyzed the data, and contributed to writing. All authors contributed to the article and approved the submitted version.

Funding

We thank Natural Sciences and Engineering Research Council of Canada (NSERC) for funding our work.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Prof. Gideon Grafi from Ben-Gurion University of the Negev, Israel, has been acknowledged for scientific input on the initial manuscript draft.

Supplementary Material

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

References

Aina, R., Sgorbati, S., Santagostino, A., Labra, M., Ghiani, A., and Citterio, S. (2004). Specific hypomethylation of DNA is induced by heavy metals in white clover and industrial hemp. Physiol. Plant. 121, 472–480.

Google Scholar

Akalin, A., Franke, V., Vlahoviček, K., Mason, C. E., and Schübeler, D. (2014). Genomation: a toolkit to summarize, annotate and visualize genomic intervals. Bioinformatics 31, 1127–1129. doi: 10.1093/bioinformatics/btu775

PubMed Abstract | CrossRef Full Text | Google Scholar

Akalin, A., Kormaksson, M., Li, S., Garrett-Bakelman, F. E., Figueroa, M. E., Melnick, A., et al. (2012). methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 13:R87. doi: 10.1186/gb-2012-13-10-r87

PubMed Abstract | CrossRef Full Text | Google Scholar

Al Mamun, A. A. M., Lombardo, M. J., Shee, C., Lisewski, A. M., Gonzalez, C., Lin, D., et al. (2012). Identity and function of a large gene network underlying mutagenic repair of DNA breaks. Science 338, 1344–1348. doi: 10.1126/science.1226683

PubMed Abstract | CrossRef Full Text | Google Scholar

Atighi, M. R., Verstraeten, B., De Meyer, T., and Kyndt, T. (2020). Genome-wide DNA hypomethylation shapes nematode pattern-triggered immunity in plants. New Phytol. 227, 545–558. doi: 10.1111/nph.16532

PubMed Abstract | CrossRef Full Text | Google Scholar

Becker, C., Hagmann, J., Müller, J., Koenig, D., Stegle, O., Borgwardt, K., et al. (2011). Spontaneous epigenetic variation in the Arabidopsis thaliana methylome. Nature 480, 245–249. doi: 10.1038/nature10555

PubMed Abstract | CrossRef Full Text | Google Scholar

Belfield, E. J., Brown, C., Ding, Z. J., Chapman, L., Luo, M., Hinde, E., et al. (2020). Thermal stress accelerates Arabidopsis thaliana mutation rate. Genome Res. 31, 40–50. doi: 10.1101/gr.259853.119

PubMed Abstract | CrossRef Full Text | Google Scholar

Bewick, A. J., and Schmitz, R. J. (2017). Gene body DNA methylation in plants. Curr. Opin. Plant Biol. 36, 103–110.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Bilichak, A., Ilnystkyy, Y., Hollunder, J., and Kovalchuk, I. (2012). The progeny of Arabidopsis thaliana plants exposed to salt exhibit changes in DNA methylation, histone modifications and gene expression. PLoS One 7:e30515. doi: 10.1371/journal.pone.0030515

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyko, A., and Kovalchuk, I. (2010). Transgenerational response to stress in Arabidopsis thaliana. Plant Signal. Behav. 5, 995–998.

Google Scholar

Boyko, A., and Kovalchuk, I. (2011). Genome instability and epigenetic modification—heritable responses to environmental stress? Curr. Opin. Plant Biol. 14, 260–266. doi: 10.1016/j.pbi.2011.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyko, A., Kathiria, P., Zemp, F. J., Yao, Y., Pogribny, I., and Kovalchuk, I. (2007). Transgenerational changes in the genome stability and methylation in pathogen-infected plants: (Virus-induced plant genome instability). Nucleic Acids Res. 35, 1714–1725. doi: 10.1093/nar/gkm029

PubMed Abstract | CrossRef Full Text | Google Scholar

Bray, E. A. Bailey-Serres, J., and Weretilnyk, E. (2000). “Responses to abiotic stress,” in Biochemistry & Molecular Biology of Plants, eds W. Gruissem and R. Jones (Rockville, MD: American Society of Plant Physiologists), 1158–1203.

Google Scholar

Chinnusamy, V., Schumaker, K., and Zhu, J. K. (2004). Molecular genetic perspectives on cross-talk and specificity in abiotic stress signalling in plants. J. Exp. Bot. 55, 225–236. doi: 10.1093/jxb/erh005

PubMed Abstract | CrossRef Full Text | Google Scholar

Cingolani, P., Platts, A., Wang, L. L., Coon, M., Nguyen, T., Wang, L., et al. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly 6, 80–92. doi: 10.4161/fly.19695

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Dalakouras, A., and Vlachostergios, D. (2021). Epigenetic approaches for crop breeding: current status and perspectives. J. Exp. Bot. 72, 5356–5371. doi: 10.1093/jxb/erab227

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Duempelmann, L., Skribbe, M., and Bühler, M. (2020). Small RNAs in the transgenerational inheritance of epigenetic information. Trends Genet. 36, 203–214. doi: 10.1016/j.tig.2019.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Eichten, S. R., and Springer, N. M. (2015). Minimal evidence for consistent changes in maize DNA methylation patterns following environmental stress. Front. Plant Sci. 6:308. doi: 10.3389/fpls.2015.00308

PubMed Abstract | CrossRef Full Text | Google Scholar

Eichten, S. R., Schmitz, R. J., and Springer, N. M. (2014). Epigenetics: beyond chromatin modifications and complex genetic regulation. Plant Physiol. 165, 933–947. doi: 10.1104/pp.113.234211

PubMed Abstract | CrossRef Full Text | Google Scholar

Erdmann, R. M., and Picard, C. L. (2020). RNA-directed DNA methylation. PLoS Genet. 16:e1009034. doi: 10.1371/journal.pgen.1009034

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, S., Jacobsen, S. E., and Reik, W. (2010). Epigenetic reprogramming in plant and animal development. Science 330, 622–627. doi: 10.1126/science.1190614

PubMed Abstract | CrossRef Full Text | Google Scholar

Filiault, D. L., Wessinger, C. A., Dinneny, J. R., Lutes, J., Borevitz, J. O., Weigel, D., et al. (2008). Amino acid polymorphisms in Arabidopsis phytochrome B cause differential responses to light. Proc. Natl. Acad. Sci. U.S.A. 105, 3157–3162. doi: 10.1073/pnas.0712174105

PubMed Abstract | CrossRef Full Text | Google Scholar

Filichkin, S. A., Priest, H. D., Givan, S. A., Shen, R., Bryant, D. W., Fox, S. E., et al. (2010). Genome-wide mapping of alternative splicing in Arabidopsis thaliana. Genome Res. 20, 45–58. doi: 10.1101/gr.093302.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganguly, D., Crisp, P. A., Eichten, S. R., and Pogson, B. J. (2017). The Arabidopsis DNA methylome is stable under transgenerational drought stress. Plant Physiol. 175, 1893–1912. doi: 10.1104/pp.17.00744

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, G., Li, J., Li, H., Li, F., Xu, K., Yan, G., et al. (2014). Comparison of the heat stress induced variations in DNA methylation between heat-tolerant and heat-sensitive rapeseed seedlings. Breed. Sci. 64, 125–133. doi: 10.1270/jsbbs.64.125

PubMed Abstract | CrossRef Full Text | Google Scholar

Gill, S. S., Anjum, N. A., Gill, R., Jha, M., and Tuteja, N. (2015). DNA damage and repair in plants under ultraviolet and ionizing radiations. Sci. World J. 2015:250158. doi: 10.1155/2015/250158

PubMed Abstract | CrossRef Full Text | Google Scholar

Gogolev, Y. V., Ahmar, S., Akpinar, B. A., Budak, H., Kiryushkin, A. S., Gorshkov, V. Y., et al. (2021). OMICs, epigenetics, and genome editing techniques for food and nutritional security. Plants 10:1423. doi: 10.3390/plants10071423

PubMed Abstract | CrossRef Full Text | Google Scholar

Groot, M. P., Kooke, R., Knoben, N., Vergeer, P., Keurentjes, J. J., Ouborg, N. J., et al. (2016). Effects of multi-generational stress exposure and offspring environment on the expression and persistence of transgenerational effects in Arabidopsis thaliana. PLoS One 11:e0151566. doi: 10.1371/journal.pone.0151566

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutzat, R., and Scheid, O. M. (2012). Epigenetic responses to stress: triple defense? Curr. Opin. Plant Biol. 15, 568–573. doi: 10.1016/j.pbi.2012.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Guyon-Debast, A., Lécureuil, A., Bonhomme, S., Guerche, P., and Gallois, J. L. (2010). A SNP associated with alternative splicing of RPT5b causes unequal redundancy between RPT5a and RPT5b among Arabidopsis thaliananatural variation. BMC Plant Biol. 10:158. doi: 10.1186/1471-2229-10-158

PubMed Abstract | CrossRef Full Text | Google Scholar

Haque, M. M., Nilsson, E. E., Holder, L. B., and Skinner, M. K. (2016). Genomic Clustering of differential DNA methylated regions (epimutations) associated with the epigenetic transgenerational inheritance of disease and phenotypic variation. BMC Genomics 17:418. doi: 10.1186/s12864-016-2748-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Herman, J. J., and Sultan, S. E. (2011). Adaptive transgenerational plasticity in plants: case studies, mechanisms, and implications for natural populations. Front. Plant Sci. 2:102. doi: 10.3389/fpls.2011.00102

PubMed Abstract | CrossRef Full Text | Google Scholar

Hossain, M. S., Kawakatsu, T., Kim, K. D., Zhang, N., Nguyen, C. T., Khan, S. M., et al. (2017). Divergent cytosine DNA methylation patterns in single-cell, soybean root hairs. New Phytol. 214, 808–819. doi: 10.1111/nph.14421

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, G. T., Ma, S. L., Bai, L. P., Zhang, L., Ma, H., Jia, P., et al. (2012). Signal transduction during cold, salt, and drought stresses in plants. Mol. Biol. Rep. 39, 969–987. doi: 10.1007/s11033-011-0823-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ilnytskyy, Y., Boyko, A., and Kovalchuk, I. (2004). Luciferase-based transgenic recombination assay is more sensitive than β-glucoronidase-based. Mutat. Res. Genet. Toxicol. Environ. Mutagen. 559, 189–197. doi: 10.1016/j.mrgentox.2004.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Jain, M., Moharana, K. C., Shankar, R., Kumari, R., and Garg, R. (2014). Genomewide discovery of DNA polymorphisms in rice cultivars with contrasting drought and salinity stress response and their functional relevance. Plant Biotechnol. J. 12, 253–264. doi: 10.1111/pbi.12133

PubMed Abstract | CrossRef Full Text | Google Scholar

Jang, S. J., Sato, M., Sato, K., Jitsuyama, Y., Fujino, K., Mori, H., et al. (2015). A single-nucleotide polymorphism in an endo-1, 4-β-glucanase gene controls seed coat permeability in soybean. PLoS One 10:e0128527. doi: 10.1371/journal.pone.0128527

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeremias, G., Barbosa, J., Marques, S. M., De Schamphelaere, K. A., Van Nieuwerburgh, F., Deforce, D., et al. (2018). Transgenerational inheritance of DNA hypomethylation in Daphnia magna in response to salinity stress. Environ. Sci. Technol. 52, 10114–10123. doi: 10.1021/acs.est.8b03225

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, C., Mithani, A., Belfield, E. J., Mott, R., Hurst, L. D., and Harberd, N. P. (2014). Environmentally responsive genome-wide accumulation of de novo Arabidopsis thaliana mutations and epimutations. Genome Res. 24, 1821–1829. doi: 10.1101/gr.177659.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Kathiria, P., Sidler, C., Golubov, A., Kalischuk, M., Kawchuk, L. M., and Kovalchuk, I. (2010). Tobacco mosaic virus infection results in an increase in recombination frequency and resistance to viral, bacterial, and fungal pathogens in the progeny of infected tobacco plants. Plant Physiol. 153, 1859–1870. doi: 10.1104/pp.110.157263

PubMed Abstract | CrossRef Full Text | Google Scholar

Kotak, S., Larkindale, J., Lee, U., von Koskull-Döring, P., Vierling, E., and Scharf, K.-D. (2007). Complexity of the heat stress response in plants. Curr. Opin. Plant Biol. 10, 310–316. doi: 10.1016/j.pbi.2007.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Kroymann, J., Donnerhacke, S., Schnabelrauch, D., and Mitchell-Olds, T. (2003). Evolutionary dynamics of an Arabidopsis insect resistance quantitative trait locus. Proc. Natl. Acad. Sci. U.S.A. 100(Suppl. 2), 14587–14592. doi: 10.1073/pnas.1734046100

PubMed Abstract | CrossRef Full Text | Google Scholar

Krueger, F., and Andrews, S. R. (2011). Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics 27, 1571–1572. doi: 10.1093/bioinformatics/btr167

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Law, J. A., and Jacobsen, S. E. (2010). Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat. Rev. Genet. 11, 204–220. doi: 10.1038/nrg2719

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Huang, Q., Sun, M., Zhang, T., Li, H., Chen, B., et al. (2016). Global DNA methylation variations after short-term heat shock treatment in cultured microspores of Brassica napus cv. Topas. Sci. Rep. 6:38401. doi: 10.1038/srep38401

PubMed Abstract | CrossRef Full Text | Google Scholar

Liégard, B., Baillet, V., Etcheverry, M., Joseph, E., Lariagon, C., Lemoine, J., et al. (2019). Quantitative resistance to clubroot infection mediated by transgenerational epigenetic variation in Arabidopsis. New Phytol. 222, 468–479. doi: 10.1111/nph.15579

PubMed Abstract | CrossRef Full Text | Google Scholar

Lind, M. I., and Spagopoulou, F. (2018). Evolutionary consequences ofepigenetic inheritance. Heredity 121, 205–209.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, L., Ou, X., Liu, J., Lin, X., Sheng, L., and Liu, B. (2009). The spaceflight environment can induce transpositional activation of multiple endogenous transposable elements in a genotype-dependent manner in rice. J. Plant Physiol. 166, 2035–2045. doi: 10.1016/j.jplph.2009.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Y., Min, L., Wang, M., Wang, C., Zhao, Y., Li, Y., et al. (2018). Disrupted genome methylation in response to high temperature has distinct affects on microspore abortion and anther indehiscence. Plant Cell 30, 1387–1403. doi: 10.1105/tpc.18.00074

PubMed Abstract | CrossRef Full Text | Google Scholar

Maloof, J. N., Borevitz, J. O., Dabi, T., Lutes, J., Nehring, R. B., Redfern, J. L., et al. (2001). Natural variation in light sensitivity of Arabidopsis. Nat. Genet. 29, 441–446.

Google Scholar

Manikkam, M., Guerrero-Bosagna, C., Tracey, R., Haque, M. M., and Skinner, M. K. (2012). Transgenerational actions of environmental compounds on reproductive disease and identification of epigenetic biomarkers of ancestral exposures. PLoS One 7:e31901. doi: 10.1371/journal.pone.0031901

PubMed Abstract | CrossRef Full Text | Google Scholar

Migicovsky, Z., and Kovalchuk, I. (2013). Changes to DNA methylation and homologous recombination frequency in the progeny of stressed plants. Biochem. Cell Biol. 91, 1–5. doi: 10.1139/bcb-2012-0046

PubMed Abstract | CrossRef Full Text | Google Scholar

Migicovsky, Z., and Kovalchuk, I. (2014). Transgenerational changes in plant physiology and in transposon expression in response to UV-C stress in Arabidopsis thaliana. Plant Signal. Behav. 9:e976490. doi: 10.4161/15592324.2014.976490

PubMed Abstract | CrossRef Full Text | Google Scholar

Migicovsky, Z., Yao, Y., and Kovalchuk, I. (2014). Transgenerational phenotypic and epigenetic changes in response to heat stress in Arabidopsis thaliana. Plant Signal. Behav. 9:e27971. doi: 10.4161/psb.27971

PubMed Abstract | CrossRef Full Text | Google Scholar

Min, L., Li, Y., Hu, Q., Zhu, L., Gao, W., Wu, Y., et al. (2014). Sugar and auxin signaling pathways respond to high-temperature stress during anther development as revealed by transcript profiling analysis in cotton. Plant Physiol. 164, 1293–1308. doi: 10.1104/pp.113.232314

PubMed Abstract | CrossRef Full Text | Google Scholar

Mindrinos, M., Katagiri, F., Yu, G. L., and Ausubel, F. M. (1994). The Arabidopsis thaliana disease resistance gene RPS2 encodes a protein containing a nucleotide-binding site and leucine-rich repeats. Cell 78, 1089–1099. doi: 10.1016/0092-8674(94)90282-8

CrossRef Full Text | Google Scholar

Mirouze, M., and Paszkowski, J. (2011). Epigenetic contribution to stress adaptation in plants. Curr. Opin. Plant Biol. 14, 267–274. doi: 10.1016/j.pbi.2011.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Miryeganeh, M., and Saze, H. (2020). Epigenetic inheritance and plant evolution. Popul. Ecol. 62, 17–27.

Google Scholar

Miura, K., Agetsuma, M., Kitano, H., Yoshimura, A., Matsuoka, M., Jacobsen, S. E., et al. (2009). A metastable DWARF1 epigenetic mutant affecting plant stature in rice. Proc. Natl. Acad. Sci. U.S.A. 106, 11218–11223. doi: 10.1073/pnas.0901942106

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicotra, A. B., Atkin, O. K., Bonser, S. P., Davidson, A. M., Finnegan, E., Mathesius, U., et al. (2010). Plant phenotypic plasticity in a changing climate. Trends Plant Sci. 15, 684–692.

Google Scholar

Paszkowski, J., and Grossniklaus, U. (2011). Selected aspects of transgenerational epigenetic inheritance and resetting in plants. Curr. Opin. Plant Biol. 14, 195–203. doi: 10.1016/j.pbi.2011.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Pecinka, A., Rosa, M., Schikora, A., Berlinger, M., Hirt, H., Luschnig, C., et al. (2009). Transgenerational stress memory is not a general response in Arabidopsis. PLoS One 4:e5202. doi: 10.1371/journal.pone.0005202

PubMed Abstract | CrossRef Full Text | Google Scholar

Perez, M. F., and Lehner, B. (2019). Intergenerational and transgenerational epigenetic inheritance in animals. Nat. Cell Biol. 21, 143–151. doi: 10.1038/s41556-018-0242-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Popova, O. V., Dinh, H. Q., Aufsatz, W., and Jonak, C. (2013). The rddm pathway is required for basal heat tolerance in Arabidopsis. Mol. Plant 6, 396–410. doi: 10.1093/mp/sst023

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, J. R. (2005). Rising temperatures are likely to reduce crop yields. Nature 436:174. doi: 10.1038/436174b

PubMed Abstract | CrossRef Full Text | Google Scholar

Provart, N., and Zhu, T. (2003). A browser-based functional classification SuperViewer for Arabidopsis genomics. Curr. Opin. Plant Biol. 2003, 271–272.

Google Scholar

Rahavi, M. R., Migicovsky, Z. D., Titov, V., and Kovalchuk, I. (2011). Transgenerational adaptation to heavy metal salts in Arabidopsis. Front. Plant Sci. 2:91. doi: 10.3389/fpls.2011.00091

PubMed Abstract | CrossRef Full Text | Google Scholar

Rahavi, S. M., and Kovalchuk, I. (2013). Changes in homologous recombination frequency in Arabidopsis thaliana plants exposed to stress depend on time of exposure during development and on duration of stress exposure. Physiol. Mol. Biol. Plants 19, 479–488. doi: 10.1007/s12298-013-0197-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramírez-Carrasco, G., Martínez-Aguilar, K., and Alvarez-Venegas, R. (2017). Transgenerational defense priming for crop protection against plant pathogens: a hypothesis. Front. Plant Sci. 8:696. doi: 10.3389/fpls.2017.00696

PubMed Abstract | CrossRef Full Text | Google Scholar

Remy, J. J. (2010). Stable inheritance of an acquired behavior in Caenorhabditis elegans. Curr. Biol. 20, R877–R878. doi: 10.1016/j.cub.2010.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Schultz, M. D., Schmitz, R. J., and Ecker, J. R. (2012). ‘Leveling’ the playing field for analyses of single-base resolution DNA methylomes. Trends Genet. 28, 583–585. doi: 10.1016/j.tig.2012.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Shastry, B. S. (2009). “SNPs: impact on gene function and phenotype,” in Single Nucleotide Polymorphisms, ed. A. Komar (Totowa, NJ: Humana Press), 3–22. doi: 10.1007/978-1-60327-411-1_1

CrossRef Full Text | Google Scholar

Shor, E., Fox, C. A., and Broach, J. R. (2013). The yeast environmental stress response regulates mutagenesis induced by proteotoxic stress. PLoS Genet 9:e1003680. doi: 10.1371/journal.pgen.1003680

PubMed Abstract | CrossRef Full Text | Google Scholar

Stroud, H., Greenberg, M. V., Feng, S., Bernatavichute, Y. V., and Jacobsen, S. E. (2013). Comprehensive analysis of silencing mutants reveals complex regulation of the Arabidopsis methylome. Cell 152, 352–364. doi: 10.1016/j.cell.2012.10.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, D.-Y., Kaplan, F., Lee, K.-J., and Guy, C. L. (2003). Acquired tolerance to temperature extremes. Trends Plant Sci. 8, 179–187. doi: 10.1016/S1360-1385(03)00047-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Susan, J. C., Harrison, J., Paul, C. L., and Frommer, M. (1994). High sensitivity mapping of methylated cytosines. Nucleic Acids Res. 22, 2990–2997. doi: 10.1093/nar/22.15.2990

PubMed Abstract | CrossRef Full Text | Google Scholar

Suter, L., and Widmer, A. (2013a). Phenotypic effects of salt and heat stress over three generations in Arabidopsis thaliana. PLoS One 8:e80819. doi: 10.1371/journal.pone.0080819

PubMed Abstract | CrossRef Full Text | Google Scholar

Suter, L., and Widmer, A. (2013b). Environmental heat and salt stress induce transgenerational phenotypic changes in Arabidopsis thaliana. PLoS One 8:e60364. doi: 10.1371/journal.pone.0060364

PubMed Abstract | CrossRef Full Text | Google Scholar

Thiebaut, F., Hemerly, A. S., and Ferreira, P. C. G. (2019). A role for epigenetic regulation in the adaptation and stress responses of non-model plants. Front. Plant Sci. 10:246. doi: 10.3389/fpls.2019.00246

PubMed Abstract | CrossRef Full Text | Google Scholar

Uthup, T. K., Ravindran, M., Bini, K., and Thakurdas, S. (2011). Divergent DNA methylation patterns associated with abiotic stress in Hevea brasiliensis. Mol. Plant 4, 996–1013. doi: 10.1093/mp/ssr039

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidal, R. O., do Nascimento, L. C., Mondego, J. M. C., Pereira, G. A. G., and Carazzolle, M. F. (2012). Identification of SNPs in RNA-seq data of two cultivars of Glycine max (soybean) differing in drought resistance. Genet. Mol. Biol. 35(1 (Suppl)), 331–334. doi: 10.1590/S1415-47572012000200014

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Qin, Q., Sun, F., Wang, Y., Xu, D., Li, Z., et al. (2016). Genome-wide differences in DNA methylation changes in two contrasting rice genotypes in response to drought conditions. Front. Plant Sci. 7:1675. doi: 10.3389/fpls.2016.01675

PubMed Abstract | CrossRef Full Text | Google Scholar

Wassenegger, M., and Dalakouras, A. (2021). Viroids as a tool to study RNA-directed DNA methylation in plants. Cells 10:1187. doi: 10.3390/cells10051187

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, T., Wang, Y., and Liu, J. H. (2020). Comparative transcriptome analysis reveals synergistic and disparate defense pathways in the leaves and roots of trifoliate orange (Poncirus trifoliata) autotetraploids with enhanced salt tolerance. Horticult. Res. 7, 1–14. doi: 10.1038/s41438-020-0311-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wibowo, A., Becker, C., Marconi, G., Durr, J., Price, J., Hagmann, J., et al. (2016). Hyperosmotic stress memory in Arabidopsis is mediated by distinct epigenetically labile sites in the genome and is restricted in the male germline by DNA glycosylase activity. eLife 5:e13546.

Google Scholar

Xu, J., Zhou, S., Gong, X., Song, Y., van Nocker, S., Ma, F., et al. (2018). Single-base methylome analysis reveals dynamic epigenomic differences associated with water deficit in apple. Plant Biotechnol. J. 16, 672–687. doi: 10.1111/pbi.12820

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadav, N. S., Shukla, P. S., Jha, A., Agarwal, P. K., and Jha, B. (2012). The SbSOS1 gene from the extreme halophyte Salicornia brachiata enhances Na+ loading in xylem and confers salt tolerance in transgenic tobacco. BMC plant biology 12:188. doi: 10.1186/1471-2229-12-188

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Sanchez, R., Kundariya, H., Maher, T., Dopp, I., Schwegel, R., et al. (2020). Segregation of an MSH1 RNAi transgene produces heritable non-genetic memory in association with methylome reprogramming. Nat. Commun. 11, 1–17. doi: 10.1038/s41467-020-16036-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, Y., and Kovalchuk, I. (2011). Abiotic stress leads to somatic and heritable changes in homologous recombination frequency, point mutation frequency and microsatellite stability in Arabidopsis plants. Mutat. Res. Fundament. Mol. Mech. Mutagen. 707, 61–66. doi: 10.1016/j.mrfmmm.2010.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, Y., Bilichak, A., Golubov, A., and Kovalchuk, I. (2012). ddm1 plants are sensitive to methyl methane sulfonate and NaCl stresses and are deficient in DNA repair. Plant Cell Rep. 31, 1549–1561. doi: 10.1007/s00299-012-1269-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., Tieman, D. M., Jiao, C., Xu, Y., Chen, K., Fei, Z., et al. (2016). Chilling-induced tomato flavor loss is associated with altered volatile synthesis and transient changes in DNA methylation. Proc. Natl. Acad. Sci. U.S.A. 113, 12580–12585.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., Chen, L., Li, M., Lou, Q., Xia, H., Wang, P., et al. (2013). Transgenerational variations in DNA methylation induced by drought stress in two rice varieties with distinguished difference to drought resistance. PLoS One 8:e80253. doi: 10.1371/journal.pone.0080253

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., Chen, L., Xia, H., Wei, H., Lou, Q., Li, M., et al. (2017). Transgenerational epimutations induced by multi-generation drought imposition mediate rice plant’s adaptation to drought condition. Sci. Rep. 7:39843. doi: 10.1038/srep39843

PubMed Abstract | CrossRef Full Text | Google Scholar