Gengyun Li1,2
Yuewan Zhao2
Fei Liu2Minnuo Shi2
Yabin Guan2
Ticao Zhang3
Fangqing Zhao2,4*
Qin Qiao1,5*
Yupeng Geng2*- 1College of Horticulture and Landscape, Yunnan Agricultural University, Kunming, China
- 2Yunnan Key Laboratory of Plant Reproductive Adaptation and Evolutionary Ecology and Centre for Invasion Biology, Institute of Biodiversity, School of Ecology and Environmental Science, Yunnan University, Kunming, China
- 3College of Chinese Material Medica, Yunnan University of Chinese Medicine, Kunming, China
- 4Beijing Institute of Life Sciences, Chinese Academy of Sciences, Beijing, China
- 5School of Agriculture, Yunnan University, Kunming, China
Transgenerational plasticity (TGP) occurs when maternal environments influence the expression of traits in offspring, and in some cases may increase fitness of offspring and have evolutionary significance. However, little is known about the extent of maternal environment influence on gene expression of offspring, and its relationship with trait variations across generations. In this study, we examined TGP in the traits and gene expression of field pennycress (Thlaspi arvense) in response to cadmium (Cd) stress. In the first generation, along with the increase of soil Cd concentration, the total biomass, individual height, and number of seeds significantly decreased, whereas time to flowering, superoxide dismutase (SOD) activity, and content of reduced glutathione significantly increased. Among these traits, only SOD activity showed a significant effect of TGP; the offspring of Cd-treated individuals maintained high SOD activity in the absence of Cd stress. According to the results of RNA sequencing and bioinformatic analysis, 10,028 transcripts were identified as Cd-responsive genes. Among them, only 401 were identified as transcriptional memory genes (TMGs) that maintained the same expression pattern under normal conditions in the second generation as in Cd-treated parents in the first generation. These genes mainly participated in Cd tolerance-related processes such as response to oxidative stress, cell wall biogenesis, and the abscisic acid signaling pathways. The results of weighted correlation network analysis showed that modules correlated with SOD activity recruited more TMGs than modules correlated with other traits. The SOD-coding gene CSD2 was found in one of the modules correlated with SOD activity. Furthermore, several TMGs co-expressed with CSD2 were hub genes that were highly connected to other nodes and critical to the network’s topology; therefore, recruitment of TMGs in offspring was potentially related to TGP. These findings indicated that, across generations, transcriptional memory of gene expression played an important role in TGP. Moreover, these results provided new insights into the trait evolution processes mediated by phenotypic plasticity.
Introduction
Phenotypic plasticity is the ability of organisms to change their phenotypes in response to environmental stimuli (Schlichting, 1986; Sultan, 2000; Pigliucci, 2005). Within-generation plasticity (WGP) allows organisms to match their traits to the local environment conditions (Colicchio and Herman, 2020). Moreover, plasticity may also occur between generations, termed as maternal effects or transgenerational plasticity (TGP), which accounts for the influences of parental environment (PE) on offspring development and phenotype (Agrawal et al., 1999; Galloway and Etterson, 2007; Bell and Hellmann, 2019). As the ecological and evolutionary significance of the phenotypic plasticity was realized and gradually elucidated (Kelly et al., 2012; Forsman, 2015), more attention has been paid to adaptive TGP because it can transfer potential ecologically and evolutionarily meaningful variation from parent to offspring in organisms (Herman and Sultan, 2011). When growth conditions are stable, the sessile lifestyle of plants and spatially limited dispersal of propagules usually make the PE a good predictor of the offspring’s environment (Sandner et al., 2018; Colicchio and Herman, 2020). As a result, the offspring’s fitness can be enhanced through adaptive TGP by shaping their phenotypes much earlier during development (Galloway and Etterson, 2007; Vu et al., 2015; Colicchio and Herman, 2020). Several studies demonstrated that TGP provides a pathway of adaptation that can influence the dynamics of selection and promote ecological breadth by maintaining populations under stressful environments (Sultan et al., 2009; Dyer et al., 2010; Herman et al., 2012). Therefore, TGP is widely believed to play an important role in plant adaptation and evolution in a challenging world (Lachmann and Jablonka, 1996; Kuijper et al., 2014; Leimar and McNamara, 2015; Bonduriansky, 2021).
Although increasing evidence has shown that TGP is common in nature (reviewed in Herman and Sultan, 2011; Donelson et al., 2017; Bonduriansky, 2021), the mechanisms responsible for the effects of PEs on offspring phenotypes remain elusive (Salinas et al., 2013; Herman and Sultan, 2016). Typically, phenotype is a result of genotype-by-environment interaction (G × E) (Des Marais et al., 2013), and the ability of phenotypic plasticity may be influenced by epigenetic regulation (Geng et al., 2013; Kooke et al., 2015). The cross-generational transmission of environmental cues may be mediated by epigenetic mechanisms and/or molecules that are packaged into seeds and transmitted from parents to offspring (Herman and Sultan, 2011). Epigenetic mechanisms, including small RNAs (sRNAs), DNA methylation, and histone modification, could regulate gene expression in offspring at transcriptional or post-transcriptional levels (Lang-Mladek et al., 2010; Pieterse, 2012; Quadrana and Colot, 2016). For example, DNA methylation induced by maternal environment for a part of sites in genome can be inherited stably by successive generations through sexual and/or asexual reproduction (Quadrana and Colot, 2016; Shi et al., 2019), thus affecting gene expression over multiple generations (Herman and Sultan, 2011; Fitz-James and Cavalli, 2022). Moreover, other maternally derived molecules, including mRNAs, proteins, phytohormones, and some primary and secondary metabolites, may also have effects on offspring traits (Verhoeven et al., 2010; Herman and Sultan, 2011). Although these mechanisms can independently or jointly influence the occurrence of TGP, all of them are ultimately functionalized by regulation of gene expression during offspring development (Bonduriansky, 2021). As a result, expression or transcriptional memory of selective genes could be one of the most important determinants of plasticity of gene expression. As gene expression connects environmental cues and trait expression, detection of changes in gene expression across generations can be used to identify biological processes and pathways underlying TGP, and to elucidate the molecular basis of TGP (Laitinen and Nikoloski, 2019; Sommer, 2020).
Recent studies mainly focused on inheritance of phenotypic trait changes and epigenetic memories in TGP (e.g., Herman and Sultan, 2016; Colicchio and Herman, 2020; Sobral et al., 2021). Although some researchers have begun to explore genomic reaction norms represented by plasticity of gene expression underlying the TGP process, knowledge about the extent to which gene expression is affected by PE is still limited (Rivera et al., 2021). Colicchio et al. (2015) first reported effects of PE on gene expression plasticity in offspring by discovering patterns of gene expression plasticity in response to parental leaf damage in Mimulus guttatus. In the same year, another study discovered the effects of parental-experienced salinity stress on expression of stored seed transcripts in Medicago truncatula (Vu et al., 2015). Some additional studies also revealed patterns of selective gene expression and transcriptional memory underlying the TGP process; however, most of these studies focused on animals (e.g., Clark et al., 2019; Shi et al., 2020). More evidence is needed to show a link between TGP and the plasticity or transcriptional memory of gene expression, especially in some widely distributed plants that are excellent at adapting to heterogeneous environments. Furthermore, it is well accepted that different genes usually work together as a complicated module and are simultaneously recruited in development to achieve a phenotype (Uller et al., 2018). However, knowledge about how co-regulated genes are recruited in the TGP process is still limited.
Field pennycress (Thlaspi arvense L.) is a diploid plant (2n = 2 × = 14) of Brassicaceae and widely distributes in temperate regions of the northern hemisphere (Best and McIntyre, 1975; Warwick et al., 2002; An et al., 2015). Field pennycress is commonly used as a biofuel to extract oil from seeds or as a cover crop to reduce soil erosion (Mitich, 1996; Sedbrook et al., 2014; Dorn et al., 2015; McGinn et al., 2019). It can also be eaten as a vegetable or used in traditional Chinese medicine. The genome of field pennycress was sequenced and assembled in our previous study, and was determined to have a genome of 527.3 Mb encoding 31,596 genes (Geng et al., 2021). Moreover, the results of an experiment that focused on transgenerational changes and dynamics of DNA methylation in response to abiotic stress, such as salinity in field pennycress, showed that the epigenetic diversity at the population level significantly increased, and some of the epi-alleles showed that transient epigenetic memory lasted at least two generations in offspring in non-stressed environment (Geng et al., 2020). It suggested that field pennycress might have some degree of TGP under environmental stress. Moreover, field pennycress mainly reproduces by self-pollination and produces seeds in less than 10 weeks (Warwick et al., 2002; Geng et al., 2021). These features make field pennycress an ideal for exploring the molecular basis underlying TGP as it is operable to control the genetic background and generate different generations in a short period.
Cadmium (Cd) is a metallic element with high pro-oxidant potential that is toxic to nearly all organisms (Kinraide, 2009). Last decades have witnessed serious Cd contamination in agricultural systems caused by increasing anthropogenic activities and industrialization (Cai et al., 2019). Some studies revealed that field pennycress is sensitive to Cd stress, and a series of plastic trait changes were identified in field pennycress response to Cd stress, including reduction in root length and shoot fresh weight, and decrease in leaf osmolality and membrane permeability (Monteiro and Soares, 2021). However, an ecotype that exhibited considerable Cd tolerance was reported in the Cd-polluted urban area of Jena, Germany, with higher activities of antioxidant enzymes and efficiency in excluding Cd compared with the normal ecotype (Martin et al., 2012). It can be inferred that TGP may participate in evolution of normal ecotypes of field pennycress in response to persistent Cd stress. Adaptive TGP is widely believed to be a fitness-enhancing strategy in similar environments across generations, thus it might play an important role in ecological adaptation and evolutionary processes (Bonduriansky, 2021). Analyses focusing on TGP and corresponding transcriptional memory of gene expression in response to Cd stress across different generations are useful for determining the role TGP plays in processes of plant adaptation and evolution under the persistence of environmental selection pressure.
This study aims to examine the pattern of trait reaction norms and its relationship with the corresponding co-regulation pattern of gene expression underlying TGP of field pennycress in response to Cd stress. We measured morphological traits (e.g., individual height, and biomass), reproductive traits (e.g., flowering time, number, and hundred grain weight of seeds), and Cd tolerance related physiological traits [e.g., content of reduced glutathione (GSH) and activity of superoxide dismutase (SOD), malondialdehyde (MDA) and photosystem II activity (Fv/Fm)] under control conditions and Cd treatment in parent individuals and offspring. The gene expression profiles in samples of different generations were also examined at the same time by RNA sequencing. This design allowed us to: (1) detect and compare reaction norms that represent WGP and TGP in response to increased soil Cd; (2) identify genes and pathways that play important roles in Cd stress adaptation in T. arvense; and (3) discover expression patterns of key genes that showed transcriptional memory in offspring, and assess their topological role in the recruitment of co-regulated gene modules underpinning TGP of field pennycress in response to Cd stress.
Materials and methods
Plant material and cadmium treatment
Two generations of field pennycress were used in this study. For generation 1 (G1), seedlings were germinated from seeds of individuals collected from a wild population in Kunming (KM, alt. 1,910 m, N 30.313°, E 99.358°). The maternal plant was maintained in the greenhouse at Yunnan University, Kunming, for more than 2 years to eliminate wild environmental effects on seedlings used in experiment. After the first two new leaves appeared, the seedlings were individually transplanted into pots (9.5 cm in diameter, 7.5 cm in depth), and each pot contained 200 g of nutrition soil. CdCl2 solution was mixed with nutritious soil before transplant to conduct Cd treatment, with concentrations of 0 (CK), 25, 50, 75 (Cd75), 125, and 185 mg/kg (Cd185), for six Cd-treatment gradients (Figure 1). The G1 individuals were grown during July and September in 2020 with approximately 12 h of daytime and harvested in the middle of September.
Figure 1. Experimental design. In G1, seedlings from maternal plant were treated with a gradient of Cd from 0 to 185 mg/kg. Seeds collected from CK and Cd185 were used to conduct experiment in G2 that treated offspring seedlings with CK and Cd185. Several morphological traits, reproductive traits and physiological traits of all individuals were measured, and samples marked by red star were further used for RNA sequencing. Transcriptional memory genes (TMGs) were identified by screened genes that showed the same expression pattern (both upregulated or downregulated) between Cd185 vs. CK, CK-S vs. CK-CK, and S-CK vs. CK-CK simultaneously.
To determine TGP, seeds collected from Cd 185 and CK individuals in G1 were used to germinate generation 2 (G2) seedlings. A total of four treatments were set for G2 individuals (Figure 1): (1) CK-CK: offspring of individuals in G1-CK treatment, continued growth in pots with the CK environment; (2) CK-S: offspring of individuals in G1-CK treatment, growth in pots with Cd185 treatment; (3) S-CK: offspring of individuals in G1-Cd185 treatment transplanted to pots with the CK environment; (4) S-S: offspring of individuals in G1-Cd185 treatment, continued growth in pots with Cd185 treatment. Individuals of G2 were grown from November to the following January in 2021, and the ambient temperature and illumination conditions were similar to those of G1 using compensation equipment for heating and lighting.
Trait measurement
The following traits were measured for G1 and G2 individuals for at least six replications. For morphological traits: (1) The total biomass was determined based on the dry weight of the plant. (2) The individual height was measured by the length from the top to the end of stem. For reproductive traits: (3) The flowering time was defined as the number of days it took for the first flower to bloom. (4) The seed traits included the number of seeds and hundred-grain weight of seeds. For physiological traits: (5) The photosystem II activity was measured as Fv/Fm, which is the ratio between light-induced variable and maximum fluorescence of chlorophyll, using a pulse-amplitude modulated fluorometer (Junior-PAM, Walz GmbH, Effeltrich, Germany). (6) The SOD activity in leaves was determined by the SOD reagent box (A001-1, Nanjing Jiancheng, Nanjing, China). (7) The content of reduced GSH was determined using a CheKine reduced GSH colorimetric assay kit (Abbkine, Wuhan, China). (8) The level of lipid peroxidation (measured as MDA) in leaves was determined using a lipid peroxidation MDA assay kit (Abbkine, Wuhan, China). All measurements of physiological traits were determined according to the manufacturers’ instructions. Tukey’s HSD was used to test for significant differences between different treatments. A two-way ANOVA was used to test the effects of G1 PE, G2 offspring environment (OE), and their interactions.
Transcriptome sequencing and data analysis
According to the results of traits measurement under different concentrations of Cd stress, we found most of these traits showed significant differences under 185 mg/kg of Cd (Figures 2A,B). Thus we selected 185 mg/kg as an effective and intense concentration of Cd stress to induce morphological trait changes, and 125 mg/kg, 75 mg/kg were considered as a moderate and slight concentration, respectively. Finally, samples of CK, Cd75, and Cd185 in G1, and CK-CK, CK-S, S-CK, S-S in G2 were used to conduct RNA sequencing. For each treatment, three independent biological replications of each treatment were randomly selected. Total RNA of young and mature leaves was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and then equally mixed for each sample. RNA quality assessment was conducted using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). High-quality RNAs were then used for RNA sequencing by Novogene Bioinformatics Technology Co., Ltd, Beijing, China. Briefly, a standard library construction process, including fragmentation, first and second strand cDNA synthesis, adaptor ligation, PCR, and quality assessment, was conducted using purified mRNA from total RNA. After clustering process, a total of 21 libraries [(3 G1 treatments × 3 biological replicates) + (4 G2 treatments × 3 biological replicates)] were then sequenced on an Illumina Novaseq platform to generate 150-bp paired-end reads.
Figure 2. Trait expression of field pennycress in response to different soil Cd concentrations. (A) Biomass; (B) individual height; (C) flowering days; (D) number of seeds; (E) hundred-grain weight; (F) Fv/Fm; (G) SOD; (H) GSH; and (I) MDA. Quoted values are means ± SE (N = 6). Different letters indicate significant differences among groups determined by Tukey’s HSD.
Clean data were obtained by removing reads containing adapter, ploy-N, and low-quality reads from raw data. Paired-end clean reads were aligned to the field pennycress reference genome (Geng et al., 2021) using Hisat2 v2.0.5. The mapped reads of each sample were assembled by StringTie (v1.3.3b), and featureCounts v1.5.0-p3 was used to count the number of reads mapped to each gene to calculate FPKM. In addition to genome annotation, transcripts were annotated to the TAIR10 database by BLAST.
Differential expression analysis was conducted based on FPKM using DESeq2 R package (1.20.0). Different standards for absolute log2 fold change for differentially expressed genes (DEGs) varied between studies, normally in a range of 0.3–2.0 (e.g., Wang et al., 2017, Montaldo et al., 2020). Genes with an adjusted absolute p-value < 0.05 in comparisons were identified, and log2 fold change was used to evaluate the degree of expression difference, with an absolute value above 0.5 assigned as significantly differentially expressed (Liu and Sun, 2021). According to this criterion, the Cd-responsive genes were identified by Cd75 vs. CK, Cd185 vs. Cd75, and CK-S vs. CK-CK. Some Cd-responsive genes enhanced their expression with increased Cd concentration, and these genes were identified by comparing Cd185 vs. Cd75.
To determine the transcriptional memory across generations, Cd-responsive genes that showed across generational persistent expression even if growth condition changed to normal in G2, were defined as “transcriptional memory genes” (TMGs). TMGs were identified by screened genes that showed the same expression pattern simultaneously in Cd185 vs. CK, CK-S vs. CK-CK, and S-CK vs. CK-CK. Briefly, we first screened genes with the similar expression pattern (both upregulated or downregulated) between Cd185 vs. CK and CK-S vs. CK-CK. This step was used to find robust Cd-responsive genes that can be induced or repressed by Cd stress across generations. The results were then cross-referenced with S-CK vs. CK-CK to select robust Cd-responsive genes with the same expression pattern in absence of Cd stress in G2. As a result, the intersection part was considered as TMGs (Figure 1). We further examined the expression enhancement of TMGs under the condition of consistent Cd stress across generations, and this was done by comparing expression pattern of TMGs with that under S-S vs. CK-S. The Gene ontology enrichment for Cd responsive genes was compared and visualized using WEGO 2.0 (Ye et al., 2018) and GO terms exhibiting a corrected p-value of ≤0.05 were considered significantly enriched.
Gene co-expression network analysis
Gene co-expression network analysis was conducted using the WGCNA R package (Langfelder and Horvath, 2008). Briefly, FPKM data of Cd-responsive genes were normalized using the limma R package, and genes with too many missing entries were filtered using goodSamplesGenes function. The function pickSoftThreshold was then used to calculate scale free topology for multiple soft thresholding powers. According to the result, the unsigned network was constructed by setting a power of 16, minModuleSize of 30, and a mergeCutHeight of 0.25. Eigengenes for each module were calculated by using moduleEigengenes function with default parameters. Pearson’s correlation coefficient (PCC) of the module eigengenes to trait data was then assessed to identify key modules involved in Cd response and plasticity of transgenerational traits. After network construction, the top 10% of weighted edges were analyzed because millions of edges were generated. The module hub genes were identified by nodes with the top 5% degree in each module.
In order to examine whether the modules were robust and reproducible, the module preservation statistics were computed using the modulePreservation function implemented in WGCNA (Langfelder et al., 2011). The multiExpr was set by G1 and G2 samples to test module preservation between generations, with G1 as reference, and the number of permutations was set to 200. The computation generated a Zsummary and a medianRank for each module. For values of Zsummary, 2 > Z represented no preservation, 2 < Z < 10 represented weak to moderate preservation, and Z > 10 represented strong preservation. Module with lower median rank exhibit higher observed preservation statistics than a module with a higher median rank (Sharma et al., 2018).
Quantitative real-time PCR verification of transcriptional memory genes
To validate the expression pattern of TMGs between generations, 1.0 μg of total RNA of each sample was reverse transcribed into cDNA using a PrimeScript™ RT reagent kit (TaKaRa, Dalian, China). The quantitative real-time PCR (qRT-PCR) was performed using TB Green™ Premix Ex Taq™ II (TaKaRa, Dalian, China) and the Roche LightCycler 96 system (Roche, Penzberg, Germany) following the manufacturers’ instructions. UBQ10 was selected as the internal reference gene. All reactions were performed in triplicate for each sample, and the relative changes for comparison between different treatments were quantified using the 2–ΔΔCt method as described in a previous study (Li et al., 2017). Primers are listed in Supplementary Table 1.
Results
Plastic responses of field pennycress to cadmium treatment in generation 1
For morphological traits measured in G1, total biomass and individual height were similar in Cd treatment of 0–125 mg/kg but significantly decreased in 185 mg/kg Cd treatment (p < 0.01, Figures 2A,B). The flowering time was significantly increased in 0–50 mg/kg Cd treatment (p < 0.01), and this was maintained in a higher concentration of soil Cd (Figure 2C). Seeds production decreased along with the increase of Cd concentration, and individuals exposed to Cd185 treatment produced fewer seeds compared with those exposed to 0 (p < 0.01) and 25 mg/kg (p = 0.03) (Figure 2D). For the physiological traits, total SOD and GSH significantly increased with the increase of soil Cd concentration (p = 0.036 and p < 0.01 for SOD and GSH, respectively; Figures 2G,H).
Effects of generation 1 parental environment on generation 2 offspring traits
Several traits in G2 showed similar changes in direction along with the change of environmental conditions (Figure 3). For example, biomass, flowering days, Number of seeds, hundred-grain weight, Fv/Fm, GSH and MDA didn’t show significant difference for CK-CK vs. S-CK and CK-S vs. S-S with Tukey’s test (p-value > 0.05) (Supplementary Table 7), and presented similar reaction norm in all offsprings (Figures 3A,C–F,H,I), indicating that PE in G1 did not affect the performance of offspring in G2 under different OEs. The individual height showed values that close to significant difference between comparison of CK-CK and S-CK (p = 0.056 in Tukey’s HSD test, Figure 3B). The SOD is the only trait that significantly differed between CK-CK and S-CK (p = 0.02 in Tukey’s HSD test, Figure 3G), ANOVA showed the significant effect of PE-by-OE interaction can only be detected on SOD (p < 0.01, Table 1).
Figure 3. Reaction norm of trait expression in G2 individuals exposed to soil Cd concentrations of 0 mg/kg (CK) and 185 mg/kg (S). The red and black lines show seedlings that originated from offspring of Cd185 (S) and CK treatment in G1, respectively. (A) Biomass; (B) individual height; (C) flowering days; (D) number of seeds; (E) hundred-grain weight; (F) Fv/Fm; (G) SOD; (H) GSH; and (I) MDA. Quoted values are means ± SE (N = 6). Different letters indicate significant differences among groups determined by Tukey’s HSD.
Table 1. Effects of the parental environment (PE), offspring environment (OE), and their interactions on traits of field pennycress in generation 2.
RNA sequencing of field pennycress in response to cadmium treatment
For each library, an average of 6.55 G of clean data was generated. The average error rate was 0.02, the average Q30 was 94.88%, and the GC contents ranged from 45.62 to 47.99% (Supplementary Table 2). A total of 95.49 to 98.26% of clean reads were mapped to the field pennycress reference genome, which had been deposited in the National Center for Biotechnology Information (NCBI) with BioProject accession number PRJNA715950 (Geng et al., 2021), with the unique mapping rate ranging from 90.67 to 95.90% for each library (Supplementary Table 3). A total of 86.76 to 90.14% of clean reads were mapped to exon regions, and approximately 2.31 to 2.91% of clean reads were mapped to intron regions (Supplementary Table 4). Different biological replicates showed a high correlation with average squared PCC of 0.93. Finally, 33061 transcripts were identified according to the mapping result. The high quality of sequencing data allowed subsequent analysis.
Genomic reaction norm in response to cadmium treatment
The Cd-responsive genes were identified by three comparison pairs of Cd-treated samples vs. control, including Cd75 vs. CK, Cd185 vs. CK, and CK-S vs. CK-CK. A total of 3,975, 6,269, and 5,356 Cd-responsive genes were identified in the Cd75, Cd185, and CK-S libraries, respectively. The results showed that high concentrations (Cd185, and CK-S) of Cd treatment led to more responsive genes than low concentrations (Cd75). After merging results, a total of 10,028 Cd-responsive genes were obtained (Figure 4A). Among them, only 762 (7.50%) genes were exclusively included in Cd75 vs. CK, indicating that these genes only participate in slight (Cd75) Cd stress, while the rest of them can be detected in both slight and intense (Cd185) Cd stress (3213 genes, 31.7%), or only participate in intense Cd stress (6,164 genes, 60.8%) (Figure 4A). However, there are only 2,248 genes identified to be differentially expressed in both G1 (Cd185 vs. CK) and G2 (CK-S vs. CK-CK) (Figure 4A), indicating that these genes can be induced differential expression robustly by Cd-stress in both G1 and G2, and differential expression of other genes (4,021 in G1 and 3,108 in G2) maybe caused by batch effects or some other stochastic factors rather than Cd-stress. The numbers of upregulated genes were 2,058, 3,331, and 2,815 for Cd75 vs. CK, Cd185 vs. CK, and CK-S vs. CK-CK, respectively (Figure 4B). The expression patterns of Cd-responsive DEGs are shown in Figure 4C.
Figure 4. Expression and function of Cd-responsive genes. (A) Comparison and distribution of Cd-responsive genes in G1 (Cd75, Cd185) and G2 (CK-S) using CK (G1) and CK-CK (G2) as controls. (B) Numbers of upregulated and downregulated genes. (C) Differential expression patterns of Cd-responsive genes. The color represents log2 fold change values; negative values (blue) represent downregulation and positive values (red) represent upregulation. (D) Gene Ontology (GO) classifications of the DEGs. (E,F) Expression and biological processes of genes with enhanced upregulation (E) or downregulation (F) along with the increase of soil Cd concentration in G1. For the heatmap, color was adjusted to log2 fold change of comparison pairs. For GO map results, the color and size were adjusted to p-value and log10 p-value, respectively.
The GO enrichment analysis showed that these DEGs participated in biological processes such as response to environmental stimulus, regulation of growth and developmental processes, cellular detoxification, and reproductive processes (Figure 4D). Furthermore, 131 and 118 genes had enhanced upregulated or downregulation along with the increased concentration of soil Cd (Figures 4E,F). These upregulated genes were enriched in cellular response to salicylic acid and abscisic stimulus, unsaturated fatty acid biosynthetic process, and protein phosphorylation etc.; while the downregulated genes were enriched in histone H3-K36 demethylation, cellular response to sulfur starvation, and oxylipin metabolic process etc. (Figures 4E,F) in the GO enrichment analysis.
Cadmium responsive transcriptional memory in generation 2 individuals
The TMGs were identified by screening genes that showed the same expression pattern (both upregulated and downregulated) in intersection of three comparison pairs (Cd185 vs. CK, CK-S vs. CK-CK, and S-CK vs. CK-CK). A total of 315 upregulated and 86 downregulated TMGs were identified (Figure 5A and Supplementary Table 5), with the expression pattern shown in Figure 5B. The GO enrichment analysis showed that TMGs participated in some important biological processes. For example, the downregulated TMGs were involved in biological processes such as response to abscisic acid (ABA), the salicylic acid-mediated signaling pathway, and regulation of multicellular organism growth (Figure 5C), whereas the upregulated TMGs participated in the biological processes of cell wall biogenesis, response to light stimulus, the hydrogen peroxide catabolic process, and DNA methylation of cytosine within a CHH sequence (Figure 5D).
Figure 5. Identification and characterization of Cd-responsive TMGs. (A) Comparison and distribution of upregulated and downregulated genes in Cd185 (G1), CK-S (G2), and S-CK (G2) using CK (G1) and CK-CK (G2) as controls. Genes with intersection were identified as TMGs. (B) Heatmap showing the expression profile of TMGs. The color scale represents log2 fold change between different comparison pairs. (C,D) Biological processes of downregulated (C) and upregulated (D) TMGs. Representative biological processes are labeled; the color and size were adjusted to p-value and log10 p-value, respectively.
According to expression pattern, the TMGs were further compared with the S-S vs. CK-S results to identify the TMGs that may enhance expression under the S-S environment when G1 and G2 both experienced Cd treatment. As a result, a total of 17 TMGs (7 downregulated and 10 upregulated) that showed enhancement of expression under the S-S environment were identified, and the detailed information was provided in Supplementary Table 6.
Gene co-expression modules related to cadmium stress adaptation and measured traits in leaves
A total of 18 modules were identified based on pairwise correlations of gene expression across all samples (Supplementary Figure 1). The results of module preservation statistics showed that four modules (gray, grey60, tan, and salmon) were not preserved as their Zsummary below 2 and medianRank above 12, five modules (red, green, yellow, midnight blue, and lightcyan) showed a weak to moderate preservation as their Zsummary range from 3.4 to 9.8 and medianRank range from 1 to 17. The rest of modules showed a good module preservation (Supplementary Figure 2).
Because only leaf samples were used in RNA sequencing, traits of other organs, including individual height, number of seeds, and hundred-grain weight, were not considered in the module-trait associations analysis. As flowers were modified leaves, flower primordia may be included when collecting young leaf samples; thus, flowering days were included in the analysis. Results of module-trait associations analysis showed that modules related to Cd-treatment also tend to associate with GSH and SOD (Figure 6A). According to the criterion of p-value < 0.05, four modules (purple, turquoise, greenyellow and red) were positively related to Cd-treatment, GSH and SOD, five modules (cyan, tan, green, black and brown) were positively related to biomass, and three modules (gray 60, greenyellow and yellow) were positively related to flowering days (Figure 6A).
Figure 6. Correlation of traits and modules, and distribution of TMGs in different modules. (A) Module–trait associations. The module eigengenes and traits are represented in rows and columns, respectively. The corresponding correlation (represented by value and color) and p-value are contained in each cell. Number beside module name indicate count of genes in this module. (B) A global view of TMGs distributed in each module. Modules are represented by different colors, and key hub genes are represented by larger circles. The TMGs are surrounded on the outside. (C) Proportion of TMGs in trait-correlated modules (PCC > 0.6 and p < 0.05). More TMGs were found in SOD-related modules than those of other traits.
Transcriptional memory genes were found to be dispersed into different modules. Among them, most TMGs were assigned to the turquoise (228 TMGs out of 1821 nodes), brown (29 out of 702), black (27 out of 251), and greenyellow modules (26 out of 118), and the remaining modules included fewer than 20 TMGs (Figure 6B). The percentage of TMGs in trait-related modules with PCC > 0.6 was calculated, and the results showed that more TMGs were found in SOD-related modules than in other traits (Figure 6C).
To discover the influences of TMGs on topological structure and overall function of gene co-expression modules in response to Cd treatment, four modules (purple, turquoise, greenyellow, and red) that positively correlated with Cd treatment, GSH, and SOD were focused on. Several upregulated TMGs and other functional related genes were included in these modules. Functional analysis showed module genes participate in biological processes such as cell wall organization or biogenesis, lipid biosynthesis and metabolism, response to Cd, and other abiotic stimuli (Figure 7A). Among them, some TMGs were identified as hub genes. For example, harbinger transposon-derived protein 2 (HDP2) (Ta.Chr.1.1011) is a transcription factor of the trihelix family. According to the annotation, this gene showed a function in prevention of DNA hypermethylation and epigenetic silencing. Moreover, HDP2 is co-expressed with the most genes (a total of 1,365 nodes) in the turquoise module, indicating its important role in the topological construction of this key module. By contrast, analysis of flowering day-related modules showed that, although the module had a large number of genes associated with flower development, no TMGS existed in this module (Figure 7B). This result may somehow be correlated with no transgenerational plastic memory in flowering day traits.
Figure 7. The TMGs in trait-related modules. (A) TMGs that participate in Cd stress adaptation and are co-expressed with other important genes in a sub-network of four modules (turquoise, red, greenyellow, and purple) that are correlated with Cd treatment, GSH, and SOD. Pictures showing the expression pattern of TMGs and its co-expression relationship with other genes. Function and corresponding gene numbers are listed at the bottom. (B) Gene co-expression relationship, expression pattern, and function in the yellow module, which is highly correlated with flowering time. No TMGs were identified in this module. (C) Sub-network formed by ACO3, CSD2, and their neighbors. Different nodes are classified according to gene function. For simplicity, neighbors of the other functions are not shown. Key hub genes are represented by larger circles, TMGs are represented by large squares, and modules that the gene belongs to are represented by different colors.
Although transgenerational memory was detected in SOD activity, we failed to detect transcriptional memory in copper/zinc superoxide dismutase 2 (CSD2) (Ta.Chr.3.2310), the unique superoxide dismutase coding gene in the turquoise module. However, ACONITASE 3 (ACO3) (Ta.Chr.5.1524), which interplays with CSD2, was found to be one of the TMGs in the brown module. We isolated ACO3, CSD2, and their neighbors from the global network; as a result, a sub-network was formed (Figure 7C). Like the sub-network in Figure 7A, this sub-network also consisted of several TMGs and other functions-related genes that participate in biological processes including developmental growth, response to light stimulus, lipid biosynthesis, response to ABA, cell wall biogenesis, response to oxidative stress, and other processes (Figure 7C). Furthermore, the sub-network of ACO3-CSD2-neighbors included a total of 110 TMGs, of which 57 (51.82%) were hub genes according to the criteria of the top 5% of degree in each module.
Validation of the expression pattern of transcriptional memory genes by quantitative real-time PCR
For validation of RNA sequencing accuracy, the expression of six TMGs was analyzed through qRT-PCR, as shown in Figure 8. The relative expressions of Cd185 vs. CK, CK-S vs. CK-CK, and S-CK vs. CK-CK were calculated using the qRT-PCR results. The results showed a similar pattern to RNA-seq. For example, for the cell wall organization gene powdery mildew resistant 6 (PMR6) (Ta.Chr.2.4195), the results of qRT-PCR showed the same pattern as fold change calculated by RNA-seq data (Log2 fold change value of 1.49, 2.47 and 1.18 for Cd185 vs. CK, CK-S vs. CK-CK, and S-CK vs. CK-CK, respectively), both of them showed an upregulated pattern (Figure 8B). These results demonstrate the accuracy of RNA-seq in this study.
Figure 8. Validation of expression pattern of TMGs by qRT-PCR. (A–F) qRT-PCR results for six randomly selected TMGs. The y-axis represents the log2 ratios between comparison pairs along the x-axis. Relative expression quantification was carried out using the 2– ΔΔCt method, with UBQ10 as the internal reference gene. Error bars represent mean ± SE (n = 3).
Discussion
Within-generation and transgenerational plasticity in field pennycress
Phenotypic plasticity is common in nature and essential for organisms surviving in a changing world (Kelly et al., 2012; Forsman, 2015), especially for plants with a sessile lifestyle. In this study, we examined within-generation and TGP in field pennycress in response to different concentrations of soil Cd. In the parental generation (G1), a series of phenotypic changes were observed along with increasing Cd concentration, including a decrease in individual height and biomass, longer time to flowering, and enhancement of SOD activity and content of GSH (Figure 2). These phenotypic variations are not necessarily adaptive. For example, the decreases in individual height and biomass in Cd-treated individuals are the same as the results of another study that conducted similar Cd treatment on field pennycress (Monteiro and Soares, 2021); this indicates that growth inhibition by Cd is universal between different genotypes. However, growth inhibition was not observed under the same conditions in the congener Thlaspi praecox, which is a hyperaccumulator that is highly tolerant to Cd stress (Tolra et al., 2006; Monteiro and Soares, 2021).
In contrast, enhancement of SOD activity and content of GSH may be useful to increase fitness under long-term Cd stress. The main mechanism of Cd phytotoxicity is the induction of oxidative stress (Loix et al., 2017). When plants deal with oxidative stress, SOD constitutes the first line of defense against ROS because it is involved in scavenging and detoxification of ROS (Alscher et al., 2002; Mishra and Sharma, 2019). Moreover, GSH participates in metal homeostasis, antioxidative defense, and signal transduction under metal-induced oxidative stress (Jozefczak et al., 2012). It was clear that SOD activity and content of GSH were increased in response to high concentrations of soil Cd, these increases may be beneficial to field pennycress to maintain normal physiological function and protect cells from injury during oxidative stress.
Furthermore, TGP was found for SOD activity in offspring. A significant parental-by-offspring interaction (PE × OE) was identified in SOD activity (Table 1). The offspring derived from G1-S-treated plants maintained higher SOD activity, even when grown under normal conditions in G2 (S-CK), compared with offspring derived from growth under normal conditions in G1 and G2 (CK-CK) (Figure 3G). For species growing in sufficiently stable environments, many lines of evidence have confirmed that effects of TGP can be adaptive when parent and progeny environments match (Walsh et al., 2016; Colicchio and Herman, 2020; Halali and Saastamoinen, 2022). In other words, inheritance of adaptive trait changes through TGP can lead to an increase of fitness in offspring when parent and progeny environments match or can be expected to be similar (Sandner et al., 2018; Colicchio and Herman, 2020). Limited ability for seed dispersal in field pennycress makes growth conditions highly predictable across generations. Therefore, transgenerational maintenance of SOD activity may be considered as adaptive for offspring growth under Cd-polluted environment in which future stress could be predicted according to the environment of maternal plants.
On the contrary, according to the theory of limit and cost of phenotypic plasticity, the expression of plastic traits may lead to impairment in fitness if there is a mismatch between final phenotype and environment (DeWitt et al., 1998; Relyea, 2002; Murren et al., 2015). TGP of SOD activity may not be adaptive for individuals whose habitat is converted to a normal condition, which is inaccurately predicted by the maternal environment. However, the consequence of adaptive TGP on SOD activity was not obvious in S-S compared with CK-S individuals. It can be inferred that differences in fitness may be detected in individuals derived from maternal plant growth under Cd stress for several generations rather than using those that were only treated for one generation.
Transcriptional memory genes played important roles in adapting to cadmium stress in Thlaspi arvense
As a sort of phenotype at the molecular level, genomic reaction norm, which is represented by gene expression under different conditions, is useful for understanding mechanisms of phenotypic plasticity (Aubin-Horth and Renn, 2009). According to the results of Cd-responsive genes in G1, approximately 1.57 times as many DEGs were induced by a high concentration of soil Cd (Cd185) compared with a low concentration (Cd75), and some Cd-responsive genes showed a pattern of enhanced expression with the increase of Cd concentration (Figures 4E,F). These results indicated that genomic reaction norm and phenotypic reaction norm present non-random changes in direction and slope along with environmental gradient. Moreover, according to function analysis of Cd-responsive genes, some biological processes, including response to oxidative stress, cell wall biosynthesis, sulfur utilization, and ABA signal transduction, were also found to be enriched in other species in response to Cd stress, such as tall fescue (Festuca arundinacea), pakchoi (Brassica chinensis L.), and rice (Oryza sativa L.) (Zhou et al., 2016; Zhu et al., 2018; Sun et al., 2019). The results indicated that these processes might be universal in Cd responses of plants.
The screening criteria for TMGs is rigorous and based on intra- and intergenerational biological replicates. As a result, only approximately 4% of Cd-responsive genes were identified as TMGs. Some TMGs have been proven to play an important role in response to Cd stress in plants. For example, several upregulated TMGs, such as transcripts of NADPH-dependent thioredoxin reductase C (NTRC) (Ta.Chr.3.760) and 2-Cys peroxiredoxin A (2CPA) (Ta.Chr.5.3838), were found to be involved in the biological processes of superoxide radical removal and hydrogen peroxide catabolism, which are useful for ROS homeostasis and chloroplast protection against oxidative damage (Pérez-Ruiz et al., 2006; Zheng et al., 2020). The ABA signaling pathway plays an important role in Cd stress tolerance (Leng et al., 2021). Several genes involved in regulation of the ABA signaling pathway showed transcriptional memory of upregulation in Cd-treated maternal plants and offspring grown under normal condition, such as transcripts of receptor dead kinase 1 (RDK1) (Ta.Chr.1.701), tetratricopeptide-repeat thioredoxin-like 1 (TTL1) (Ta.Chr.4.2161), and LYSM-CONTAINING RECEPTOR-LIKE KINASE 3 (LYK3) (Ta.Chr.4.2285) (Rosado et al., 2006; Paparella et al., 2014; Kumar et al., 2017). Additionally, the cell wall is a barrier to prevent Cd from entering and damaging the protoplast. Several upregulated TMGs were involved in cell wall organization or biogenesis, such as xyloglucan endotransglucosylase/hydrolase 9 (XTH9) (Ta.Chr.0.3619), GPI-anchored protein (GPI) (Ta.Chr.5.3158), and PLAT DOMAIN PROTEIN 2 (PLAT2) (Ta.Chr.3.4335) (Kushwah et al., 2020; Depuydt and Vandepoele, 2021).
The expression patterns of some of these genes were examined by qRT-PCR (Figure 8). Furthermore, a DNA methyltransferase associated with maintenance of genome-wide DNA methylation, chromomethylase 2 (CMT2) (Ta.Chr.1.2289) (He et al., 2021), also showed transcriptional memory of upregulation in Cd-treated maternal plants and offspring grown under normal conditions. The transcription memory of these key Cd stress-responsive genes may influence the fitness of offspring under long-term Cd stress. Our previous study showed that some of the DNA methylation modification of epi-alleles caused by salinity stress can be inherited by at least two generations of offspring in non-stressed environments (Geng et al., 2020); thus, it can be inferred that epigenetic memory may contribute to the transcriptional memory in T. arvense. In addition to epigenetics, many other mechanisms contribute to transcriptional memory, such as transmission of mRNAs, proteins, and hormones from parent to offspring (Herman and Sultan, 2011). More experiments are needed to examine the causes of transcriptional memory in field pennycress in response to Cd stress.
Recruitment of transcriptional memory genes in co-expression gene modules correlated with transgenerational plasticity
The recruitment of most TMGs by a limited number of modules reveals that co-regulation of TMGs may play an important role in the TGP of T. arvense. The superoxide dismutase coding gene CSD2, which was recruited in the turquoise module, was not identified as a TMG according to our criteria. However, we found that CSD2 was co-expressed with many other oxidative stress-responsive genes, and several of these genes showed transcriptional memory (Figure 7C). ACO3 encodes aconitase and is co-expressed with CSD2, and its expression pattern was validated by qRT-PCR (Figure 8F). Evidence showed that aconitase acts as a negative regulator that is specifically bound to 5′ UTR of the Arabidopsis CSD2 in vitro, an increased level of CSD2 transcript was observed in aconitase-knockout plants, leading to more tolerance to oxidative stress in aconitase-knockout plants (Moeder et al., 2007).
In this study, ACO3 was one of the TMGs that was downregulated in both G1 and G2 Cd-treated individuals, and also showed transcriptional memory in G2 S-CK individuals. Transcriptional memory of ACO3 may indirectly affect SOD activity in G2 S-CK individuals, thus contributing to TGP of SOD activity. The sub-network constructed by CSD2, ACO3, and their neighbors included a considerable amount of hub TMGs, indicating recruitment of TMGs in offspring potentially related to TGP. In addition to ACO3, the other TMGs involved in biological processes, including lipid biosynthesis, cell wall biogenesis, response to ABA, developmental growth, and response to light stimulus, were co-expressed as neighbors of ACO3 and CSD2 (Figure 7C), some of them were proven to play important roles in Cd stress tolerance (Abd_Allah et al., 2015; Loix et al., 2017; Leng et al., 2021). This result indicates that these pathways may be associated with TGP of SOD activity in T. arvense. In contrast, there were no TMGs found in the yellow module, which was highly related to flowering time. Many genes that participate in flower development and transition from vegetative to reproductive phases were recruited in the yellow module; however, none of them were identified as TMGs or co-expressed with TMGs (Figure 7B). This may be one of the reasons why TGP was not observed for this trait.
The increase of SOD activity induced by TGP in offspring may be recognized as the first step in the process of genetic assimilation. Genetic assimilation is a phenomenon that describes the decrease of environmental sensitivity of gene expression, which may lead to fixed or constitutive expression of environmentally induced traits in the absence of a recurrent environmental cue (Waddington, 1961; Ehrenreich and Pfennig, 2016; Nijhout et al., 2021). It can be inferred from the results that high SOD activity maintained in offspring may ultimately be fixed if a high soil Cd concentration is maintained for many generations, just as in the Cd-tolerant field pennycress ecotype reported in the Cd-polluted urban area of Jena, Germany (Martin et al., 2012). Transcriptional memory, which is often neglected in research on TGP, may play an important role in the processes of trait evolution mediated by phenotypic plasticity.
Data availability statement
The data presented in this study are deposited in the SRA database of NCBI repository, accession numbers: SRR18883184 – SRR18883204.
Author contributions
GL performed the data analysis and drafted the manuscript. YZ, FL, MS, TZ, and YGu performed the wet lab work. FZ, QQ, and YGe participated in data analysis and design of the study. YGe conceived the idea, participated in the design of the study, and finalized the manuscript. All authors have read and approved the final manuscript.
Funding
This study was supported by grants from the National Natural Science Foundation of China (32060059 to YGe, 32060085 to QQ, and 32060237 to TZ) and doctoral research initial funding of Yunnan Agricultural University (KY2022-65 to GL).
Acknowledgments
We thank Jinjing Jian, from Fudan University, for beautifying figures, and Mallory Eckstut, Ph.D, from Liwen Bianji (Edanz) (www.liwenbianji.cn), for editing the English text of a draft of this manuscript.
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.953794/full#supplementary-material
Supplementary Figure 1 | Gene co-expression network construction. (A) Clustering dendrogram of genes with the original module colors and assigned merged module colors. (B) Soft threshold of network construction. Power = 16 was selected.
Supplementary Figure 2 | The medianRank preservation statistics (A) and the Zsummary preservation statistics (B) of the modules. The dashed blue and green lines indicate the thresholds of 2 < Zsummary < 10.
References
Abd_Allah, E. F., Abeer, H., Alqarawi, A. A., and Hend, A. A. (2015). Alleviation of adverse impact of cadmium stress in sunflower (Helianthus annuus L.) by arbuscular mycorrhizal fungi. Pak. J. Bot. 47, 785–795.
Agrawal, A. A., Laforsch, C., and Tollrian, R. (1999). Transgenerational induction of defences in animals and plants. Nature 401, 60–63. doi: 10.1038/43425
Alscher, R. G., Erturk, N., and Heath, L. S. (2002). Role of superoxide dismutases (SODs) in controlling oxidative stress in plants. J. Exp. Bot. 53, 1331–1341. doi: 10.1093/jexbot/53.372.1331
An, M., Zeng, L., Zhang, T., and Zhong, Y. (2015). Phylogeography of Thlaspi arvense (Brassicaceae) in China inferred from chloroplast and nuclear DNA sequences and ecological niche modeling. Int. J. Mol. Sci. 16, 13339–13355. doi: 10.3390/ijms160613339
Aubin-Horth, N. A. D. I. A., and Renn, S. C. (2009). Genomic reaction norms: Using integrative biology to understand molecular mechanisms of phenotypic plasticity. Mol. Ecol. 18, 3763–3780. doi: 10.1111/j.1365-294X.2009.04313.x
Bell, A. M., and Hellmann, J. K. (2019). An integrative framework for understanding the mechanisms and multigenerational consequences of transgenerational plasticity. Annu. Rev. Ecol. Evol. Syst. 50, 97–118. doi: 10.1146/annurev-ecolsys-110218-024613
Best, K. F., and McIntyre, G. I. (1975). The biology of Canadian weeds: 9. Thlaspi arvense L. Can. J. Plant. Sci. 55, 279–292. doi: 10.4141/cjps75-039
Bonduriansky, R. (2021). “Plasticity across generations,” in Phenotypic Plasticity & Evolution: Causes, Consequences, Controversies, ed. D. W. Pfennig (Taylor and Francis: CRC Press), 327–348.
Cai, L. M., Wang, Q. S., Wen, H. H., Luo, J., and Wang, S. (2019). Heavy metals in agricultural soils from a typical township in guangdong province, China: Occurrences and spatial distribution. Ecotox. Environ. Safe. 168, 184–191. doi: 10.1016/j.ecoenv.2018.10.092
Clark, M. S., Suckling, C. C., Cavallo, A., Mackenzie, C. L., Thorne, M. A., Davies, A. J., et al. (2019). Molecular mechanisms underpinning transgenerational plasticity in the green sea urchin Psammechinus miliaris. Sci. Rep. 9:952. doi: 10.1038/s41598-018-37255-6
Colicchio, J. M., and Herman, J. (2020). Empirical patterns of environmental variation favor adaptive transgenerational plasticity. Ecol. Evol. 10, 1648–1665. doi: 10.1002/ece3.6022
Colicchio, J. M., Monnahan, P. J., Kelly, J. K., and Hileman, L. C. (2015). Gene expression plasticity resulting from parental leaf damage in Mimulus guttatus. New Phytol. 205, 894–906. doi: 10.1111/nph.13081
Depuydt, T., and Vandepoele, K. (2021). Multi-omics network-based functional annotation of unknown Arabidopsis genes. Plant J. 108, 1193–1212. doi: 10.1111/tpj.15507
Des Marais, D. L., Hernandez, K. M., and Juenger, T. E. (2013). Genotype-by-environment interaction and plasticity: Exploring genomic responses of plants to the abiotic environment. Annu. Rev. Ecol. Evol. Syst. 44, 5–29. doi: 10.1146/annurev-ecolsys-110512-135806
DeWitt, T. J., Sih, A., and Wilson, D. S. (1998). Costs and limits of phenotypic plasticity. Trends. Ecol. Evol. 13, 77–81. doi: 10.1016/S0169-5347(97)01274-3
Donelson, J. M., Salinas, S., Munday, P. L., and Shama, L. N. (2017). Transgenerational plasticity and climate change experiments: Where do we go from here? Global Change Biol. 24, 13–34. doi: 10.1111/gcb.13903
Dorn, K. M., Fankhauser, J. D., Wyse, D. L., and Marks, M. D. (2015). A draft genome of field pennycress (Thlaspi arvense) provides tools for the domestication of a new winter biofuel crop. DNA Res. 22, 121–131. doi: 10.1093/dnares/dsu045
Dyer, A. R., Brown, C. S., Espeland, E. K., McKay, J. K., Meimberg, H., and Rice, K. J. (2010). SYNTHESIS: The role of adaptive trans-generational plasticity in biological invasions of plants. Evol. Appl. 3, 179–192. doi: 10.1111/j.1752-4571.2010.00118.x
Ehrenreich, I. M., and Pfennig, D. W. (2016). Genetic assimilation: A review of its potential proximate causes and evolutionary consequences. Ann. Bot. 117, 769–779. doi: 10.1093/aob/mcv130
Fitz-James, M. H., and Cavalli, G. (2022). Molecular mechanisms of transgenerational epigenetic inheritance. Nat. Rev. Genet. 23, 325–341. doi: 10.1038/s41576-021-00438-5
Forsman, A. (2015). Rethinking phenotypic plasticity and its consequences for individuals, populations and species. Heredity 115, 276–284. doi: 10.1038/hdy.2014.92
Galloway, L. F., and Etterson, J. R. (2007). Transgenerational plasticity is adaptive in the wild. Science 318, 1134–1136. doi: 10.1126/science.1148766
Geng, Y., Chang, N., Zhao, Y., Qin, X., Lu, S., Crabbe, M. J. C., et al. (2020). Increased epigenetic diversity and transient epigenetic memory in response to salinity stress in Thlaspi arvense. Ecol. Evol. 10, 11622–11630. doi: 10.1002/ece3.6795
Geng, Y., Gao, L., and Yang, J. (2013). “Epigenetic flexibility underlying phenotypic plasticity,” in Progress in Botany, eds U. Lüttge, W. Beyschlag, D. Francis, and J. Cushman (Berlin: Springer), 153–163.
Geng, Y., Guan, Y., Qiong, L., Lu, S., An, M., Crabbe, M. J. C., et al. (2021). Genomic analysis of field pennycress (Thlaspi arvense) provides insights into mechanisms of adaptation to high elevation. BMC Biol. 19:143. doi: 10.1186/s12915-021-01079-0
Halali, S., and Saastamoinen, M. (2022). Exploring links between climatic predictability and the evolution of within-and transgenerational plasticity. bioRxiv [preprint]. doi: 10.1101/2022.03.20.484890
He, L., Zhao, C., Zhang, Q., Zinta, G., Wang, D., Lozano-Durán, R., et al. (2021). Pathway conversion enables a double-lock mechanism to maintain DNA methylation and genome stability. Proc. Natl Acad. Sci. U. S. A 118:e2107320118. doi: 10.1073/pnas.2107320118
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
Herman, J. J., and Sultan, S. E. (2016). DNA methylation mediates genetic variation for adaptive transgenerational plasticity. Proc. Roy. Soc. B-Biol. Sci. 283:20160988. doi: 10.1098/rspb.2016.0988
Herman, J. J., Sultan, S. E., Horgan-Kobelski, T., and Riggs, C. (2012). Adaptive transgenerational plasticity in an annual plant: Grandparental and parental drought stress enhance performance of seedlings in dry soil. Integr. Comp. Biol. 52, 77–88. doi: 10.1093/icb/ics041
Jozefczak, M., Remans, T., Vangronsveld, J., and Cuypers, A. (2012). Glutathione is a key player in metal-induced oxidative stress defenses. Int. J. Mol. Sci. 13, 3145–3175. doi: 10.3390/ijms13033145
Kelly, S. A., Panhuis, T. M., and Stoehr, A. M. (2012). Phenotypic plasticity: Molecular mechanisms and adaptive significance. Compr. Physiol. 2, 1417–1439. doi: 10.1002/cphy.c110008
Kinraide, T. B. (2009). Improved scales for metal ion softness and toxicity. Environ. Toxicol. Chem. 28, 525–533. doi: 10.1897/08-208.1
Kooke, R., Johannes, F., Wardenaar, R., Becker, F., Etcheverry, M., Colot, V., et al. (2015). Epigenetic basis of morphological variation and phenotypic plasticity in Arabidopsis thaliana. Plant Cell 27, 337–348. doi: 10.1105/tpc.114.133025
Kuijper, B., Johnstone, R. A., and Townley, S. (2014). The evolution of multivariate maternal effects. PLoS Comput. Biol. 10:e1003550. doi: 10.1371/journal.pcbi.1003550
Kumar, D., Kumar, R., Baek, D., Hyun, T. K., Chung, W. S., Yun, D. J., et al. (2017). Arabidopsis thaliana receptor dead kinase1 functions as a positive regulator in plant responses to ABA. Mol. Plant 10, 223–243. doi: 10.1016/j.molp.2016.11.011
Kushwah, S., Banasiak, A., Nishikubo, N., Derba-Maceluch, M., Majda, M., Endo, S., et al. (2020). Arabidopsis XTH4 and XTH9 contribute to wood cell expansion and secondary wall formation. Plant Physiol. 182, 1946–1965. doi: 10.1104/pp.19.01529
Lachmann, M., and Jablonka, E. (1996). The inheritance of phenotypes: An adaptation to fluctuating environments. J. Theor. Biol. 181, 1–9. doi: 10.1006/jtbi.1996.0109
Laitinen, R. A., and Nikoloski, Z. (2019). Genetic basis of plasticity in plants. J. Exp. Bot. 70, 739–745. doi: 10.1093/jxb/ery404
Langfelder, P., and Horvath, S. (2008). WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 9:559. doi: 10.1186/1471-2105-9-559
Langfelder, P., Luo, R., Oldham, M. C., and Horvath, S. (2011). Is my network module preserved and reproducible? PLoS Comput. Biol. 7:e1001057. doi: 10.1371/journal.pcbi.1001057
Lang-Mladek, C., Popova, O., Kiok, K., Berlinger, M., Rakic, B., Aufsatz, W., et al. (2010). Transgenerational inheritance and resetting of stress-induced loss of epigenetic gene silencing in Arabidopsis. Mol. Plant 3, 594–602. doi: 10.1093/mp/ssq014
Leimar, O., and McNamara, J. M. (2015). The evolution of transgenerational integration of information in heterogeneous environments. Am. Nat. 185, E55–E69. doi: 10.1086/679575
Leng, Y., Li, Y., Ma, Y. H., He, L. F., and Li, S. W. (2021). Abscisic acid modulates differential physiological and biochemical responses of roots, stems, and leaves in mung bean seedlings to cadmium stress. Environ. Sci. Pollut. R. 28, 6030–6043. doi: 10.1007/s11356-020-10843-8
Li, G., Deng, Y., Geng, Y., Zhou, C., Wang, Y., Zhang, W., et al. (2017). Differentially expressed microRNAs and target genes associated with plastic internode elongation in Alternanthera philoxeroides in contrasting hydrological habitats. Front. Plant Sci. 8:2078. doi: 10.3389/fpls.2017.02078
Liu, K., and Sun, Q. (2021). Intragenic tRNA-promoted R-loops orchestrate transcription interference for plant oxidative stress responses. Plant Cell 33, 3574–3591. doi: 10.1093/plcell/koab220
Loix, C., Huybrechts, M., Vangronsveld, J., Gielen, M., Keunen, E., and Cuypers, A. (2017). Reciprocal interactions between cadmium-induced cell wall responses and oxidative stress in plants. Front. Plant Sci. 8:1867. doi: 10.3389/fpls.2017.01867
Martin, S. R., Llugany, M., Barceló, J., and Poschenrieder, C. (2012). Cadmium exclusion a key factor in differential Cd-resistance in Thlaspi arvense ecotypes. Biol. Plantarum. 56, 729–734. doi: 10.1007/s10535-012-0056-8
McGinn, M., Phippen, W. B., Chopra, R., Bansal, S., Jarvis, B. A., Phippen, M. E., et al. (2019). Molecular tools enabling pennycress (Thlaspi arvense) as a model plant and oilseed cash cover crop. Plant Biotechnol. J. 17, 776–788. doi: 10.1111/pbi.13014
Mishra, P., and Sharma, P. (2019). “Superoxide Dismutases (SODs) and their role in regulating abiotic stress induced oxidative stress in plants,” in Reactive oxygen, nitrogen and sulfur species in plants: Production, metabolism, signaling and defense mechanisms, eds M. Hasanuzzaman, V. Fotopoulos, K. Nahar, and M. Fujita (Hoboken: John Wiley and Sons Ltd), 53–88.
Mitich, L. W. (1996). Field pennycress (Thlaspi arvense L.)—the stinkweed. Weed technol. 10, 675–678. doi: 10.1017/S0890037X00040604
Moeder, W., Del Pozo, O., Navarre, D. A., Martin, G. B., and Klessig, D. F. (2007). Aconitase plays a role in regulating resistance to oxidative stress and cell death in Arabidopsis and Nicotiana benthamiana. Plant Mol. Biol. 63, 273–287. doi: 10.1007/s11103-006-9087-x
Montaldo, P., Cunnington, A., Oliveira, V., Swamy, R., Bandya, P., Pant, S., et al. (2020). Transcriptomic profile of adverse neurodevelopmental outcomes after neonatal encephalopathy. Sci. Rep. 10:13100. doi: 10.1038/s41598-020-70131-w
Monteiro, M. S., and Soares, A. M. (2021). Physiological and Biochemical Effects of Cd Stress in Thlaspi arvense L—A Non-Accumulator of Metals. Arch. Environ. Con. Tox. 81, 285–292. doi: 10.1007/s00244-021-00873-9
Murren, C. J., Auld, J. R., Callahan, H., Ghalambor, C. K., Handelsman, C. A., Heskel, M. A., et al. (2015). Constraints on the evolution of phenotypic plasticity: Limits and costs of phenotype and plasticity. Heredity 115, 293–301. doi: 10.1038/hdy.2015.8
Nijhout, H. F., Kudla, A. M., and Hazelwood, C. C. (2021). Genetic assimilation and accommodation: Models and mechanisms. Curr. Top. Dev. Biol. 141, 337–369. doi: 10.1016/bs.ctdb.2020.11.006
Paparella, C., Savatin, D. V., Marti, L., De Lorenzo, G., and Ferrari, S. (2014). The Arabidopsis LYSIN MOTIF-CONTAINING RECEPTOR-LIKE KINASE3 regulates the cross talk between immunity and abscisic acid responses. Plant Physiol. 165, 262–276. doi: 10.1104/pp.113.233759
Pérez-Ruiz, J. M., Spinola, M. C., Kirchsteiger, K., Moreno, J., Sahrawy, M., and Cejudo, F. J. (2006). Rice NTRC is a high-efficiency redox system for chloroplast protection against oxidative damage. Plant Cell 18, 2356–2368. doi: 10.1105/tpc.106.041541
Pieterse, C. M. (2012). Prime time for transgenerational defense. Plant Physiol. 158, 545–545. doi: 10.1104/pp.112.900430
Pigliucci, M. (2005). Evolution of phenotypic plasticity: Where are we going now? Trends Ecol. Evol. 20, 481–486. doi: 10.1016/j.tree.2005.06.001
Quadrana, L., and Colot, V. (2016). Plant transgenerational epigenetics. Annu. Rev. Genet. 50, 467–491. doi: 10.1146/annurev-genet-120215-035254
Rivera, H. E., Aichelman, H. E., Fifer, J. E., Kriefall, N. G., Wuitchik, D. M., Wuitchik, S. J., et al. (2021). A framework for understanding gene expression plasticity and its influence on stress tolerance. Mol Ecol. 30, 1381–1397. doi: 10.1111/mec.15820
Rosado, A., Schapire, A. L., Bressan, R. A., Harfouche, A. L., Hasegawa, P. M., Valpuesta, V., et al. (2006). The Arabidopsis tetratricopeptide repeat-containing protein TTL1 is required for osmotic stress responses and abscisic acid sensitivity. Plant Physiol. 142, 1113–1126. doi: 10.1104/pp.106.085191
Salinas, S., Brown, S. C., Mangel, M., and Munch, S. B. (2013). Non-genetic inheritance and changing environments. Non-Genetic Inherit. 1, 38–50. doi: 10.2478/ngi-2013-0005
Sandner, T. M., van Braak, J. L., and Matthies, D. (2018). Transgenerational plasticity in Silene vulgaris in response to three types of stress. Plant Biol. 20, 751–758. doi: 10.1111/plb.12721
Schlichting, C. D. (1986). The Evolution of Phenotypic Plasticity in Plants. Ann. Rev. Ecol. Evol. S. 17, 667–693. doi: 10.1146/annurev.es.17.110186.003315
Sedbrook, J. C., Phippen, W. B., and Marks, M. D. (2014). New approaches to facilitate rapid domestication of a wild plant to an oilseed crop: Example pennycress (Thlaspi arvense L.). Plant Sci. 227, 122–132. doi: 10.1016/j.plantsci.2014.07.008
Sharma, R., Singh, G., Bhattacharya, S., and Singh, A. (2018). Comparative transcriptome meta-analysis of Arabidopsis thaliana under drought and cold stress. PLoS One 13:e0203266. doi: 10.1371/journal.pone.0203266
Shi, D., Zhao, C., Chen, Y., Ding, J., Zhang, L., and Chang, Y. (2020). Transcriptomes shed light on transgenerational and developmental effects of ocean warming on embryos of the sea urchin Strongylocentrotus intermedius. Sci. Rep. 10:7931. doi: 10.1038/s41598-020-64872-x
Shi, W., Chen, X., Gao, L., Xu, C. Y., Ou, X., Bossdorf, O, et al. (2019). Transient stability of epigenetic population differentiation in a clonal invader. Front. Plant Sci. 9:1851. doi: 10.3389/fpls.2018.01851
Sobral, M., Neylan, I. P., Narbona, E., and Dirzo, R. (2021). Transgenerational plasticity in flower color induced by caterpillars. Front. Plant Sci. 12:617815. doi: 10.3389/fpls.2021.617815
Sommer, R. J. (2020). Phenotypic plasticity: From theory and genetics to current and future challenges. Genetics 215, 1–13. doi: 10.1534/genetics.120.303163
Sultan, S. E. (2000). Phenotypic plasticity for plant development, function and life history. Trends Plant Sci. 5, 537–542. doi: 10.1016/S1360-1385(00)01797-0
Sultan, S. E., Barton, K., and Wilczek, A. M. (2009). Contrasting patterns of transgenerational plasticity in ecologically distinct congeners. Ecology 90, 1831–1839. doi: 10.1890/08-1064.1
Sun, L., Wang, J., Song, K., Sun, Y., Qin, Q., and Xue, Y. (2019). Transcriptome analysis of rice (Oryza sativa L.) shoots responsive to cadmium stress. Sci. Rep. 9:10177. doi: 10.1038/s41598-019-46684-w
Tolra, R., Pongrac, P., Poschenrieder, C., Vogel-Mikuš, K., Regvar, M., and Barcelo, J. (2006). Distinctive effects of cadmium on glucosinolate profiles in Cd hyperaccumulator Thlaspi praecox and non-hyperaccumulator Thlaspi arvense. Plant Soil 288, 333–341. doi: 10.1007/s11104-006-9124-1
Uller, T., Moczek, A. P., Watson, R. A., Brakefield, P. M., and Laland, K. N. (2018). Developmental bias and evolution: A regulatory network perspective. Genetics 209, 949–966. doi: 10.1534/genetics.118.300995
Verhoeven, K. J., Jansen, J. J., Van Dijk, P. J., and Biere, A. (2010). Stress-induced DNA methylation changes and their heritability in asexual dandelions. New Phytol. 185, 1108–1118. doi: 10.1111/j.1469-8137.2009.03121.x
Vu, W. T., Chang, P. L., Moriuchi, K. S., and Friesen, M. L. (2015). Genetic variation of transgenerational plasticity of offspring germination in response to salinity stress and the seed transcriptome of Medicago truncatula. BMC Evol. Biol. 15:59. doi: 10.1186/s12862-015-0322-4
Waddington, C. H. (1961). Genetic assimilation. Adv. Genet. 10, 257–293. doi: 10.1016/S0065-2660(08)60119-4
Walsh, M. R., Castoe, T., Holmes, J., Packer, M., Biles, K., Walsh, M., et al. (2016). Local adaptation in transgenerational responses to predators. P. Roy. Soc. B Biol. Sci. 283:20152271. doi: 10.1098/rspb.2015.2271
Wang, Q., Lu, F., and Lan, R. (2017). RNA-sequencing dissects the transcriptome of polyploid cancer cells that are resistant to combined treatments of cisplatin with paclitaxel and docetaxel. Mol. Biosyst. 13, 2125–2134. doi: 10.1039/C7MB00334J
Warwick, S. I., Francis, A., and Susko, D. J. (2002). The biology of Canadian weeds. 9. Thlaspi arvense L.(updated). Can. J. Plant. Sci. 82, 803–823.
Ye, J., Zhang, Y., Cui, H., Liu, J., Wu, Y., Cheng, Y., et al. (2018). WEGO 2.0: A web tool for analyzing and plotting GO annotations, 2018 update. Nucleic Acids Res. 46, 71–75. doi: 10.1093/nar/gky400
Zheng, M., Zhu, C., Yang, T., Qian, J., and Hsu, Y. F. (2020). GSM2, a transaldolase, contributes to reactive oxygen species homeostasis in Arabidopsis. Plant Mol. Biol. 104, 39–53. doi: 10.1007/s11103-020-01022-x
Zhou, Q., Guo, J. J., He, C. T., Shen, C., Huang, Y. Y., Chen, J. X., et al. (2016). Comparative transcriptome analysis between low-and high-cadmium-accumulating genotypes of pakchoi (Brassica chinensis L.) in response to cadmium stress. Environ. Sci. Technol. 50, 6485–6494. doi: 10.1021/acs.est.5b06326
Keywords: field pennycress, phenotypic plasticity, maternal effect, transgenerational plasticity, transcriptional memory
Citation: Li G, Zhao Y, Liu F, Shi M, Guan Y, Zhang T, Zhao F, Qiao Q and Geng Y (2022) Transcriptional memory of gene expression across generations participates in transgenerational plasticity of field pennycress in response to cadmium stress. Front. Plant Sci. 13:953794. doi: 10.3389/fpls.2022.953794
Received: 26 May 2022; Accepted: 29 August 2022;
Published: 30 September 2022.
Edited by:
Bi-Cheng Dong, Beijing Forestry University, ChinaReviewed by:
Mozhu Wang, Institute of Botany (CAS), ChinaYuepeng Song, Beijing Forestry University, China
Copyright © 2022 Li, Zhao, Liu, Shi, Guan, Zhang, Zhao, Qiao and Geng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Fangqing Zhao, zhfq@biols.ac.cn; Qin Qiao, qiaoqin@ynu.edu.cn; Yupeng Geng, ypgeng@ynu.edu.cn