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

HAL
open science

# Plasticity across levels: Relating epigenomic, transcriptomic, and phenotypic responses to osmotic stress in a halotolerant microalga

Christelle Leung, Daphné Grulois, Luis-miguel Chevin

► To cite this version:

Christelle Leung, Daphné Grulois, Luis-miguel Chevin. Plasticity across levels: Relating epigenomic, transcriptomic, and phenotypic responses to osmotic stress in a halotolerant microalga. Molecular Ecology, 2022, 31 (18), pp.4672-4687. 10.1111/mec.16542. hal-04913451

HAL Id: hal-04913451

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

Submitted on 5 Mar 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Plasticity across levels: Relating epigenomic, transcriptomic, and phenotypic responses to osmotic stress in a halotolerant microalga

Christelle Leung ID |

Daphné Grulois ID |

Luis-Miguel Chevin ID |

CEFE, Université de Montpellier, CNRS, EPHE, IRD, Montpellier, France

**Correspondence**

Christelle Leung, CEFE, Université de Montpellier, CNRS, EPHE, IRD, Montpellier, France.

Email: christelle.leung@umontreal.ca

**Funding information**

Fonds de Recherche du Québec - Nature et Technologies; H2020 European Research Council, Grant/Award Number: Grant 678140-FluctEvol; Agence Nationale de la Recherche, Grant/Award Number: ANR-21-ERCC-0001-01 InterAdapt

Handling Editor: Mitchell Cruzan

**Abstract**

Phenotypic plasticity, the ability of a given genotype to produce alternative phenotypes in response to its environment of development, is an important mechanism for coping with variable environments. While the mechanisms underlying phenotypic plasticity are diverse, their relative contributions need to be investigated quantitatively to better understand the evolvability of plasticity across biological levels. This requires relating plastic responses of the epigenome, transcriptome, and organismal phenotype, and investigating how they vary with the genotype. Here we carried out this approach for responses to osmotic stress in *Dunaliella salina*, a green microalga that is a model organism for salinity tolerance. We compared two strains that show markedly different demographic responses to osmotic stress, and showed that these phenotypic responses involve strain- and environment-specific variation in gene expression levels, but a relative low-albeit significant-effect of strain × environment interaction. We also found an important genotype effect on the genome-wide methylation pattern, but little contribution from environmental conditions to the latter. However, we did detect a significant marginal effect of epigenetic variation on gene expression, beyond the influence of genetic differences on epigenetic state, and we showed that hypomethylated regions are correlated with higher gene expression. Our results indicate that epigenetic mechanisms are either not involved in the rapid plastic response to environmental change in this species, or involve only few changes in *trans* that are sufficient to trigger concerted changes in the expression of many genes, and phenotypic responses by multiple traits.

**KEYWORDS**

DNA methylation, gene expression, phenotypic plasticity, population dynamics, programmed cell death, RNA-sequencing, whole-genome bisulphite sequencing

## 1 | INTRODUCTION

Phenotypic plasticity, the ability of a given genotype to produce alternative phenotypes in response to its environment of development,

can evolve as an adaptation to a variable environment that is sufficiently predictable (de Jong, 1999; Gavrilets & Scheiner, 1993; Levins, 1963; Tufto, 2015; Via & Lande, 1985), but it may also reflect non- or maladaptive phenotypic responses to the environment,

This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2022 The Authors. *Molecular Ecology* published by John Wiley & Sons Ltd.
caused by various constraints (Ghalambor et al., 2007, 2015). To advance our understanding of the mechanisms and evolution of plasticity, we can rely on organisms that are able to tolerate a particularly broad range of a challenging environment, because phenotypic responses to these environments are likely to have played a prominent role in their evolutionary history.

An interesting example is provided by halotolerant organisms, such as the microalgae *Dunaliella salina*. This unicellular organism has the broadest salinity tolerance of all known eukaryotes, ranging in the laboratory from near freshwater to salt saturation (Ben-Amotz et al., 2009). This huge salinity range is afforded notably by outstanding physiological mechanisms of osmotic regulation (primarily through glycerol metabolism), great morphological flexibility, and unique specificities of its cell membrane (Ben-Amotz et al., 2009; Ginzburg, 1988; Pick, 2002). For these reasons, *Dunaliella* has been proposed as a model organism for investigating salinity tolerance in plants (Cowan et al., 1992). More generally, this also makes it a good model for studying the mechanisms of adaptive phenotypic plasticity, as its phenotypic responses to salt have played a key ecological role in its evolutionary history, and were even shown recently to evolve in the laboratory in response to the predictability of environmental fluctuations, as predicted by theory (Leung et al., 2020).

Most existing theoretical predictions about the evolution of plasticity come from theoretical models that rely on simplifying assumptions about the mechanisms and inheritance of plasticity, often assuming quantitative genetic inheritance without explicit loci (de Jong, 1990; Gavrilets & Scheiner, 1993; Via & Lande, 1985). These assumptions were made for mathematical convenience, but also out of lack of empirical evidence on the mechanisms of plasticity. For instance, the classic debate of the 1990s about the putative existence of plasticity genes, and the relative importance of allelic sensitivity versus gene regulation in plasticity (Scheiner, 1993; Schlichting & Pigliucci, 1993; Via, 1993a, 1993b), was eventually settled because these alternatives were essentially undistinguishable based on standard quantitative genetics (both theoretical and empirical), which only rests on analysis of phenotype distributions (de Jong, 1995; Via et al., 1995). However, when explicitly considering the loci that contribute to the evolution of plasticity, the details of its genetic architecture do matter. For instance, the outcome of evolution may differ depending on whether the genes that influence variation in plasticity also cause variation in the "nonplastic component of the trait", or more precisely the phenotypic value in a reference environment (de Jong & Gavrilets, 2000; Scheiner & Holt, 2012). More generally, the way molecular mechanisms influence the expression and inheritance of plasticity at different levels of the organism are probably key to how plasticity can evolve at these levels, and should thus be investigated more systematically.

Recent advances have demonstrated that the molecular mechanisms contributing to phenotypic plasticity are diverse, but that many of them rely on variation in gene expression with the environment (Beldade et al., 2011; Gibert et al., 2016; Monteiro et al., 2015). Such variation in gene expression may in turn be regulated by hormones and/or epigenetic processes, that is, a set of enzyme-mediated modifications resulting in the alteration of gene expression without any change in DNA sequences. Environmentally induced epigenetic

variation has been proposed as a molecular mechanism underlying phenotypic plasticity (Angers et al., 2010; Beldade et al., 2011; Bollati & Baccarelli, 2010). However, epigenetic variation can also result from other sources, most importantly genetic variation (Angers et al., 2020; Leung et al., 2016), and it is unclear to what extent plasticity in gene expression results from environmentally induced epigenetic variation.

To address these questions, we here investigated rapid responses to osmotic stress of two closely related strains of the microalga *Dunaliella salina*, at three biological levels: the epigenome, the transcriptome, and macroscopic phenotypes. Not only does the outstanding ecology of this species make it ideal for studying plasticity, but it also has a reference genome since 2017 (Polle et al., 2017). This has recently led to work in comparative genomics aiming at identifying gene families involved in adaptation to salt (Polle et al., 2020), as well as investigations of plastic transcriptomic responses to salinity (Fang et al., 2017; He et al., 2020; Zhao et al., 2013), and specific epigenetic mechanisms such as small RNA (Lou et al., 2020). However, there has been no attempt to connect plastic responses to salinity at different levels (epigenome, transcriptome, and macroscopic phenotypes) across different genotypes, to decipher the underlying mechanisms of plasticity and their putative genetic variation. Our goal is thus two-fold: (1) to further our understanding of the molecular mechanisms of plasticity in a model organism for salinity tolerance, including by shedding light on the (yet little investigated) contribution of DNA methylation; and (2) to more generally investigate how plastic and genetic variation are related at multiple levels of the organisms.

Our approach rests on comparing DNA methylation levels, gene expression levels, and demographic phenotypes, across environments and genotypes. Detecting a marginal effect of the environment on methylation patterns would confirm the environment as a source of epigenetic variation. Furthermore, considering epigenetic processes as an intermediate step between the genotype and the phenotype, we expected to detect a substantial contribution of the genotype to both epigenetic and phenotypic variation. And under the hypothesis that epigenetic processes play a role in fine-tuning gene expression, we expected to detect a correlation between the DNA methylation level and gene expression level. If epigenetic states and gene expression levels jointly change in response to environment, then this implies a role of epigenetics in gene-expression plasticity, which is thought to underlie the plasticity of most higher-level organismal traits (Beldade et al., 2011; Gibert et al., 2016; Monteiro et al., 2015). Finally, a significant genotype × environment interaction on gene expression and epigenetics would indicate an evolutionary potential for plasticity at these basal phenotypic levels of the organism.

