# Phenotypic plasticity evolves at multiple biological levels in response to environmental predictability in a long-term experiment with a halotolerant microalga

Christelle Leung, Daphné Grulois, Leandro Quadrana, Luis-Miguel Chevin

► To cite this version:

Christelle Leung, Daphné Grulois, Leandro Quadrana, Luis-Miguel Chevin. Phenotypic plasticity evolves at multiple biological levels in response to environmental predictability in a long-term experiment with a halotolerant microalga. PLoS Biology, 2023, 21 (3), pp.e3001895. 10.1371/journal.pbio.3001895. hal-04285689

HAL Id: hal-04285689

https://hal.science/hal-04285689v1

Submitted on 24 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

# Phenotypic plasticity evolves at multiple biological levels in response to environmental predictability in a long-term experiment with a halotolerant microalga

Christelle Leung<sup>1*</sup>, Daphné Grulois<sup>1</sup>, Leandro Quadrana<sup>2</sup>, Luis-Miguel Chevin<sup>1*</sup>

1 CEFE, Université de Montpellier, CNRS, EPHE, IRD, Montpellier, France, 2 Institut of Plant Science Paris-Saclay (IPS2), UMR Université Paris Saclay—CNRS 9213—INRAE 1403, Université Evry Val d'Essonne—Université Paris Diderot, Paris, France

* christelle.leung@umontreal.ca (CL); luis-miguel.chevin@cefe.cnrs.fr (L-MC)

<figure><img src="image_2.png" /></figure>

## Abstract

Phenotypic plasticity, the change in the phenotype of a given genotype in response to its environment of development, is a ubiquitous feature of life, enabling organisms to cope with variation in their environment. Theoretical studies predict that, under stationary environmental variation, the level of plasticity should evolve to match the predictability of selection at the timing of development. However, the extent to which patterns of evolution of plasticity for more integrated traits are mirrored by their underlying molecular mechanisms remains unclear, especially in response to well-characterized selective pressures exerted by environmental predictability. Here, we used experimental evolution with the microalgae *Dunaliella salina* under controlled environmental fluctuations, to test whether the evolution of phenotypic plasticity in responses to environmental predictability (as measured by the squared autocorrelation $\rho^2$) occurred across biological levels, going from DNA methylation to gene expression to cell morphology. Transcriptomic analysis indicates clear effects of salinity and $\rho^2$ × salinity interaction on gene expression, thus identifying sets of genes involved in plasticity and its evolution. These transcriptomic effects were independent of DNA methylation changes in *cis*. However, we did find $\rho^2$-specific responses of DNA methylation to salinity change, albeit weaker than for gene expression. Overall, we found consistent evolution of reduced plasticity in less predictable environments for DNA methylation, gene expression, and cell morphology. Our results provide the first clear empirical signature of plasticity evolution at multiple levels in response to environmental predictability, and highlight the importance of experimental evolution to address predictions from evolutionary theory, as well as investigate the molecular basis of plasticity evolution.

## 1. Introduction

Phenotypic plasticity, the ability of a given genotype to produce different phenotypes depending on environmental conditions, is an important mechanism enabling organisms to cope

<figure><img src="image_3.png" /></figure>

### OPEN ACCESS

**Citation:** Leung C, Grulois D, Quadrana L, Chevin L-M (2023) Phenotypic plasticity evolves at multiple biological levels in response to environmental predictability in a long-term experiment with a halotolerant microalga. PLoS Biol 21(3): e3001895. https://doi.org/10.1371/journal.pbio.3001895

**Academic Editor:** Csaba Pál, Biological Research Centre, HUNGARY

**Received:** October 19, 2022

**Accepted:** February 20, 2023

**Published:** March 24, 2023

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pbio.3001895

**Copyright:** © 2023 Leung et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** Raw sequence data (RNA-seq and WGB-seq) used in this study are deposited in the NCBI's Sequence Read Archive
(SRA) database under BioProject PRJNA736997 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA736997). The specific samples used in this study are under the BioSample accessions listed in S1 and S2 Tables. DNA methylation 100 bp tiles percentage values, raw RNA-seq read counts (i.e. transcript expression values), structural and functional gene annotations, and cell morphology data have been deposited in Figshare https://doi.org/10.6084/m9 figshare.21905670.

**Funding:** This work was supported by the European Research Council (Grant 678140-FluctEvol) to LMC, a Fonds de Recherche du Québec - Nature et Technologies (FRQNT) postdoctoral fellowship to CL, and a travel grant for collaboration provided by the GDR Plasticité Phénotypique (GDR 3715) from CNRS (CL). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

**Abbreviations:** BH, Benjamini-Hochberg; DE, differentially expressed; DMR, differentially methylated region; FSC, Forward Scatter; GO, Gene Ontology; LRT, likelihood ratio test; RDA, redundancy analyses; RNA-seq, RNA-sequencing; SSC, Side Scatter; TSS, transcription start site; WGBS, whole-genome bisulfite sequencing.

with variation in their environment. Understanding what drives the evolution of plasticity, from selective causes to underlying mechanisms, is thus important not only for basic research, but also for on our ability to predict the fate of populations, especially under global change [1,2]. Despite this interest, we still lack critical empirical information on the extent to which selection on phenotypic plasticity propagates across hierarchical levels of the organisms, from more integrated phenotypes that are directly exposed to selection, to their underlying molecular basis.

Theoretical and empirical studies have demonstrated that the adaptiveness of phenotypic plasticity arises from its interplay with environmental variation in selection, when the latter is partly predictable [3-7]. A high degree of plasticity is expected to be favored when the environmental cue causing induction of a given phenotype is a reliable predictor of selection acting on this phenotype. Conversely, under unpredictable environmental changes (or unreliable cues), plasticity leads to phenotypes that are often mismatched with their optimum, such that lower degrees of plasticity are selected [3-7]. In the limit of completely unpredictable environments, bet-hedging strategies that do not rely on environmental cues may even be selected [8,9].

In terms of mechanisms and pathways, the increased availability of *-omics* methods has greatly facilitated our understanding of the relationship between molecular changes and changes in phenotypic expression. In particular, many studies have shown that phenotypic plasticity often involves at its most basic level variation in gene expression with the environment [10-14]. This suggests that the molecular basis of phenotypic plasticity should also encompass mechanisms that regulate gene expression, as clearly demonstrated in a number of cases [12-14]. In particular epigenetic processes, i.e., modifications of chromatin beyond DNA sequence, transmitted through mitosis (and possibly meiosis), and potentially influencing gene expression, were shown to be important molecular mechanisms for phenotypic plasticity [15-18]. As both gene expression and epigenetic marks (such as DNA methylations) are known to be under genetic control [19-21], selection on phenotypic plasticity for more integrated traits is thus expected to cascade down to cause evolution of plasticity for its underlying molecular mechanisms (e.g., [22]), but this process has seldom been investigated experimentally. More specifically, the extent to which classical theoretical predictions about evolution of plasticity in response to environmental predictability hold from molecular phenotypes to more integrated ones remains unexplored.

Recently, we have shown that environmental variation in CpG methylation and gene expression in the halotolerant microalga *Dunaliella salina* contributes to phenotypic plasticity in this species [23]. In addition, we confirmed the contribution of the genotype to salinity responses at the levels of DNA methylation, gene expression, and population growth rate, highlighting the evolutionary potential of phenotypic plasticity at multiple levels [23]. Importantly, long-term experimental evolution in this species has demonstrated that reduced plasticity in cell shape and content evolved in populations confronted to less predictable environments [5]. Here, we used experimental evolution under controlled environmental fluctuations to assess whether the evolution of plasticity in responses to environmental predictability permeates across biological levels, from DNA methylation to gene expression to more integrated phenotypes.

