Environment-induced heritable variations are common in Arabidopsis thaliana

Abstract

Parental or ancestral environments can induce heritable phenotypic changes, but whether such environment-induced heritable changes are a common phenomenon remains unexplored. Here, we subject 14 genotypes of Arabidopsis thaliana to 10 different environmental treatments and observe phenotypic and genome-wide gene expression changes over four successive generations. We find that all treatments caused heritable phenotypic and gene expression changes, with a substantial proportion stably transmitted over all observed generations. Intriguingly, the susceptibility of a genotype to environmental inductions could be predicted based on the transposon abundance in the genome. Our study thus challenges the classic view that the environment only participates in the selection of heritable variation and suggests that the environment can play a significant role in generating of heritable variations.

Ancestral environments can induce heritable phenotypic changes, but whether it is a common phenomenon remains unknown. This long-term experiment in Arabidopsis thaliana reveals that environment-induced heritable changes that are common, reproducible, and predictable.

Introduction

Heritable variation is the foundation of biodiversity and adaptation to the changing climate^1,2, ultimately generated through mutations and moulded by forces of environmental selection—a tenet well established in textbooks³. Several pioneering studies suggest that heritable variation is also directly introduced by parental or ancestral environments^{4, 5, 6–7}, causing stably transmitted phenotypic and gene expression changes that potentially enhance offspring resistance to environmental stress^4,7,8. However, whether such environment-induced heritable changes are a common phenomenon remains unexplored, and it is also unknown whether the phenomenon is an occasional exception or a fundamental principle. Therefore, exploring the prevalence of environment-induced heritable changes will challenge the textbook knowledge of the origin of heritable variation and provide critical insights for predicting the adaptive response to global environmental changes.

The phenomenon in which the parental or ancestral environment induces offspring phenotypic changes is conventionally referred to as the parental or transgenerational effect⁴, and the major challenge in understanding the prevalence of this phenomenon could be the scale of the experiments⁹. Expanding the scale of these experiments to include different genetic backgrounds has revealed that such environment-induced changes are genotype-specific^10,11 and predictable based on the climate of the genotype’s origin^12,13, while this predictability decreases when sacrificing the number of genotypes to include various environmental treatments⁹. Nevertheless, these studies span only one or two offspring generations, and it is unclear whether environment-induced changes are heritably stable. A large-scale experiment encompassing multiple genotypes, environments, and generations is necessary to reveal the prevalence and predictability of environment-induced heritable changes.

In this work, to address the prevalence and predictability of environment-induced heritable changes, we subject 14 natural accessions (genotypes) of Arabidopsis thaliana to the Control and ten environmental treatments to establish the ancestral generation. Before this generation, we plant the seeds of these accessions in the control environment for two generations to remove potential influence due to seed collection or storage. Then, we cultivate the offspring of the ancestral generation in the control environment to assess whether those environmental treatments induce phenotypic changes in offspring and whether the induced changes are heritable over four offspring generations (Test I, Fig. 1, Supplementary Fig. 1a, Supplementary Table 1). To explore the reproducibility of environment-induced changes, we plant all four offspring generations together under uniform conditions to eliminate any heterogeneity between generations in the growth conditions (Test II, Fig. 1). Besides performing phenotypic observations, we analyse the genome-wide gene expression of samples from Test II to investigate whether heritable phenotypic changes are predicted by heritable gene expression changes and whether the induced expression changes were nonrandom with respect to environments. We find that all ancestral environmental treatments caused heritable phenotypic and gene expression changes. The induced phenotypic changes are quantitatively reproducible between offspring generations and tests. The susceptibility of a genotype to environmental inductions could be predicted based on the transposon abundance in the genome. Therefore, our study reveals that environment-induced heritable changes are common, reproducible, and predictable, and highlights the significance of the environment in generating heritable variations.

Fig. 1 Design of the experiment. [Images not available. See PDF.]

The experiment subjected 14 genotypes of Arabidopsis thaliana to one control and ten environmental treatments and lasted seven generations. Plants were grown for two generations before launching the “ancestral generation” to eliminate variation due to seed collection or storage. The ten ancestral environmental treatments included jasmonic acid addition (JA), control for JA, nutrient addition, nutrient depletion, leaf removal, drought, low cadmium, high cadmium, low salinity and high salinity. The control and each treatment group consisted of 12 and six replicates, resulting in a total of 1008 lines. To test the ancestral environmental effect on offspring, we established Test I to grow all these lines in the control environment over four successive generations (F1 to F4), following the single-seed descendant approach. To explore the reproducibility and assess the gene expression changes, we established Test II, where we selected ten genotypes and six environments and planted all four offspring generations together in the control environment.

Results

Prevalence of environment-induced heritable phenotypic changes

In the ancestral generation, we subjected A. thaliana genotypes to the control environment and 10 treatment environments: jasmonic acid addition (JA), control for JA, nutrient addition, nutrient depletion, leaf removal, drought, low and high concentrations of cadmium, and low and high concentrations of salinity. The statistics showed that the genotype, environmental treatment and their interaction significantly affected all phenotypes measured (except for the interaction effect on flowering time, Supplementary Table 2). Among the specific treatments, nutrient addition significantly increased aboveground biomass, and jasmonic acid addition slightly enhanced fruit production, representing benign environments (Supplementary Fig. 2a). Other treatments significantly decreased fruit production and thus represented stressful environments (Supplementary Fig. 2a). High concentrations of cadmium and salinity generated more stressful conditions than the corresponding low-concentration treatments, resulting in either smaller plants with fewer fruits or shorter plants with fewer fruits (Supplementary Fig. 2a). Therefore, our experiment provided a gradient of benign and stressful environments to assess their effects on offspring phenotypes.

In the offspring generations, we found little evidence of significant main effects of the ancestral environments on phenotypes but significant genotype-by-ancestral_environment interactions (G × A.E) on all phenotypes measured except for fruit number (Supplementary Fig. 2b, Supplementary Table 2). More importantly, the G × A.E interactions remained significant in the four successive offspring generations. These results demonstrated that the effects of the ancestral environments on offspring phenotypes (i.e., transgenerational effects) depended on the genotypes and persisted for four offspring generations. To understand the stable transmission of such genotype-specific changes caused by ancestral environments, we estimated the effect sizes of changes for different genotypes and environments and plotted them in the order of generations (Fig. 2a). The plot showed considerable variation in the magnitude and direction of the phenotypic changes. In some notable cases, such as the response of genotype Abd-0 to high cadmium and the response of TRE-1 to high salinity, we found that the environment-induced phenotypic changes were stably transmitted over four generations of offspring (Fig. 2a). Therefore, our study demonstrates that environment-induced phenotypic changes are genotype-specific and can be stably transmitted (i.e., are heritable) for at least four offspring generations.

Fig. 2 Prevalence of environment-induced heritable phenotypic changes. [Images not available. See PDF.]

a The effect sizes of environment-induced phenotypic changes measured in Test I. The effect sizes are extracted from linear mixed-effects models (LMM) and two-sided tests are applied. b The occurrence of significant and marginally significant environment-induced changes and this occurrence summed over phenotypes. Environment-induced heritable changes (significant changes in at least two generations) marked with dash lines are found in each environment. Source data including the statistics and exact P  values are provided as a Source Data file.

While several cases of environment-induced heritable changes have been identified, a fundamental question remains: how prevalent are these changes across different environments? To answer this question, we transformed the quantitative estimates of effect sizes into qualitative assessments, defining an environment-induced change as occurring when the effect size significantly deviated from zero. The occurrence data revealed environment-induced phenotypic changes in each genotype and environment (Fig. 2b). Importantly, a substantial proportion (36.8% 10.4%, and 7.5%) of these changes occurred in two, three, and four offspring generations, suggesting that these environment-induced changes were heritable. Furthermore, we observed heritable changes (changes occurring in at least two generations) in all environments and persistent changes (changes occurring in three generations and remaining significant in the fourth generation) in seven environments (Fig. 2b). Consequently, these findings demonstrate that environment-induced heritable changes commonly occur across various environments.

Predictability and reproducibility of environment-induced phenotypic changes

