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

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.

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 .
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, 1993aVia, , 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.

| 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 . Our standard growth conditions are 50 ml 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-500 ml), for a total 25 ml (including the inoculate), incubated at a constant temperature 24°C, and a 12:12 h light/dark cycle with a 200 μmol.m −2 .s −1 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.

| 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. 140 g/L, but is considered as intermediate for our model species). We then inoculated c. 2 × 10 4 cells/m l 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 3 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 subsample 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/50 nm) and Yellow-B (583/26 nm) 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 , 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.

| 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 wholegenome 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 (250 ml) than for the demographic assays, and at a density of c. 1 × 10 5 cells/ml so as to ensure enough material for high-throughput sequencing. After 24 h following the salinity changes, the microalgae cells were harvested by centrifugation at 5000 rpm 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-24 h for D. salina exposed to salinity change (Chen & Jiang, 2009;Fang et al., 2017), we have chosen 24 h 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 construc-
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.

| 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 Finally, for each locus that diverged between the strains (i.e., with alleles nearly fixed within strain but different between strains, H S ≤0.1 and G ST ≥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.

| 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 10×, and maximum tolerated variant fraction of 0.95). Cytosine methylation levels were determined for each CpG, CHG and CHH context. The high bisulphite 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.

| 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.

| 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 (p-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 nonoverlapping 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

| 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 2 , and tested the significance of each R 2 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.

| Genetic differences between strains
RNA-sequencing generated 1.92 × 10 8 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 7 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 −4 per base pair between the two strains, consistent with previous observations from ITS sequences that led to placing these strains closeby in the Dunaliella 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 betweenspecies 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 S = 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 T = 0.439, 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 S ≤0.1 and G ST ≥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.   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  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).

| 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 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.

TA B L E 1 Strain and salinity effects on population growth
to the salinity treatment (second PC axis, 17% of variance, Figure 3a).   (Table S1).

| DNA methylation variation is structured by genomic contexts
After the data filtering, we performed our methylation analyses on an average of 5.19 × 10 8 (s.d. 1.08 × 10 8 ) cytosines per samples, and observed that DNA methylation is not randomly distributed along the genome (Figure 4a). Cytosine in CpG context displayed the highest methylation level, as compared to CHG and CHH contexts. Furthermore, the nuclear genome had the most methylated cytosine, as compared to the mitochondria and chloroplast genomes Redundancy analyses revealed a significant effect of strain (R 2 = 58.69%; p = .005) on total DNA methylation at the CpG context. However, contrary to gene expression, we did not detect any significant marginal effect of salinity (R 2 = 6.47%; p = .602) on epigenetic state at the genomic level, as also observed on the PCA plot ( Figure 5a). On a finer scale, considering non-overlapping 100 bp windows instead of the overall methylation profile, we detected 932 DMRs between the two strains, but only 27 DMR between salinities when all strains were pooled together, and no significant strain × salinity interaction (R 2 = 5.83%; p = .566). Nevertheless, we did detect a few strain-specific DMRs between salinities for each strain (n = 40 and 54, for CCAP 19/12 and 19/15, respectively; Figure 5b-e).

| Gene expression depends on both the genotype and the environment
The genotype, the environment, and the CpG methylation level jointly significantly explained an important part of total transcriptomic variation (adjusted R 2 = 88.11%; p = .005; Figure 6a). More precisely, the strong correlation between the genotype and DNA methylation reported above (R 2 = 58.69%; p = .005) resulted in a great confounding effect (R 2 = 62.73%), whereby the relative contri-

| DISCUSS ION
We investigated how genetic variation relates to plastic responses at multiple biological levels, by comparing DNA methylation patterns, gene expression levels, and demographic phenotypes across salinities, in two strains of a halotolerant eukaryote microbe known to be predominantly under selection for salinity tolerance in its natural habitat (Ben-Amotz et al., 2009;Kirst, 1990;Oren, 2005). Our results shed light on the role of DNA methylation in the rapid response to osmotic shock in Dunaliella salina, and the importance of the genotype in such response.

| Variation of reaction norms between genotypes
We first showed that the two strains used in our study displayed distinct phenotypic response to environmental change, both in 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 withinstrain 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 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.

| 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 in the nuclear genome. More specifically, we observed the highest 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 context. Interestingly, this general pattern of methylation variation according 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 reinhardi 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 methylation 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 represses 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, cytosine methylation, observed in the chloroplast genome during gametogenesis, was suggested to be involved in the uniparental inheritance of mating type (Nishiyama et al., 2002(Nishiyama et al., , 2004Sager & Grabowy, 1983;Umen & Goodenough, 2001). Here, we could associate only few differentially methylated regions with differentially expressed transcripts. Nonetheless, we showed a negative correlation between CpG methylation level for TSS-proximal cytosine and gene expression, suggesting that DNA methylation at the promoters 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).

| 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 methylations 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 sa- (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)  (log2FC) virginica 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 genotypephenotype 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. 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.

CO N FLI C T O F I NTE R E S T S
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
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).

B EN EFIT-S H A R I N G
Benefits from this research accrue from the sharing of our data and results on public databases as described above.