Adaptive divergence generates distinct plastic responses in two closely related Senecio species

Organisms rely on plasticity to track environmental variation within their native range. However, it remains unclear how adaptation and plasticity interact, and how adaptive divergence affects the evolution of plasticity. To test for variation in plastic responses among two closely related but ecologically divergent ragwort species (Senecio, Asteraceae), we sampled c.40 genotypes of each species from natural populations. We then transplanted multiple clones of each genotype into four field sites along an elevational gradient representing each species’ native range, the edge of their range, and conditions outside their native range. At each transplant site, we quantified survival, growth, leaf investment, leaf morphology, chlorophyll fluorescence and gene expression. Both species performed better at their home sites, but the high elevation species showed lower tolerance to conditions outside its range than the low elevation species, suggesting stronger specialisation to the high elevation habitat. The two species also differed substantially in the direction of phenotypic and gene expression change across elevation, suggesting that distinct plastic responses have rapidly evolved in these two species. Adaptive divergence has led to the evolution of distinct plastic responses to environmental variation with distinct genomic architectures, despite these two species having shared a recent common ancestry.


38
The resilience of natural populations and communities to novel or changing environments relies on the ability 39 of genotypes to adjust their phenotype to track local conditions (Chevin et al. 2010). Phenotypic plasticity 40 generates different phenotypes from the same genotype depending on the environment to which it is exposed 41 (Via et al. 1995;Ghalambor et al. 2007; Charmantier et al. 2008). The ability for plasticity to track 42 environmental variation is shaped by selection within environments routinely experienced (Ghalambor et al. 43 2007). Plasticity therefore evolves to buffer populations in response only to particular environmental regimes 44 (Bradshaw 1965;Schlichting 1986; Baythavong and Stanton 2010), which can increase ecological 45 specialisation when populations adapt to highly predictable environments (Poisot et al. 2011). Characterising 46 plasticity in closely related, but ecologically divergent species can test to what extent adaptive divergence 47 and ecological speciation is also associated with divergence in phenotypic plasticity. 48 The effect of adaptation on the nature of plastic responses will depend on how plasticity and selection and amount of variation in plastic responses. We would predict that contrasting environments select for 52 different forms and different magnitudes of plasticity, but it is not known to what extent a shared common 53 ancestry will constrain such divergence in plasticity. We would also predict that plasticity can maintain 54 fitness within a narrower range of environmental variation in species that are more specialised to a particular 55 environment, reflecting a narrower range of ecological tolerances in such species (Lortie and Aarssen 1996; 56 Debat and David 2001). However, few studies have assessed whether adaptation to contrasting environments 57 causes not only adaptive divergence, but also divergence in plastic responses when exposed to the same 58 environmental variation. 59 expression profiles for closely related but ecologically divergent species can therefore reveal the effect of 78 adaptive divergence on the genomic basis of plasticity. 79 In this study, we first identify physiological differences in a common garden experiment between two closely 80 related species of Senecio that inhabit contrasting elevations on Mt. Etna, Sicily. We then quantify variation 81 in environmental sensitivity and phenotypic plasticity across four field transplant sites along the elevational 82 range of both species. Senecio chrysanthemifolius (Fig. 1a) is a short-lived perennial with highly dissected 83 leaves that occupies disturbed habitats (e.g., vineyards, abandoned land and roadsides) in the foothills of Mt. 84 Etna c. 400-1,000m a.s.l (above sea level) and at similar elevations throughout Sicily. By contrast, S. 85 aethnensis (Fig. 1b) is a perennial with entire, glaucous leaves and is endemic to lava flows c. 2,000-2,600m  Given that S. aethnensis is endemic to high elevations on Mt. Etna, while S. chrysanthemifolius is found in a 92 variety of lower-elevation habitats across Sicily, we predicted that: (1) Their occupation of highly contrasting 93 habitats would be reflected by differences in physiology in a common garden, suggesting functional trait 94 divergence.
(2) As evidence of adaptive divergence, species would perform better at their native elevation, 95 and better than the foreign species. In addition, S. aethnensis would be less tolerant of novel environments, 96 providing evidence that it is more specialised to its native habitat than S. chrysanthemifolius. (3) Despite 97 clear patterns of adaptive divergence, their recent common ancestry would mean that both species will show 98 similar patterns of plasticity, and similar patterns of G×E underlying plasticity. By contrast, if ecological 99 speciation is associated with rapid changes in plasticity, then functional traits associated with adaptive 100 divergence would also show large differences in plasticity between the two species, and different levels of 101 transplanted multiple cuttings of each genotype to four transplant sites across an elevational range that 104 included the home range of each species, and two intermediate elevations. We quantified survival, growth, 105 leaf gene expression, photosynthetic activity and leaf investment. Given that rapid divergence in leaf shape 106 and dissection is key feature of speciation in this system, as well as in other Senecio species (Walter et al. 107 2018), we also quantified leaf morphology traits, which are likely to be associated with fitness at different 108 elevations.