## 2. Results

We analyzed 9 populations of the halotolerant microalga *D. salina* (strain CCAP 19/15) that have evolved under regimes of randomly fluctuating environments, with controlled and variable predictability [5,24]. During experimental evolution, lines derived from a single ancestral population were exposed to randomly fluctuating salinity, with changes every 3 or 4 days (with
<figure><img src="image_4.png" /><figcaption>**Fig 1. Experimental design.** The experiment included 3 steps: (i) long-term experimental evolution (left); (ii) acclimation at a constant, intermediate salinity (middle); and (iii) plasticity assays at high versus low salinity (right). Each colored time series on the left represents an actual realization of salinity fluctuations for one of the populations used in this study, with the color denoting the treatment of stationary (i.e., expected long-term) temporal autocorrelation ($\rho$). At the end of the experiment, cells were harvested for DNA methylation, gene expression, and cell morphology analyses.</figcaption></figure>

https://doi.org/10.1371/journal.pbio.3001895.g001

about one generation per day), for a total of ca. 500 generations. Salinity had a normal distribution over time, with the same mean (2.4 M NaCl) and standard deviation (1 M NaCl) across treatments, but variable autocorrelation, and, hence, variable predictability [5,24]. Previous morphological analysis of 32 of these lines revealed that reduced morphological plasticity has evolved in lines that experienced less predictable environments [5]. To further characterize the molecular basis of this evolution, we set out to determine whether the plasticity of DNA methylation and gene expression levels evolved for a subset of these lines that experienced 3 different predictability treatments (3 lines per treatment): low ($\rho^2$ = 0), intermediate ($\rho^2$ = 0.25), and high ($\rho^2$ = 0.81) predictability, where $\rho$ is the stationary (long-term) temporal autocorrelation of salinity time series. All populations were subsequently subjected to a 10-day acclimation step at intermediate salinity ([NaCl] = 2.4 M), to ensure they had similar physiological states and population densities before the phenotypic and molecular plasticity assays. They were then placed for 24 h at 2 salinities near the extremes of their historical range ([NaCl] = 0.8 M and 4.0 M), to assess their degree of plasticity in DNA methylation, gene expression, and individual cell morphology (Fig 1). At any of these levels, an effect of salinity is indicative of phenotypic plasticity, an effect of the evolutionary treatment $\rho^2$ denotes evolution, while $\rho^2$ × salinity interaction indicates evolution of plasticity.

## 2.1 DNA methylation and gene expression plasticity evolved in response to environmental predictability

DNA methylation can contribute to phenotypic differentiation by influencing gene expression regulation [18, 25–28]. To assess whether experimental evolution of *D. salina* may lead to epigenetic differentiation, we performed whole-genome bisulfite sequencing (WGBS) for all samples, yielding a total of 1.23 × 10⁹ 150 bp paired-end raw reads, and an estimated average depth
<figure><img src="image_5.png" /><figcaption>**Fig 2. Evolution and plasticity of DNA methylation.** (A) Variation of DNA methylation patterns across evolved populations. RDA plot performed on DNA methylation patterns according to evolution conditions (colors) and salinity during plasticity assay (shapes). (B) Distribution of DMRs. Venn diagram describing the number of DMRs (q-value < 0.05 and |diff-Methylation| > 20%) across evolutionary conditions (Predictability, light gray), and between salinities within each of the 3 evolutionary conditions: low (green), intermediate (blue), and high (red) predictability of environmental changes. The raw data underlying this figure are available in the Figshare repository https://doi.org/10.6084/m9 figshare.21905670. DMR, differentially methylated region; RDA, redundancy analysis.</figcaption></figure>

https://doi.org/10.1371/journal.pbio.3001895.g002

of coverage of 43.76 × (s.d. 3.14 ×) per sample (S1 Table). After data filtering, we carried out our methylation analyses on an average of 7.56 × 10⁷ (s.d. 8.28 × 10⁶) cytosines per samples at the CpG context (S1 Table), where methylations are most prevalent in this species [23]. Redundancy analyses (RDA) based on overall CpG methylation revealed a significant effect of evolutionary treatments ρ² (R²adj = 4.66%; P = 0.005) on DNA methylation (Fig 2A). Fig 2A, which is a constrained ordination that maximizes variation in DNA methylation that can be explained by salinity and environmental predictability, suggests that lines from highly predictable environments displayed higher salinity differences than those from lowly predictable environments. However, we did not detect a significant marginal effect of salinity (R²adj = 0.79%; P = 0.356) or ρ² × salinity interaction (R²adj = 0.30%; P = 0.562) on the overall DNA methylation pattern, suggesting the absence of overall DNA methylation changes in response to salinity or a lack of power to detect small variation being explained.
That we did not detect significant differentiation at the whole epigenome level does not pre-
clude the existence of more localized epigenetic differences along the genome, so we also inves-
tigated regional changes in DNA methylation at a finer scale, by considering nonoverlapping
100 bp windows, hereafter denoted differentially methylated regions (DMRs). We detected
227 DMRs (with FDR < 0.05 under Benjamini-Hochberg (BH) adjustment of *P* values and |
diff-Methylation| > 20%) among the evolutionary treatments, as summarized by their envi-
ronmental predictability *ρ*² (Fig 2B). We also detected ρ²-specific DMRs between salinities
within each evolutionary treatment, among the 14,357 total 100 bp regions (Fig 2B). Interest-
ingly, populations that evolved under less predictable environmental fluctuations displayed the
least number of DMRs between salinities (*n* = 3), as compared to populations from intermedi-
ate (*n* = 78) or high (*n* = 29) environmental predictability (Fig 2B), indicative of reduced epige-
netic plasticity. We then assessed whether changes in DNA methylation patterns in response
to a given environmental challenge involved similar genomic regions in the different evolved
lines. Comparison of the list of DMRs between salinities revealed no overlap across evolution-
ary conditions (Fig 2B), suggesting that evolution of plastic epigenetic responses involved
modifications of methylations in distinct genomic regions in different treatments.

We next investigated variation in gene expression, by analyzing 32,718 transcripts through
RNA-sequencing (RNA-seq). We obtained 5.62 × 10⁸ 150 bp paired-end raw reads in total (S2
Table). As with DNA methylation, we detected a significant effect of the evolutionary treat-
ments ρ² (R²<sub>adj</sub> = 9.79%; P < 0.001) on gene expression at the whole-transcriptome level. How-
ever, the assay salinity now explained the greatest part of variation in gene expression (R²<sub>adj</sub> =
32.91%; P < 0.001), indicating pervasive and highly significant plasticity. The ρ² × salinity
interaction was also significant (R²<sub>adj</sub> = 4.30%; P = 0.029) (Fig 3A), indicating evolution of
transcriptional plasticity. At a more local level, analyses of differentially expressed (DE) tran-
scripts yielded similar results: The numbers of transcripts that were significantly DE (Likeli-
hood ratio test, FDR < 0.05 and |Log₂FC| > 1) were highest for contrasts between salinities (*n*
= 4,283), followed by evolutionary treatments ρ² (*n* = 1,315), and, finally, ρ² × salinity interac-
tion (*n* = 837). As for DMRs, populations that evolved in less predictable environments dis-
played fewer DE transcripts than populations that evolved in highly predictable environments
(Fig 3B, *n* = 2,638, 2,844, and 4,086 for ρ² = low, intermediate, and high, respectively). How-
ever, in contrast to what we found for DNA methylation, we observed substantial overlap
among DE transcripts identified between salinities for the different evolved lines (Fig 3B),
indicating that the plastic response to salinity largely involved transcriptional regulation of a
common pool of genes, regardless of their evolutionary trajectory. Nonetheless, we still
detected some ρ²-specific DE transcripts between salinities (Fig 3B).