The qualitative data revealed that the genotype, ancestral treatment, offspring generation, and phenotype significantly affected the occurrence of environment-induced phenotypic changes (Fig. 3a). The occurrence frequency estimated on the fitness proxy ( = number of siliques per plant or fruit number) was the lowest among different phenotypes, marginally significantly lower than those of flowering time and plant height. This frequency also varied across different offspring generations, but no consistent decline was observed over successive generations. In general, the magnitude of the ancestral treatment effect did not predict the frequency of occurrence. For instance, the leaf removal treatment, affecting only fruit number, appeared to be much less stressful than the low-salinity treatment (affecting four phenotypes), but the occurrence frequencies associated with these two treatments were similar (Supplementary Fig. 2a, Fig. 3a). Only within the same treatment type was a higher intensity of treatment associated with an increase in the occurrence frequency of induced phenotypic changes (high vs. low concentrations of cadmium/salinity, Fig. 3a).

Fig. 3 Predictability and reproducibility of environment-induced phenotypic changes. [Images not available. See PDF.]

a The occurrence frequency (mean ± SE) of significant environment-induced changes varied between genotypes, ancestral environmental treatments, offspring generations and phenotypes. The occurrence frequencies are estimated from six independent biological replicates (n = 6) per genotype/treatment/generation. Significance is tested using generalized linear mixed-effects models (GLMM) with two-sided chi-square statistics. The P values for multiple comparisons are adjusted with the “single-step” method in the function glht (R package “multcomp”). b The occurrence frequency of different genotypes could not be predicted from the climate at origin or methylation, but from the abundance of transposons in certain superfamilies. c The effect sizes of induced changes show limited correlations with ancestor response. d The effect sizes are quantitatively reproducible between offspring generations. e The effect sizes are quantitatively reproducible between Test I and Test II. Solid circles represent significant effect sizes in either the generation or test, and transparent circles represent nonsignificant effect sizes in both generations and tests. FT Flowering time, PH Plant height, AB Aboveground biomass, RD Rosette diameter, and FN Fruit number. Source data are provided as a Source Data file.

Among the factors affecting the occurrence of environment-induced changes, genotype emerged as the primary explanatory variable. This factor accounted for the largest proportion of the variance in occurrence frequency (Fig. 3a), and removing the three genotypes with the highest frequencies (TRE-1, Abd-0, Sap-0) rendered the effect of ancestral environments on occurrence insignificant (Supplementary Fig. 3a). To explore which factor predicted the susceptibility of a genotype to environmental induction, we regressed the occurrence frequency of genotypes against different explanatory variables (climate of origin, methylation levels, and transposon abundance) and found that the abundance of transposons in the Gypsy (r = −0.840) and MuDR (r = 0.815) superfamilies explained a significant proportion of the variation (Fig. 3b, Supplementary Fig. 3b, Supplementary Table 1). Therefore, these findings suggest that genotype may be a primary predictor of the occurrence of environment-induced changes, and this property of a genotype can be predicted by the number of transposons in its genome.

Although the occurrence of environment-induced changes is predictable to some extent, predicting the effect sizes of these changes proved challenging. Consistent with a previous study⁹, the responses of the ancestors showed little, if any, predictability related to the effect sizes of offspring (Fig. 3c). Nevertheless, it was encouraging to find that the effect sizes of offspring could exceed the range of ancestral responses, as observed for flowering time (Fig. 3c). For this trait, the fixed G × A.E effect also explained a greater proportion of the total variation than G × E (Supplementary Fig. 3c). Compared to ancestor-offspring correlations, the coefficients estimated between offspring generations were significantly greater and consistently significantly deviated from zero, except for fruit number, for which the coefficient reached 0.682 (median = 0.474) (Fig. 3d, Supplementary Fig. 4a, b). The coefficient increased to 0.857 (median = 0.614) when including only significant effect sizes in the analysis and declined substantially and insignificantly deviated from zero for most cases when excluding significant effect sizes from the analysis (Supplementary Fig. 4c). A similar pattern of correlation coefficients was also observed when comparing effect sizes between Test I and Test II (Fig. 3e, Supplementary Fig. 4a, please also see below for Test II). Intriguingly, when we scaled down the experiment to include only one genotype (but multiple treatments) or one treatment (but multiple genotypes), we found that the chances of obtaining significantly correlated effect size correlations between generations declined substantially (Supplementary Fig. 4d). Therefore, our study demonstrates that the effect sizes of environment-induced phenotypic changes are reproducible between offspring generations and tests, and conducting large-scale experiments to obtain abundant significant effect sizes is crucial for achieving this reproducibility.

Despite the findings regarding the predictability and reproducibility of environment-induced changes, some unresolved questions remained: why did some significant phenotypic changes disappear in one generation but reappear in the next (Fig. 2b), and why did the occurrence frequency of induced changes vary between offspring generations (Fig. 3a)? A potential explanation is that the growth conditions differed between offspring generations, as evidenced by significant phenotypic variation between generations in the control group (Supplementary Fig. 2c, Supplementary Table 3). Such variation may also contribute to the inconsistent effect sizes observed between generations (significant G × A.E × F, Supplementary Table 4). To eliminate the influence of variable growth conditions, we conducted a second test of ancestral treatment effects (Test II), where we selected a subset of 10 genotypes and five treatments and grew F1 to F4 plants together under the same conditions (Fig. 1). Consistent with the results of Test I, Test II revealed that the effect of ancestral environments depended on the genotypes (Supplementary Fig. 5a, Supplementary Table 5); that environment-induced heritable changes were common across environments (Supplementary Fig. 6); that the occurrence frequency of these changes varied significantly between genotypes (Supplementary Fig. 7a); and that the variation among genotypes could be predicted by the abundance of several transposon superfamilies, including HAT (r = 0.893, P = 0.041), Mariner (r = 0.873, P = 0.053), and MuDR (r = 0.846, P = 0.071) (Supplementary Fig. 7b, c).

Surprisingly, the repeatability of the effect sizes between generations in Test II did not increase compared to that in Test I (Fig. 3d, Supplementary Fig. 7d). We still found significant phenotypic variations between generations in the control group (Supplementary Fig. 5b, Supplementary Table 6), and the genotype-specific ancestral treatment effects varied significantly between generations for flowering time and aboveground biomass (significant G × A.E × F, Supplementary Table 7). As all plants in Test II were grown under the same conditions, the observed generation effects could only be introduced by different growth conditions of previous generations in Test I. Taken together, these results suggest that both the current conditions and conditions of previous generations can affect the measurements of environment-induced changes, and conducting large-scale experiments under well-controlled conditions is the key to achieving high reproducibility and an in-depth understanding of environment-induced changes.

Environment-induced heritable gene expression changes

In Test II, we conducted transcriptomic analysis of offspring generations (F1-F4) to investigate genome-wide gene expression changes induced by ancestral environments (Fig. 1). The patterns of genome-wide gene expression were primarily determined by genotypes (Supplementary Fig. 8a), with noticeable treatment-specific differentiation observed only in the high-cadmium treatment of Abd-0 and drought treatment of TRE-1 (Supplementary Fig. 8b). The number of differentially expressed genes (DEGs) induced by the ancestral treatment ranged from 146 to 3694 (Fig. 4a), which varied significantly between genotypes but not between environmental treatments or generations (Fig. 4b). Importantly, a significant proportion (15.5%, 3.43%, and 1.54%) of these genes exhibited differential expression in two, three, and four offspring generations (Fig. 4c), indicating that the gene expression changes were heritable. We considered those genes that were differentially expressed for at least two generations to be heritable DEGs, and the number of heritable DEGs ranged from 69 to 2668. Importantly, we identified heritable DEGs in each genotype and environment; heritable DEGs that remained differentially expressed over four offspring generations were also found in each environment (Fig. 4c). The number of heritable DEGs for each genotype and treatment was significantly related to the occurrence frequency of induced phenotypic changes (Supplementary Fig. 9a) and was also significantly explained by the number of transposons in several superfamilies (Fig. 4d, Supplementary Fig. 9b). Therefore, these results demonstrate that environment-induced heritable gene expression changes were common across different environments, and the number of genes showing heritable expression changes correlates with the abundance of transposons in the genomes.