116
Sampling natural populations 117 We sampled seeds, and took cuttings from naturally growing individuals of both species after plants started 118 flowering. This was conducted in May-June 2017 for S. chrysanthemifolius and July 2017 for S. aethnensis 119 because S. aethnensis grows more slowly and flowers later than S. chrysanthemifolius, given its high 120 elevation habitat. For S. chrysanthemifolius, we sampled from 88 individuals from five sites, each a 121 geographically separated patch of individuals representing potentially discrete populations, all below 122 800m.a.s.l (Fig 1c, Table S1). For S. aethnensis, we sampled from 87 individuals at four different elevations 123 (2,600m, 2,500m, 2,400m and 2,300m.a.s.l) on both the North and South slopes of Mt. Etna (Fig 1c, Table   124 S1). To minimise the risk of sampling close relatives, most plants sampled were more than 10m apart.

125
Physiological differences between species under common garden conditions 126 To assess differences in physiology between the species, we grew plants from seeds in a growth cabinet (see In the glasshouse (Giarre, Italy), cuttings from all individuals sampled from natural populations (hereafter, 141 genotypes) were cut into 5cm stem segments, each possessing 2-3 leaf nodes. Each smaller cutting was then 142 dipped in rooting plant growth regulator for softwood cuttings (Germon® Bew., Der. NAA 0.5%, L. Gobbi, 143 Italy) and placed in a compressed mix of coconut coir and perlite (1:1) in one cell of an 84-cell tray. All 144 cuttings (i.e. clones) from each genotype were kept together in one half of a tray, with tray positions 145 randomised regularly to prevent systematic environmental or positional effects. Trays were kept moist and 146 checked regularly for cuttings that successfully produced roots (roots extending out of the bottom of tray).

147
For each genotype, rooted cuttings were randomised into experimental blocks and transplanted at four field 148 sites. From the initial genotypes sampled, we transplanted 37 S. chrysanthemifolius genotypes and 42 S. 149 aethnensis genotypes that produced enough cuttings with roots.

150
Field transplant sites were at four elevations (500m, 1,000m, 1,500m and 2,000m a.s.l) along a transect on 151 the south-eastern side of Mt. Etna (Fig. 1c). The 500m site was located in a garden among fruit trees and 152 grape vines, the 1,000m site on an abandoned vineyard among Quercus ilex, the 1,500m site among an apple 153 and pear orchard, and the 2,000m site surrounded by pine trees on a lava flow from 1983. Both the native 154 elevations (500m for S. chrysanthemifolius and 2,000m for S. aethnensis) were located less than 1km from 155 natural populations. Furthermore, plants of both species were often observed at intermediate elevations 156 (including close to the 1,000m and 1,500m transplant sites), but were never observed within the native range 157 of the other species. Soil is characterised as a silty sand at elevations between 500m and 1,500m, but changes 158 to volcanic sand at 2,000m. At each transplant site we deployed four data loggers (Tinytag Plus, Gemini Data 159 Loggers, UK), which measured temperature hourly. We also took three soil samples at each transplant site, 160 which were analysed for 36 variables that included nutrients, salts and ions (Nucleo Chimico Mediterraneo 161 laboratories, Catania, Italy). To analyse the soil data, we used Multi-Dimensional Scaling (MDS) to calculate 162 the scaled distance between the soil samples taken at all transplant sites. 163 We transplanted multiple cuttings of each genotype into three experimental blocks at each transplant site.   and II (see Methods S1).

203
Statistical analyses of survival, growth and plasticity 204 We first tested for significant differences in survival across elevation, and quantified plasticity as the change , (1)  aethnensis). All leaves for a cutting were placed in the same Eppendorf tube and stored in RNAlater at -20ºC.