To confirm that differences in gene expression (and, to a lesser extent, DNA methylation) among salinities were due to plasticity, rather than resulting from putative strong selection taking place in the polymorphic populations during the short duration of the assay, we also sequenced isogenic populations—founded from a single, presumably haploid cell from the evolved populations—subjected to different salinities. Despite the genetic homogeneity of these population, we confirmed the salinity effect on both DNA methylation and gene expression (S1 Fig). Furthermore, we did not detect any significant differences in either DNA methylation (P = 0.423) or gene expression (P = 0.520) between isogenic lines and the experimental populations they originated from (S2 Fig).

We next performed Gene Ontology (GO) enrichment analysis, to assess gene functions
involved in *D. salina* response to salinity changes and its evolution, using the functional
genome annotation constructed by Leung and colleagues [23]. Transcript functional annota-
tion analysis revealed that DE transcripts between salinities mostly involved genes associated
with cellular components and molecular functions involving the chloroplast and protein
<figure><img src="image_6.png" /><figcaption>**Fig 3. Evolution and plasticity of gene expression.** (A) Variation of gene expression levels across evolved populations. RDA plot performed on gene expression levels according to evolution conditions (colors) and salinity during plasticity assay (shapes). (B) Distribution of DE transcripts. Venn diagram describing the numbers of DE transcripts across evolutionary conditions (Predictability, light gray) identified by performing LRT (FDR < 0.05) as implemented in *DESeq2*, and between salinities (Wald test, FDR < 0.05 after BH adjustment and |log₂FC| > 1) for 3 evolutionary conditions: low (green), intermediate (blue), and high (red) predictability of environmental changes. The raw data underlying this figure are available in the Figshare repository https://doi.org/10.6084/m9 figshare.21905670. DE, differentially expressed; LRT, likelihood-ratio test; RDA, redundancy analysis.</figcaption></figure>

https://doi.org/10.1371/journal.pbio.3001895.g003
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Gene functions involved in salinity-induced plasticity and its evolution. Enriched GO terms of identified DE transcripts between salinities (Salinity) and ρ²-specific transcriptional response (ρ² × Salinity). GO categories included MF (purple), CC (orange), and BP (green) and were sorted by decreasing order of evidence within each category based on GO enrichment test *P* value (for FDR ≤ 0.1). No GO term enrichment with FDR ≤ 0.1 was detected for DE transcripts among evolutionary treatments (ρ²). The raw data underlying this figure are available in the Figshare repository https://doi.org/10.6084/m9 figshare.21905670. BP, biological process; CC, cellular component; DE, differentially expressed; GO, Gene Ontology; MF, molecular function.</figcaption></figure>

https://doi.org/10.1371/journal.pbio.3001895.g004

transport (Fig 4). While no GO term enrichment was detected for DE transcripts among evolutionary treatments (for FDR < 0.1), we found that ρ² × salinity interaction effects on gene expression essentially involved gene functions associated with different metabolic processes within the GO category “biological process” (Fig 4).

## 2.2 DNA methylation in *cis* was not associated with gene expression

We then asked whether DNA methylation could influence the expression of neighboring genes in *cis*. For each methylated cytosine, we searched for the nearest transcription start site (TSS). For this analysis, all cytosines associated to the same transcript were merged into a common gene-associated methylation region, following the same method used by Leung and colleagues [23]. Among the 3,064 methylated regions associated to a transcript, we detected 18 regions that were significantly differentially methylated among evolutionary treatments (including all comparisons of ρ² pairs), and 21 regions among salinities (for each evolutionary treatments ρ²). However, only 4 of these 39 transcript-associated regions were associated with significant DE transcripts for the same comparison, and only one involved a comparison between salinities. This suggested that, should DNA methylation have an effect on gene expression plasticity, this effect must be acting on regulation in *trans*.

## 2.3 The degree of plasticity evolved consistently across levels

We detected significant effect of salinity ($R^2_{adj} = 2.51\%$; *P* < 0.001), evolutionary treatment ρ² ($R^2_{adj} = 4.13\%$; *P* < 0.001), and their interaction ($R^2_{adj} = 1.06\%$; *P* < 0.001) on cell morphology,
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Multilevel evolution of plasticity in response to environmental predictability. The degree of plasticity of populations, as measured by the Euclidean distance between low and high salinities for (A) DNA methylation patterns at CpG context, (B) gene expression levels, and (C) cell morphology, is plotted against the predictability ρ² of the environmental fluctuations that these populations have experienced during experimental evolution. The dashed line is the regression slope, and the gray area represents the 95% confidence interval of the linear regression. The raw data underlying this figure are available in the Figshare repository https://doi.org/10.6084/m9 figshare.21905670.</figcaption></figure>

https://doi.org/10.1371/journal.pbio.3001895.g005

thus providing independent replication of the results from Leung and colleagues [5]. To investigate the consistency in the direction of experimental evolution of plasticity across levels, we then quantified the overall degree of plasticity for DNA methylation, gene expression, and cell morphology, measured as the Euclidean distance between the multivariate means at low versus high salinities. When regressing the magnitude of plasticity against the environmental predictability of the experimental evolution treatment, we consistently found a positive relationship, with weak to moderate evidence (sensu [29]) for DNA methylation (P = 0.057; Fig 5A), but strong evidence for gene expression (P = 0.018; Fig 5B) and cell morphology (P = 0.004; Fig 5C). Hence, populations that have experienced less predictable environments during experimental evolution have evolved reduced plasticity, not only for cell morphology (as shown by Leung and colleagues [5]), but also for DNA methylation and gene expression, two important molecular mechanisms potentially underlying phenotypic plasticity.

**3. Discussion**

We aimed to understand whether the evolution of plasticity in response to environmental predictability is a phenomenon that can be observed at different hierarchical levels. We therefore exposed experimental populations of the microalga *D. salina* to randomly fluctuating salinity with controlled predictability, before assessing how their plasticity evolved at 3 levels: DNA methylation, gene expression, and cell morphology. Our results highlighted important aspects of the molecular mechanisms of plasticity and its evolution.

**3.1 Gene expression plasticity and its regulation**

A wealth of studies have investigated the molecular underpinnings of phenotypic plasticity, and most of them have identified gene expression as a key mechanism of phenotypic changes
[10–14]. A crucial step in gene expression is transcription to mRNA, which explains a large fraction of variation in protein abundance [30] and, hence, in more integrated phenotype. Here, we found that salinity-induced transcriptional plasticity of *D. salina* involved largely overlapping DE transcripts among different populations, indicating that this species has specific genes to cope with salinity changes in its environments, consistent with its halotolerant ecology [31,32]. In addition, we detected significant $\rho^2 \times$ salinity interaction for gene expression, indicative of evolution of transcriptional plasticity with respect to environmental predictability. While only ca. 30% of the total transcripts were successfully assigned to at least one GO term [23], we were able to assess some gene functions involved in *D. salina* response to salinity. On the one hand, common salinity responses mostly involved genes associated with chloroplast and membrane-associated functions, like protein transport. On the other hand, $\rho^2$-specific salinity responses entailed biological processes necessary for the production of energy for various functions of the organisms [33–35]. From a functional standpoint, these results suggest that *D. salina* preserves its complex molecular machineries for coping consistently with osmotic stress, while evolutionary differentiation in salinity-induced plasticity may involve potential energy allocation trade-offs among different functions.