Fig. 4 Environment-induced heritable gene expression changes are found in each ancestral environmental treatment, and the predictability of environment-induced gene expression changes. [Images not available. See PDF.]

a The number of differentially expressed genes (DEGs) for each genotype, ancestral environment treatment and offspring generation. b The number of DEGs (mean ± SE) varied significantly between different genotypes, but not between treatments and generations. The number of DEGs are estimated from three independent biological replicates (n = 3) per genotype/treatment/generation. Significance is tested using linear models with two-sided F statistics. c The number of heritable DEGs (genes significantly differentially expressed in two, three or four offspring generations) for each genotype and treatment. d The number of heritable DEGs for different genotypes could be predicted from the abundance of transposons in several superfamilies. Source data are provided as a Source Data file.

An important unanswered question is whether the induced gene expression changes are directional or unidirectional. If these changes are directional, we expect the DEGs induced by different environments to have distinct functions. To test this hypothesis, we conducted GO and KEGG functional enrichment analyses of heritable DEGs and found that 80.2% (73 out of 91) of enriched GO terms and 58% (27 out of 46) of enriched KEGG pathways were shared between treatments (Fig. 5a, b, Supplementary Fig. 9c, d). The overlapping functions appeared to be primarily associated with the stress response (Fig. 5a) and metabolic pathways (Supplementary Fig. 9c). These results indicated that, contrary to our expectations, heritable DEGs induced by different environmental treatments were functionally similar, which may be associated with these treatments (drought and low and high concentrations of cadmium) generating similar stresses (i.e., inducing hypoxia in the root system).

Fig. 5 Heritable gene expression changes induced by various environments were functionally similar and explained the induced phenotype changes. [Images not available. See PDF.]

a Gene ontology (GO) enrichment results of heritable differentially expressed genes (DEGs) show that induced gene expression changes are functionally similar between environmental treatments and genotypes. Significance was calculated by hypergeometric distribution by R package clusterProfiler. b Venn plot shows the numbers of shared terms between environmental treatments and genotypes. c Venn plot shows the number of heritable DEGs within the shared terms. d The heritable DEGs explained the environment-induced heritable changes in flowering time. Source data are provided as a Source Data file.

Intriguingly, in addition to being shared between different environmental treatments, a substantial proportion of terms/pathways were also shared between genotypes (Fig. 5b, Supplementary Fig. 9d). The shared terms/pathways likely originated from heritable DEGs that were shared between treatments and/or genotypes (Fig. 5b, c, Supplementary Fig. 9d, e). As the heritable expression changes in these genes were repeatedly induced in two or more treatments or genotypes, these genes could be considered susceptible to environmental induction. We hypothesize that the susceptibility of these genes is explained by certain genomic features, such as the abundance of transposons in their promoter and gene body regions. To test this hypothesis, we identified all heritable DEGs shared between genotypes or treatments (Supplementary Fig. 10a) and found that the abundance of members of several superfamilies, including Gypsy, HAT, Mariner, and Helitron, upstream and body of these susceptible genes was significantly higher than that of background genes (Supplementary Fig. 10b). Therefore, these results suggest that environment-induced changes may be nonrandom in the sense of genomic features.

To investigate whether the heritable phenotypic changes were associated with heritable gene expression changes, we focused on the phenotype of flowering time, which has a comparatively well-characterized molecular basis^{14, 15–16}. Among the 79 genes underlying flowering time changes identified in the literature, 16 were differentially expressed in at least one genotype or treatment, and seven were heritable DEGs (Supplementary Fig. 10c). The heritable gene expression changes of six genes co-occurred with heritable flowering time changes, including the early flowering of Abd-0 in response to high cadmium and drought-induced late flowering in TRE-1 in response to drought (Fig. 5d, Supplementary Fig. 6). In Abd-0, ancestral treatment with high cadmium caused downregulation of FLC and MAF4 and upregulation of FT and SOC1, promoting early flowering (Fig. 5d). In TRE-1, ancestal drought treatment caused the upregulation of MAF3 and MAF4 and downregulation of FT and SOC1, which was associated with delayed flowering time (Fig. 5d). FLC is the primary determinant of the vernalization requirement and negatively regulates the expression of FT^17,18. MAFs are paralogues of FLC, which act redundantly with FLC as floral repressors^17,19. Therefore, these results suggest that heritable expression changes in key genes underlie environmentally induced heritable phenotypic changes.

Discussion

Here, we conducted a large-scale multigenerational experiment, and the main finding of our study was that environment-induced heritable changes were caused by each ancestral environmental treatment. These results extend the understanding of the phenomenon of environment-induced heritable changes beyond occasional findings under specific environmental conditions^4,5,12 and suggest that this phenomenon could be a fundamental principle of genetics. Our work challenges the textbook understanding that environments are only involved in the selection of heritable variation³ and suggests that environments are broadly involved in the generation of heritable variation.

The scale of our experiment allowed us to explore the stochasticity of environment-induced heritable changes and reveal genotypes with a higher abundance of transposons that are more susceptible to environmental inductions. Furthermore, our study revealed that genes with more transposons upstream are also prone to environmental induction. These results indicate that transposons play a central role in mediating environment-induced changes, which could occur through epigenetic modifications of transposons (i.e., DNA methylation^20,21 and RNA interference²²) that regulate the expression of downstream genes or through de novo transposon insertions^23,24. Transposons are mobile elements that constitute a significant proportion of eukaryotic genomes²⁵. While these elements were once considered junk DNA, it became evident that they can contribute to the transcriptional binding region²⁶, methylation loci²¹, transcriptional isoform sites²⁷, and even coding sequences²⁸ to the host genomes. Barbara McClintock, who discovered transposons, proposed the idea of “genome shock,” suggesting that stress can activate the movement of transposons, causing insertional mutations and structural variants²⁹. This viewpoint has been increasingly supported by evidence that stress increases transposition rate^{30, 31, 32–33}, with the phenotypic effects of de novo transposon insertions including enhanced resistance of insects to pesticides^{34, 35–36} and decreased sensitivity of plants to abiotic stress^{36, 37, 38–39}. The stress-induced transposition may be mediated by Hsp70 chaperone, a response factor to various biotic and abiotic stressors in plants and animals^40,41, whose inducible expression could disturb the biogenesis of transposon-silencing piRNA⁴². These studies, along with our own, bridge and revive these previous ideas and explanations, suggesting that the environment is not only the selector but can also act as an inducer of genetic variation^36,43.

The finding that environment-induced heritable changes are common across environments also indicates that, besides standing genetic variation within populations^44,45, de novo heritable variations caused by environmental changes can serve as another primary source for rapid adaptation. These de novo heritable variations may arise from the mechanisms discussed earlier, including epigenetic changes (DNA methylation⁴⁶, histone modifications⁴⁷, and small interfering RNAs⁶) or de novo transposon insertion^23,24,48. Some notable examples of rapid evolution based on de novo transposon insertion include the industrial melanism mutation of the peppered moth²³ and the rapid increase in resistance to DDT in Drosophila melanogaster²⁴. While many transposon insertions induced by environmental factors are deleterious, they enable the creation of alleles with significant, positively selected effects^37,49,50. If the environment can act as both a selector and an inducer of heritable variation, then the evolutionary speed constraints imposed by standing genetic variations will be liberated^43,48. Therefore, exploring environmentally induced heritable variations and the molecular mechanisms related to transposon regulation is not only significant for fundamental research in genetics and evolution, but also crucial for assessing the threats of global changes and species responses.

Methods

Plant materials and the experimental design