244
Three genotypes of each species were then selected at random, and three clones sampled from each transplant 245 site (36 samples in total per species). We extracted RNA from each sample using QIAgen RNeasy kits.

246
Library preparation and RNA sequencing was performed at the Oxford Genomics Centre on an Illumina 247 Hiseq4000 platform, producing 75bp paired-end reads.

248
Quantifying differential expression across transplant sites, genotypes, and species 249 A reference transcriptome was assembled for each species (Methods S1).  The transplant sites experience contrasting climatic conditions associated with elevation, with extreme heat 301 (regularly exceeding 40ºC) at 500m and 1,000m during summer, and extreme cold (regularly below 0ºC) at 302 1,500m and 2,000m during winter (Fig. 3a). Soil profiles separated the four transplant sites in a linear 303 fashion along the first axis (MDS1), which represented a transition in soil type and reduction in nutrients 304 (amount of organic material, total nitrogen, cation exchange capacity and exchangeable ions) at higher 305 elevations (Fig. 3b). The second axis (MDS2) characterised differences between the 1,000m site and the 306 other sites, associated with greater concentrations of various salts at 1,000m (soluble nitrates, calcium and 307 magnesium). remained at 500m and 1,000m, respectively, as compared to 79% and 39% for S. chrysanthemifolius (Fig.   325   S1). This consistency in patterns of mortality suggests that the 2017 experiment represented typical patterns 326 of mortality. confidence intervals for the estimate of the mean and letters denote significant differences (full statistical summaries are located in 332 Table S2).  95% confidence intervals for the estimate of the mean of each species at each elevation. Letters denote significant differences 360 between transplant sites calculated using pairwise tests conducted within each species and adjusted for multiple comparisons (full 361 statistical summaries are located in Table S2). Grey points represent the mean of all cuttings for each genotype, within species.

362
Senecio aethnensis showed significant decreases in the number of indents with elevation, while S. chrysanthemifolius showed a 363 decrease in leaf complexity at higher elevations. Leaf area changed similarly across elevation for both species, while specific leaf 364 area was much greater for S. aethnensis at lower elevations. Senecio aethnensis also showed lower total photosynthetic 365 performance (PI total ) at lower elevations, while S. chrysanthemifolius did not change. 366