## 2 | MATERIALS AND METHODS

### 2.1 | Study strains and culture conditions

We compared two genetically related strains of *Dunaliella salina* (CCAP 19/12 and CCAP 19/15) originating from the Culture Collection of Algae and Protozoan (UK). For each of these strains, we obtained 10 different lines that were propagated independently
in our laboratory for over two years (Leung et al., 2020; Rescan et al., 2020). Our standard growth conditions are 50ml suspension flasks (CELLSTAR; VWR 392-0016) containing artificial seawater with additional NaCl, complemented with 2% Guillard's F/2 marine water enrichment solution (Sigma; G0154-500ml), for a total 25ml (including the inoculate), incubated at a constant temperature 24°C, and a 12:12h 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.

## 2.2 | Population dynamics under osmotic stresses

To quantify salinity tolerance, we assessed the demographic responses of *Dunaliella salina* to three levels of salinity. To ensure similar physiological states and densities among all populations at the beginning of the population dynamics assays, we first performed an acclimation step during 10 days, by diluting all populations at 1:125 in fresh medium at intermediate salinity ([NaCl] = 2.4 M; note that this corresponds in absolute to a high salinity of c. 140g/L, but is considered as intermediate for our model species). We then inoculated c. 2 × 10⁴ cells/ml of each populations into low (0.8 M), intermediate (2.4 M) or high (4.0 M) salinity, and tracked population density for the next 10 days under our standard growth conditions. We also followed the population dynamics of six randomly chosen populations of CCAP 19/15 strain starting at a lower density of 5 × 10³ cells/ml, to assess the effect of initial density on population growth rate in hyperosmotic condition ([NaCl] = 4.0 M).

To measure population growth rates under the different salt concentrations, we assessed population densities by passing a sub-sample of 150µl of each populations through a Guava EasyCyte HT flow cytometer (Luminex Corporation), following the protocol described in Leung et al. (2020). Discrimination between alive and dead algae was possible thanks to chlorophyll auto-fluorescence detected through a cytogram of emissions at Red-B (695/50nm) and Yellow-B (583/26nm) band pass filters (Papageorgiou, 2004), and the particle size was assessed through the forward scatter (FSC) and side scatter (SSC) parameters (Adan et al., 2017). Note that, even though osmotic stress causes immediate changes in cell volume in *D. salina*, live cell size could still be discriminated from other particles unambiguously (Leung et al., 2020), and population density could thus be assessed in all conditions as the number of cells per volume of assayed medium. To estimate population dynamics, cell counts were performed for live algae at 11 time points: end of the acclimation step, 4 h after the transfer to fresh media, and once per day for the following nine days.

## 2.3 | 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 two strains (CCAP 19/12 and CCAP 19/15), in two environmental conditions (hypo- and hyperosmotic). At the end of an acclimation step as described above, we transferred two biological replicates per strain to low ([NaCl] = 0.8 M) and high ([NaCl] = 4.0 M) salinities, in a greater volume (250ml) than for the demographic assays, and at a density of c. 1×10⁵ cells/ml so as to ensure enough material for high-throughput sequencing. After 24h following the salinity changes, the microalgae cells were harvested by centrifugation at 5000rpm for 15 min at room temperature, and cell pellets were stored at -80°C until nucleic acid extraction. As it has been shown that physiological regulations involving changes in gene expression usually start within 12-24h for *D. salina* exposed to salinity change (Chen & Jiang, 2009; Fang et al., 2017), we have chosen 24h as the time point for nucleic acid extraction.

Total RNA extraction and purification of the eight samples (2 strains ×2 salinities ×2 replicates) 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 et al. (1989). 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×150bp, Illumina HiSeq) were performed by Genewiz. We then performed all the bioinformatic preprocessing analyses with publicly available software implemented in the European UseGalaxy server (Afgan et al., 2018).

### 2.3.1 | Gene expression analyses

The RNA-seq raw reads were checked for quality using FastQC version 0.72 (Andrews, 2010) and subjected to adapter trimming and quality filtering using Trim Galore! version 0.4.3.1 (Krueger, 2015). 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 aligned the trimmed reads on the *D. salina* CCAP 19/18 reference nuclear (Dunsal1 v. 2, GenBank accession: GCA_002284615.2), chloroplastic (GenBank accession: GQ250046) and mitochondrial (GenBank accession: GQ250045) genomes, using HISAT2 version 2.1.0 (Kim, Langmead, et al., 2015) with default parameters for PE reads and spliced alignment option. We used Stringtie version 2.1.1 (Pertea et al., 2015) to predict transcript structures of each library based on the aligned reads, and performed de novo transcriptome assembly using the Stringtie merge tool, thus generating a unified and nonredundant set of transcripts across the different RNA-seq samples. We finally quantified the number of reads per transcript with FeatureCounts version 2.0.1 (Liao et al., 2014) using the alignment files from HISAT2 and the transcript annotation file from Stringtie.
## 2.3.2 | Genetic variant calling and genetic diversity analysis

We assessed the genetic differences between strains based on the sequences from the transcriptomic data. We first calculated the genomic range from which variant calling was performed as the number of bases in all exons, after merging any overlapping exons from different transcripts. BAM files from previous read alignment analyses with HISAT2 were submitted to freebayes (Garrison & Marth, 2012) for variant identification, with the following parameters: joint variant calling of the eight samples simultaneously, minimum alignment quality of three, ploidy set to one, and samples assumed to result from pooled sequencing. After left-alignment and normalization of indels using bcftools (Li et al., 2009), we filtered the obtained freebayes multisample VCF file for strand bias (SPR and SAP > 20), placement bias (EPP > 20), variant quality (QUAL > 30), and depth of coverage (DP > 20). To assess genetic diversity within and differentiation between strains, we calculated the gene diversity for each library ($H_E$), mean gene diversity for each strain ($H_s$), total gene diversity ($H_T$), and genetic differentiations ($G_{ST}$) among samples within each strain, and also between the two strains for all variants. We then represented the Euclidean genetic distances among samples with a dendrogram computed with ape R package (Paradis & Schliep, 2019).

Finally, for each locus that diverged between the strains (i.e., with alleles nearly fixed within strain but different between strains, $H_S \le 0.1$ and $G_{ST} \ge 0.9$), we assessed the synonymous versus nonsynonymous status of the substitution relative to the reference genome by using SnpEff version 4.3 (Cingolani et al., 2012). Sites where one of the focal strains had the reference allele while the other strain exhibited a synonymous mutation relative to the reference genome represented synonymous substitutions between our focal strains. As we detected only a single locus displaying a nonsynonymous mutation in both strains relative to the references, we did not consider it in the synonymous versus nonsynonymous mutation comparison between the strains.

## 2.3.3 | Methylation calling

The WGB-seq raw reads were also 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 and 5 bp were trimmed at the 5' and 3' extremity, respectively, to remove the tail added during library preparation, 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 (Krueger & Andrews, 2011). Only uniquely mapping reads were retained, 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 10x, and maximum tolerated variant fraction of 0.95). Cytosine methylation levels were determined for each CpG, CHG and CHH context. The high bisulfite conversion rate (>99%) was assessed by Genewiz by spiking in unmethylated lambda DNA in three randomly chosen libraries, and it was also confirmed by our analyses by estimating the number of methylated cytosine calls in the organelle genomes.

## 2.4 | Statistical analysis

### 2.4.1 | Population dynamics

To investigate how the population dynamics of our strains varied in response to salinity, we performed general linear models (GLMs) using cells count data as responses variables, with a negative binomial distribution and logarithm link function. We performed GLMs on population size using strain, salinity (treated as categorical variables), day (as continuous variable), and their interactions as fixed effect. In such models with logarithmic link functions, any effect of time (here day) translates into a rate of exponential growth (or decline if negative) per unit time (here per day), and interactions of time with other factors estimate effects on population growth, which are our main interest. We performed such models on initial growth rate (Day 0 to Day 1: GLM no. 1) and growth rate in the exponential phase (Day 1 to Day 4: GLM no. 2). We also investigated whether difference in growth rates could be explained by a release of density-dependent competition. To do so, we tested the effect of initial density on population growth rate of one of the strains in hyperosmotic condition, using populations starting at 20,000 or 5000 cells/ml (Day 1 to Day 4 at [NaCl] = 4.0 M: GLM no. 3 in Table S2). Statistical analyses were conducted using the statistical environment R version 4.0.3 (R Core Team, 2020) with the MASS package (Venables & Ripley, 2002) for the GLMs.

