Genetic assimilation of ancestral plasticity during parallel adaptation to zinc contamination in Silene uniflora

Phenotypic plasticity in ancestral populations is hypothesized to facilitate adaptation, but evidence is piecemeal and often contradictory. Further, whether ancestral plasticity increases the probability of parallel adaptive changes has not been explored. The most general finding is that ancestral responses to a new environment are reversed following adaptation (known as reversion). We investigated the contribution of ancestral plasticity to adaptive evolution of gene expression in two independently evolved lineages of zinc-tolerant Silene uniflora. We found that the general pattern of reversion is driven by the absence of a widespread stress response in zinc-adapted plants compared with zinc-sensitive plants. We show that ancestral plasticity that moves expression closer to the optimum value in the new environment influences the evolution of gene expression among genes that are likely to be involved in adaptation and increases the chance that genes are recruited repeatedly during adaptation. However, despite convergence in gene expression levels between independently adapted lineages, ancestral plasticity does not influence how similar expression values of adaptive genes become. Surprisingly, we also observed that ancestral plasticity that increases fitness often becomes genetically determined and fixed, that is, genetically assimilated. These results emphasize the important role of ancestral plasticity in parallel adaptation.

Phenotypic plasticity in ancestral populations is hypothesized to facilitate adaptation, but evidence is piecemeal and often contradictory. Further, whether ancestral plasticity increases the probability of parallel adaptive changes has not been explored. The most general finding is that ancestral responses to a new environment are reversed following adaptation (known as reversion). We investigated the contribution of ancestral plasticity to adaptive evolution of gene expression in two independently evolved lineages of zinc-tolerant Silene uniflora. We found that the general pattern of reversion is driven by the absence of a widespread stress response in zinc-adapted plants compared with zinc-sensitive plants. We show that ancestral plasticity that moves expression closer to the optimum value in the new environment influences the evolution of gene expression among genes that are likely to be involved in adaptation and increases the chance that genes are recruited repeatedly during adaptation. However, despite convergence in gene expression levels between independently adapted lineages, ancestral plasticity does not influence how similar expression values of adaptive genes become. Surprisingly, we also observed that ancestral plasticity that increases fitness often becomes genetically determined and fixed, that is, genetically assimilated. These results emphasize the important role of ancestral plasticity in parallel adaptation.
The contributions of determinism and contingency in shaping evolution are hotly debated [1][2][3] . Whether repeated adaptation to the same environment results in similar changes at the molecular level is key to understanding this balance 1,4-6 , as well as the predictability of future responses to environmental change 7 . Adaptation to novel environments often involves gene expression changes, but previous studies have found varying degrees of parallelism during repeated adaptation [8][9][10][11] .
These changes occur at various levels, including in the overlap of shared differentially expressed genes, fold changes of these genes or final expression levels 9,12 . Understanding the mechanisms that influence the extent of parallelism is an important step in predicting evolutionary responses to new environmental challenges 6,7,13 .
Phenotypic plasticity in ancestral populations (ancestral plasticity) is suspected to play a role in facilitating adaptation to new Article https://doi.org/10.1038/s41559-022-01975-w (5) the relationship between ancestral plasticity and convergent gene expression changes. In so doing, we establish the extent to which rapid adaptation is shaped by constraint and plasticity, disentangling the influence of general stress responses versus adaptive responses on patterns of reversion and reinforcement.
Heavy metals are highly phytotoxic and high concentrations of zinc have a considerable impact on growth and fitness of coastal populations of S. uniflora [38][39][40] . Transcriptome-wide ancestral plasticity (that is, the response to zinc in sensitive populations) was dominated by a general and widespread stress response. In total, 51.1% of the transcriptome (14,327 genes) was differentially expressed in both sensitive populations between treatments, with an overwhelming majority being shared across populations (Extended Data Fig. 1a). Shared upregulated genes were enriched for 15 gene ontology (GO) terms related to stress (Supplementary File 1). Further, the major difference in expression between susceptible and tolerant populations was the lack of this extreme response to zinc stress in tolerant populations. Only 223 genes were differentially expressed between treatments in both tolerant populations. In the zinc treatment, 9,549 genes were differentially expressed between tolerant and sensitive populations in both pairs (Extended Data Fig. 1b), which were enriched for 12 stress-related GO terms (Supplementary File 2). Of these genes, 87.0% were ancestrally plastic (that is, also differentially expressed between treatments in both sensitive populations), but only 1.4% showed derived plasticity (DP; that is, were also differentially expressed between treatments in both tolerant populations; Extended Data Fig. 2c). This reveals a substantial disruption to transcription in sensitive plants, consistent with the broad impact of zinc toxicity on cellular processes 41 . It also indicates that, in general, greater transcriptomic perturbations in ancestral populations exposed to new environments may be driven by general stress responses 20,30,36,37 .