Arabidopsis thaliana is a model plant species for molecular biology research and is increasingly used to address questions in ecology and evolution⁵¹. It is a self-pollinating annual weed that originates from North Africa⁵² and is broadly distributed across Europe, Asia, and North America. Seeds collected from different geographic sites, known as natural accessions or ecotypes, are preserved at the Arabidopsis Biological Resource Centre (ABRC). As this species is predominantly selfing, different accessions are also taken as different genotypes^10,12. Here, we adopted 14 A. thaliana genotypes ordered from the ABRC, including Abd-0 (from Aberdeen, UK), Ang-0 (Namur, Belgium), Col-0 (Columbia, US), Dja-1 (Chui, Kyrgyzstan), En-1 (Frankfurt, Germany), Fei-0 (Aveiro, Portugal), Kar-1 (Suusamyr, Kyrgyzstan), Ler-0 (Munich, Germany), Olympia-2 (Olympia, Greece), Sap-0 (Hokkaido, Japan), Tol-0 (Ohio, USA), TRE-1 (Marne, France), Tri-0 (Sevilla, Spain) and Ws-2 (Belarus, Belarus) (Supplementary Fig. 1a, Supplementary Table 1). Four seeds per genotype were planted in a common environment for two generations to eliminate the variation due to seed collection or storage. In the second generation, seeds of the four resulting plants were pooled for each genotype to establish the ancestral generation.

To investigate whether environmental changes caused heritable changes in offspring (i.e., transgenerational effects), we subjected the ancestral generation to different environmental treatments and assessed the offspring phenotypic changes over four generations in a common environment. The ancestral environments included the control and ten treatments: jasmonic acid addition (JA), control for JA, nutrient addition, nutrient depletion, leaf removal, drought, low cadmium, high cadmium, low salinity and high salinity (Fig. 1). Each treatment included six replicate plants per genotype, and the control included 12 plants per genotype, resulting in 1008 plants grown or lines established in the ancestral generation. Seeds were collected and reserved for each line, and in the offspring generation, we planted one descendant per line, following the single-seed descendant approach. In Test I, all 1008 lines were planted generation-by-generation in the control environment for four generations (F1-F4), and together 4032 plants were planted. For each offspring generation, seeds were collected and reserved for each line and generation separately.

The growth room contained four shelves, each with four layers (Supplementary Fig. 1b). In each layer, we organized the plants into six rows and 14 columns, and this organization allowed one replicate per treatment/genotype to be planted in two layers. For one replicate of the ancestral generation, treatments were randomly assigned to rows in two layers, and genotypes were distributed randomly among columns in a row (Supplementary Fig. 1c). For the offspring of this replicate, we assigned them to two random layers and distributed the plants randomly among different positions (Supplementary Fig. 1c). All randomization was conducted using the “RAND” function in Excel.

To eliminate the potential influence of growth conditions that varied between generations, we planted all four offspring generations together in the control environment (Test II). Due to limited space in the growth room, Test II included a subset of ten genotypes (Abd-0, Col-0, Ang-0, Dja-1, Fei-0, Ler-0, Sapporo-0, Tol-0, TRE-1 and Tri-0), six different environments (the control, drought, low cadmium, high cadmium, low salinity and high salinity), and thus 360 lines ( = 10 genotypes * 6 environments * 6 replicates). For each line, we regenerated one plant from reserved seeds of the previous generation, and 1440 plants were grown to establish F1-F4. We adopted a completely random design, assigning the 1440 plants randomly into positions among the 16 layers, with six rows and 15 columns each (Supplementary Fig. 1c). Due to the insufficient storage of seeds, especially for the F1 generation, Test II eventually involved 1252 plants. The randomization was conducted using the “RAND” function in Excel.

Growth conditions, phenotype measurements, and sample collection

The experiment was conducted by four graduate students who assisted each other (Supplementary Fig. 1d). Prior to sowing, we placed 10–15 seeds in a 2 mL centrifuge tube, added 1 mL of sterilized distilled water, and kept the tubes in a 4 °C freezer in the dark for 7 days. The soaked seeds were then evenly distributed into pots (5.5 × 5.5 × 6.0 cm) filled with a mixture of 1:1 sterilized soil and vermiculite. The plants were subjected to a photoperiod of 16 h light (7:00 am to 11:00 pm) and 8 h dark (11:00 pm to 7:00 am), with temperatures ranging between 22 and 24 °C. On the tenth day, one A. thaliana seedling was randomly selected and retained. Lines that failed to germinate were transplanted by seedling of a randomly selected replicate from the same genotype and treatment. For the control environment, we watered the plants from the bottom every four days and substituted water with half Hoagland solution (Hope Bio-Technology Co., Ltd) every 16 days. For the jasmonic acid addition (JA) treatment, we sprayed the plants with 100 μmol/L jasmonic acid solution (Sigma-Aldrich) dissolved in 0.5% ethanol every four days, and in the control for JA treatment, we sprayed 0.5% ethanol (sinopharm Chemical Reagent Co.,Ltd) solution to plants. For the nutrient addition treatment, we watered plants with half Hoagland solution every four days and never gave Hoagland solution to plants in the nutrient depletion treatment. For the leaf removal treatment, we removed half of the total number of leaves 30 days after germination. For the drought treatment, plants were watered every eight days. We dissolved CdCl₂ (Shanghai Macklin Biochemical Technology Co., Ltd.) in water or half Hoagland solution to form 500 μmol/L and 1500 μmol/L CdCl₂ solutions and watered the plants with these solutions in low- and high-cadmium treatments. For the low- and high-salinity treatments, we watered the plants with 50 mmol/L NaCl solutions (Sinopharm Chemical Reagent Co.,Ltd.). These treatments started 25 days after sowing and lasted for 32 days. Seeds for the next generation were collected 50–60 days after sowing. In the ancestral generation, 26 out of 1008 (2.58%) lines died before seed production, and their offspring generation was replaced with seeds of a randomly selected replicate from the same genotype and treatment.

We collected phenotypic data from the ancestral and offspring generations in Test I and Test II. Flowering time was checked daily and recorded as the days from sowing to the opening of the first flower. Fruit number was the number of siliques for a plant, used as a proxy for plant fitness. Rosette diameter was measured as the length of the most extended rosette cross centre. Plant height was the height of a plant when pulled straight. Aboveground biomass was the weight of the aboveground parts after drying in an oven at 80 °C for 48 h.

Statistical analysis of phenotypic data

To assess the effects of the environmental treatment and genotype on the ancestral generation, we applied a linear mixed-effect model to each phenotype (flowering time, plant height, aboveground biomass, rosette leaf diameter and fruit number). The model included genotype (G), environment (E), and their interaction (G × E) as fixed effects and the layer (A.Layer) and row (A.Row) as random effects, which was coded as Y_Anc ~ G + E + G:E + (1│A.Layer) + (1│A.Row) (1). Note that the row number was unique and differed between layers.

To assess the effects of ancestral treatment and genotype on offspring generations, we applied a linear mixed-effect model to each phenotype measured in Tests I and II. The model included offspring layer (O. Layer), genotype (G), ancestral environment (A.E), and their interaction (G × A.E) as fixed effects and ancestor layer (A. Layer) and row (A. Row) as random effects and was coded as Y_Off ~ O. Layer + G + A.E + G:A.E + (1│A. Layer) + (1│A. Row) (2). We also integrated all four offspring generations in a model. The model included offspring layer (O. Layer), genotype (G), ancestral environment (A.E), and offspring generation (F) and their two and three interactions as fixed effects and line (Line), ancestral layer (A.Layer) and row (A.Row) as random effects and was coded as Y_Off ~ O.Layer + G + A.E + F + G:A.E + G:F + A.E: F + G:A.E:F + (1│Line) + (1│A.Layer) + (1│A.Row) (3). In Test I, generation was nested in the offspring layer, removing the effect of which also removed that of generation, and Model (3) ~.- F was applied to Test I. In addition, we applied a generalized linear mixed-effect model to the subsets of data, including only the offspring of the control environment. This model included the genotype (G) and offspring generation (F) and their interaction as fixed effects and was coded as Y ~ G + F + G:F+ (1│A. Layer) + (1│A. Row) (4).

These models were formulated with function lmer in the R ver.3.6.2 package “lme4” v.1.1.30⁵³. Statistics of fixed effects were estimated with the Anova function in the “car” package v3.1.1⁵⁴, and those of random effects were estimated with the ranova function in the “lmerTest” package v3.1.3⁵⁵. The effect sizes of ancestral responses to environmental treatments for different genotypes were extracted from model (1), and those of induced changes in offspring generations for different genotypes were extracted from model (2) with the summary function following Groot et al. ¹². For the JA treatment, the effect sizes were estimated by applying the models to a subset of data consisting of JA and JA control, and for other treatments, the effect sizes were estimated by applying the models to the respective treatment and the control. The significance of effect sizes was obtained with the summary function in the “lmerTest” package v3.1.3. According to these significances, we converted the effect sizes to qualitative occurrence data, where significant effect sizes were assigned a value of one, and nonsignificant effect sizes were assigned a value of zero.