Regulation of gene expression involves different mechanisms, including epigenetic processes [36,37]. Here, we investigated one such epigenetic process, namely DNA methylation. We detected DMRs between salinities in every evolutionary treatment, confirming that DNA methylation is a process that could be involved in phenotypic plasticity in this species [23]. However, we failed to detect a correlation between DMRs and the expression of proximal genes, suggesting that methylation-based regulation of gene expression plasticity for this species may be driven by *trans*, instead of *cis*, regulatory changes [38,39]. Alternatively, mechanisms of gene expression regulation could be highly stress- and organism-dependent [40]. The time lapse (24 h) between osmotic changes and DNA methylation measurements in our study, despite being long enough to allow for epigenetic marks to be modified during DNA replication, could still be too short to enable detection of significant correlations between DNA methylation and gene expression in *D. salina* (e.g., [41]). The first responses to osmotic shock in *D. salina* may implicate homeostatic regulation, relying on mechanisms other than DNA methylation [42].

Exploring the genetic basis of phenotypic plasticity and its evolution requires the identification of genes whose expression varies across environments, as well as loci contributing to this variable degree of plasticity among genotypes. Interestingly, the $\rho^2$-specific DMRs for the 3 categories of evolutionary treatments did not display any overlap, suggesting that different regions in the genome were involved in the control of salinity-induced plasticity. This could be the result of epistatic interactions, where different loci could contribute to fine-tuning the phenotype toward its optimum [43–48].

## 3.2 Multilevel evolution of plasticity

To what extent do we expect evolution to take the same course across different biological levels, along the hierarchy of traits that relates the genotype to more integrated phenotypes, and fitness? This question has received increasing attention in the evolutionary literature, where it has mostly been expressed with respect to parallel/convergent evolution in response to a constant environmental challenge [49,50]. Theoretical and empirical work has established that the amount of parallelism at each biological level depends on the degree of redundancy in the mapping from one level to the next [46–48]. This, in turn, depends on the extent to which each item at a given level is connected to one or many items at the level above, e.g., how many traits a mutation can modify (pleiotropy) [51–53]. Because mappings along hierarchies of traits
are often redundant, we do not necessarily expect to find the same degree of parallelism at the
level of genes as we find at the level of functional categories, organs, or highly integrated traits
[49,50]. And reciprocally, the degree of parallelism across levels conveys information about the
redundancies in mappings along hierarchies of traits.

Here, we applied a similar approach to a classic topic from theoretical evolutionary ecology: the evolution of phenotypic plasticity in response to environmental predictability [3–7]. Just like for evolutionary parallelism, we do not necessarily expect the evolution of plasticity in response to environmental predictability to be consistent across levels. Indeed, it was recently highlighted that the comparison of phenotypic plasticity (and its evolution) across levels can provide important information about constraints in the genotype-(environment)-phenotype map, including—but not restricted to—the abovementioned redundancies. For example, this can yield insights about the extent to which mutations and environments have similar influences on phenotypic variation, and which one is likely to drive evolution of the other [54,55]. Here, the significant DMRs and DE transcripts detected across evolutionary treatments can be interpreted as parallel evolution of DNA methylation and gene expression, respectively. Moreover, we not only found substantial plasticity at the levels of DNA methylation, gene expression, and cell morphology, but also that plasticity evolved in a consistent direction—and consistent with theoretical predictions—in response to environmental predictability across these levels. This response was, however, less pronounced for DNA methylation, which only displayed marginally significant relationship between total plasticity and predictability (Fig 5A), but with substantial variation in the number of salinity DMRs among evolutionary treatments (Fig 2B).

Such consistent evolution of plasticity across levels may indicate that many unmeasured
integrated traits exhibit similar evolutionary responses to environmental predictability as
those we have measured here (cell morphology), such that their transcriptomic (and, to some
extent, epigenetic) basis is less subject to redundancies. Interestingly, in a study focusing on
other lines from this same experiment [56], we found that plasticity of intracellular glycerol
concentration—a major mechanism for osmoregulation in this species [57,58]—also evolved
in response to environmental predictability, but in the opposite direction (higher plasticity in
populations from unpredictable environments). However, this was mostly explained by the
higher glycerol levels maintained by populations from unpredictable treatments (relative to
those from predictable treatments) at intermediate and high salinities (but not at low salinity),
which can be interpreted as a benefit of maintaining hyperoptimal phenotypes when fitness
functions are highly asymmetrical [59], as shown in this species with respect to salinity transi
tions [24,56]. Alternatively, consistent evolution of plasticity across levels may indicate that
selection operates, to some extent, independently at each of these hierarchical levels. For exam-
ple, expression plasticity of a given gene may partly influences fitness by itself, rather than only
via its effect on the plasticity of more integrated phenotypes. Investigating these questions
would require mechanistic studies beyond the scope of this work. Nevertheless, our finding
that plasticity can experimentally evolve at multiple levels in response to environmental
predictability, and in a consistent direction, brings much needed empirical insights into the
mechanisms underlying plasticity evolution, and highlights the usefulness of experimental evo-
lution for investigating ambitious questions at the forefront of evolutionary biology.

**4. Materials and methods**

**4.1 Experimental evolution conditions and plasticity assays**

To investigate the evolution of plasticity in responses to environmental predictability at the
molecular level, we analyzed different populations of the *D. salina* strain CCAP 19/15 that
evolved in a fluctuating environment for ca. 500 generations, from the experiment described
in Rescan and colleagues [24] and Leung and colleagues [5]. Briefly, different populations
starting from the same ancestor were exposed to randomly fluctuating salinity, with changes
every 3 or 4 generations (i.e., twice a week, assuming one generation per day [31]). Each of
these populations were subjected to independent random time series of salinity changes, over
a continuous range. The time series were characterized by the same stationary mean (μ = 2.4
M [NaCl]) and variance (σ = 1) but differed in how salinity at a given time depends on the
salinity prior the transfer (i.e., temporal autocorrelation ρ of salinity). The predictability of
salinity changes for a given time series was assessed by ρ², the proportion of temporal variance
in salinity explained by the previous salinity, prior to the latest transfer (see Rescan and col-
leagues [24] and Leung and colleagues [5] for detailed protocol). We specifically analyzed 9
populations from 3 different target autocorrelation treatments (ρ = 0, -0.5, and 0.9; 3 popula-
tions per autocorrelation treatments). These autocorrelation treatments correspond to low (ρ²
= 0), intermediate (ρ² = 0.25), and high (ρ² = 0.81) predictability of environmental changes.
While similar levels of morphological plasticity were previously found between populations
that evolved in intermediate predictability changes but with different autocorrelation treat-
ments (ρ = -0.5 and ρ = 0.5; Leung and colleagues [5]), here we chose the negative autocorrela-
tion to include a larger range of ρ.