### 2.4.2 | Differential gene expression and DNA methylation analyses

The differential expression analyses were performed with the Bioconductor's package DESeq2 version 1.30.1 (Love et al., 2014). We first normalized the count matrix using DESeq2 regularized log transformed (rlog). We applied a redundancy analyses (RDA: Borcard et al., 1992), computed with the function rda() available in the vegan R package (Oksanen et al., 2020), to quantify the proportions of the total gene expression variation that are significantly explained by strain, salinity or the strain × salinity interaction, with population identity as covariates to account for paired samples between salinities. We also identified differentially expressed transcripts among the same three factors, by building a general linear model as implemented in DESeq2 and using Wald significance tests (Love et al., 2014). Transcripts with FDR < 0.05 (*-values after Benjamini-Hochberg [BH] adjustment) and $|log_{2}FC| > 1$ were
considered as differentially expressed. Significance of the strain × salinity interaction term was assessed using likelihood ratio tests (LRT) comparing models with and without the interaction term (Love et al., 2014).

Similarly, we performed a RDA to quantify the proportions of the total variation in DNA methylation level that are significantly explained by strain, salinity or strain × salinity interaction. We then used Bioconductor's methylKit package (Akalin et al., 2012) to identify differentially methylated regions (DMRs), corresponding to non-overlapping 100 bp windows between strains, or between salinities for each strain. The significance of calculated differences was determined using Fisher's exact tests. We used the Benjamini-Hochberg (BH) adjustment of *p*-values (FDR < 0.05) and methylation difference cutoffs of 20%.

For both RNA-seq and WGB-seq data, we performed principal component analyses (PCA) to represent the total variation among samples along its major axes. We also illustrated the expression level of DE transcripts and methylation levels of identified DMRs of all samples with heat-maps.

### 2.4.3 | Correlation between DNA methylation and gene expression

Since gene expression is a crucial step in the mapping from genotypes to phenotypes, and is thought to be a key mechanism underlying phenotypic plasticity, we wished to quantify to what extent gene expression levels can be predicted by key covariates in our data set. We first evaluated the different sources of variation in gene expression by achieving a global partitioning of variation, assessing the influence of (i) the genotype, (ii) the environment, and (iii) epigenetic variation on the total gene expression. We performed a RDA using the rlog transformed transcript count table as the response variable, and strain identity, salinity and epigenetic variation as explanatory variables. To account for paired samples among salinities, partial RDA was performed by removing the effect of population identity prior to assessing the effect of the environment or epigenetic variation on gene expression variation. Prior to RDA, we performed a PCA on DNA methylation levels and used the principal component (PCs) factors explaining at least 10% of the variation as a multivariate summary of epigenetic variation. We quantified the contributions to the total gene expression variation using the adjusted *R*², and tested the significance of each *R*² by ANOVA-like permutation tests using 999 randomization of the data (Legendre & Legendre, 1998).

To further quantify the role of DNA methylation in gene expression, we correlated differences in local DNA methylation to fold-changes in expression of the corresponding transcript. In the absence of well-annotated *D. salina* genome, we associated each cytosine to a given transcript based on its distance to the nearest transcription start site (TSS). Using the *genomation* package (Akalin et al., 2015), we first calculated TSS coordinates using the

gene structure file from the de novo transcriptome assembly of the RNA-seq preprocessing analyses. We then obtained the distance to nearest TSS and associated transcript identity for each cytosine. All cytosines associated to the same transcript were merged into a common gene-associated methylation region. We identified methylation differences between strains, and between salinities for each strain, using the same criterion for these gene-associated methylation regions as for DMRs from sliding windows above (i.e., FDR < 0.05 and |ΔmCG| > 20%). We finally compared the mean expression fold change associated to hypo- (<-20%) or hyper- (>20%) methylated DMRs for a given comparison. We restricted our correlation analysis to genomic regions that were both detected as significantly methylated between strains or salinities, and associated to transcripts that also showed significant differential expression between strains or salinities.

## 3 | RESULTS

### 3.1 | Genetic differences between strains

RNA-sequencing generated 1.92×10⁸ 150 bp paired-end raw reads from eight samples (Table S1). The genomic range from which variant calling were processed represented a total of 4.043×10⁷ nucleotides. Across these loci, we identified a total of 6201 (0.015%) variants among all samples. Among these polymorphic loci, 5500 loci displayed synonymous substitutions between the strains. This represents a moderate synonymous divergence of c. 10⁻⁴ per base pair between the two strains, consistent with previous observations from ITS sequences that led to placing these strains closeby in the *Dundaliella* phylogeny (Assunção et al., 2012; Emami et al., 2015). Where this has been investigated (in animals), such a level of synonymous divergence is consistent with within- rather than between-species variation (Roux et al., 2016), even though this notion becomes less clear for highly clonal micro-organisms.

Despite this low divergence, genetic variation was highly structured between the strains. Across the 6201 variants, 6170 (99.50%) displayed only two alleles across all the populations, and we measured a very low mean genetic diversity for these variants within each sample of a given strain (*H*<sub>S</sub> = 0.03, sd = 0.097 and 0.039, sd = 0.102 for CCAP 19/12 and 19/15, respectively, Figures S1-a,b), but a very high total genetic diversity across all samples (*H*<sub>T</sub> = 0.439, sd = 0.313, Figures S1-c), indicating highly structured genetic variation across strains (Figure 1). Specifically, 5589 (89.97%) variants displayed a fixed allele within a given strain, but distinct alleles between the strains (*H*<sub>S</sub> ≤0.1 and *G*<sub>ST</sub> ≥0.9). While we confirmed the genetic differences between the two strains, we also observed that samples of the same strain are quite genetically similar (Figure 1 and Figures S1-d,e). As a consequence, in subsequent analyses we used strain identity as a qualitative explanatory factor for genetic effects, instead of quantitative values based on genetic distances among samples.
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** Genetic differences between the strains. Genetic distance refers to proportion of genetic differences among samples, across all 6201 variants detected in the transcriptome. Low and High refer to salinity and numbers to samples</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**FIGURE 2** Population dynamics of two *Duniella salina* strains under different osmotic stresses. (a) Mean population growth curve in three different salinities. For each strain (CCAP 19/12 and CCAP 19/15 in red and black, respectively), mean cell density and standard error were calculated for the experiments starting with an initial density of 20,000 cells × ml⁻¹ (circle) for both strains (10 distinct populations for each), or 5000 cells × ml⁻¹ (cross) for CCAP 19/15 only (6 populations). Densities of individual populations also appear as smaller and lighter symbols. (b) Relationship between the initial and exponential growth rates following different osmotic stresses. Rug plots illustrate the distribution of the initial (days 0 to 1) and exponential (days 1 to 4) growth rates on their respective axes</figcaption></figure>
## 3.2 | Strains had markedly different demographic responses to osmotic stress

The two strains displayed very different population dynamics in response to osmotic stress. On the first day following transfer to a new salinity, both strains showed significantly positive growth in [NaCl] = 0.8 M and 2.4 M, but significant decline at [NaCl] = 4.0 M (Figure 2, GLMs no. 1_19/12 and no. 1_19/15 in Table S2). However, this decline in hyperosmotic conditions was much more pronounced for CCAP 19/12 compared to CCAP 19/15, as evidenced by the highly significant positive Day:Salinity_4.0 M:Strain_19/15 term (Table 1, GLM no. 1). While CCAP 19/12 declined by almost 80% (exp[-1.537] = 0.215 in GLM no. 1_19/12 in Table S2), CCAP 19/15 declined only by about 35% (exp[-0.435] = 0.647 in GLM no. 1_19/15 in Table S2). Strikingly, both strains then recovered from this initial decline at high salinity, and started growing again after Day 1, reaching a phase of exponential growth until c. Day 4 (Figure 2). But in this phase also, the dynamics markedly differed between the strains. Notably, CCAP 19/15 grew significantly slower than CCAP 19/12 at [NaCl] = 4.0 M (highly significant negative Day:Salinity_4.0 M:Strain_19/15 term in Table 1, GLM no. 2). In fact, the exponential growth rate of CCAP 19/12 did not significantly differ among salinities (GLM no. 2_19/12, Table S2), while that of CCAP