367
To identify how multivariate phenotype changed across elevation, we quantified elevational changes in the 368 first two principal components (Fig. 6a). Both species showed significant changes in principal component 369 scores across the four transplant sites, but patterns differed between the two species. The first principal 370 component described differences between species as well as phenotypic differences associated with elevation 371 for S. aethnensis ( Fig. 6b; PC1 species×elevation  (Fig. 6a). By contrast, 375 S. chrysanthemifolius outside its range shifted its phenotype towards that expressed by S. aethnensis at high 376 elevations.  G×E can be created either by a change in the magnitude of differences among genotypes across sites (i.e., a 390 change in among-genotype variance), or by genotype-specific differences in their response to the 391 environment, where genotypes change rank in different environments (i.e., where genotype reaction norms 392 cross each other). We found that genotypes of both species varied significantly in their response to elevation 393 (Fig. 7a). Generally, this meant that the magnitude of the variance among genotypes (within species) 394 changed significantly between elevations, whereby >90% of G×E effects are characterised by changes in 395 among-genotype variance between transplant sites (Fig. 7a), and visualised as a greater magnitude of 396 differences among genotypes at 500m compared to 2,000m (Fig. 7b-c and Fig S2). However, genotypes did 397 not change in rank between elevations. Furthermore, S. chrysanthemifolius showed no change in mean for 398 PC1 across elevation, but a large change in PC2 (Fig. 6), and this was reflected by the opposite pattern in 399 G×E: strong G×E in PC1 and no G×E in PC2 (Fig. 7). By contrast, S. aethnensis showed changes in mean 400 (Fig. 6) and significant G×E for both principal components (Fig. 7). The two methods of estimating differential expression (DESeq2 vs limma/voom) showed similar patterns of 414 gene expression variation, although limma/voom detected fewer strongly differentially expressed genes in 415 each species (Fig. S3). Patterns of gene expression between species revealed a high proportion of 416 differentially expressed genes at each transplant site, with the most at 2,000m (1,677 genes) and decreasing at 417 1,500m (1,451 genes), 1,000m (1,079 genes) and 500m (1,051 genes) (Fig. 8a). In total 383 genes were  differentially expressed across all elevations. Functional enrichment of differentially expressed genes at each 419 transplant site reflected differences in photoprotection, pigmentation and water use efficiency (Fig. 8a). 420 Within each species, more genes were differentially expressed as genotypes were moved further from their 421 native elevations (Fig. S3). Genes that were differentially expressed between the home site and the most 422 novel environment (i.e. the elevational extremes) for each species indicated little overlap between the two 423 species, with just 5.5% and 6.5% of overexpressed and underexpressed genes shared between species (Fig.   424   8b). For the ten genes of each species with the largest change in overexpression and underexpression 425 between 2,000m and 500m, we observed contrasting patterns between the two species: strong overexpression 426 or underexpression in one species but a relatively unchanged expression profile in the other species (Fig. 9). 427 Functional enrichment analyses of differentially expressed genes between the elevational extremes revealed   (Table S4). 434 Network reconstruction resulted in 32 network modules ranging in size from 196 to 3,108, with 170 genes 435 unassigned (Table S5). 17 modules showed a significant correlation with elevation (Fig. 8c). 7 modules 436 correlated with elevation uniquely in S. aethnensis and were annotated with terms including phyllome 437 development, translational initiation, protein targeting to chloroplasts, response to heat and immune response.   Given the two Senecio species on Mt. Etna are the result of recent ecological speciation, we predicted that (1) 457 the two species would exhibit differences in physiology when grown in common garden conditions, and (2) 458 that species would perform well in their native habitat, but poorly in the habitat of the other species. In 459 addition, we predicted that, given S. aethnensis is endemic to high elevations on Mt. Etna, it would maintain 460 high fitness across a narrower range of elevations than S. chrysanthemifolius.
(3) We also predicted that due 461 to their recent common ancestry, both species would display similar patterns of plasticity, and similar 462 patterns of G×E underlying plasticity.

463
Consistent with our first prediction, genotypes of both species showed distinct behaviour in the common 464 garden experiment (Fig. 2), indicating that adaptive divergence in high and low elevation habitats has 465 generated substantial shifts in physiology. In support of our second prediction, S. aethnensis showed reduced 466 survival at all elevations outside its native elevation, while S. chrysanthemifolius only showed reduced 467 survival at the elevation furthest from its native elevation (Fig. 4a). Furthermore, S. aethnensis suffered 468 greater mortality at the start of the transplant experiment (i.e. during summer), which occurred early in their 469 flowering season and likely reduced flower and seed production (Fig. 4b). By contrast, high mortality for S. 470 chrysanthemifolius was only recorded at the highest elevation after an extensive flowering season (i.e. during 471 winter; Fig. 4c), suggesting that fecundity would likely be higher for S. chrysanthemifolius across the entire 472 elevation gradient during the experiment. Reduced photosynthetic activity at lower elevations was also 473 associated with lower survival for S. aethnensis, which contrasted with S. chrysanthemifolius that maintained 474 a lower, but constant photosynthetic activity across elevation (Fig. 5e). These results support our prediction 475 that S. chrysanthemifolius can tolerate a wider range of environmental variation compared to S. aethnensis, 476 reflecting its broader distribution across a wider range of (lowland) habitats across Sicily.

477
However, contrary to our third prediction, we found different patterns of plasticity in phenotype (Figs. 5-6) 478 and gene expression (Figs. 8-9) between the two species. Senecio chrysanthemifolius showed substantial 479 reductions in leaf complexity at higher elevations but showed little change in the number of leaf indents. By the two species were also associated with differences in trait mean, suggesting an important link between 483 adaptive divergence and plasticity. Both species produced larger, thinner leaves at lower elevations, but S. 484 aethnensis produced significantly thinner leaves than S. chrysanthemifolius at the lowest elevation. Taken 485 together, phenotypic plasticity in S. chrysanthemifolius for these traits increased its resemblance to the 486 multivariate phenotype of S. aethnensis at high elevation. By contrast, phenotypic plasticity in S. aethnensis 487 reduced its resemblance to S. chrysanthemifolius at lower elevations. We also detected significant G×E 488 interactions underlying plasticity in leaf traits (Fig. 7), suggesting standing variation that could allow 489 plasticity to evolve in response to novel environmental conditions. However, contrary to our prediction that 490 both species would show similar levels of G×E interactions, we found evidence that patterns of G×E 491 interactions were stronger for genotypes of S. aethnensis (Fig. 7). 492 Consistent with the distinct forms of phenotypic plasticity observed, but also contrary to our third prediction, 493 gene expression changes across elevation involved surprisingly distinct gene networks for the two species.

494
The number and nature of differentially expressed genes changed with elevation (Fig 8a), with only 1.6% of these data reflect clear differences in plasticity between species for genotypes growing in natural populations.

506
Putative function of differentially expressed genes 507 We identified c.600 loci per species that showed differential expression between the extreme elevations.

514
At each transplant site, the genes differentially expressed between species were associated with functions that 515 reflected the differences observed in physiological measures between species grown under common garden 516 conditions. For example, genes relating to water use were differentially expressed between the species at 517 lower elevations, which were reflected by greater water use efficiency in the leaves of S. 518 chrysanthemifolius in the common garden experiment. Likewise, genes differentially expressed at higher 519 elevations were associated with light responses, photosynthesis and cuticle composition, which could reflect 520 the differences in pigment composition observed between the leaves of the species grown in the common 521 garden. This functional significance of the genes that vary most in expression, and differ between these 522 species, provides further support that adaptive divergence has shaped the evolution of distinct forms of 523 plasticity in these species.

524
Adaptive divergence between the two Sicilian Senecio species 525 Our data indicate that the plastic responses of either species are not sufficient for them to maintain fitness at 526 the native elevations of the other species, and that this is especially the case for S. aethnensis (Fig. 4a). 527 Adaptation of these species to the contrasting habitats on and around Mt. Etna is likely to result in selection 528 against genotypes of either species outside their native habitat, which is likely to reduce any gene flow 529 between the species (Ross 2010). Results from the common garden experiment suggest that such ecological 530 divergence has resulted in differentiation in functional traits related to light defence and leaf pigments.

531
Furthermore, the field transplants reveal that divergence of functional traits at the gene network and 532 phenotypic level was associated with adaptation to the contrasting elevations. Therefore, reductions in 533 survival outside each species native ranges seems likely to be due to divergence in traits that contribute to 534 adaptation in their native habitat (Ross 2010).

536
Although the parental forms can form hybrids, there is evidence of intrinsic reproductive isolation that results 537 in low fitness of the F2 and F3 generations (Hegarty et al. 2009). There is also evidence that parental 538 individuals rarely meet on Mt. Etna, instead the hybrid zone is almost entirely populated by hybrid 539 phenotypes that have only occasional contact with parental individuals due to the separation between the 540 hybrid zone and parental populations ). However, it remains possible that some mutually 541 beneficial alleles are able to spread into the parental genomes of either species, making it even more 542 remarkable that we detected divergence in plasticity between the two species.