To estimate the variance explained by different variables in mixed-effects models (1) and (2), we converted these models into full random effect models and applied the function VarCorr (package “lme4”) to the models. Model (1) was converted to Y_Anc ~ 1 + (1│G) + (1│E) + (1│G × E) + (1│A.Layer) + (1│A.Row), and Model (2) was converted to Y_Off ~ 1 + O.Layer + (1│G) + (1│A.E) + (1│G × A.E) + (1│A.Layer) + (1│A.Row). The estimated variance was visualized with the function geom_scatterpie in the R package “scatterpie” v0.1.8 [https://github.com/GuangchuangYu/scatterpie].

Predictability and reproducibility analysis

To estimate whether the genotype, treatment, generation, and phenotype significantly affected the occurrence of transgenerational effects, we applied generalized linear mixed-effect models to the qualitative occurrence data of Test I and Test II. The model included the genotype, treatment, generation and trait as fixed effects and an additional factor group as a random effect. As effect sizes of different traits were estimated from the same group of six replicate plants per genotype/treatment/generation, we defined the six replicates as one group, which was assigned a group ID. The model was formulated with function glmer in the R package “lme4” v.1.1.30 and was coded as Y ~ G + A.E + F + Trait + (1│group.id) (5). Statistics of fixed effects were estimated with the Anova function in the “car” package. Furthermore, we used the glht function in the package “multcomp”⁵⁶ v1.4.20 to further investigate the differences in occurrence frequency among different ancestral treatments, generations, and phenotypes. The adjusted P values were calculated with the “single-step” method based on the joint normal or t distribution of the linear function.

To explore which factors could predict the genotypic variation in the occurrence of transgenerational effects, we related the occurrence probability of genotypes to climates at the origin sites, the DNA methylation level, and the number of transposons in the genomes. We extracted climate data (Bio1-Bio19) from WorldClim [https://worldclim.org], including monthly climate data between 1901 and 2000. The data were extracted with the coordinates of genotypes from a 30-s spatial resolution ( ~ 1 km²) by the R package “raster” v3.6.20 (Supplementary Table 1). The raw methylation data of 14 genotypes were obtained from published sources⁵⁷. Quality control was performed using Trimmomatic v0.36⁵⁸ with the following parameters: ILLUMINACLIP: TruSeq3-PE.fa:2:30:10” LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:50. We aligned the reads against the TAIR10 reference with bwa-meth⁵⁹ v0.2.5 with the default parameters. The aligned reads with mapping quality higher than 30 were obtained with SAMtools v1.7⁶⁰ with the -q 30 option, and duplicates were marked and removed by Picard v2.18.10 [http://broadinstitute.github.io/picard] and SAMtools (-F 1024). The methylation level of each cytosine for different contexts (CG, CHG, and CHH) was extracted with MethylDackel v0.6.1 [https://github.com/dpryan79/MethylDackel]. We excluded cytosines with a read depth lower than 5, and the average methylation level per cytosine was estimated for each context (Supplementary Table 1).

We extracted the number of transposable elements from published sources⁶¹. The original file contained copy number estimation based on read coverage for the 317 TE families analysed across 211 A. thaliana accessions. We summarized the total number of TEs according to superfamilies and then used log-transformed counts for further analysis. These data were available for seven genotypes, including Abd-0, Ang-0, Col-0, Dja-1, Sap-0, Tol-0, and Ws-2 (Supplementary Table 1). The transposon data of Col-0 were not included in the subsequent analysis because Col-0 is the reference genome and shows biased (much more detailed) annotation of transposons compared to other genotypes. We applied arcsine square root transformation to the occurrence frequency and log transformation to the number of transposons in different superfamilies, regressed the transformed occurrence frequency per genotype against the climate data, average methylation levels and transformed transposon number, and extracted the correlation coefficient and P value from the regression.

To investigate whether the effect sizes of environment-induced changes can be quantitatively predicted from the ancestor’s responses, we calculated the pairwise Pearson correlation coefficients between the ancestor’s responses and the effect sizes of different offspring generations. To assess the reproducibility of effect sizes between offspring generations and tests, we also computed pairwise Pearson correlation coefficients between offspring generations for each test and between the tests. These correlation coefficients were estimated for subsets of data, including significant effects in either generation only or excluding these significant effects. In addition, for Test I, we used subsets of data, including the effect sizes of only one genotype or only one treatment, to calculate correlation coefficients between generations. Visualization was conducted using the ggplot function in the R package “ggplot2” v3.4.2⁶².

Genome-wide gene expression analyses

To assess patterns of gene expression changes caused by ancestral treatments, we collected fresh leaf samples from plants in Test II and conducted transcriptomic analysis. This analysis included four genotypes (Abd-0, Ang-0, Tol-0, TRE-1), four environments (the control, drought, low cadmium, and high cadmium), four generations (F1-F4), each with three replicates, and a total of 144 samples. These genotypes and environments were carefully selected to include genotypes and treatments representing a high/low occurrence probability of transgenerational effects. For instance, Abd-0 and TRE-1 represent genotypes sensitive to environmental induction, whereas Ang-0 and Tol-0 represent genotypes resistant to environmental induction. We sampled two to three rosette leaves per plant, immediately froze the samples in liquid nitrogen, and stored them at −80 °C before RNA extraction. We extracted total RNA with TRIzol reagent (Invitrogen)⁶³. The samples of total RNA were subjected to quality assessment with an Agilent 2100 Bioanalyzer. The transcriptome sequencing library was prepared with the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA). The library was sequenced with the Illumina NovaSeq platform to produce 150 bp paired-end reads. Quantity assessment, library preparation and sequencing were conducted by Novogene (Beijing, China). Sequencing yielded 19.5 to 29.4 million pairs of reads (average: 22.6) per sample (Supplementary Data 1).

To improve the quality of reads, we trimmed and filtered the reads with Trimmomatic v0.36 with the following parameters: “ILLUMINACLIP: TruSeq3-PE.fa:2:30:10” LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:50. Reads with improved quality were aligned to the TAIR10 reference genome using STAR v2.7.9⁶⁴ with the default parameters and a TAIR10 gff-converted gtf-format file (parameter ‘-sjdbGTFfile’) (Supplementary Data 1). The conversion was conducted with gffread v.0.12.7⁶⁵. We extracted the gene expression matrix using featureCounts v1.22.2⁶⁶, with the ‘-p’ option to count the meta-features of genes using the TAIR10 gtf file. Genes were included in subsequent analysis if the sum CPM (counts per million) over all samples for the treatment and control per genotype was greater than one. The retained gene and the expression were subjected to cluster analysis of different samples with the plotPCA function in the “DESeq2” v1.36.0 package⁶⁷, and cluster analysis results were visualized using the R package “ggplot2”. We subjected the retained genes and their expression to differential expression analysis with the R package DESeq2 v1.36.0 (these genes were also used for random sampling, please see below). For each genotype, the difference in gene expression was estimated by comparing the expression of a particular gene between the ancestral treatment and the control. The expression matrix was normalized by function vst in package “DESeq2”. The difference in gene expression was estimated with function DESeq, and the differentially expressed genes (DEGs) were filtered by∣log2FoldChange∣>1.5 and adjusted P < 0.05. A gene was considered to be a heritable DEG if it was significantly differentially expressed for at least two generations. To explore the factors that explained the number of heritable DEGs, we regressed the log-transformed number per genotype and treatment against the arcsine-transformed occurrence frequency per genotype and treatment as well as the log-transformed number of transposons per genotype and extracted the correlation coefficient and P value from the regression.

To investigate the function of DEGs, we focused on genes that were differentially expressed in at least two offspring generations (heritable DEGs) and conducted KEGG and Gene Ontology (GO) enrichment analysis. The enrichment analysis was conducted with functions enrichKEGG and enrichGO in the package “clusterProfiler” v4.4.4⁶⁸. Enrichment analysis was conducted separately for four genotypes and three environments, which resulted in 12 sets of enrichment results (Supplementary Data 2). For the KEGG results, we took the union of significantly enriched pathways of the 12 sets, clustered the pathways according to the maps and visualized the results. For the GO results, we removed duplicated terms for each set with function simplify in the package (Supplementary Data 3). Then, we filtered out terms with less than five genes enriched and selected the top 15 terms with the smallest P values for each set, took their union, and visualized the results. KEGG and GO visualization was conducted with the R package “pheatmap” v1.0.12 [https://github.com/raivokolde/pheatmap]. To facilitate visualization, we characterized the terms into functions of response to stimulus, biosynthetic and metabolic processes, cellular process, signalling, biological regulation, rhythmic process, and developmental process. Then, we extracted the mapped heritable DEGs in pathways and terms and summarized their overlapping patterns between genotypes and environments.

We found 2669 heritable DEGs shared between genotypes or environments, which are considered genes sensitive to environmental induction. We hypothesized that their sensitivity was related to the enrichment of transposons in their regulatory and gene body regions. To test this hypothesis, we randomly sampled 2669 genes from the genes subjected to the differential expression analysis and estimated the sum of the number of transposons per superfamily overlapping the 2 kb upstream and gene body regions. Through random sampling 10000 times, we obtained the null distribution of transposon abundance upstream of genes for different superfamilies. Then, we compared the observed transposon abundance upstream of the 2669 heritable DEGs against this distribution. The P value was estimated as the proportion of randomly sampled values greater than the observed value, and significant P values indicate significant enrichment of transposons upstream of heritable DEGs compared with the null expectation. We adopted the transposon annotation of TAIR10 and performed the analysis using a customized script coded with R package “GenomicRanges” v.1.42.0⁶⁹ [available at https://github.com/Xiaohe-Lin/Heritable_variation/]. Notably, our study adopted Col-0 as the reference, and genotype-specific genome assembly is necessary for a comprehensive understanding of this nonrandom phenomenon.

To investigate whether the heritable phenotypic changes were caused by heritable gene expression changes, we focused on flowering time and identified 79 flowering time-related genes from the literature^{14, 15–16} (Supplementary Data 4). We visualized the expression patterns of a subset of genes that were DEGs for at least one generation using the R package “pheatmap” v1.0.12.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Acknowledgements

We thank Wenjing Hong, Gang Wang, and Rong Wang for their constructive suggestions. We thank Ming Zhou, Xianwen Ma, Jiayi Ding, and Xinwei Wang for their help in the growth room and Xiuxiu Wang and Haidong Qu for their technical assistance. This study was supported by the National Natural Science Foundation of China (32371558 and 32071485 to Y.-Y.Z.) and the Fundamental Research Funds for the Central Universities of China (20720210075 to Y.-Y.Z.).

Author contributions

Y.-Y.Z., O.B., and V.L. conceived the idea of this project; Y.-Y.Z. J.J.Y. and X.L designed the experiments; Y-Y.Z. supervised and Q.Q.L. co-supervised the project; J.J.Y., Y.W., J.Y, and X.L. carried out the experiments; X.L, JJ.Y., and Y.W. analysed the data; X.L wrote the first draft of this manuscript, and all authors contributed to the revision of the manuscript.

Peer review

Peer review information

Nature Communications thanks Igor Kovalchuk and Sergio Pimpinelli for their contribution to the peer review of this work. A peer review file is available.

Data availability

The phenotypic data of the ancestral generation, Test I and Test II are available on GitHub [https://github.com/Xiaohe-Lin/Heritable_variation]. Raw RNA-seq data are available from the NCBI SRA database with accession ID PRJNA997595. Source data are provided as a Source Data file. The geographical and climatic data sourced from WorldClim [https://worldclim.org], and information of DNA methylation and transposon abundance were obtained from published literature (Kawakatsu. et al. Cell, 2016, 166:492-505; Quadrana. et al. eLife, 2016, 5:e15716). The Genome, gene and Transposon annotations of Arabidopsis thaliana are from TAIR website [https://www.arabidopsis.org].  are provided with this paper.

Code availability

Codes for data processing are available on GitHub [https://github.com/Xiaohe-Lin/Heritable_variation] and Zenodo [https://zenodo.org/records/11020990]

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41467-024-49024-3.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Hoffmann, AA; Sgrò, CM. Climate change and evolutionary adaptation. Nature; 2011; 470, pp. 479-485.2011Natur.470.479H1:CAS:528:DC%2BC3MXisVags7o%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21350480][DOI: https://dx.doi.org/10.1038/nature09670]

2. Scheffers, BR et al. The broad footprint of climate change from genes to biomes to people. Science; 2016; 354, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27846577][DOI: https://dx.doi.org/10.1126/science.aaf7671] aaf7671.

3. Futuyma, D. J. & Kirkpatrick, M. Evolution. (Oxford University Press, 2017).

4. Agrawal, AA; Laforsch, C; Tollrian, R. Transgenerational induction of defences in animals and plants. Nature; 1999; 401, pp. 60-63.1999Natur.401..60A1:CAS:528:DyaK1MXlvVGmsrw%3D [DOI: https://dx.doi.org/10.1038/43425]

5. Klosin, A; Casas, E; Hidalgo-Carcedo, C; Vavouri, T; Lehner, B. Transgenerational transmission of environmental information in C. elegans. Science; 2017; 356, pp. 320-323.2017Sci..356.320K1:CAS:528:DC%2BC2sXmtFCgurg%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28428426][DOI: https://dx.doi.org/10.1126/science.aah6412]

6. Rechavi, O et al. Starvation-induced transgenerational inheritance of small RNAs in C. elegans. Cell; 2014; 158, pp. 277-287.1:CAS:528:DC%2BC2cXhtFChsrzL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25018105][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4377509][DOI: https://dx.doi.org/10.1016/j.cell.2014.06.020]

7. Yin, J; Zhou, M; Lin, Z; Li, QQ; Zhang, Y-Y. Transgenerational effects benefit offspring across diverse environments: a meta-analysis in plants and animals. Ecol. Lett.; 2019; 22, pp. 1976-1986. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31436014][DOI: https://dx.doi.org/10.1111/ele.13373]

8. Galloway, LF. Parental environmental effects on life history in the herbaceous plant Campanula americana. Ecology; 2001; 82, pp. 2781-2789. [DOI: https://dx.doi.org/10.2307/2679960]

9. Latzel, V et al. Parental environmental effects are common and strong, but unpredictable, in Arabidopsis thaliana. New Phytol.; 2023; 237, pp. 1014-1023.1:CAS:528:DC%2BB38XjtVyls77J [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36319609][DOI: https://dx.doi.org/10.1111/nph.18591]

10. Alvarez, M; Bleich, A; Donohue, K. Genotypic variation in the persistence of transgenerational responses to seasonal cues. Evolution; 2020; 74, pp. 2265-2280.1:CAS:528:DC%2BB3cXitV2rsLrE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32383475][DOI: https://dx.doi.org/10.1111/evo.13996]

11. Colicchio, J. Transgenerational effects alter plant defence and resistance in nature. J. Evol. Biol.; 2017; 30, pp. 664-680.1:STN:280:DC%2BC1c7nt12jsw%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28102915][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5382043][DOI: https://dx.doi.org/10.1111/jeb.13042]

12. Groot, MP et al. Transgenerational effects of mild heat in Arabidopsis thaliana show strong genotype specificity that is explained by climate at origin. New Phytol.; 2017; 215, pp. 1221-1234. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28590553][DOI: https://dx.doi.org/10.1111/nph.14642]

13. Munch, SB et al. A latitudinal gradient in thermal transgenerational plasticity and a test of theory. Proc. Royal Soc. B; 2021; 288, 20210797. [DOI: https://dx.doi.org/10.1098/rspb.2021.0797]

14. Gaudinier, A; Blackman, BK. Evolutionary processes from the perspective of flowering time diversity. New Phytol.; 2020; 225, pp. 1883-1898. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31536639][DOI: https://dx.doi.org/10.1111/nph.16205]

15. He, L et al. DNA methylation-free Arabidopsis reveals crucial roles of DNA methylation in regulating gene expression and development. Nat. Commun.; 2022; 13, 2022NatCo.13.1335H1:CAS:528:DC%2BB38XntVamtLY%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35288562][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8921224][DOI: https://dx.doi.org/10.1038/s41467-022-28940-2] 1335.

16. Roux, F; Touzet, P; Cuguen, J; Le Corre, V. How to be early flowering: an evolutionary perspective. Trends Plant Sci.; 2006; 11, pp. 375-381.1:CAS:528:DC%2BD28Xns1Kjurw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16843035][DOI: https://dx.doi.org/10.1016/j.tplants.2006.06.006]

17. Kim, DH; Doyle, MR; Sung, S; Amasino, RM. Vernalization: winter and the timing of flowering in plants. Annu. Rev. Cell Dev. Biol.; 2009; 25, pp. 277-299.1:CAS:528:DC%2BD1MXhsVKisbvN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19575660][DOI: https://dx.doi.org/10.1146/annurev.cellbio.042308.113411]

18. Searle, I et al. The transcription factor FLC confers a flowering response to vernalization by repressing meristem competence and systemic signaling in Arabidopsis. Genes Dev.; 2006; 20, pp. 898-912.1:CAS:528:DC%2BD28XjsFGjs7o%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16600915][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1472290][DOI: https://dx.doi.org/10.1101/gad.373506]

19. Ratcliffe, OJ; Nadzan, GC; Reuber, TL; Riechmann, JL. Regulation of flowering in Arabidopsis by an FLC homologue. Plant Physiol.; 2001; 126, pp. 122-132.1:CAS:528:DC%2BD3MXjslWjsLo%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11351076][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC102287][DOI: https://dx.doi.org/10.1104/pp.126.1.122]

20. Morgan, HD; Sutherland, HGE; Martin, DIK; Whitelaw, E. Epigenetic inheritance at the agouti locus in the mouse. Nat. Genet.; 1999; 23, pp. 314-318.1:CAS:528:DyaK1MXnt1Gns7w%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10545949][DOI: https://dx.doi.org/10.1038/15490]

21. Ong-Abdullah, M et al. Loss of Karma transposon methylation underlies the mantled somaclonal variant of oil palm. Nature; 2015; 525, pp. 533-537.2015Natur.525.533O1:CAS:528:DC%2BC2MXhsVOls7zM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26352475][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4857894][DOI: https://dx.doi.org/10.1038/nature15365]

22. Ito, H et al. An siRNA pathway prevents transgenerational retrotransposition in plants subjected to stress. Nature; 2011; 472, pp. 115-119.2011Natur.472.115I1:CAS:528:DC%2BC3MXjtVGmt74%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21399627][DOI: https://dx.doi.org/10.1038/nature09861]

23. Van’t Hof, AE et al. The industrial melanism mutation in British peppered moths is a transposable element. Nature; 2016; 534, pp. 102-105.2016Natur.534.102H [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27251284][DOI: https://dx.doi.org/10.1038/nature17951]

24. Schmidt, JM et al. Copy number variation and transposable elements feature in recent, ongoing adaptation at the Cyp6g1 locus. PLoS Genet.; 2010; 6, e1000998. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20585622][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2891717][DOI: https://dx.doi.org/10.1371/journal.pgen.1000998]

25. Bourque, G et al. Ten things you should know about transposable elements. Genome Biol.; 2018; 19, 1:CAS:528:DC%2BC1MXhtVKksbrK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30454069][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6240941][DOI: https://dx.doi.org/10.1186/s13059-018-1577-z] 199.

26. Modzelewski, AJ et al. A mouse-specific retrotransposon drives a conserved Cdk2ap1 isoform essential for development. Cell; 2021; 184, pp. 5541-5558.e5522.1:CAS:528:DC%2BB3MXit1alu7nE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34644528][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8787082][DOI: https://dx.doi.org/10.1016/j.cell.2021.09.021]

27. Kapusta, A et al. Transposable elements are major contributors to the origin, fiversification, and regulation of vertebrate long noncoding RNAs. PLoS Genet.; 2013; 9, e1003470.1:CAS:528:DC%2BC3sXnt1aqs7o%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23637635][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3636048][DOI: https://dx.doi.org/10.1371/journal.pgen.1003470]

28. Pastuzyn, ED et al. The neuronal gene Arc encodes a repurposed retrotransposon Gag protein that mediates intercellular RNA transfer. Cell; 2018; 172, pp. 275-288.e218.1:CAS:528:DC%2BC1cXhtVegtLo%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29328916][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5884693][DOI: https://dx.doi.org/10.1016/j.cell.2017.12.024]

29. McClintock, B. The significance of responses of the genome to challenge. Science; 1984; 226, pp. 792-801.1984Sci..226.792M1:STN:280:DC%2BD2M7gvV2juw%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15739260][DOI: https://dx.doi.org/10.1126/science.15739260]

30. Baduel, P et al. Genetic and environmental modulation of transposition shapes the evolutionary potential of Arabidopsis thaliana. Genome Biol.; 2021; 22, 1:CAS:528:DC%2BB3MXht1Krt7jO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33957946][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8101250][DOI: https://dx.doi.org/10.1186/s13059-021-02348-5] 138.

31. Ho, EKH et al. Engines of change: transposable element mutation rates are high and variable within Daphnia magna. PLoS Genet.; 2021; 17, e1009827.1:CAS:528:DC%2BB3MXisFCjurrP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34723969][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8594854][DOI: https://dx.doi.org/10.1371/journal.pgen.1009827]

32. Ratner, VA; Zabanov, SA; Kolesnikova, OV; Vasilyeva, LA. Induction of the mobile genetic element Dm-412 transpositions in the Drosophila genome by heat shock treatment. Proc. Natl. Acad. Sci. USA.; 1992; 89, pp. 5650-5654.1992PNAS..89.5650R1:CAS:528:DyaK38Xks1Kqtrg%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/1319068][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC49350][DOI: https://dx.doi.org/10.1073/pnas.89.12.5650]

33. Strand, DJ; McDonald, JF. Copia is transcriptionally responsive to environmental stress. Nucleic Acids Res.; 1985; 13, pp. 4401-4410.1:CAS:528:DyaL2MXkvF2ltr4%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/2409535][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC321795][DOI: https://dx.doi.org/10.1093/nar/13.12.4401]

34. Aminetzach, YT; Macpherson, JM; Petrov, DA. Pesticide resistance via transposition-mediated adaptive gene truncation in Drosophila. Science; 2005; 309, pp. 764-767.2005Sci..309.764A1:CAS:528:DC%2BD2MXmsFWmsL4%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16051794][DOI: https://dx.doi.org/10.1126/science.1112699]

35. Chung, H et al. Cis-regulatory elements in the Accord retrotransposon result in tissue-specific expression of the Drosophila melanogaster insecticide resistance gene Cyp6g1. Genetics; 2007; 175, pp. 1071-1077.1:CAS:528:DC%2BD2sXltFKhsrw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17179088][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1840086][DOI: https://dx.doi.org/10.1534/genetics.106.066597]

36. Pimpinelli, S; Piacentini, L. Environmental change and the evolution of genomes: transposable elements as translators of phenotypic plasticity into genotypic variability. Funct. Ecol.; 2020; 34, pp. 428-441. [DOI: https://dx.doi.org/10.1111/1365-2435.13497]

37. Kanazawa, A; Liu, B; Kong, F; Arase, S; Abe, J. Adaptive evolution involving gene duplication and insertion of a novel Ty1/copia-like retrotransposon in soybean. J. Mol. Evol.; 2009; 69, pp. 164-175.2009JMolE.69.164K1:CAS:528:DC%2BD1MXpsFelsrs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19629571][DOI: https://dx.doi.org/10.1007/s00239-009-9262-1]

38. Liu, B et al. Genetic redundancy in soybean photoresponses associated with duplication of the phytochrome A gene. Genetics; 2008; 180, pp. 995-1007.1:CAS:528:DC%2BD1cXhtlWnurrE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18780733][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2567397][DOI: https://dx.doi.org/10.1534/genetics.108.092742]

39. Thieme, M et al. Experimentally heat-induced transposition increases drought tolerance in Arabidopsis thaliana. New Phytol.; 2022; 236, pp. 182-194.1:CAS:528:DC%2BB38Xit1ens7fJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35715973][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9544478][DOI: https://dx.doi.org/10.1111/nph.18322]

40. Yu, A et al. Roles of Hsp70s in stress responses of microorganisms, plants, and animals. BioMed Res. Int.; 2015; 2015, 510319. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26649306][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4663327][DOI: https://dx.doi.org/10.1155/2015/510319]

41. Berka, M; Kopecká, R; Berková, V; Brzobohatý, B; Černý, M. Regulation of heat shock proteins 70 and their role in plant immunity. J. Exp. Bot.; 2022; 73, pp. 1894-1909.1:CAS:528:DC%2BB38XisFGltbfE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35022724][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8982422][DOI: https://dx.doi.org/10.1093/jxb/erab549]

42. Cappucci, U et al. The Hsp70 chaperone is a major player in stress-induced transposable element activation. Proc. Natl. Acad. Sci. USA.; 2019; 116, pp. 17943-17950.2019PNAS.11617943C1:CAS:528:DC%2BC1MXhslWktLbF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31399546][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6731680][DOI: https://dx.doi.org/10.1073/pnas.1903936116]

43. Piacentini, L et al. Transposons, environmental changes, and heritable induced phenotypic variability. Chromosoma; 2014; 123, pp. 345-354. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24752783][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4107273][DOI: https://dx.doi.org/10.1007/s00412-014-0464-y]

44. Barrett, RDH; Schluter, D. Adaptation from standing genetic variation. Trends Ecol. Evol.; 2008; 23, pp. 38-44. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18006185][DOI: https://dx.doi.org/10.1016/j.tree.2007.09.008]

45. Bitter, MC; Kapsenberg, L; Gattuso, JP; Pfister, CA. Standing genetic variation fuels rapid adaptation to ocean acidification. Nat. Commun.; 2019; 10, 2019NatCo.10.5821B1:CAS:528:DC%2BC1MXisVyrtrbO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31862880][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6925106][DOI: https://dx.doi.org/10.1038/s41467-019-13767-1] 5821.

46. Weaver, ICG et al. Epigenetic programming by maternal behavior. Nat. Neurosci.; 2004; 7, pp. 847-854.1:CAS:528:DC%2BD2cXmtVamt7c%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15220929][DOI: https://dx.doi.org/10.1038/nn1276]

47. Wan, Q-L et al. Histone H3K4me3 modification is a transgenerational epigenetic signal for lipid metabolism in Caenorhabditis elegans. Nat. Commun.; 2022; 13, 2022NatCo.13.768W1:CAS:528:DC%2BB38Xjt1Cmu7k%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35140229][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8828817][DOI: https://dx.doi.org/10.1038/s41467-022-28469-4] 768.

48. Fanti, L; Piacentini, L; Cappucci, U; Casale, AM; Pimpinelli, S. Canalization by selection of de Novo induced mutations. Genetics; 2017; 206, pp. 1995-2006.1:CAS:528:DC%2BC1cXitVShtbfE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28576865][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5560803][DOI: https://dx.doi.org/10.1534/genetics.117.201079]

49. Feiner, N. Accumulation of transposable elements in Hox gene clusters during adaptive radiation of Anolis lizards. Proc. Biol. Sci. 283, 20161555 (2016).

50. Schrader, L et al. Transposable element islands facilitate adaptation to novel environments in an invasive species. Nat. Commun.; 2014; 5, 2014NatCo..5.5495S1:CAS:528:DC%2BC2MXjsVagsr8%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25510865][DOI: https://dx.doi.org/10.1038/ncomms6495] 5495.

51. Pigliucci, M. Ecology and evolutionary biology of Arabidopsis. Arabidopsis Book; 2002; 1, e0003. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22303188][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3243336][DOI: https://dx.doi.org/10.1199/tab.0003]

52. Durvasula, A et al. African genomes illuminate the early history and transition to selfing in Arabidopsis thaliana. Proc. Natl. Acad. Sci. USA.; 2017; 114, pp. 5213-5218.2017PNAS.114.5213D1:CAS:528:DC%2BC2sXntVOrsr8%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28473417][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5441814][DOI: https://dx.doi.org/10.1073/pnas.1616736114]

53. Bates, D; Mächler, M; Bolker, B; Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw.; 2015; 67, pp. 1-48. [DOI: https://dx.doi.org/10.18637/jss.v067.i01]

54. Fox, J. W., S. An R companion to applied regression. (SAGE, 2019).

55. Kuznetsova, A; Brockhoff, PB; Christensen, RHB. lmerTest package: tests in linear mixed effects models. J. Stat. Softw.; 2017; 82, pp. 1-26. [DOI: https://dx.doi.org/10.18637/jss.v082.i13]

56. Hothorn, T; Bretz, F; Westfall, P. Simultaneous inference in general parametric models. Biom. J.; 2008; 50, pp. 346-363.2521547 [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18481363][DOI: https://dx.doi.org/10.1002/bimj.200810425]

57. Kawakatsu, T et al. Epigenomic diversity in a global collection of Arabidopsis thaliana accessions. Cell; 2016; 166, pp. 492-505.1:CAS:528:DC%2BC28XhtFynu7bO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27419873][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5172462][DOI: https://dx.doi.org/10.1016/j.cell.2016.06.044]

58. Bolger, AM; Lohse, M; Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics; 2014; 30, pp. 2114-2120.1:CAS:528:DC%2BC2cXht1Sqt7nP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24695404][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4103590][DOI: https://dx.doi.org/10.1093/bioinformatics/btu170]

59. Pedersen, B. S., Eyring, K. R., De, S., Yang, I. V. & Schwartz, D. A. Fast and accurate alignment of long bisulfite-seq reads. Preprint at arXivhttps://doi.org/10.48550/arXiv.1401.1129 (2014).

60. Li, H et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics; 2009; 25, pp. 2078-2079. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19505943][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2723002][DOI: https://dx.doi.org/10.1093/bioinformatics/btp352]

61. Quadrana, L et al. The Arabidopsis thaliana mobilome and its impact at the species level. eLife; 2016; 5, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27258693][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4917339][DOI: https://dx.doi.org/10.7554/eLife.15716] e15716.

62. Hadley, W. ggplot2: elegant graphics for data analysis. (Springer Cham, 2016).

63. Rio, DC; Ares, M, Jr.; Hannon, GJ; Nilsen, TW. Purification of RNA using TRIzol (TRI reagent). Cold Spring Harb. Protoc.; 2010; 2010, pdb.prot5439. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20516177][DOI: https://dx.doi.org/10.1101/pdb.prot5439]

64. Dobin, A et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics; 2013; 29, pp. 15-21.1:CAS:528:DC%2BC38XhvV2gsbnF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23104886][DOI: https://dx.doi.org/10.1093/bioinformatics/bts635]

65. Pertea, G; Pertea, M. GFF Utilities: GffRead and GffCompare. F1000Res.; 2020; 9, 304. [DOI: https://dx.doi.org/10.12688/f1000research.23297.1] 304

66. Liao, Y; Smyth, GK; Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics; 2014; 30, pp. 923-930.1:CAS:528:DC%2BC2cXltFGqu7c%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24227677][DOI: https://dx.doi.org/10.1093/bioinformatics/btt656]

67. Love, MI; Huber, W; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.; 2014; 15, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25516281][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4302049][DOI: https://dx.doi.org/10.1186/s13059-014-0550-8] 550.

68. Yu, GC; Wang, LG; Han, YY; He, QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics; 2012; 16, pp. 284-287.1:CAS:528:DC%2BC38XmsFarsLw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22455463][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339379][DOI: https://dx.doi.org/10.1089/omi.2011.0118]

69. Lawrence, M et al. Software for computing and annotating genomic ranges. PLoS Comput. Biol.; 2013; 9, e1003118.1:CAS:528:DC%2BC3sXhsVensr%2FK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23950696][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3738458][DOI: https://dx.doi.org/10.1371/journal.pcbi.1003118]