19/15 was significantly lower in [NaCl] = 2.4 M and 4.0 M as compared to 0.8 M (GLM no. 2_19/15, Table S2). These dynamics led to a negative relationship between initial growth rate and later exponential growth rate at high salinity, with the two strains occupying different regions along this relationship (rapid decline and growth for CCAP 19/12, slow decline and growth for CCAP 19/15), while no such pattern was observed at lower salinities (Figure 2b). We further showed that the faster growth rate of CCAP 19/12 during the exponential phase was not explained by released density-dependent competition as a result of its drastic initial decline (from GLM no. 2_19/12 and GLM no. 2_19/15, and GLM no. 3, Table S2), and that the dynamics of decline and rebound was not compatible with evolutionary rescue (Bell & Gonzalez, 2009; Gomulkiewicz & Holt, 1995), as it also occurred in isogenic populations (Figure S2 and Table S3).

## 3.3 | Gene expression responses to osmotic stresses

Structural gene annotation resulted in 31,926 putative transcripts present across the eight samples. PCA on the total gene expression variation revealed that samples first clustered according to the strain of origin (first PC axis, 61% of variance, Figure 3a), and then according

<table><caption>TABLE 1 Strain and salinity effects on population growth</caption><thead><tr><th></th><th>Estimate</th><th>SE</th><th>z value</th><th>Pr(&gt;|z|)</th><th></th></tr></thead><tbody><tr><td colspan="6"><em>GLM #1: Initial population response (days 0 to 1), starting from same initial density</em></td></tr><tr><td>Intercept (Strain_19/12 at Salinity 0.8 M and Day 0)</td><td>9.925</td><td>0.022</td><td>458.649</td><td>&lt;0.001</td><td>***</td></tr><tr><td>Day</td><td>0.197</td><td>0.072</td><td>2.723</td><td>0.006</td><td>**</td></tr><tr><td>Day:Salinity_2.4 M</td><td>-0.024</td><td>0.096</td><td>-0.254</td><td>0.800</td><td></td></tr><tr><td>Day:Salinity_4.0 M</td><td>-1.797</td><td>0.096</td><td>-18.673</td><td>&lt;0.001</td><td>***</td></tr><tr><td>Day:Strain_19/15</td><td>0.096</td><td>0.096</td><td>0.994</td><td>0.320</td><td></td></tr><tr><td>Day:Salinity_2.4 M:Strain_19/15</td><td>-0.048</td><td>0.136</td><td>-0.349</td><td>0.727</td><td></td></tr><tr><td>Day:Salinity_4.0 M:Strain_19/15</td><td>1.135</td><td>0.136</td><td>8.340</td><td>&lt;0.001</td><td>***</td></tr><tr><td colspan="6"><em>GLM #2: Population growth in exponential phase (days 1 to 4)</em></td></tr><tr><td>Intercept (Strain_19/12 at Salinity 0.8 M and Day 1)</td><td>10.317</td><td>0.071</td><td>145.898</td><td>&lt;0.001</td><td>***</td></tr><tr><td>Day</td><td>0.541</td><td>0.038</td><td>14.305</td><td>&lt;0.001</td><td>***</td></tr><tr><td>Salinity_2.4 M</td><td>-0.134</td><td>0.100</td><td>-1.337</td><td>0.181</td><td></td></tr><tr><td>Salinity_4.0 M</td><td>-1.904</td><td>0.100</td><td>-19.032</td><td>&lt;0.001</td><td>***</td></tr><tr><td>Strain_19/15</td><td>0.056</td><td>0.100</td><td>0.565</td><td>0.572</td><td></td></tr><tr><td>Day:Salinity_2.4 M</td><td>0.018</td><td>0.053</td><td>0.340</td><td>0.734</td><td></td></tr><tr><td>Day:Salinity_4.0 M</td><td>0.095</td><td>0.053</td><td>1.784</td><td>0.075</td><td></td></tr><tr><td>Day:Strain_19/15</td><td>0.056</td><td>0.053</td><td>1.043</td><td>0.297</td><td></td></tr><tr><td>Salinity_2.4 M:Strain_19/15</td><td>-0.021</td><td>0.141</td><td>-0.151</td><td>0.880</td><td></td></tr><tr><td>Salinity_4.0 M:Strain_19/15</td><td>1.095</td><td>0.141</td><td>7.741</td><td>&lt;0.001</td><td>***</td></tr><tr><td>Day:Salinity_2.4 M:Strain_19/15</td><td>-0.118</td><td>0.076</td><td>-1.566</td><td>0.117</td><td></td></tr><tr><td>Day:Salinity_4.0 M:Strain_19/15</td><td>-0.411</td><td>0.076</td><td>-5.436</td><td>&lt;0.001</td><td>***</td></tr></tbody></table>

Notes: General linear models (GLMs) with a negative binomial distribution were performed on cells count data, where the interaction of time (Day) with salinity or strain estimates effects of the latter on population exponential growth (or decline) rates. ** *p* < .01 and *** *p* < .001.
<figure><img src="image_5.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Strain</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) Salinity</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) Strain:Salinity interaction</figcaption></figure>

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

FIGURE 3 Gene expression response to osmotic stresses of two *D. salina* strains. (a) Principal component analysis (PCA) of gene expression level for two *D. salina* strains (colour) facing hypo- and hyperosmotic stress (shape). (b) Venn diagram showing the numbers of transcript that are significantly differentially expressed (DE; FDR<0.05 and |Log~2~FC|>1) between strains (light grey), salinities for CCAP 19/12 (red), or for CCAP 19/15 (dark grey). (c-e) Heat-maps of RNA-Seq transcriptome analyses for significant DE transcripts between strains (c; *n* = 3122), exclusively between salinities and common to the two strains (d; *n* = 284), and for strain × salinity interaction, identified by performing likelihood-ratio test (LRT, FDR<0.05) as implemented in DESeq2 (e; *n* = 2199). Each row and column represent a transcript and a biological replicate, respectively. Relative expression intensities among replicates vary from blue (under-expressed) to red (overexpressed), as shown on the right-hand side of the heat-maps. Dendrograms on the top resulted from a hierarchical clustering analysis using the Euclidean distance of the relative transcript expression level among replicates

term of population dynamics and gene expression level. We observed strain-specific growth rates within the first 24 h following an osmotic stress, and later during exponential growth. Specifically, strain CCAP 19/12 displayed a greater initial population decline as compared to CCAP 19/15 under hyperosmotic stress. This is possibly caused by a variable intensity of programmed cell death, a phenomenon that is widely present in unicellular microalgae, and can be triggered by a variety of environmental stresses, including an increase in salinity (Bidle, 2016; Durand & Ramsey, 2019; Zuppini et al., 2010). Interestingly, the exponential growth rate following initial decline was also strain-specific, rather than only resulting from variable competition levels following different initial declines among strains. Indeed at high salinity, CCAP 19/12 experienced faster growth than CCAP 19/15, even when the latter strain started at low density.

At the cellular level, the response of *D. salina* to osmotic stress is characterized by immediate changes in cell volume and intracellular ions, followed by slower changes in gene expression, starting within

4 to 24 h under the osmotic stress (Chen & Jiang, 2009), and regulating notably glycerol metabolism (Zhao et al., 2013). In this study, we detected differentially expressed transcripts between salinities 24 h after the osmotic stress, confirming gene-expression plasticity in response to salinity in *D. salina*. As we measured a very low within-strain genetic variation, we are confident that observed phenotypic differences among salinities, both in term of population dynamics and gene expression, resulted from phenotypic plasticity, rather than from selection in a heterogeneous cell population. This conclusion was also confirmed by the observation of the same population dynamics in isogenic populations.

We showed that responses to osmotic stress partly involved the regulation of strain-specific genes, which thus probably underlie the observed strain-specific population dynamics. We were able to successfully assign only c. 30% of the total transcripts to at least one gene ontology term (Figure S3). This transcript annotation analysis revealed that differentially expressed genes are mostly related to chloroplast structure and activities (Figure S3), but located in the
<figure><img src="image_10.png" /><figcaption>FIGURE 4 Methylation level and pattern in *Dunaliella salina*. (a) Distribution of the methylation level according to cytosine contexts. Histogram (grey) and density (red line) plots were calculated for CpG (n = 46,800), CHG (n = 69,039) and CHH (n = 227,531) context. (b) Distribution of methylation level for the nuclear (n = 32,409, 58,391 and 173,653 for CpG, CHG and CHH, respectively), mitochondria (n = 963, 929 and 4381 for CpG, CHG and CHH, respectively) and chloroplast (n = 13,428, 9719 and 49,497 for CpG, CHG and CHH, respectively) genomes, for each cytosine context</figcaption></figure>