Plasticity assays were performed following the protocol described in Leung and colleagues [5]. We first acclimatized all evolved lines during 10 days in the same environmental conditions ([NaCl] = 2.4 M) to ensure that all cells were in similar physiological states and at similar population densities at the beginning of the phenotypic assays. Cells grew in suspension flasks containing artificial seawater with additional NaCl to reach the required salinity, complemented with 2% Guillard's F/2 marine water enrichment solution (Sigma; G0154-500 ML), and incubated at a constant temperature 24°C with a 12:12 h light/dark cycle with a 200-µmol m⁻² s⁻¹ light intensity. Target salinity was achieved by mixing the required volumes of hypo- ([NaCl] = 0 M) and hyper- ([NaCl] = 4.8 M) saline media, accounting for the salinity of the inoculate. At the end of the acclimation step, we transferred ca. 1 × 10⁵ cells mL⁻¹ of each populations to low ([NaCl] = 0.8 M) and high ([NaCl] = 4.0 M) salinities, for a total volume of 250 mL. After 24 h, following the salinity changes, we harvested the cells by centrifugation at 5,000 rpm for 15 min at room temperature, and cell pellets were stored at -80°C until acid nucleic extraction.

We also measured intrinsic structural parameters of cells by passing a subsample of 150 µL of each populations through a Guava EasyCyte HT flow cytometer (Luminex Corporation, TX, USA), also following the protocol described in Leung and colleagues [5]. We specifically assessed the environment-specific cell morphology using the Forward Scatter (FSC) and Side Scatter (SSC) as proxies for cells size and complexity (cytoplasmic contents) [60], respectively, and fluorescence emission at 695/50 nm band pass filter (Red-B) values for chlorophyll content [61]. The cell morphology matrix consisted of values for these 3 parameters (FSC, SSC, and Red-B), for 150 randomly sampled cells identified as alive *D. salina* for each populations.

## 4.2 Sample preparation, sequencing, and bioinformatic preprocessing

To investigate the molecular mechanisms involved in osmotic stress responses, we performed
whole-transcriptome shotgun sequencing (RNA-seq) for the comparison of gene expression
levels, and whole-genome bisulphite sequencing (WGB-seq) for the comparison of DNA
methylation variation among the 9 evolved lines and 2 environmental conditions (hypo- and
hyperosmotic). The different lines started from potentially genetically diverse population and
could thus be polymorphic at the end of experimental evolution. To assess whether genetic
variation within population affected the observed plasticity levels for gene expression or intra-
population DNA methylation variation, we also founded 3 isogenic populations (one per auto-
correlation treatment) started from a single cell using cells-sorting flow cytometry (BD
FACSAria IIu; Biosciences-US). As *D. salina* is haploid, we expected all derived cells of a given
population to be genetically identical. Each of these isogenic populations were also subjected
to the plasticity assay described above.

Total RNA extraction and purification of 24 samples ((9 lines + 3 isogenic populations) × 2 salinities) was carried out using Nucleozol, following Macherey Nagel's protocol, and whole genomic DNA was isolated according to the phenol-chloroform purification and ethanol precipitation method of Sambrook and colleagues [62]. Library construction (TruSeq RNA Library Preparation kit for RNA-seq and Swift Bioscience Accel-NGS Methyl-Seq DNA library Kit for WGB-seq) and high-throughput sequencing steps (Paired-End (PE) 2 × 150 bp, Illumina HiSeq) were performed by Genewiz (Leipzig, Germany). We performed all the bioinformatic preprocessing analyses using publicly available software implemented in the European UseGalaxy server [63].

**Methylation calling.** The WGB-seq raw reads were checked for quality using *FastQC*. Adapter and low-quality sequences were then trimmed using *Trim Galore!* Version 0.4.3.1. As specified by the Accel-NGS Methyl-seq Kit manual, additional 15 bp and 5 bp were also trimmed at the 5' and 3' extremity, respectively, to remove the tail added during library preparation and thus avoiding nonquality-related bias. Mapping was performed on the same references genomes as in RNA-seq analyses, using *Bismark Mapper* version 0.22.1 [64]. Only uniquely mapping reads were retained and PCR duplicates were removed using *Bismark Deduplicate* tool. We then extracted the methylation status from the resulting alignment files using *MethylDackel* (Galaxy Version 0.3.0.1), where only cytosines covered by a minimum of 10 reads in each library were considered, and with the option of excluding likely variant sites (i.e., minimum depth for variant avoidance of 10×, and maximum tolerated variant fraction of 0.95). In a previous study, we showed that cytosines at CpG context displayed the highest methylation levels (and variation thereof) in *D. salina* ([23]; S2 Table). Furthermore, CpG-methylations have been suggested to play a role in gene regulations [65–67] and proposed as a molecular mechanisms underlying phenotypic plasticity [15,17]. We thus investigated the genomic DNA methylation patterns in response to salinity across the different evolutionary treatments for cytosines at the CpG context. The high bisulfite conversion rate (>99%) was assessed by Genewiz, by spiking in unmethylated lambda DNA in 3 randomly chosen libraries.

**Gene expression analyses.** The RNA-seq raw reads were checked for quality using *FastQC* version 0.72 [68] and subjected to adapter trimming and quality filtering using *Trim Galore!* version 0.4.3.1 [69]. Additional 12 bp and 3 bp were also removed at the 5' and 3' extremity, respectively, to avoid bias not directly related to adapter sequences or basecall quality according to *FastQC* outputs, and only reads with a minimum length of 50 bp were retained. We used the reference nuclear (Dunsall 1 v. 2, GenBank accession: GCA_002284615.2), chloroplastic (GenBank accession: GQ250046), and mitochondrial (GenBank accession: GQ250045) genomes of *D. salina* strain CCAP 19/18 (closely related to CCAP 19/15 used here) for trimmed reads alignment, using *HISAT2* version 2.1.0 [70] with default parameters for PE reads and spliced alignment option. We finally quantified the number of reads per transcript with *FeatureCounts* version 2.0.1 [71] using the alignment files from *HISAT2* and the de novo transcript annotation produced for this species from Leung and colleagues [23] using *StringTie* v. 2.1.1 [72] and including the evolved lines of the current study.
## 4.3 Statistical analyses

**Differential DNA methylation and gene expression.** We used Bioconductor's *methylKit* package [73] to identify DMRs, i.e., nonoverlapping 100 bp windows with methylation levels that varied significantly among evolutionary treatments or between salinities within each treatments. The significance of calculated methylation differences was determined using Fisher's exact tests. We used the BH adjustment of *P* values (FDR < 0.05) and methylation difference cutoffs of 20%. Similarly, the differential gene expression analyses were performed with the Bioconductor's package *DESeq2* version 1.30.1 [74]. We identified DE transcripts among the evolutionary treatments, salinity, and their interaction, by building a general linear model as implemented in *DESeq2*. We used the Wald test when comparing 2 salinities, and transcripts with FDR < 0.05 (*P* values after BH adjustment) and |log₂FC| > 1 were considered as DE. Significance of evolution (predictability *ρ*²) and evolution of plasticity (*ρ*² × salinity interaction) were assessed using likelihood ratio tests (LRTs) comparing models with and without the corresponding terms [74].

We then used the GO assignments from Leung and colleagues [23] to classify the functions of *D. salina* transcripts and to functionally annotate the identified DE transcripts. Enriched GO terms of the DE transcripts were identified using the ‘classic’ algorithm from the *topGO* R package [75] and based on *p*-value generated using Fisher’s exact method, for the 3 GO categories, i.e., molecular function, cellular component, and biological process. GO terms were then sorted by decreasing order of evidence within each category, based on the GO enrichment test *P* value after BH adjustment, and we showed the most probable gene function candidates with a threshold of FDR ≤ 0.1.