Rapid evolution of highly parallel gene expression changes
Silene uniflora has independently colonized mines and evolved tolerance to the very high levels of zinc (2,400-48,100 ppm) in the contaminated soils [38][39][40] . Given that this phenotype has evolved in parallel due to strong selection, we also expected a component of the transcription profiles to show parallel changes in tolerant populations. In the control treatment, principal component analysis (PCA) of transcriptome-wide gene expression levels separated populations by zinc tolerance (that is, tolerant versus sensitive) on PC1 and by geographic origin (that is, T1 and S1 versus T2 and S2) on PC2 ( Fig. 1b and Extended Data Fig. 3a,b). Within-population variation was low relative to between populations/treatments (Extended Data Fig. 3c). In these benign conditions, the trajectories of whole transcriptome evolution were divergent and almost orthogonal, rather than parallel (sensu ref. 6 ).
In total, 2,119 and 2,884 genes were differentially expressed in control conditions between T1 and S1, and T2 and S2, respectively, of which 532 were shared (Extended Data Fig. 2d). We categorized 400 of these shared genes as displaying parallel constitutive evolutionary changes of expression (CEC genes); these were differentially expressed in both tolerant-sensitive pairs 'and' had expression differences in the same direction (that is, increased or decreased expression in both T1 versus S1 and T2 versus S2; Extended Data Figs. 3 and 4a). Genes with expression shifts in the same direction are more likely to be the result of parallel adaptation across the mines [8][9][10][11] and 400 genes represents a greater overlap than expected by chance (one-sided Fisher's Exact test, odds ratio = 2.2, P < 2.2 × 10 −16 ). The degree of similarity in gene expression levels between populations can be quantified by comparing the absolute per-gene log 2 -transformed shrunken fold changes (FC; see Methods for rationale). For a set of genes, a small median |FC| indicates high expression similarity between a pair of populations. In control conditions, transcriptome-wide expression levels of tolerant populations were less similar than the coastal populations were to each other (|FC| S1-S2 = 0.056 versus |FC| T1-T2 = 0.12; two-sided paired Wilcoxon environments [14][15][16] . In addition to generally preserving the genetic variability of a colonizing population 17 , plastic responses to new environments could provide the basis for adaptation by moving the trait values in some individuals closer to the new local optimum 18 . Beneficial plasticity of this kind could be retained in locally adapted populations or genetically assimilated and canalized into constitutive expression differences 19 . Alternatively, ancestral plasticity that takes expression levels further away from the new optimum is potentially maladaptive and could hinder adaptation to the novel environment 20,21 .
Current evidence suggests a variety of possible impacts of ancestral plasticity on adaptation 12,[21][22][23][24] , but the relationship between plasticity and evolutionary parallelism has received limited attention 6,25 . Other properties of gene expression in ancestral populations, such as ancestral expression level or tissue expression location, are associated with increased co-option and potentially parallelism 26,27 . If phenotypic plasticity substantially facilitates the repurposing of traits during adaptation 28 , then beneficial plasticity may result in greater parallelism than when plasticity is maladaptive.
Previous studies have generally found that most ancestral plasticity across transcriptomes is reversed in derived populations, that is, it takes expression values further from the new optimum 22,[29][30][31] (although there are exceptions 32,33 ). However, there are examples of ancestral plasticity in particular genes or traits facilitating subsequent adaptation 18,21,34,35 . Most expression studies on the topic examine transcriptome-wide patterns in ancestrally plastic genes, rarely considering whether genes involved in evolutionary adaptation to the new environment are more likely to have possessed beneficial ancestral plasticity, when compared with the whole transcriptome 20,22,30,31,36,37 . Transcriptome-wide assessments include changes that may not directly contribute to adaptation (in the evolutionary sense), such as those stemming from general stress responses. As a result, estimates of the contribution of ancestral plasticity to adaptation may be distorted in whole transcriptome analysis.
Here we investigate the relationship between ancestral plasticity, adaptation and parallelism using independently evolved lineages of zinc-tolerant Silene uniflora from contaminated metal mines and local zinc-sensitive coastal populations 38 . In this species, ancestral coastal populations have repeatedly colonized contaminated mine soils throughout Great Britain and Ireland over the past 250 years 39 , producing locally adapted populations that can grow at high concentrations of zinc [38][39][40] . As a result, expression differences between closely related mine-coast pairs should resemble the expression differences between the current mine populations and their coastal ancestors. Common changes across replicates are likely to represent adaptive changes rather than drift 6 . Extant coastal populations also provide an approximation of the ancestral plastic response to zinc. This system provides an ideal opportunity to investigate the role of ancestral plasticity in adaptation across multiple evolutionary replicates.

Results and discussion
Zinc-tolerant populations of S. uniflora largely exclude zinc from their shoots, preferentially accumulating zinc in their roots 39,40 . We quantified gene expression in the roots of two independently derived, zinc-tolerant populations from geographically distant, derelict mines (T1, England; T2, Wales) and their nearest and most closely related zinc-sensitive coastal populations (S1 and S2; Fig. 1a). Extant zinc-sensitive coastal populations acted as proxies for ancestral expression. We exposed clones of the same individuals to two treatment conditions (control or zinc-contaminated) and collected RNA-seq data from the roots of the experimental plants. Our experimental design allowed us to quantify: (1) the ancestral plastic response to zinc contamination; (2) the extent of convergent gene expression changes during rapid parallel adaptation; (3) the evolutionary response to ancestral plasticity at a transcriptome-wide level; (4) whether the evolutionary response differs for genes plausibly involved in adaptation; and  The CEC genes had similar expression values in sensitive populations (CEC|FC| S1-S2 = 0.077), but expression was also highly similar in tolerant populations (CEC|FC| T1-T2 = 0.12), despite substantial expression divergence and genome-wide genetic differentiation from the nearest coastal populations ( Fig. 1c; mean F ST = 0.36 between susceptible and tolerant populations 38 ). In other words, for the 400 CEC genes, parallel evolution in mine populations produced expression similarity comparable to that observed between sensitive populations, which is the product of shared ancestry, gene flow, drift and selection. Unlike in the control treatment, there was a higher degree of parallelism in the evolved response to zinc treatment across the whole transcriptome (solid black arrows, Fig. 1d). However, this was largely driven by the widespread transcriptomic response of the sensitive plants, with a less dramatic shift in expression of tolerant populations in zinc versus control treatments. Genes with significant expression responses to zinc in both tolerant populations and that, in both tolerant populations, show expression differences from susceptible populations in either the control or the zinc treatment (or both; Extended Data Fig. 3), are likely to play some role in zinc tolerance (hereafter called DP genes). Of the 245 and 653 genes with this expression pattern in T1 and T2, respectively, 137 were shared. This is a greater overlap than expected by chance (one-sided Fisher's Exact Test, odds ratio = 66.8, P < 2.2 × 10 −16 ). This level of parallelism is high compared with other systems, such as repeated adaptation to elevation 42,43 . This difference may be due to the strength and specificity of selection that metal toxicity imposes, rather than more multifarious selection along elevation gradients. These shared genes had highly correlated expression shifts (log 2 fold changes between treatments; linear model slope = 0.84, t = 26, P < 2.2 x−10 −16 , adjusted R 2 = 0.83; Fig. 1e). Many DP genes (83%) were also differentially expressed between treatments in both susceptible populations and may constitute a stress response that is partially inherited from their coastal ancestors; indeed, 'response to stress' was the most highly enriched GO term for DP genes with ancestral plasticity ( Supplementary Files 3 and 4). Nevertheless, there were also convergent changes in expression levels in these genes between tolerant populations. Expression profiles for DP genes were similar in the control treatment, but when exposed to zinc, evolved responses were almost perfectly parallel in tolerant populations ( Fig. 1f and Extended Data Fig. 5), consistent with previous studies indicating that phenotypic plasticity can result in increased phenotypic parallelism 25 .
There were three times as many genes with constitutive differences between the sensitive ecotype and the tolerant ecotype (CEC genes) as genes with DP. In the literature, there is considerable variability across taxa in the ratios of constitutive to plastic differences associated with local adaptation 19,30,31,[44][45][46][47][48] . This may be a function of the degree to which a stressor varies in strength temporally and spatially within a habitat 39,49,50 . However, CEC genes that do not respond to zinc could be involved in zinc tolerance and/or adaptation to other aspects of the mine environment (for example, exposure, water availability and so on), which may also explain this difference. Overall, these results suggest that highly parallel patterns of differential gene expression across evolutionary replicates can be acquired very early in adaptation and over very short timescales. This is true for both the identity of the genes and the magnitude of expression shifts. Papadopulos et al. 38 identified both shared and non-shared genetic changes across mine-adapted populations and concluded that there may be a highly polygenic basis to adaptation. These evolved expression shifts could be caused by the same or different underlying genetic variants. The responsible variants may be cis-or trans-acting, may have arisen via gene duplications, and may either directly affect gene expression or target a few upstream regulators-we are unable to assess this from transcriptomic data alone. Regardless of the nature of the genetic changes that have occurred, they have produced remarkably similar gene expression across independent mine colonization events. Previous experimental evolution studies in Drosophila, Tribolium and Ipomoea have demonstrated the evolution of gene expression plasticity in response to heterogeneous environments within 22-130 generations 20,33,36,51 . We demonstrate that this can also occur in wild plant populations in comparable timeframes and is repeatable between independent colonizations of a novel habitat.

Convergent zinc tolerance pathways
Examining sets of shared genes with expression patterns consistent with a role in adaptation sheds light on the mechanisms underlying zinc tolerance. CEC genes were enriched for 222 GO terms, including terms associated with metal tolerance (for example, zinc ion transport; see Supplementary File 5). This included homologues of A. thaliana Zinc Transporter 1 ZIP1, which encodes a protein that mediates the uptake of zinc from the rhizosphere 52 , Heavy Metal Atpase 2 (HMA2, encoding a plasma membrane protein that transports zinc from cells 53,54 ) and Metal Tolerance Protein 1 (ZAT, encoding a protein that sequesters zinc into vacuoles and controls zinc accumulation in roots 55,56 ). These are upregulated in zinc hyperaccumulators such as Arabidopsis halleri 57 and when overexpressed confer increased metal accumulation and tolerance 55,58,59 . The function of these genes is consistent with increased zinc accumulation in the roots of zinc-tolerant S. uniflora populations 38,40 .

Fig. 2 | Conceptual overview of evolutionary responses to ancestral plasticity.
When an ancestral population reaches a novel environment, an immediate PC moves the trait from an initial value of L o in the old environment to L p in the new environment. As populations adapt over time, a further EC shifts L p to a new value of L a . a, The evolutionary response to ancestral plasticity can be divided into three categories depending on the values of PC and EC. b-e, Cartoon representations of scenarios. Dashed line represents transition from ancestral to novel environment and associated trait shift, PC. b, Reinforcement occurs when the subsequent EC is in the same direction as PC. c, Overshooting occurs when PC has moved the trait value closer to the new optimum (that is, L a is closer to L p than L o ). In this scenario, EC is in the opposite direction to PC, but |EC| < 0.5 × |PC|. d,e, Reversion occurs when the optimum in the new habitat is nearer to the value of the unstressed ancestor in its home environment than the ancestor's response (that is, L a is closer to L o than L p ), so EC is in the opposite direction to PC, but |EC| < 0.5 × |PC|. Reversion can include the restoration of the ancestral state in the old environment (|EC| = |PC|) (e) or move beyond this value in the opposite direction (|EC| >|PC|) (d). Reinforcement and overshooting suggest that ancestral plasticity was adaptive, whereas reversion indicates it was maladaptive.  60 ). Overexpression of GSTs results in enhanced zinc and cadmium tolerance 61,62 . These GSTs were also differentially expressed between conditions in susceptible populations, further hinting at a role of ancestral plasticity in adaptation.
These results indicate that genes that have been repeatedly recruited for a role in zinc tolerance across multiple species 41 have also undergone repeated gene expression changes in zinc-tolerant populations over a few hundred generations.

Ancestral plasticity is generally reversed during adaptation
To understand the relationship between ancestral plasticity and adaptation, an established approach is to investigate mean differences in about the role of ancestral plasticity during adaptation, we can compare the direction and magnitude of the initial plastic response of an ancestral population when it is exposed to a new environment (ancestral plasticity/plastic change, PC = L p −L o ) with the subsequent change in expression between the ancestral population and an adapted population in the new environment 12,22 (evolutionary change, EC = L a −L p ). The relationship between PC and EC (that is, the evolutionary response to ancestral plasticity) can be characterized in three ways: (1) 'reinforcement', where the initial PC and subsequent EC both move expression in the same direction towards the new optimum (Fig. 2a,b); (2) 'overshooting', where PC takes expression beyond the new optimum and EC then adjusts expression in the opposite direction (Fig. 2a,c); and (3) 'reversions', where the new optimum is closer to the level of the ancestor in its home environment, so EC largely counteracts the change observed in PC (Fig. 2a,d,e). During both reinforcement and overshooting, the ancestral PC moves expression closer to the new optimum, so both can be interpreted as ancestral plasticity facilitating adaptation to the new environment. Conversely, reversions are likely to be the outcome when ancestral plasticity is maladaptive. We evaluated the degree of reversion, reinforcement and overshooting in our transcriptome dataset. To avoid spurious assignment to these categories resulting from very small expression changes, only genes showing substantial changes in PC and EC (those that were differentially expressed: (1) in susceptible plants between conditions (PC); and (2) between mine and coast plants in zinc (EC); see Methods) were placed into these three categories (12,679 genes in total; Fig. 2a). To establish the general pattern of evolutionary responses to ancestral plasticity, we first considered these patterns transcriptome-wide. Across the entire transcriptome, 95.2% of genes showed reversion, with only 1.1% showing reinforcement and 3.7% showing overshooting. Therefore, in the vast majority of cases, ancestral plasticity does not move expression closer to the new optimum ( Fig. 3a and Extended Data Fig. 6a). Our transcriptome-wide results are consistent with previous studies in animals and microorganisms, which generally find that reversion is dominant 22,30,36 .
The majority of genes displaying substantial PC and EC across the transcriptome undergo high stress responses in sensitive plants in zinc and remain at unstressed levels in tolerant populations in the zinc treatment. As such, most of the transcriptome is not directly involved in adaptation. Examining the evolutionary response to ancestral plasticity across the predominantly non-adaptive transcriptome provides an indication of the probability that an ancestral plastic response moves expression closer to the new optimum in the zinc-contaminated environment. The large number of subsequent evolutionary reversions indicates that this probability is low (Fig. 3a). Whether this probability increases for genes directly involved in adaptation is more informative for understanding the role of plasticity in adaptation. The DP and CEC genes plausibly have a role in repeated adaptation to zinc as they are consistently recruited across parallel replicates [8][9][10][11] , but they account for only 1.8% of the transcriptome. DP and CEC genes are unlikely to include all genes that are involved in adaptation to zinc contamination-some may only be important in a single population or not detected under the framework applied here. Nevertheless, the DP and CEC sets are likely to be enriched for genes involved in adaptation, making them informative as to whether adaptive genes have different responses to ancestral plasticity versus the largely non-adaptive background transcriptomic response.

Ancestral plasticity is less likely to be reversed in adaptive genes
To understand whether ancestral plasticity facilitates adaptive evolution, we considered the proportion of genes undergoing reversion, reinforcement and overshooting in the DP and CEC gene sets. Among DP genes with substantial PC and EC (82.5% of the total), 79.6% underwent reversion, 3.5% reinforcement and 16.8% overshooting (Fig. 3b and Extended Data Fig. 6b). The higher proportion of overshooting in DP genes relative to the whole transcriptome (16.8% versus 3.7%; P = 2.48 × 10 −7 ; binomial two-sided test) suggests that DP genes may carry a fitness cost for being expressed at an inappropriate/inaccurate level for a given concentration of zinc and have fine-tuned the ancestral level of plasticity.
Adaptation to zinc contamination has also produced constitutive gene expression differences between tolerant and sensitive populations in the absence of zinc (CEC genes). Ancestral plasticity may facilitate the evolution of differences by moving expression closer to the new optimum, which could then lead to constitutive adaptive changes 19 . Among CEC genes showing substantial PC and EC (56.2% of the total), only 68.4% show signs of reversion, with 28.0% undergoing reinforcement and 3.6% overshooting (Fig. 3c and Extended Data Fig. 6c). This is significantly higher than in either DP genes (28.0% versus 3.5%, P < 2.2 × 10 −16 ; binomial two-sided test) or transcriptome-wide (28.0% versus 1.1%, P < 2.2 × 10 −16 ; binomial two-sided test). We conducted parametric bootstrapping as recommended in ref. 29 to reduce bias stemming from the presence of L p in calculations of PC and EC (Fig. 2). Bootstrapping (see Methods) generally increased the proportion of reversions and reduced the proportion of overshooting, but substantial enrichment of reinforcement in the CEC genes remained (Supplementary Table 1). The increase in reinforcement among CEC genes suggests that ancestral plastic responses make an important contribution during adaptation and may be genetically assimilated in the process.
Here we define genetic assimilation as: when a trait with an environmentally induced response that increases fitness becomes genetically determined and canalized (that is, there is a loss of plasticity) 15,63-65 . Of the 400 CEC genes, 310 are not zinc-responsive in either tolerant population but display substantial PC; these have been repeatedly canalized. Other definitions of genetic assimilation only include cases where the derived trait value is similar to the ancestral value in the new environment 30,66 . Of the 310 canalized genes, 114 do not display substantial EC, that is, the ancestral response was close to the new optimum. These included HMA2 and ZAT (see Zinc tolerance pathways section). For an additional 69 of these canalized CEC genes, the ancestral response took expression closer to the new optimum (that is, overshooting or reinforcement). Altogether 183 genes have undergone genetic assimilation (46% of CEC genes, 0.7% of the transcriptome), emphasizing the importance of ancestral plasticity during rapid adaptation to new environments.
Other studies have looked for a role for ancestral plasticity in producing constitutive expression differences by establishing a positive correlation between ancestral plasticity (which they define as L p /L o ) and evolutionary change in control conditions (defined as L c /L o , where L c is the level of the adapted population in the ancestral environment 20,37 ). However, the common denominator of L o in both variables would tend to produce a positive correlation 67 , potentially making these results unreliable. Most constitutive differences have been found to have evolutionary changes in the opposite direction to ancestral plasticity (reversion and overshooting were not distinguished) 12 , but whether there was an increase compared with the transcriptome-wide pattern was not assessed. Here, we demonstrated that although most ancestral plasticity is maladaptive, ancestral plasticity that can move expression closer to the new optimum contributes to adaptation.

Ancestral plasticity is not necessary for substantial gene expression convergence
Given this evidence of ancestral plasticity contributing to adaptation, the question of its importance for parallelism in adaptation arises. Plasticity may also increase the propensity of genes to be repeatedly recruited during adaptation. Unlike the shared CEC genes, which had relatively low rates of reversion (68.4%), genes differentially expressed in the control in only one population pair were more likely to show reversion (74.8% in T1/S1 and 80.5% in T2/S2, Supplementary Table 2). Article https://doi.org/10.1038/s41559-022-01975-w In other words, genes repeatedly recruited during adaptation are more likely to have had ancestral plasticity that moved expression closer to the new optimum, than those that were only recruited in one event.
In addition to affecting gene recruitment, ancestral plasticity may also affect the degree of expression convergence in repeatedly recruited genes. Comparisons of expression levels for DP or CEC genes that had ancestral plasticity versus those without returned no significant differences (Fig. 3d; DP genes, |FC| NOPLAST = 0.33, |FC| PLAST = 0.20, two-sided Wilcoxon signed-rank test, W = 815, P = 0.18; CEC genes |FC| NOPLAST = 0.15, |FC| PLAST = 0.19, W = 1.2 × 10 4 , P = 0.86). Genes lacking ancestral plasticity can rapidly evolve plastic responses with comparable expression convergence to ancestrally plastic genes. In summary, ancestral plasticity facilitates the repeated recruitment of genes but does not necessarily lead to greater convergence in expression levels during adaptation.

Experimental considerations
Despite our modest sample size, we controlled for between-treatment expression variation that might stem from genetic differences between individuals within a population by using clones paired across treatments. Further, within-population relative to between-population variability is very low (Fig. 1b and Extended Data Fig. 2c). We acknowledge that the between-residuals effects among cuttings from the same individual may not be zero, but these are likely to be very small given the common starting conditions and identical genotype. Removing the genotype term that pairs individuals across treatments did not alter the observed patterns (Supplementary Table 5). We also could not directly observe expression in the mine populations' ancestors, but the very recent colonization from coasts means responses in these extant coastal populations are likely to be very similar to the ancestral plastic response. Although some expression shifts might have taken place in the coastal populations since the mine populations diverged, these differences are unlikely to explain the patterns we observed consistently across the replicated events. Additionally, the design limits, but does not eliminate, maternal effects on expression. As such, it is possible that residual maternal effects might have affected some individual genes; however, this would not account for the patterns in large groups of genes and the relationship between putatively adaptive genes and ancestral plasticity. Finally, our experiment only considers gene expression responses; other forms of gene regulation, or mutational effects besides transcription (for example, coding sequence change) could also be important in zinc tolerance evolution.

Conclusions
Highly parallel gene expression phenotypes have evolved in S. uniflora during the repeated colonization of zinc-contaminated mines, despite the short timescales involved and a lack of gene flow between the tolerant populations 38 . By using coastal relatives to approximate the ancestral state, we show that genes displaying beneficial patterns of ancestral plasticity are overrepresented in these highly parallel gene sets, indicating that ancestral plasticity facilitates repeated adaptation to novel environments. The results of our experiment and others confirm that most ancestral plasticity is non-adaptive 22,31,36 . Nevertheless, the considerable proportion of fixed adaptive differences that co-opt ancestral plastic responses suggests that it is a major force in rapid adaptation. Despite a role for ancestral plasticity in enhancing the recruitment of genes, it does not result in an increased level of phenotypic convergence at the level of gene expression compared with genes showing no significant ancestral plasticity. In other words, ancestral plasticity only facilitates parallel evolution at certain levels of biological organization. Overall, our results indicate that genetic assimilation and modification of ancestral plastic responses play an important role in adaptation to novel environments and may be partially responsible for parallelism in gene expression during local adaptation.

Plant materials and experimental procedure
Populations T1, S1, T2 and S2 correspond to WWA-M, WWA-C, ENG-M and ENG-C in ref. 38 ; seeds were collected as described in that study. Three seeds per population, collected from different mothers, were germinated and cuttings propagated at 10 weeks (see Supplementary Methods for conditions). Cuttings were transferred to six deep water culture tanks containing dilute Hoagland's solution. Susceptible and tolerant populations grow normally in these benign conditions [38][39][40] . Cuttings from each individual were included in each tank and there was approximately equal representation of populations per tank. The use of cuttings should reduce any maternal effects from differences in resource allocation to seeds between populations. After 1 week of acclimation, the hydroponic solution was replaced with fresh solution in three tanks (control treatment) and the solution adjusted to 600 µM ZnSO 4 solution in the remaining three tanks (zinc treatment). Eight days later, roots from each individual cutting were flash frozen in liquid nitrogen and stored at −80 °C. For each individual within a treatment, roots of one cutting per tank (three in total) were pooled, homogenized and RNA extracted using a Qiagen RNeasy plant mini kit (see Supplementary Methods for full experimental and extraction conditions). RNA-seq libraries were sequenced at the Beijing Genomics Institute in Hong Kong on a BGISEQ500 with 100 bp paired-end reads (mean insert size 161 bp), producing 25.1-26.0 M read pairs per sample.

Transcriptome assembly and transcript quantification
After quality control and trimming of sequencing reads (see Supplementary Methods for details), de novo transcriptome assembly was performed using Trinity v2.10.0 68 using data from one individual per population per treatment. Completeness was assessed using the Eudicots dataset in BUSCO 69 v.4.0.5: 75% complete (72.2% single copy, 2.8% duplicated), 8.4% fragmented, 16.6% missing. After filtering (see Supplementary Methods for details), 27,970 genes were retained for downstream analysis. Transcripts were annotated using hmmer 70

Differential gene expression
Abundance estimates for transcripts were summarized at the gene level using tximport 72 v1.4.2. Gene expression analysis was performed using DESeq2 73 v1.26.0. Genes with low counts (<10) across all samples were removed. Variance-stabilizing transformed counts for 27,970 genes across all conditions were calculated and used in downstream analysis. This transformation reduces the dependence of the variance on mean expression values, making it more suitable for visualizing between-sample differences 73,74 . Principal components analysis of these counts for (1) all genes in control conditions (Fig. 1b), (2) all genes across all conditions (Fig. 1c) and (3) for DP genes (Fig. 1f) were calculated using the R prcomp function.
Genes differentially expressed between two populations within a treatment (control or zinc) were identified using DESeq2's in-built models with a single combined factor for population + condition (adjusted P = 0.05). Differentially expressed genes between T1 and S1, and T2 and S2 were identified in (1) control and (2) zinc treatments separately using contrasts (see Supplementary Methods section 5 for more details on models and contrasts used for all sets of differentially expressed genes). CEC genes were defined as those differentially expressed between both T1 and S1 in the control, and T2 and S2 in the control in the same direction (that is, both increasing or decreasing in T1 relative to S1 and T2 relative to S2). For between-treatment, within-population comparisons, a model with terms '~ Population + Population:Individual + Population:Condition' was fitted to account for individual-specific variation, which could be accounted for due to pools of clones from each individual being represented in both treatments. Genes differentially expressed between control and zinc treatment were identified for S1, S2, T1 and T2 using individual contrasts Article https://doi.org/10.1038/s41559-022-01975-w (see Supplementary Methods section 5). DP genes were defined as those differentially expressed between conditions in both T1 and T2 in the same direction (that is, both increasing or both decreasing from control to zinc treatment), and were differentially expressed between tolerant and susceptible populations in the control or zinc (or both). The significance of overlaps between sets of differentially expressed genes was determined using a one-sided Fisher's Exact test. GO enrichment analysis of gene sets was performed using GOseq 75 v1.38.0 with a false discovery rate of 0.05.
Quantification of fold changes of genes between populations and/or treatments used empirical Bayes shrinkage, calculated with the lfcShrink() function in DESeq2 76 . Values of |FC| were calculated for each gene as the absolute log 2 fold change between pairs of population/treatment groups (for example, T1 and T2 in the zinc) for a given set of genes. The sign of the log 2 fold change depends on the order of comparisons being made (for example, a value of +1 between T1 and T2 is equivalent of −1 between T2 and T1); the absolute value must be taken to meaningfully summarize the difference in expression levels (for example, the mean of −2 and +2 would be lower than that of 0.5 and 0.6). The median was used to summarize the values of |FC| as their distribution is highly skewed. Pairwise Wilcoxon signed-rank tests with Benjamini-Hochberg correction were used to detect significant differences in the distributions of |FC| between different pairs of population/ treatment groups.

Classifying responses to ancestral plasticity
To classify evolutionary responses to ancestral plasticity in the transcriptome-wide, DP and CEC gene sets, the following parameters were calculated for each gene: L o , mean expression value across S1 and S2 in control; L p , mean expression value across S1 and S2 in zinc; L a , mean expression value across T1 and T2 in zinc. These were used to calculate the initial plastic change (PC = L p −L o ) and subsequent evolutionary change (EC = L a −L p ) for each gene, as in ref. 12 . Only genes with substantial plastic and evolutionary change were assigned as undergoing reversion, reinforcement or plasticity (very small values of EC or PC due to measurement error would lead to spurious assignment of genes to categories 22 ). Genes were defined as having substantial (1) PC if they were differentially expressed between conditions in susceptible populations, combining data across S1 and S2 (using model ~Ecotype + Ecotype:Individual_plant + Ecotype:Condition, and contrast EcotypeS.CondZ, where ecotype (S1, S2) = S and (T1, T2) = T) and (2) EC if they were differentially expressed between tolerant and susceptible populations in zinc, combining data across both population pairs (using model ~Eco_Cond; a combined term of ecotype and condition, and the contrast SZ versus TZ). Data across ecotypes were combined to gain maximum power to detect small shifts in expression; alternative approaches outlined in Supplementary Methods gave similar results. Genes were assigned to one of three categories of evolutionary response to ancestral plasticity 36 : (1) Reinforcement, if EC×PC > 0; (2) Overshooting, if EC×PC < 0 and |EC | < 0.5×|PC|; or (3) Reversion, if EC×PC < 0 and |EC| > 0.5×|PC|. Significant differences in the relative proportions of these categories between sets of genes (for example, CEC genes compared to the transcriptome as a whole) were assessed using a two-tailed binomial test. Parametric bootstrapping of gene assignment to these categories following ref. 29 was implemented in R and repeated 100 times per gene (see Supplementary Methods); classification of genes passing this threshold is reported in Supplementary Table 1. For genes showing DP/CEC expression patterns but in T1/S1 or T2/S2 only, values of L o , L p , L a , EC and PC were only calculated using the samples from T1/S1 and T2/S2 separately (Supplementary Table 2) and categorized on the basis of these values. Assignment of categories for transcriptome-wide, CEC and DP genes was also calculated using T1/S1 and T2/S2 separately; these did not differ substantially between evolutionary replicates or the combined calculations (Supplementary Table 2).

Genotyping
For genotyping, cleaned reads were mapped to the transcriptome using HISAT2 77 v2.2.1. Genotypes were called using bcftools and a phylogenetic tree was constructed on the basis of 15,285 single-nucleotide polymorphisms (SNPs) using SNPhylo 78 v20180901 (see Supplementary Methods for details).

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

Data availability
RNA-seq data are deposited on the NCBI databases under Bioproject PRJNA706929. Source data are provided with this paper. Fig. 1 | Overlaps between sets of differentially expressed genes. Venn diagrams outlining the numbers of significantly differentially expressed genes for a given comparison (separated by a 'v'). T1/S1/T2/S2 correspond to populations, (C) corresponds to control conditions, (Z) to zinc conditions: for example, T1(C) v S1(C) refers to genes differentially expressed between T1 in the control, and S1 in the control. A '+' corresponds to genes that fulfil both sets of criteria: for example, T1 (C) v T1 (Z) + T2 (C) v T2 (Z) corresponds to the set of genes that are differentially expressed between both (i) T1 in the control vs. T1 in the zinc, and (ii) T2 in the control vs. T2 in the zinc. Panels illustrate (A) between treatment expression changes within populations, and their overlaps, (B) between ecotype-changes in the zinc in the two geographic groups, and their overlap, (C) the overlap between genes differentially expressed in both tolerant populations between conditions, both susceptible populations between conditions, and between both ecotype pairs within the zinc treatment, and (D) between-ecotype changes in the control in the two geographic groups, and their overlap.

Extended Data
Corresponding author(s): Alexander S.T. Papadopulos Last updated by author(s): Nov 23, 2022 Reporting Summary Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our Editorial Policies and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted 1 (mapping reads to transcriptome for genotyping), samtools v1.7 (removing reads with low mapping quality), SNPhylo v2014071 (phylogenetic tree construction), DESeq2 v1.2.6 (differential gene expression), tximport v1.14.2 (differential gene expression), R v 3.3.6 (statistical analysis), R v3. 63 For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

March 2021
Data Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A description of any restrictions on data availability -For clinical datasets or third party data, please ensure that the statement adheres to our policy RNA-seq data is deposited on the NCBI databases under Bioproject PRJNA706929.

Field-specific reporting
Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. All studies must disclose on these points even when the disclosure is negative.

Study description
This study aimed to identify whether zinc-tolerant populations of Silene uniflora have evolved parallel gene expression changes, and whether these gene expression changes were facilitated by ancestrally plastic responses to zinc. For each of 2 mine populations, and their nearest zinc-sensitive coastal populations (as identified in Papadopulos et al. 2021), 3 individuals were grown from wild collected seed and 6 cloned cuttings were propagated per individual; half of which were exposed to a neutral hydroponic solution, and half to a zinc-contaminated solution. Root tissue of cloned cuttings from a single individual were pooled in each treatment, and RNA extracted and sequenced. Gene expression profiles were calculated and compared between populations to identify the extent of shared gene expression changes in mine populations, and the extent to which these were facilitated by ancestral plastic responses to zinc (determined from the coastal zinc sensitive populations).

Research sample
Seeds were sampled previously at two mine sites and coastal sites where S. uniflora was known to occur. The nearest coastal populations to mine populations were identified using the BSBI database.

Sampling strategy
For the four sites, metal tolerance has previously been conducted and the evolutionary relationships of the populations established using reduced representation sequencing and phylogenetic/population genetic methods. At each site, we collected seeds from a minimum of twelve individuals, which were dried and stored separately with silica gel. Each individual grown from the wild seed was collected from a different individual at each site. Sample sizes were chosen to provide an accurate estimate gene expression changes between populations and treatments.

Data collection
Plants were grown by JAH and DPW. Hydroponics experiments and RNA extraction was performed by DPW, library preparation and sequencing were performed externally by the Beijing Genomics Institute.
Timing and spatial scale Seeds were collected between 2015-2017. Plants were grown from seed in October 2019; 10 weeks later cuttings were taken and grown for 5.5 weeks (timings allowed plants/cuttings to grow to sufficient size for experimentation). Rooted cuttings were transferred to 6 48x39x20cm 8L hydroponics tanks in a greenhouse in a neutral solution for 1 week then exposed either to refreshed neutral solution or zinc contaminated solution for 8 days before root tissue was harvested.

Data exclusions
No samples were excluded. Reproducibility 3 replicate hydroponic tanks per treatment each contained one cutting per individual; these were later pooled and analysed together. The whole experiment was not repeated.

Randomization
Cuttings of approximately equal size were distributed evenly between replicate hydroponics tanks, with each individual tank having the same proportions of cuttings from each population. Cuttings were randomly distributed within each tank.

Blinding
Blinding was not relevant to the study; data generation did not involve subjective human judgements.
Did the study involve field work?