nuclear genome. Interestingly, this lends support to the hypothesis above that the rapid decline followed by rebound that we observed under hyperosmotic stress may involve chloroplast-mediated programmed cell death (PCD), as studies have underlined the potential roles of cytochrome f, responsible for oxygenic photosynthesis, and thylakoid membrane complexes, in PCD for both plants and unicellular photosynthetic organisms (Ambastha al., 2015; Bidle, 2016; Murik et al., 2014; Thamatrakoln et al., 2013; Zuppini et al., 2009).

Our aim was to better understand how a given phenotype is connected to different levels of the genotype-phenotype (GP) map, by empirically relating gene expression to the expression of a higher phenotype. However, explanatory power in the GP map is necessarily limited by the phenotypic level with the lowest dimensionality. Here, we only analysed a few higher-level phenotypes chosen for their ecological meaningfulness. Achieving satisfying resolution in the GP map would probably require characterizing the phenotype more finely, as well as comparing additional genotypes and using more environmental contrasts (Chevin et al., 2021). For example, quantifying different chloroplast-related phenotypes (e.g., reactive oxygen species production, photosystem activity,

etc.) over a gradient of salinity would help break down correlations in the plastic responses, as well as elucidate the sequence of molecular events associated with chloroplasts influencing the rate of cellular death.

## 4.2 | DNA methylation patterns and gene silencing in *Dunaliella salina*

The levels and patterns of DNA methylation are known to vary drastically among organisms, including within microalgae (Feng et al., 2010; Zemach et al., 2010), such that characterizing these epigenetic traits and their potential roles in gene regulation remains an important goal, especially for organisms with special ecological or biological interest. Here, we present the first study of whole-genome DNA methylation in *D. salina*, a model organism for salinity tolerance, and a biotechnologically important species for carotene production (Ben-Amotz et al., 2009). We found that cytosine methylation is nearly absent in both the chloroplast and the mitochondrial genomes, and depends on the genomic context
<figure><img src="image_11.png" /><figcaption>FIGURE 5 DNA methylation response to osmotic stresses of two *D. salina* strains. (a) Principal component analysis (PCA) of DNA methylation level for two *D. salina* strains (colour) facing hypo- and hyperosmotic stress (shape). (b) Venn diagram showing the numbers of differentially methylated regions (DMR; q-value <0.05 and |diff-Methylation| >20%) between strains (light grey), salinities for CCAP 19/12 (red), or for CCAP 19/15 (dark grey). (c–e) Heat-maps of WGB-Seq analysis for DMRs between strains (c; *n* = 911), or between salinities for CCAP 19/12 (d; *n* = 40) and CCAP 19/15 (e; *n* = 54). Each row and column represent a DMRs and a biological replicate, respectively. DNA methylation levels vary from blue (unmethylated) to red (methylated), as shown on the right-hand side of the heat-maps. Dendrograms on the top resulted from a hierarchical clustering analysis using the Euclidean distance of DNA methylation level among replicates</figcaption></figure>

in the nuclear genome. More specifically, we observed the high-
est methylation levels of nuclear cytosines in the CpG context, but
with a highly bimodal distribution, including many low methylated
sites. We detected a lower degree of methylation (c. 10%) in the
CHG context, and quasi-absence of methylation in the CHH con-
text. Interestingly, this general pattern of methylation variation ac-
cording to C contexts is very similar to that observed in *Arabidopsis
thaliana* (Cokus et al., 2008; Feng et al., 2010; Zhang et al., 2020),
but quite different from those in more close related green algae
such as *Chlamydomonas reinhardtii* or *Volvox carteri*, which display
very low methylation levels, or *Chlorella variabilis*, where genes
are universally methylated (Feng et al., 2010; Zemach et al., 2010).
This confirms that the type and extent of DNA methylation varies
at relatively low phylogenetic scales (Alonso et al., 2019; Bewick
et al., 2017).

Beyond levels and patterns, the roles of DNA cytosine methyl-
ation in microalgae remains poorly understood, and could be very
variable. For example, in *C. variabilis* methylation level in the CpG

context in promoters was shown to be inversely correlated to gene
expression, suggesting that promoter-proximal methylation re-
presses transcription (Zemach et al., 2010). At the opposite, only a
weak negative correlation between promoter methylation and gene
transcription was observed in *V. carteri*, while CpG methylation is
enriched in transposons (Zemach et al., 2010). In *C. reinhardtii*, cy-
tosine methylation, observed in the chloroplast genome during
gametogenesis, was suggested to be involved in the uniparental
inheritance of mating type (Nishiyama et al., 2002, 2004; Sager &
Grabowy, 1983; Umen & Goodenough, 2001). Here, we could asso-
ciate only few differentially methylated regions with differentially
expressed transcripts. Nonetheless, we showed a negative correla-
tion between CpG methylation level for TSS-proximal cytosine and
gene expression, suggesting that DNA methylation at the promot-
ers represses transcription. Methylation at CHG context remained
however unclear, but could be involved in the silencing of repetitive
sequences and transposon, as in land plants or *C. variabilis* (Saze &
Kakutani, 2011; Kim, Ma, et al., 2015).
<figure><img src="image_12.png" /></figure>

FIGURE 6 Correlation between gene expression and epigenetic variation. (a) Transcriptomic variation partitioning, showing the proportions of the total transcriptomic variation explained by the genotype (strains CCAP 19/12 vs CCAP 19/15), the environment (low vs high salinities, controlling for paired samples) and cytosine methylation on CpG context, based on the adjusted $R^2$ of RDA analyses. Bold number outside circles refer to total percentage of variation (adjusted $R^2$) explained by the genotype, the environment and CpG methylation (top, full model), or by only the genotype, the environment or CpG methylation (reduced models). Percentages inside circles indicate pure contributions to gene expression level resulting from partial RDAs that isolate the effect of each single explanatory variable, but taking into account others as covariates. Percentages within intersections indicate shared contribution across different variables to gene expression variation. ANOVA-like permutation tests were calculated on total and pure contributions. NS, nonsignificant; \*p < .05; \*\*p < .01. (b) Differential transcript expression against differential DNA methylation. Boxplot showing the level of expression differentiation (Log2 fold change) against the difference in methylation level (ΔmCG). Boxplot and mean comparison analysis were performed for only DMRs detected as significantly differentially methylated between strains or salinities, and associated to a given transcript that also showed significant differential expression between strains or salinities. Differences between strain (n = 27), between salinities for CCAP 19/12 (n = 1), between salinities for CCAP 19/15 (n = 2)

## 4.3 | Sources of phenotypic variation and plasticity

One of the aims of our study was to assess the role of the genotype and non-genetic (environmental and epigenetic) factors on gene expression, by contrasting two strains with divergent phenotypic responses to osmotic shocks. While our results are necessarily limited by the use of only two strains, they still allowed refining our understanding of the role of CpG methylation in gene expression regulation for *D. salina*. First, we detected a substantial contribution of the genotype to both the epigenetic and phenotypic variation. Different studies have indeed suggested a genetic control of DNA methylation variation (Dubin et al., 2015; Hagmann et al., 2015; Richards, 2006). We also observed a negative correlation between levels of DMRs and DE transcripts, at least for few loci, suggesting that DNA methylation at CpG context is an intermediate step between the genotype and the gene-expression phenotype.

It has been shown through experimental evolution that reducing the amount of epigenetic variation, in terms of DNA methylation and histone acetylation, can diminish the plastic responses and reduce the rate of adaptation in response to some environmental stresses (Kronholm et al., 2017). However, the absence of a significant marginal effect of salinity or strain × salinity interaction on epigenetic state in our study questioned the role of overall methylation profile in plasticity. In plants and vertebrates, cytosine methylation responsible for gene expression regulation is mostly organized in specific regions, such as CpG islands, frequently in the gene promoters, or