**Evolution of the degree plasticity at multiple levels.** Since phenotypic plasticity can evolve as an adaptation to fluctuating and predictable environment [3–7], we wished to quantify to what extent environmental predictability during experimental evolution contributed to key mechanisms underlying phenotypic plasticity. We first applied a RDA [76], computed with the function *rda*() from the *vegan* R package [77], to quantify the proportion of the total epigenetic and gene expression variation that are significantly explained by the evolutionary treatments (predictability *ρ*²) and assay environment (salinity). We performed a RDA using the table of DNA methylation levels of 100 bp windows regions, or the table of *rlog* transformed transcript count, as response variables, and the evolutionary treatments (predictability *ρ*²), salinity, and the *ρ*² × salinity interaction, as explanatory variables, with population identity as covariates to account for paired samples between salinities. We treated the target *ρ*² (low, intermediate, or high) and salinity (low or high) as categorical variables. We quantified the contribution to the total genetic variation using the adjusted *R*² and tested its significance by ANOVA-like permutation tests using 999 randomizations of the data [78]. For both RNA-seq and WGB-seq data, we presented the results of the RDA ordinations as a biplot, to visualize the variation among samples along its major axes. We performed all multivariate statistical analyses with the *vegan* R package [77].

To investigate how the overall degree of plasticity at multiple levels evolved in our experiment, we followed the protocol described in Leung and colleagues [5], which we here applied to DNA methylation level, gene expression, and cells morphological variation. (Note that the morphological measurements were replicated here, rather than just reproducing values from Leung and colleagues [5].) First, we assessed the degree of plasticity of each experimental population by computing the Euclidean distance between the multivariate means from the CpG-methylation level table, the *rlog* transformed transcript count table, or the raw cells morphological data measured, at low ([NaCl] = 0.8 M) versus high ([NaCl] = 4.0 M) salinity. The realized autocorrelation of a given time series can vary, to some extent, from its long-term
stationary expectation, because of the randomness of the stochastic process in finite time. For a more quantitative relationship between plasticity and environmental predictability, we thus used the realized (rather than the target) $\rho^2$ as index of environmental predictability (as also done in [5]). For each time series, we thus calculated the realized environmental autocorrelation $\rho$ as the correlation between salinities at 2 subsequent transfers. We then tested whether the degree of plasticity evolved according to environmental predictability, by regressing the Euclidean distance of plastic changes against the realized $\rho^2$. The linear regression t test was applied to determine whether the slope of the regression line differs significantly from zero.

**Supporting information**

**S1 Table. WGB-seq mapping statistics.** For each sample, population ID, realized $\rho^2$, and salinity are detailed in S2 Table. We used a genome size of 350 M bp to estimate depth of coverage. Total C's analyzed were after reads deduplication and number of methylated C's and percentage of methylation are given by cytosine context.
(XLSX)

**S2 Table. RNA-seq mapping statistics.**
(XLSX)

**S1 Fig. Salinity effect on DNA methylation and gene expression in isogenic populations.**
We founded 3 populations from single isolated cells from 3 evolved populations, following the protocol in Leung and colleagues [5]. As *D. salina* is haploid, a population founded from a single cell is expected to be isogenic. (A) Heat-maps of WGB-seq analysis for DMRs between salinities (*n* = 27). Each row represents a DMR, and column names are the population identity. Relative DNA methylation levels vary from blue (under-methylated) to red (over-methylated), as shown on the right-hand side of the heat-maps. Dendrograms on the top result from a hierarchical clustering analysis using the Euclidean distance of DNA methylation level among populations. (B) Volcano plot illustrating significant (for FDR < 0.05 and |Log₂FC| > 1) and nonsignificant DE transcripts between salinities as red and gray points, respectively. Salinity effect was assessed by comparing three isogenic populations (i.e., found from a single cell). The raw data underlying this figure are available in the Figshare repository https://doi.org/10.6084/m9 figshare.21905670. DE, differentially expressed; DMR, differentially methylated region; WGB-seq, whole-genome bisulfite sequencing.
(PDF)

**S2 Fig. Comparison of isogenic and experimental populations.** Principal component analysis (PCA) of (A) DNA methylation and (B) gene expression levels among isogenic (open symbols) and experimental non-isogenic (filled symbols) populations. The raw data underlying this figure are available in the Figshare repository https://doi.org/10.6084/m9 figshare.21905670.
(PDF)

**Acknowledgments**

We thank Rachel Steward for useful feedback on previous versions of this manuscript.

**Author Contributions**

**Conceptualization:** Christelle Leung, Daphné Grulois, Luis-Miguel Chevin.

**Formal analysis:** Christelle Leung.

**Funding acquisition:** Christelle Leung, Luis-Miguel Chevin.
**Methodology:** Christelle Leung, Daphné Grulois, Leandro Quadrana.

**Supervision:** Luis-Miguel Chevin.

**Writing – original draft:** Christelle Leung, Luis-Miguel Chevin.

**Writing – review & editing:** Christelle Leung, Daphné Grulois, Leandro Quadrana, Luis-Miguel Chevin.

## References

1. Vedder O, Bouwhuis S, Sheldon BC. Quantitative assessment of the importance of phenotypic plasticity in adaptation to climate change in wild bird populations. PLoS Biol. 2013; 11:e1001605. https://doi.org/10.1371/journal.pbio.1001605 PMID: 23874152

2. Chevin L-M, Lande R, Mace GM. Adaptation, plasticity, and extinction in a changing environment: towards a predictive theory. PLoS Biol. 2010; 8:e1000357. https://doi.org/10.1371/journal.pbio.1000357 PMID: 20463950

3. Dey S, Proulx SR, Teotonio H. Adaptation to temporally fluctuating environments by the evolution of maternal effects. PLoS Biol. 2016; 14:e1002388. https://doi.org/10.1371/journal.pbio.1002388 PMID: 26910440

4. Gavrilets S, Scheiner SM. The genetics of phenotypic plasticity. V. Evolution of reaction norm shape. J Evol Biol. 1993; 6:31–48.

5. Leung C, Rescan M, Grulois D, Chevin LM. Reduced phenotypic plasticity evolves in less predictable environments. Ecol Lett. 2020; 23:1664–1672. https://doi.org/10.1111/ele.13598 PMID: 32869431

6. Levins R. Theory of fitness in a heterogeneous environment. II. Developmental flexibility and niche selection. Am Nat. 1963; 97:75–90.

7. King JG, Hadfield JD. The evolution of phenotypic plasticity when environments fluctuate in time and space. Evol Lett. 2019; 3:15–27. https://doi.org/10.1002/evl3.100 PMID: 30788139

8. Botero CA, Weissing FJ, Wright J, Rubenstein DR. Evolutionary tipping points in the capacity to adapt to environmental change. Proc Natl Acad Sci U S A. 2015; 112:184–189. https://doi.org/10.1073/pnas.1408589111 PMID: 25422451

9. Tufto J. Genetic evolution, plasticity, and bet-hedging as adaptive responses to temporally autocorrelated fluctuating selection: a quantitative genetic model. Evolution. 2015; 69:2034–2049. https://doi.org/10.1111/evo.12716 PMID: 26140293

10. Aubin-Horth N, Renn SC. Genomic reaction norms: using integrative biology to understand molecular mechanisms of phenotypic plasticity. Mol Ecol. 2009; 18:3763–3780. https://doi.org/10.1111/j.1365-294X.2009.04313.x PMID: 19732339