543
Adaptive divergence creates differences in plasticity 544 Our results suggest that adaptive divergence has a rapid effect on patterns of phenotypic plasticity and on the 545 genetic basis of plasticity, even in two closely related species. It remains possible that genetic drift during the 546 formation of these two species could cause the species-specific differences in plasticity that we observed.

547
However, given substantial divergence in leaf form and physiology, the differences in plasticity of these 548 same traits among species, and the higher fitness of each species at their native versus novel habitats, it 549 seems far more likely that adaptive divergence between S. aethnensis and S. chrysanthemifolius is 550 distinct differences in plasticity.

553
Understanding how plasticity evolves is important for understanding how species can respond to novel 554 environmental variation (Bradshaw 1965 decline in fitness than S. chrysanthemifolius across the elevational gradient (Fig. 4a), and this was associated 560 with stronger reductions in leaf investment (Fig. 5d), lower photosynthetic activity (Fig. 5e) and reduced 561 capacity to approach the native phenotype of S. chrysanthemifolius at lower elevations. Therefore, the 562 distinct patterns of plasticity observed for the two species may be caused by S. aethnensis exhibiting 563 maladaptive plasticity as a consequence of greater specialisation to the high elevation habitat. 564 We suggest two possible explanations for the distinct patterns of plasticity in these species. First, species may habitat.

576
The second explanation is that species will show reduced plasticity in the genes and traits that need to be 577 continually expressed in their native habitat. This would occur if plastic responses required to initially 578 colonise the high elevation habitat become genetically-based (via genetic assimilation) because consistently 579 strong stabilising selection for a single phenotype is favoured (Waddington 1953). This would suggest that as 580 divergent selection drives trait divergence between the two species, it changes the level of plasticity in the 581 traits that are diverging. In this scenario, plasticity would only be retained in traits that maintain fitness in 582 response to the environmental variation experienced within their native habitat. For example, S. aethnensis 583 shows a greater number of indents, as well as greater plasticity in this trait (Fig. 5a), while S. 584 chrysanthemifolius shows both greater leaf complexity than S. aethnensis, as well as greater plasticity in this 585 trait (Fig. 5b). At high elevations Senecio aethnensis is likely to experience consistently strong selection for 586 less complex leaves, leading to a loss of plasticity in this trait, perhaps via genetic assimilation. By contrast, 587 S. chrysanthemifolius likely experiences selection for more complex leaves, but the variety of habitats that 588 this species occupies could maintain plasticity in this trait.

589
Traits, such as leaf complexity or leaf indentation, that show strong divergence between species as well as  Genotype-by-environment interactions and the evolutionary potential of plasticity 597 The lower tolerance of S. aethnensis to conditions beyond its native range suggests that adaptation to the high 598 elevation environment has reduced its ability to respond to novel environmental variation at lower elevations.

599
Such limits to plasticity threaten the persistence of high elevation species in response to climate change 600 unless genotypic variation in plasticity (G×E) can promote rapid adaptation. In both species we found 601 significant genotypic variation in plastic responses to the elevational gradient, suggesting substantial 602 evolutionary potential. Senecio aethnensis showed substantial G×E in both the multivariate axis that 603 separated the two species, and in the axis associated with consistent phenotypic change across elevation for 604 both species. By contrast, S. chrysanthemifolius displayed substantial G×E in the multivariate axis that 605 separated the two species, but not in the axis associated with consistent phenotypic change across elevation.

606
This suggests little genetic variation for plasticity (and therefore adaptive potential) in response to elevational  (Fig. 7), rather than by genotype-specific (i.e., change in genotype rank between