within coding sequences (Bird, 2002; Deaton & Bird, 2011; Gallusci et al., 2016; Jaenisch & Bird, 2003; Law & Jacobsen, 2010; Zilberman et al., 2007), while the functional relevance of single CG methylation remain ambiguous (Denkena et al., 2021). Here, we detected strain-specific epigenetic differences between salinities, but only when considering 100 bp non-overlapping regions, rather than overall methylation profiles. This result highlighted that response to environmental changes involve fine-tuning of the expression of specific genes, rather than whole-methylome reprogramming. Accordingly, to better ascertain the role of epigenetic processes in plasticity, an interesting direction for future work would be to investigate methylation profiles of genes belonging to the same module of coexpression, which are likely to underlie expression of adaptive phenotypes in response to specific environmental changes.

We found that DNA methylation levels vary little across salinity for both strains, and also covary little with gene expression in *cis*. This result could suggest that environmentally induced changes in methylation are not involved in regulating gene expression's rapid responses to new environments. Alternatively, DNA methylation may play an important role in variation of gene expression across salinity, but mediated by cis-regulation at a small number of genes, such as transcription factors that cause *trans*-regulation of expression of many other genes. This hypothesis is supported by the significant marginal environmental effect that we detected on gene expression levels, and is also consistent with observations in, for example, the easter oyster *Crassostrea*
virginica (Johnson et al., 2021; Sirovy et al., 2021). Another possible explanation for the small detectable contribution of DNA methylation to gene expression plasticity in our experiment may be that the responses to osmotic stress that we observed may be too rapid to be explained by DNA methylation. DNA methyltransferases establish DNA methylations during cell divisions (Law & Jacobsen, 2010), which occur roughly once per day in *D. salina* (Ben-Amotz et al., 2009), so assessing DNA methylation 24 h after an osmotic stress as we did here should in theory be sufficient to observe epigenetic responses to salinity. However, the cell cycle is likely to depend on salinity and the strain, as reflected by their demographic dynamics (Figure 2). Other mechanisms of gene expression regulation are widespread in microalgae, such as post-translational histone modifications, and small-RNA mediated pathways (Kim, Ma, et al., 2015), which could also be involved the gene expression plasticity we observed in *D. salina*.

Although we detected that gene expression variation is mostly explained by the genotype and CpG methylation, we could not disentangle their influences, due to their strong correlation. Nevertheless, the significant genotype × environment interaction on gene expression, as well as the strain-specific epigenetic responses to environmental conditions, highlight the evolutionary potential of plasticity at different levels of the genotype-phenotype map.

In conclusion, we confirmed that the molecular mechanisms contributing to phenotypic plasticity in the halotolerant microalga *Dunaliella salina* involve variation in DNA methylation and gene expression with the environment, but found relatively weak contribution of DNA methylation to gene expression. Importantly, the large contribution of the genotype to the observed variation at multiple levels, including to demographic traits that are direct components of absolute fitness (Figure 2), highlights the evolutionary potential of phenotypic plasticity at multiple molecular levels. An interesting avenue for future research will be to use experimental evolution under controlled environmental conditions to investigate the evolution of plasticity at different levels, from DNA methylation to gene expression to higher phenotypes.

## AUTHOR CONTRIBUTIONS

Christelle Leung conceived and designed the study, with input from Daphné Grulois and Luis-Miguel Chevin. Christelle Leung and Daphné Grulois performed the experiment and collected the data. Christelle Leung analysed the data, prepared figures and tables, and wrote the original draft, Luis-Miguel Chevin reviewed and edited the draft. All the authors reviewed and approved the final draft of the manuscript.

## ACKNOWLEDGEMENTS

This work was supported by Grant 678140-FluctEvol from the European Research Council and Grant ANR-21-ERCC-0001-01 InterAdapt from the Agence Nationale de la Recherche to LMC, and a Fonds de Recherche du Québec—Nature et Technologies (FRQNT)

fellowship to CL. We thank N. Planidin, Samuel Scheiner and two anonymous reviewers for useful feedback on a previous version of this manuscript.

## CONFLICT OF INTERESTS

The authors declare no conflict of interest.

## 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 ID PRJNA736997. The specific samples used in this study are under the BioSample accessions SAMN19677492, SAMN19677495, SAMN19677496, SAMN19677497, SAMN19677507, SAMN19677510, SAMN19677511 and SAMN19677512. The de novo transcriptome assembly and growth rate data have been made available in the Dryad repository (doi:10.5061/dryad.3ffbg79m0).

## BENEFIT-SHARING

Benefits from this research accrue from the sharing of our data and results on public databases as described above.

## ORCID

Christelle Leung https://orcid.org/0000-0002-7242-3075

Daphné Grulois https://orcid.org/0000-0002-3054-1364

Luis-Miguel Chevin https://orcid.org/0000-0003-4188-4618

## REFERENCES

Adan, A., Alizada, G., Kiraz, Y., Baran, Y., & Nalbant, A. (2017). Flow cytometry: Basic principles and applications. *Critical Reviews in Biotechnology*, 37, 163–176.

Afgan, E., Baker, D., Batut, B., van den Beek, M., Bouvier, D., Cech, M., Chilton, J., Clements, D., Corao, N., Grüning, B. A., Guerler, A., Hillman-Jackson, J., Hiltemann, S., Jalili, V., Rasche, H., Soranzo, N., Goecks, J., Taylor, J., Nekrutenko, A., & Blankenberg, D. (2018). The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. *Nucleic Acids Research*, 46, W537–W544.

Akalin, A., Franke, V., Vlahoviček, K., Mason, C. E., & Schübeler, D. (2015). Genomation: A toolkit to summarize, annotate and visualize genomic intervals. *Bioinformatics*, 31, 1127–1129.

Akalin, A., Kormaksson, M., Li, S., Garrett-Bakelman, F. E., Figueroa, M. E., Melnick, A., & Mason, C. E. (2012). methylKit: A comprehensive R package for the analysis of genome-wide DNA methylation profiles. *Genome Biology*, 13, 1–9.

Alonso, C., Medrano, M., Perez, R., Canto, A., Parra-Tabla, V., & Herrera, C. M. (2019). Interspecific variation across angiosperms in global DNA methylation: Phylogeny, ecology and plant features in tropical and Mediterranean communities. *New Phytologist*, 224, 949–960.

Ambastha, V., Tripathy, B. C., & Tiwari, B. S. (2015). Programmed cell death in plants: A chloroplastic connection. *Plant Signalling & Behaviour*, 10, e989752.

Andrews, S. (2010). *FastQC: A quality control tool for high throughput sequence data*. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

Angers, B., Castonguay, E., & Massicotte, R. (2010). Environmentally induced phenotypes and DNA methylation: How to deal with unpredictable conditions until the next generation and after. *Molecular Ecology*, 19, 1283–1295.
Angers, B., Perez, M., Menicucci, T., & Leung, C. (2020). Sources of epi-
genetic variation and their applications in natural populations.
*Evolutionary Applications*, **13**, 1262–1278.

Assunção, P., Jaén-Molina, R., Caujapé-Castells, J., de la Jara, A., Carmona, L., Freijanes, K., et al. (2012). Molecular taxonomy of *Dunaliella* (Chlorophyceae), with a special focus on *D. salina*: ITS2 sequences revisited with an extensive geographical sampling. *Aquatic Biosystems*, **8**, 1–11.

Beldade, P., Mateus, A. R. A., & Keller, R. A. (2011). Evolution and mo-
lecular mechanisms of adaptive developmental plasticity. *Molecular
Ecology*, **20**, 1347–1363.

Bell, G., & Gonzalez, A. (2009). Evolutionary rescue can prevent extinction following environmental change. *Ecology Letters*, **12**, 942–948.

Ben-Amotz, A., Polle, J. E. W., & Subba, R. D. V. (2009). *The alga Dunaliella: Biodiversity, physiology, genomics and biotechnology*. Science Publishers.

Bewick, A. J., Vogel, K. J., Moore, A. J., & Schmitz, R. J. (2017). Evolution of DNA methylation across insects. *Molecular Biology and Evolution*, **34**, 654–665.

Bidle, K. D. (2016). Programmed cell death in unicellular phytoplankton. *Current Biology*, **26**, R594–R607.

Bird, A. (2002). DNA methylation patterns and epigenetic memory. *Genes and Development*, **16**, 6–21.

Bollati, V., & Baccarelli, A. (2010). Environmental epigenetics. *Heredity*, **105**, 105–112.

Borcard, D., Legendre, P., & Drapeau, P. (1992). Partialling out the spatial component of ecological variation. *Ecology*, **73**, 1045–1055.