11. Pavey SA, Collin H, Nosil P, Rogers SM. The role of gene expression in ecological speciation. Ann N Y Acad Sci. 2010; 1206:110–129. https://doi.org/10.1111/j.1749-6632.2010.05765.x PMID: 20860685

12. Beldade P, Mateus ARA, Keller RA. Evolution and molecular mechanisms of adaptive developmental plasticity. Mol Ecol. 2011; 20:1347–1363. https://doi.org/10.1111/j.1365-294X.2011.05016.x PMID: 21342300

13. Gibert J-M, Mouchel-Vielh E, De Castro S, Peronnet F. Phenotypic plasticity through transcriptional regulation of the evolutionary hotspot gene tan in *Drosophila melanogaster*. PLoS Genet. 2016; 12: e1006218.

14. Monteiro A, Tong X, Bear A, Liew SF, Bhardwaj S, Wasik BR, et al. Differential expression of ecdysone receptor leads to variation in phenotypic plasticity across serial homologs. PLoS Genet. 2015; 11: e1005529. https://doi.org/10.1371/journal.pgen.1005529 PMID: 26405828

15. Angers B, Castonguay E, Massicotte R. Environmentally induced phenotypes and DNA methylation: how to deal with unpredictable conditions until the next generation and after. Mol Ecol. 2010; 19:1283–1295. https://doi.org/10.1111/j.1365-294X.2010.04580.x PMID: 20298470

16. Holoch D, Moazed D. RNA-mediated epigenetic regulation of gene expression. Nat Rev Genet. 2015; 16:71–84. https://doi.org/10.1038/nrg3863 PMID: 25554358

17. Bollati V, Baccarelli A. Environmental epigenetics. Heredity. 2010; 105:105–112. https://doi.org/10.1038/hdy.2010.2 PMID: 20179736

18. Foret S, Kucharski R, Pellegrini M, Feng S, Jacobsen SE, Robinson GE, et al. DNA methylation dynamics, metabolic fluxes, gene splicing, and alternative phenotypes in honey bees. Proc Natl Acad Sci. 2012; 109:4968–4973. https://doi.org/10.1073/pnas.1202392109 PMID: 22416128
19. Shen X, De Jonge J, Forsberg SK, Pettersson ME, Sheng Z, Hennig L, et al. Natural CMT2 variation is associated with genome-wide methylation changes and temperature seasonality. PLoS Genet. 2014; 10:e1004842. https://doi.org/10.1371/journal.pgen.1004842 PMID: 25503602

20. Angers B, Perez M, Menicucci T, Leung C. Sources of epigenetic variation and their applications in natural populations. Evol Appl. 2020; 13:1262–1278. https://doi.org/10.1111/eva.12946 PMID: 32684958

21. Gerken AR, Eller OC, Hahn DA, Morgan TJ. Constraints, independence, and evolution of thermal plasticity: probing genetic architecture of long-and short-term thermal acclimation. Proc Natl Acad Sci. 2015; 112:4399–4404. https://doi.org/10.1073/pnas.1503456112 PMID: 25805817

22. Leung C, Breton S, Angers B. Facing environmental predictability with different sources of epigenetic variation. Ecol Evol. 2016; 6:5234–5245. https://doi.org/10.1002/ece3.2283 PMID: 27551379

23. Leung C, Grulois D, Chevin LM. Plasticity across levels: relating epigenomic, transcriptomic, and phenotypic responses to osmotic stress in a halotolerant microalga. Mol Ecol. 2022; 31:4672–4687. https://doi.org/10.1111/mec.16542 PMID: 35593517

24. Rescan M, Grulois D, Ortega-Aboud E, Chevin L-M. Phenotypic memory drives population growth and extinction risk in a noisy environment. Nature Ecology & Evolution. 2020; 4:193–201.

25. Mohandas T, Sparkes R, Shapiro L. Reactivation of an inactive human X chromosome: evidence for X inactivation by DNA methylation. Science. 1981; 211:393–396. https://doi.org/10.1126/science.6164095 PMID: 6164095

26. Jaenisch R, Bird A. Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals. Nat Genet. 2003; 33:245–254. https://doi.org/10.1038/ng1089 PMID: 12610534

27. Kucharski R, Maleszka J, Foret S, Maleszka R. Nutritional control of reproductive status in honeybees via DNA methylation. Science. 2008; 319:1827–1830. https://doi.org/10.1126/science.1153069 PMID: 18339900

28. Matsumoto Y, Buemio A, Chu R, Vafaee M, Crews D. Epigenetic Control of Gonadal Aromatase (cyp19a1) in Temperature-Dependent Sex Determination of Red-Eared Slider Turtles. PLoS ONE. 2013; 8:e63599. https://doi.org/10.1371/journal.pone.0063599 PMID: 23762231

29. Muff S, Nilsen EB, O'Hara RB, Nater CR. Rewriting results sections in the language of evidence. Trends Ecol Evol. 2021; 37:203–210. https://doi.org/10.1016/j.tree.2021.10.009 PMID: 34799145

30. Wang Z-Y, Leushkin E, Liechti A, Ovchinnikova S, Mößinger K, Brüning T, et al. Transcriptome and translatome co-evolution in mammals. Nature. 2020; 588:642–647. https://doi.org/10.1038/s41586-020-2899-z PMID: 33177713

31. Ben-Amotz A, Polle JEW, Subba RDV. The alga *Dunaliella*: Biodiversity, physiology, genomics and biotechnology. Enfield, NH: Science Publishers; 2009.

32. Oren A. A hundred years of *Dunaliella* research: 1905–2005. Saline Syst. 2005; 1:2.

33. Kirst G. Salinity tolerance of eukaryotic marine algae. Annu Rev Plant Biol. 1990; 41:21–53.

34. Müller M, Mentel M, van Hellemond JJ, Henze K, Woehle C, Gould SB, et al. Biochemistry and evolution of anaerobic energy metabolism in eukaryotes. Microbiol Mol Biol Rev. 2012; 76:444–495. https://doi.org/10.1128/MMBR.05024-11 PMID: 22688819

35. Atteia A, van Lis R, Tielens AG, Martin WF. Anaerobic energy metabolism in unicellular photosynthetic eukaryotes. Biochim Biophys Acta. 2013;210–223. https://doi.org/10.1016/j.babio.2012.08.002 PMID: 22902601

36. Zilberman D. An evolutionary case for functional gene body methylation in plants and animals. Genome Biol. 2017; 18:1–3.

37. Neri F, Rapelli S, Krepelova A, Incarnato D, Parlato C, Basile G, et al. Intragenic DNA methylation prevents spurious transcription initiation. Nature. 2017; 543:72–77. https://doi.org/10.1038/nature21373 PMID: 28225755

38. Necsulea A, Kaessmann H. Evolutionary dynamics of coding and non-coding transcriptomes. Nat Rev Genet. 2014; 15:734–748. https://doi.org/10.1038/nrg3802 PMID: 25297727

39. Albert FW, Bloom JS, Siegel J, Day L, Kruglyak L. Genetics of trans-regulatory variation in gene expression. Elife. 2018; 7:e35471. https://doi.org/10.7554/eLife.35471 PMID: 30014850

40. De Nadal E, Ammerer G, Posas F. Controlling gene expression in response to stress. Nat Rev Genet. 2011; 12:833–845. https://doi.org/10.1038/nrg3055 PMID: 22048664