Chen, H., & Jiang, J. G. (2009). Osmotic responses of *Dunaliella* to the changes of salinity. *Journal of Cellular Physiology*, **219**, 251–258.

Chen, X., Schönberger, B., Menz, J., & Ludewig, U. (2018). Plasticity of DNA methylation and gene expression under zinc deficiency in Arabidopsis roots. *Plant and Cell Physiology*, **59**, 1790–1802.

Chevin, L.-M., Leung, C., Le Rouzic, A., & Uller, T. (2021). Using phenotypic plasticity to understand the structure and evolution of the genotype-phenotype map. *Genetica*, **1-13**. https://doi.org/10.1007/s10709-021-00135-5

Cingolani, P., Platts, A., Wang le, L., Coon, M., Nguyen, T., Wang, L., Land, S. J., Lu, X., & Ruden, D. M. (2012). A program for annotating and predicting the effects of single nucleotide polymorphisms. *SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3*. Fly, 6, 80–92.

Cokus, S. J., Feng, S., Zhang, X., Chen, Z., Merriman, B., Haudenschild, C. D., Pradhan, S., Nelson, S. F., Pellegrini, M., & Jacobsen, S. E. (2008). Shotgun bisulphite sequencing of the *Arabidopsis* genome reveals DNA methylation patterning. *Nature*, **452**, 215–219.

R Core Team. (2020). *R: A language and environment for statistical computing*. R Foundation for Statistical Computing Vienna.

Cowan, A., Rose, P., & Horne, L. (1992). *Dunaliella salina*: A model system for studying the response of plant cells to stress. *Journal of Experimental Botany*, **43**, 1535–1547.

de Jong, G. (1990). Quantitative genetics of reaction norms. *Journal of Evolutionary Biology*, **3**, 447–468.

de Jong, G. (1995). Phenotypic plasticity as a product of selection in a variable environment. *American Naturalist*, **145**, 493–512.

de Jong, G. (1999). Unpredictable selection in a structured population leads to local genetic differentiation in evolved reaction norms. *Journal of Evolutionary Biology*, **12**, 839–851.

de Jong, G., & Gavrilets, S. (2000). Maintenance of genetic variation in phenotypic plasticity: The role of environmental variation. *Genetics Research*, **76**, 295–304.

de Mendoza, A., Lister, R., & Bogdanovic, O. (2020). Evolution of DNA methylome diversity in eukaryotes. *Journal of Molecular Biology*, **432**, 1687–1705.

Deaton, A. M., & Bird, A. (2011). CpG islands and the regulation of transcription. *Genes and Development*, **25**, 1010–1022.

Denkena, J., Johannes, F., & Colomé-Tatché, M. (2021). Region-level epi-
mutation rates in *Arabidopsis thaliana*. *Heredity*, **127**, 190–202.

Dubin, M. J., Zhang, P., Meng, D., Remigereau, M.-S., Osborne, E. J., Casale, F. P., et al. (2015). DNA methylation in *Arabidopsis* has a genetic basis and shows evidence of local adaptation. *Elife*, 4, e05255.

Durand, P. M., & Ramsey, G. (2019). The nature of programmed cell death. *Biological Theory*, **14**, 30–41.

Emami, K., Hack, E., Nelson, A., Brain, C. M., Lyne, F. M., Mesbahi, E., et al. (2015). Proteomic-based biotyping reveals hidden diversity within a microalgae culture collection: An example using *Dunaliella*. *Scientific Reports*, **5**, 1–15.

Fang, L., Qi, S., Xu, Z., Wang, W., He, J., Chen, X., et al. (2017). De novo transcriptomic profiling of *Dunaliella salina* reveals concordant flows of glycerol metabolic pathways upon reciprocal salinity changes. *Algal Research*, **23**, 135–149.

Feng, S., Cokus, S. J., Zhang, X., Chen, P.-Y., Bostick, M., Goll, M. G., et al. (2010). Conservation and divergence of methylation patterning in plants and animals. *Proceedings of the National Academy of Sciences*, **107**, 8689–8694.

Gallusci, P., Hodgman, C., Teyssier, E., & Seymour, G. B. (2016). DNA methylation and chromatin regulation during fleshy fruit development and ripening. *Frontiers in Plant Science*, **7**, 807.

Garrison, E. & Marth, G. (2012). Haplotype-based variant detection from short-read sequencing. *arXiv preprint arXiv:1207.3907 [q-bio.GN].

Gavrilets, S., & Scheiner, S. M. (1993). The genetics of phenotypic plasticity. *V. Evolution of reaction norm shape. Journal of Evolutionary Biology*, **6**, 31–48.

Ghalambor, C. K., Hoke, K. L., Ruell, E. W., Fischer, E. K., Reznick, D. N., & Hughes, K. A. (2015). Non-adaptive plasticity potentiates rapid adaptive evolution of gene expression in nature. *Nature*, **525**, 372–375.

Ghalambor, C. K., McKay, J. K., Carroll, S. P., & Reznick, D. N. (2007). Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments. *Functional Ecology*, **21**, 394–407.

Gibert, J.-M., Mouchel-Vielh, E., De Castro, S., & Peronnet, F. (2016). Phenotypic plasticity through transcriptional regulation of the evolutionary hotspot gene tan in *Drosophila melanogaster*. *PLOS Genetics*, **12**, e1006218.

Ginzburg, M. (1988). *Dunaliella*: A green alga adapted to salt. *Advances in Botanical Research*, **14**, 93–183.

Gomulkiewicz, R., & Holt, R. D. (1995). When does evolution by natural selection prevent extinction? *Evolution*, **49**, 201–207.

Hagmann, J., Becker, C., Müller, J., Stegle, O., Meyer, R. C., Wang, G., et al. (2015). Century-scale methylome stability in a recently diverged *Arabidopsis thaliana* lineage. *PLOS Genetics*, **11**, e1004920.

He, Q., Lin, Y., Tan, H., Zhou, Y., Wen, Y., Gan, J., et al. (2020). Transcriptomic profiles of *Dunaliella salina* in response to hyper saline stress. *BMC Genomics*, **21**, 1–17.

Jaenisch, R., & Bird, A. (2003). Epigenetic regulation of gene expression: How the genome integrates intrinsic and environmental signals. *Nature Genetics*, **33**, 245–254.

Johnson, K. M., Sirovy, K. A., & Kelly, M. W. (2021). Differential DNA methylation across environments has no effect on gene expression in the eastern oyster. *The Journal of Animal Ecology*, **1-13**.

Kim, D., Langmead, B., & Salzberg, S. L. (2015). HISAT: A fast spliced aligner with low memory requirements. *Nature Methods*, **12**, 357–360.

Kim, E.-J., Ma, X., & Cerutti, H. (2015). Gene silencing in microalgae: Mechanisms and biological roles. *Bioresource Technology*, **184**, 23–32.

Kirst, G. (1990). Salinity tolerance of eukaryotic marine algae. *The Annual Review of Plant Biology*, **41**, 21–53.

Kronholm, I., Bassett, A., Baulcombe, D., & Collins, S. (2017). Epigenetic and genetic contributions to adaptation in *Chlamydomonas*. *Molecular Biology and Evolution*, **34**, 2285–2306.

Krueger, F. (2015). *Trim Galore*. https://github.com/FelixKrueger/TrimG.alore
Krueger, F., & Andrews, S. R. (2011). Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. *Bioinformatics*, **27**, 1571-1572.

Law, J. A., & Jacobsen, S. E. (2010). Establishing, maintaining and modifying DNA methylation patterns in plants and animals. *Nature Reviews Genetics*, **11**, 204-220.

Legendre, P., & Legendre, L. (1998). *Numerical ecology* (2nd English edition ed.). Elsevier Science.

Leung, C., Breton, S., & Angers, B. (2016). Facing environmental predictability with different sources of epigenetic variation. *Ecology and Evolution*, **6**, 5234-5245.

Leung, C., Rescan, M., Grulois, D., & Chevin, L. M. (2020). Reduced phenotypic plasticity evolves in less predictable environments. *Ecology Letters*, **23**, 1664-1672.

Levins, R. (1963). Theory of fitness in a heterogeneous environment. II. Developmental flexibility and niche selection. *The American Naturalist*, **97**, 75-90.

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. *Bioinformatics*, **25**, 2078-2079.

Liao, Y., Smyth, G. K., & Shi, W. (2014). featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. *Bioinformatics*, **30**, 923-930.