41. Putnam HM, Davidson JM, Gates RD. Ocean acidification influences host DNA methylation and phenotypic plasticity in environmentally susceptible corals. Evol Appl. 2016; 9:1165–1178. https://doi.org/10.1111/eva.12408 PMID: 27695524

42. López-Maury L, Marguerat S, Bähler J. Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. Nat Rev Genet. 2008; 9:583–593. https://doi.org/10.1038/nrg2398 PMID: 18591982
43. Pickett FB, Meeks-Wagner DR. Seeing double: appreciating genetic redundancy. Plant Cell. 1995; 7:1347. https://doi.org/10.1105/tpc.7.9.1347 PMID: 8589620

44. Nordin K, Heino P, Palva ET. Separate signal pathways regulate the expression of a low-temperature-induced gene in *Arabidopsis thaliana* (L.) Heynh. Plant Mol Biol. 1991; 16:1061-1071.

45. Callahan HS, Pigliucci M, Schlichting CD. Developmental phenotypic plasticity: where ecology and evolution meet molecular biology. Bioessays. 1997; 19:519-525. https://doi.org/10.1002/bies.950190611 PMID: 9204769

46. Barghi N, Hermisson J, Schlötterer C. Polygenic adaptation: a unifying framework to understand positive selection. Nat Rev Genet. 2020; 21:769-781. https://doi.org/10.1038/s41576-020-0250-z PMID: 32601318

47. Barghi N, Tobler R, Nolte V, Jakšić AM, Mallard F, Otte KA, et al. Genetic redundancy fuels polygenic adaptation in Drosophila. PLoS Biol. 2019; 17:e3000128. https://doi.org/10.1371/journal.pbio.3000128 PMID: 30716062

48. Láruson ÁJ, Yeaman S, Lotterhos KE. The importance of genetic redundancy in evolution. Trends Ecol Evol. 2020; 35:809-822. https://doi.org/10.1016/j.tree.2020.04.009 PMID: 32439075

49. Tenaillon O, Rodríguez-Verdugo A, Gaut RL, McDonald P, Bennett AF, Long AD, et al. The molecular diversity of adaptive convergence. Science. 2012; 335:457-461. https://doi.org/10.1126/science.1212986 PMID: 22282810

50. Lenormand T, Chevin L-M, Bataillon T, Ramsey G, Pence C. Parallel evolution: what does it (not) tell us and why is it (still) interesting. Chance in evolution. 2016:196-220.

51. Chevin LM, Martin G, Lenormand T. Fisher's model and the genomics of adaptation: restricted pleiotropy, heterogenous mutation, and parallel evolution. Evolution. 2010; 64:3213-3231. https://doi.org/10.1111/j.1558-5646.2010.01058.x PMID: 20662921

52. Storz JF. Causes of molecular convergence and parallelism in protein evolution. Nat Rev Genet. 2016; 17:239-250. https://doi.org/10.1038/nrg.2016.11 PMID: 26972590

53. Rennison DJ, Peichel CL. Pleiotropy facilitates parallel adaptation in sticklebacks. Mol Ecol. 2022; 31:1476-1486. https://doi.org/10.1111/mec.16335 PMID: 34997980

54. Chevin L-M, Leung C, Le Rouzic A, Uller T. Using phenotypic plasticity to understand the structure and evolution of the genotype-phenotype map. Genetica. 2022; 150:209-221. https://doi.org/10.1007/s10709-021-00135-5 PMID: 34617196

55. Brun-Usan M, Rago A, Thies C, Uller T, Watson RA. Development and selective grain make plasticity 'take the lead' in adaptive evolution. BMC Ecol Evol. 2021; 21:1-17.

56. Rescan M, Leurs N, Grulois D, Chevin LM. Experimental evolution of environmental tolerance, acclimation, and physiological plasticity in a randomly fluctuating environment. Evol Lett. 2022; 6:522-536. https://doi.org/10.1002/evl3.306 PMID: 36579167

57. Ben-Amotz A, Avron M. The role of glycerol in the osmotic regulation of the halophilic alga *Dunaliella parva*. Plant Physiol. 1973; 51:875-878.

58. Chen H, Jiang JG. Osmotic responses of *Dunaliella* to the changes of salinity. J Cell Physiol. 2009; 219:251-258.

59. Urban MC, Bürger R, Bolnick DI. Asymmetric selection and the evolution of extraordinary defences. Nat Commun. 2013; 4:2085. https://doi.org/10.1038/ncomms3085 PMID: 23820378

60. Adan A, Alizada G, Kiraz Y, Baran Y, Nalbant A. Flow cytometry: Basic principles and applications. Crit Rev Biotechnol. 2017; 37:163-176. https://doi.org/10.3109/07388551.2015.1128876 PMID: 26767547

61. Papageorgiou GC. Fluorescence of photosynthetic pigments in vitro and in vivo. In: Papageorgiou GC, Govindjee, editors. Chlorophyll a fluorescence: A signature of photosynthesis. Advances in photosynthesis and respiration, Volume 19. Dordrecht: Springer; 2004. p. 43-63.

62. Sambrook J, Fritsch EF, Maniatis T. Molecular cloning: a laboratory manual. 2nd ed. New York: Cold Spring Harbor Laboratory Press; 1989.

63. Afgan E, Baker D, Batut B, Van Den Beek M, Bouvier D, Čech M, et al. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. Nucleic Acids Res. 2018; 46: W537-W544. https://doi.org/10.1093/nar/gky379 PMID: 29790989

64. Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011; 27:1571-1572. https://doi.org/10.1093/bioinformatics/btr167 PMID: 21493656

65. Zhang H, Lang Z, Zhu J-K. Dynamics and function of DNA methylation in plants. Nat Rev Mol Cell Biol. 2018; 19:489-506. https://doi.org/10.1038/s41580-018-0016-z PMID: 29784956

66. Chen X, Schönberger B, Menz J, Ludewig U. Plasticity of DNA methylation and gene expression under zinc deficiency in Arabidopsis roots. Plant Cell Physiol. 2018; 59:1790-1802. https://doi.org/10.1093/pcp/pcy100 PMID: 29800330
67. Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010; 11:204-220. https://doi.org/10.1038/nrg2719 PMID: 20142834

68. Andrews S. FastQC: a quality control tool for high throughput sequence data. Cambridge, United Kingdom: Babraham Bioinformatics, Babraham Institute; 2010.

69. Krueger F. Trim Galore. In A wrapper tool around Cutadapt and FastQC to consistently apply quality and adapter trimming to FastQ files; 2015.

70. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015; 12:357-360. https://doi.org/10.1038/nmeth.3317 PMID: 25751142

71. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014; 30:923-930. https://doi.org/10.1093/bioinformatics/btt656 PMID: 24227677

72. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015; 33:290-295. https://doi.org/10.1038/nbt.3122 PMID: 25690850

73. Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012; 13:1-9. https://doi.org/10.1186/gb-2012-13-10-r87 PMID: 23034086

74. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014; 15:1-21. https://doi.org/10.1186/s13059-014-0550-8 PMID: 25516281

75. Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for gene ontology. R package version 2. 2010:2010.

76. Borcard D, Legendre P, Drapeau P. Partialling out the spatial component of ecological variation. Ecol- ogy. 1992; 73:1045-1055.

77. Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: Community Ecol- ogy Package. 2020. Available from: http://CRAN.R-project.org/package=vegan [cited May 2021].

78. Legendre P, Legendre L. Numerical ecology. 2nd English ed. Amsterdam: Elsevier Science; 1998.