Lou, S., Zhu, X., Zeng, Z., Wang, H., Jia, B., Li, H., et al. (2020). Identification of microRNAs response to high light and salinity that involved in beta-carotene accumulation in microalga *Dunaliella salina*. *Algal Research*, **48**, 101925.

Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biology*, **15**, 1-21.

Monteiro, A., Tong, X., Bear, A., Liew, S. F., Bhardwaj, S., Wasik, B. R., et al. (2015). Differential expression of ecdysone receptor leads to variation in phenotypic plasticity across serial homologues. *PLOS Genetics*, **11**, e1005529.

Murik, O., Elboher, A., & Kaplan, A. (2014). Dehydroascorbate: A possible surveillance molecule of oxidative stress and programmed cell death in the green alga *Chlamydomonas reinhardtii*. *New Phytologist*, **202**, 471-484.

Nishiyama, R., Ito, M., Yamaguchi, Y., Koizumi, N., & Sano, H. (2002). A chloroplast-resident DNA methyltransferase is responsible for hypermethylation of chloroplast genes in *Chlamydomonas* maternal gametes. *Proceedings of the National Academy of Sciences*, **99**, 5925-5930.

Nishiyama, R., Wada, Y., Mibu, M., Yamaguchi, Y., Shimogawara, K., & Sano, H. (2004). Role of a nonselective de novo DNA methyltransferase in maternal inheritance of chloroplast genes in the green alga, *Chlamydomonas reinhardtii*. *Genetics*, **168**, 809-816.

Oksanen, J., Blanchet, F.G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D. et al. (2020). *vegan: Community Ecology Package*. http://CRAN.R-project.org/package=vegan

Oren, A. (2005). A hundred years of *Dunaliella* research: 1905-2005. *Saline Systems*, **1**, 2.

Papageorgiou, G. C. (2004). Fluorescence of photosynthetic pigments in vitro and in vivo. In G. C. Papageorgiou & R. Govindjee (Eds.), *Chlorophyll a fluorescence: A signature of photosynthesis*. *Advances in photosynthesis and respiration* (pp. 43-63). Springer.

Paradis, E., & Schliep, K. (2019). ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. *Bioinformatics*, **35**, 526-528.

Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T.-C., Mendell, J. T., & Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. *Nature Biotechnology*, **33**, 290-295.

Pick, U. (2002). Adaptation of the halotolerant alga *Dunaliella* to high salinity. In A. Läuchli & U. Lüttge (Eds.), *Salinity: Environment-plants-molecules* (pp. 97-112). Springer.

Polle, J. E., Barry, K., Cushman, J., Schmutz, J., Tran, D., Hathwaik, L. T., et al. (2017). Draft nuclear genome sequence of the halophilic and beta-carotene-accumulating green alga *Dunaliella salina* strain CCAP19/18. *Genome Announcements*, **5**, e01105-e01117.

Polle, J. E., Calhoun, S., McKie-Krisberg, Z., Prochnik, S., Neofotis, P., Yim, W. C., et al. (2020). Genomic adaptations of the green alga *Dunaliella salina* to life under high salinity. *Algal Research*, **5**, 101990.

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

Richards, E. J. (2006). Inherited epigenetic variation—revisiting soft inheritance. *Nature Reviews Genetics*, **7**, 395-401.

Roux, C., Fraisse, C., Romiguier, J., Anciaux, Y., Galtier, N., & Bierne, N. (2016). Shedding light on the grey zone of speciation along a continuum of genomic divergence. *PLOS Biology*, **14**, e2000234.

Sager, R., & Grabowy, C. (1983). Differential methylation of chloroplast DNA regulates maternal inheritance in a methylated mutant of *Chlamydomonas*. *Proceedings of the National Academy of Sciences*, **80**, 3025-3029.

Sambrook, J., Fritsch, E. F., & Maniatis, T. (1989). *Molecular cloning: A laboratory manual* (2nd ed.). Cold Spring Harbour Laboratory Press.

Saze, H., & Kakutani, T. (2011). Differentiation of epigenetic modifications between transposons and genes. *The Current Opinion Journals*, **14**, 81-87.

Scheiner, S. M. (1993). Plasticity as a selectable trait: Reply to Via. *The American Naturalist*, **142**, 371-373.

Scheiner, S. M., & Holt, R. D. (2012). The genetics of phenotypic plasticity. X. Variation versus uncertainty. *Ecology and Evolution*, **2**, 751-767.

Schlichting, C. D., & Pigliucci, M. (1993). Control of phenotypic plasticity via regulatory genes. *The American Naturalist*, **142**, 366-370.

Sirovy, K. A., Johnson, K. M., Casas, S. M., La Peyre, J. F., & Kelly, M. W. (2021). Lack of genotype-by-environment interaction suggests limited potential for evolutionary changes in plasticity in the eastern oyster, *Crassostrea virginica*. *Molecular Ecology*, **30**, 5721-5734.

Thamatrakoln, K., Bailleul, B., Brown, C. M., Gorbunov, M. Y., Kustka, A. B., Frada, M., et al. (2013). Death-specific protein in a marine diatom regulates photosynthetic responses to iron and light availability. *Proceedings of the National Academy of Sciences*, **110**, 20,123-20,128.

Tufto, J. (2015). Genetic evolution, plasticity, and bet-hedging as adaptive responses to temporally autocorrelated fluctuating selection: A quantitative genetic model. *Evolution*, **69**, 2034-2049.

Umen, J. G., & Goodenough, U. W. (2001). Chloroplast DNA methylation and inheritance in *Chlamydomonas*. *Genes and Development*, **15**, 2585-2597.

Venables, W. N., & Ripley, B. D. (2002). *Modern Applied Statistics with S* (4th ed.). Springer-Verlag.

Via, S. (1993a). Adaptive phenotypic plasticity: Target or by-product of selection in a variable environment? *The American Naturalist*, **142**, 352-365.

Via, S. (1993b). Regulatory genes and reaction norms. *The American Naturalist*, **142**, 374-378.

Via, S., Gomulkiewicz, R., De Jong, G., Scheiner, S. M., Schlichting, C. D., & Van Tienderen, P. H. (1995). Adaptive phenotypic plasticity: Consensus and controversy. *Trends in Ecology and Evolution*, **10**, 212-217.

Via, S., & Lande, R. (1985). Genotype-environment interaction and the evolution of phenotypic plasticity. *Evolution*, **39**, 505-522.

Zemach, A., McDaniel, I. E., Silva, P., & Zilberman, D. (2010). Genome-wide evolutionary analysis of eukaryotic DNA methylation. *Science*, **328**, 916-919.

Zhang, H., Lang, Z., & Zhu, J.-K. (2018). Dynamics and function of DNA methylation in plants. *Nature reviews Molecular cell biology*, **19**, 489-506.
Zhang, Y., Liu, C., Cheng, H., Tian, S., Liu, Y., Wang, S., et al. (2020). DNA methylation and its effects on gene expression during primary to secondary growth in poplar stems. *BMC Genomics, 21*, 1–18.

Zhao, L.n., Gong, W.f., Chen, X.w., & Chen, D.f. (2013). Characterization of genes and enzymes in *Dunaliella salina* involved in glycerol me-tabolism in response to salt changes. *Phycological Research, 61*, 37–45.

Zilberman, D., Gehring, M., Tran, R. K., Ballinger, T., & Henikoff, S. (2007). Genome-wide analysis of *Arabidopsis thaliana* DNA methylation uncovers an interdependence between methylation and transcription. *Nature Genetics, 39*, 61–69.

Zuppini, A., Gerotto, C., & Baldan, B. (2010). Programmed cell death and adaptation: Two different types of abiotic stress response in a uni-cellular chlorophyte. *Plant and Cell Physiology, 51*, 884–895.

Zuppini, A., Gerotto, C., Moscatiello, R., Bergantino, E., & Baldan, B. (2009). *Chlorella saccharophila* cytochrome f and its involvement in the heat shock response. *Journal of Experimental Botany, 60*, 4189–4200.

## SUPPORTING INFORMATION

Additional supporting information may be found in the online version of the article at the publisher's website.

**How to cite this article:** Leung, C., Grulois, D., & Chevin, L-M (2022). Plasticity across levels: Relating epigenomic, transcriptomic, and phenotypic responses to osmotic stress in a halotolerant microalga. *Molecular Ecology, 31*, 4672–4687. https://doi.org/10.1111/mec.16542