Marie Rescan

Daphné Grulois

Enrique Ortega-Aboud

Luis-Miguel Chevin
email: luis-miguel.chevin@cefe.cnrs.fr

Phenotypic memory drives population growth and extinction risk in a noisy environment

come

The demography and extinction risk of wild populations is largely determined by responses to random fluctuations in their environment [START_REF] Ovaskainen | Stochastic models of population extinction[END_REF][START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Lewontin | On population growth in a randomly varying environment[END_REF][START_REF] Lawson | Environmental variation and population responses to global change[END_REF][START_REF] Drake | Population effects of increased climate variation[END_REF][START_REF] García-Carreras | Are Changes in the Mean or Variability of Climate Signals More Important for Long-Term Stochastic Growth Rate?[END_REF][START_REF] Tuljapurkar | Population dynamics in variable environment[END_REF] , the magnitude and predictability of which are currently altered by anthropogenic activities [START_REF] Boer | Changes in interannual variability and decadal potential predictability under global warming[END_REF][START_REF] Wigley | Anthropogenic influence on the autocorrelation structure of hemispheric-mean temperatures[END_REF] . In addition to their amplitude, stochastic (i.e., random) fluctuations are characterized by their temporal autocorrelation, which determines the extent to which one time step is a good predictor for the next. Temporal autocorrelation, sometimes described as the color of noise in reference to the spectral decomposition of random fluctuations, is a ubiquitous feature of natural systems [START_REF] Vasseur | The color of environmental noise[END_REF][START_REF] Sabo | Quantifying periodic, stochastic, and catastrophic environmental variation[END_REF] . Because environmental autocorrelation characterizes the pattern and ordering of environmental noise rather than its magnitude, it was originally thought to be irrelevant to expected population dynamics in a stochastic environment, which instead were predicted to only depend on the mean and variance of environmental fluctuations in population growth [START_REF] Lewontin | On population growth in a randomly varying environment[END_REF] . However, more recent theory has shown that environmental autocorrelation may have a dramatic impact on population growth and extinction risk in a stochastic environment [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF][START_REF] Ripa | Noise colour and the risk of population extinctions[END_REF][START_REF] Wichmann | Extinction risk in periodically fluctuating environments[END_REF][START_REF] Heino | Extinction risk under coloured environmental noise[END_REF][START_REF] Ruokolainen | Ecological and evolutionary dynamics under coloured environmental variation[END_REF] . The reason is that positive temporal autocorrelation facilitates phenotypic tracking of environmental fluctuations by phenotypic plasticity [START_REF] Reed | Phenotypic plasticity and population viability: the importance of environmental predictability[END_REF][START_REF] Ashander | Predicting evolutionary rescue via evolving plasticity in stochastic environments[END_REF] , evolutionary responses to natural selection [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF] , or both, thereby modifying the mean and variance of population growth rates. However theoretical predictions are scarcer for negatively autocorrelated fluctuations, where high and low conditions alternate rapidly but predictably. More importantly, we lack empirical validation of this theory in a form that combines the realism of a continuously varying environment (such as temperature, humidity, salinity) with the high replication required to properly investigate responses to stochasticity. High replication is all the more important as environmental autocorrelation is expected to influence the variance of population size among independent realizations of stochastic population dynamics, therefore increasing extinction risk [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF] .

An increasing number of studies have inferred the level of environmental stochasticity in population growth in the wild [START_REF] Saether | Pattern of variation in avian population growth rates[END_REF][START_REF] Bjørkvoll | Stochastic population dynamics and life-history variation in marine fish species[END_REF] and the demographic consequences of environmental autocorrelation [START_REF] Engen | Estimating the effect of temporally autocorrelated environments on the demography of density-independent age-structured populations[END_REF][START_REF] Paniw | Interactive life-history traits predict sensitivity of plants and animals to temporal autocorrelation[END_REF] .

However, attributing population size fluctuations to measured environmental factors is extremely challenging in the wild, where such factors may largely covary among each other, and the origin of their variation is unknown. Furthermore, observations from natural population generally consist of a single realization of a process with a limited number of time points (usually below 50), whereas assessing the range of likely outcomes for any stochastic process in the future requires many replicates, and/or long time spans. Experimental studies in the laboratory with short-lived organisms appears as a necessary bridge between theory and nature, because they allow testing theoretical predictions with ample replication and precision, potentially yielding general insights about population responses to environmental stochasticity. However, the few experimental studies about demographic consequences of environmental autocorrelation have had contrasting findings [START_REF] Pike | The effect of autocorrelation in environmental variability on the persistence of populations: An experimental test[END_REF][START_REF] Petchey | Environmental colour affects the dynamics of single-species populations[END_REF][START_REF] Laakso | Environmental noise and population dynamics of the ciliated protozoa Tetrahymena thermophila in aquatic microcosms[END_REF][START_REF] Duncan | Stochastic environmental fluctuations drive epidemiology in experimental host-parasite metapopulations[END_REF] . To achieve significant progress on this question, we need to overcome the specificities of each individual system by establishing robust and general mechanistic links between environments and population demographic processes.

Such a link can be provided by environmental tolerance curves, which measure how fitness or its life history components (survival and fecundity) respond to abiotic environmental variables such as temperature or salinity. Tolerance curves are increasingly popular tools connecting physiology and ecology in response to climate change [START_REF] Deutsch | Impacts of climate warming on terrestrial ectotherms across latitude[END_REF][START_REF] Palumbi | Mechanisms of reef coral resistance to future climate change[END_REF][START_REF] Calosi | Thermal tolerance, acclimatory capacity and vulnerability to global climate change[END_REF] , notably because they share many commonalities across taxa and environmental factors, thus allowing for generalization. They are typically hump-shaped, with an optimum at some intermediate environmental value that coincides to some extent with the historical environment of the species 31,[START_REF] Thomas | A global pattern of thermal adaptation in marine phytoplankton[END_REF] . This non-linear relationship between the environment and absolute fitness is predicted to strongly drive effects of environmental stochasticity on population dynamics and extinction risk [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF] , but these predictions still largely lack empirical validation. Furthermore, environmental tolerance curves can potentially be influenced by the history of environments experienced by an individual [START_REF] Calosi | Thermal tolerance, acclimatory capacity and vulnerability to global climate change[END_REF][START_REF] Hoffmann | Acclimation: Increasing survival at a cost[END_REF][START_REF] Nougué | The ontogeny of tolerance curves: habitat quality vs. acclimation in a stressful environment[END_REF][START_REF] Lande | Evolution of phenotypic plasticity and environmental tolerance of a labile quantitative character in a fluctuating environment[END_REF] (phenotypic plasticity) or its ancestors [START_REF] Proulx | How differing modes of transgenerational inheritance affect population viability in fluctuating environments[END_REF] (transgenerational plasticity), but it is unclear to what extent population dynamics will be affected by such environmental memory effects. To answer these questions, we carried out a highly replicated experiment with a eukaryotic microorganism, in order to investigate how environmental autocorrelation influences population dynamics and extinction risk in a stochastic environment.

Results

We exposed six ancestral lineages (hereafter genotypes) of the facultative sexual microalgae Dunaliella salina to randomly varying salinities during 37 transfers (~100 generations). Using a liquid handling robot, we transferred each line twice a week (every 3 or 4 days), diluting 15% of the population of origin into fresh medium with controlled salinity, which differed across lines. More precisely, our treatments consisted of 164 time series with the same mean (µ = 2.4M NaCl) and variance (standard deviation σ = 1M NaCl), but 4 autocorrelation levels, from negative ( = -0.5) to highly positive ( = 0.9). In doing so, we modulated the pattern of salinity change between two successive transfers, without modifying the overall magnitude of fluctuations (Figure 1 a. and b., upper panel). Importantly, 39 (out of 41) time series within each treatment were different, independent realizations of the same stochastic process, allowing investigation of the variance introduced by environmental stochasticity, which is known to be an important driver of extinction risk [START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF] . We tracked population densities and the proportion of extinct populations at each of 37 transfers (~100 generations, over 4 months) using optical density, fluorescence and flow cytometry. We then estimated the intrinsic rate of increase between each successive transfers, assuming continuous logistic growth following the 15% dilution (Figure 1c).

In addition to three single strains (genotypes A: CCAP 19/18, B: CCAP 19/18 and C: CCAP 19/15), we used the three possible mixes of two of these strains (genotype AB, AC and BC), to test for potential advantages of genetic diversity on population dynamics. One of the single genotypes (B) grew too slowly to compensate for dilution at each transfer and got extinct rapidly (Extended Data 1) even in constant treatments, so it was discarded from further analysis.

Environmental autocorrelation strongly affects population size. Environmental stochasticity

initially caused a rapid increase of variation in population size among our different lines of a same stochasticity treatment (as materialized by the dashed confidence interval in Figure 2.a.), as expected because population growth accumulates linearly on the log scale during a phase of exponential growth [START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF] . The stochastic variation among lines then settled to an approximately constant value, characteristic of the stationary phase (days > 45, Figure 2a), where fluctuating population sizes spanned several orders of magnitude, in contrast with the small variance among replicates in constant environments (Extended Data 2). This is consistent with theory showing that stochastic environments can cause large fluctuations in population sizes around their expected carrying capacity [START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF][START_REF] Engen | A Quantitative Genetic Model of r -and K -Selection in a Fluctuating Population[END_REF] .

The distribution of log-population size was asymmetric for all treatments and genotypes, with a negative skew causing an excess of very small populations undergoing high extinction risk (Figure 2.b. and e). Remarkably, this skewness is consistent with models of a moving optimum phenotype in a stochastic environment [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF] -equivalent to models of tolerance curves with an optimum environment [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF][START_REF] Lande | Evolution of phenotypic plasticity and environmental tolerance of a labile quantitative character in a fluctuating environment[END_REF][START_REF] Chevin | Adaptation, plasticity, and extinction in a changing environment: Towards a predictive theory[END_REF] -, but not necessarily with classic phenomenological models where noise is added to the population growth rate without an explicit tolerance curve [START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Lewontin | On population growth in a randomly varying environment[END_REF][START_REF] Tuljapurkar | Population dynamics in variable environments I. Long-run growth rates and extinction[END_REF] . This suggests that fluctuations in our experimental populations are driven to some extent by the interplay of a stochastic environment with a hump-shaped tolerance curve with an optimum. However, density dependence in population growth also impacts the distribution of population size, modifying the influence of the stochastic environment [START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Diserud | A general and dynamic species abundance model, embracing the lognormal and the gamma Models[END_REF] . To control for this influence of density dependence, we analyzed the distribution of the intrinsic growth rate r, i.e. the rate of exponential growth at low population density, extracted from a model of continuous logistic growth (Figure 1.c). We found strong support for an asymmetrical distribution of r (ΔAIC = -3713 between a normal and a reverse gamma distribution), with an estimated skewness between -1.7 and -0.4 (Figure 3.d).

Using the population sizes estimations obtained from the reverse gamma model, we found important differences between the dynamics of populations exposed to different autocorrelation treatments. For each of the genotype x autocorrelation treatment (from 438 to 833 data points), we computed the distributions of the population size during the stationary phase, and analyzed the influence of environmental autocorrelation on moments of these distributions. As environmental autocorrelation increased, populations underwent fewer episodes of negative growth, and smaller fluctuations in size (Figure 1). As a consequence, the mean (log) population size was larger under larger environmental autocorrelation (Figure 2c, slope = 0.71±0.11, p = 4.0.10 -5 ), contrary to classical population dynamics predictions where the expected intrinsic growth rate only depends on the mean and variance of the environment [START_REF] Lewontin | On population growth in a randomly varying environment[END_REF][START_REF] Lawson | Environmental variation and population responses to global change[END_REF][START_REF] García-Carreras | Are Changes in the Mean or Variability of Climate Signals More Important for Long-Term Stochastic Growth Rate?[END_REF] . The variance in log population size among replicate lines decreased with increasing autocorrelation (Figure 2d, slope = -3.0±0.71, p = 1.3.10 -3 ). This contrasts with the common theoretical prediction that higher autocorrelation generates long batches of good/bad environments in different lines, thereby increasing population size variance among lines [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF][START_REF] Ripa | Noise colour and the risk of population extinctions[END_REF][START_REF] Wichmann | Extinction risk in periodically fluctuating environments[END_REF][START_REF] Heino | Extinction risk under coloured environmental noise[END_REF] , but is consistent with a contribution of environmental autocorrelation to phenotypic tracking of the environment (via phenotypic plasticity or genetic evolution), reducing the magnitude of population fluctuations [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF] . We found a general trend for lower asymmetry (negative skewness closer to 0) with increasing temporal autocorrelation, but this effect highly depended on the strain (or mix) used and was not significant (Figure 3.d, slope = 0.19±0.25). Temporal autocorrelation of the stationary population size was highly variable within genotype x autocorrelation treatment (Figure 3.e). In particular, it was never found significantly different from 0, even in the highly autocorrelated treatment (ρ = 0.9), highlighting that temporal autocorrelation of the environmental factors translates only very weakly into autocorrelation of population sizes [START_REF] Petchey | Environmental colour affects the dynamics of single-species populations[END_REF][START_REF] Ferguson | An updated perspective on the role of environmental autocorrelation in animal populations[END_REF][START_REF] García-Carreras | An empirical link between the spectral colour of climate and the spectral colour of field populations in the context of climate change[END_REF][START_REF] Knape | Effects of weather and climate on the dynamics of animal population time series[END_REF] .

Environmental memory governs population dynamics.

To reach more mechanistic insights into drivers of stochastic population dynamics in our experiment, we used the known salinity and population measurements to estimate a hump-shaped tolerance curve relating the intrinsic rate of increase r(T) following transfer T to the current salinity S T , , 1 where the subscript 'univ' is for univariate. This was compared to a bivariate tolerance curve that also includes an effect of memory of the previous environment (prior to the last transfer), via the salinity

S T -1 ,
2 where e determines the tolerance breadth with respect to the current environment, d the breadth of the environmental memory effect, and f the interaction between past and current salinity.

We fitted parameters of the tolerance curve with (equation 2) or without memory (equation 1), under the same model of logistic growth as in the first section, assuming a Gaussian distribution of the population size around its expected value (given the value at the previous time step). We found strong support for the generalized bivariate tolerance curve (ΔAIC < -10 -3 and p-values from likelihood ratio test < 10 -9 for all genotypes, see Supplementary Table 1), revealing an important impact of environmental memory on population growth. Growth rates were highest for populations transferred from high to low salinity, while populations transferred from low to high salinities declined (Figure 3.a), notably because of high mortality that could be detected by flow cytometry (Extended Data 3a, b, c).

To quantify the extent to which tolerance curves and memory effects explain the variation in growth within and across autocorrelation treatments, we derived the moments of population growth rates predicted by combining these tolerance curves with the realized salinity time series, and compared them to their estimated counterparts that do not use information about salinity (Supplementary Table 1). The tolerance curve with memory correctly predicted increases in the mean (slope = 0.069±0.009, p = 1.5.10 -5 ) and skewness (slope = 0.58±0.09, p = 2.9.10 -5 ) of population growth rates with increasing environmental autocorrelation, and a decrease in their variance (slope = -0.12±0.01, p = 2.18.10 -8 ), while all these effects were entirely missed by a memoryless tolerance curve (continuous vs dashed lines in Figure 3b,c,d). The predicted and observed responses to environmental autocorrelation were in very good quantitative agreement for the mean and variance of population growth rate (Figure 3b,c), and were qualitatively consistent for the skewness of population growth rate (Figure 3.d). Temporal autocorrelation of the growth rate displayed discrepancy between the fit from the data and prediction from the tolerance curves (Figure 3e). Nevertheless, the tolerance curve with memory allowed for negative autocorrelation of population growth rates despite positive environmental autocorrelation, as observed in the data but not possible with a memoryless tolerance curve.

A candidate mechanism: plasticity of glycerol content. Although a range of physiological mechanisms are probably involved in these population responses to salinity fluctuations, glycerol content is a key candidate phenotypic trait. Glycerol acts as an osmoregulator, with intracellular concentration responding nearly linearly to the environmental salinity [START_REF] Ben-Amotz | The role of glycerol in the osmotic regulation of the halophilic alga Dunaliella parva[END_REF] . It is the major mechanism that makes Dunaliella a model species for salinity tolerance [START_REF] Ben-Amotz | The role of glycerol in the osmotic regulation of the halophilic alga Dunaliella parva[END_REF] . To investigate the possible role of this metabolite in our population responses to salinity, we tracked the dynamics of intracellular glycerol following a change in salinity. When going from high (3.5M NaCl) to low salinity, intracellular glycerol concentration first changed almost immediately, with significant differences in glycerol contents as soon as 16 minutes after transfer in 0.5, 2 or 3.5M NaCl (Supplementary Table 1), probably through excretion [START_REF] Zidan | Loss of intracellular glycerol from Dunaliella tertiolecta after decreasing the external salinity[END_REF] , followed by a slower adjustment of glycerol level (Figure 4). In contrast, the reciprocal transfer from low (0.5M NaCl) to high salinity, which requires de novo production of glycerol [START_REF] Ben-Amotz | The role of glycerol in the osmotic regulation of the halophilic alga Dunaliella parva[END_REF] , triggered slower change in glycerol content, with significant differences between transfer in 0.5 and in 2M NaCl found only after 135 minutes (Figure 4 and Supplementary Table 1). Most cells even died before their glycerol content could be measured in the transfer to the highest salinity (Figure 4). This is consistent with our observation that transitions from low to high salinity are more detrimental to the population dynamics than the reverse (Figure 3a).

Extinction rates increase with decreasing autocorrelation.

Nearly half the lines went extinct by the end of the experiment (Figure 5a). Environmental autocorrelation had a strong impact on extinction risk. Populations from low autocorrelation treatments (ρ = -0.5 and ρ = 0) went more extinct (72 and 61% extinct by the end of the experiment, median survival time: 77 days and 83 days, respectively) than populations from positive autocorrelation treatments (33 and 34% extinct by the end of the experiment for ρ = 0.5 and ρ = 0.9, respectively; p < 1.10 -4 , log ranked test performed independently for all genotypes, see Supplementary Table 1 for details). We found no evidence for an effect of genetic variation on extinction risk in our experiment, as single-strain lines did not differ significantly from mixed-strain lines, and there was no indication that recombination contributed to population persistence in mixed-strain lines (Supplementary Information). Simulations using the tolerance curve with environmental memory combined with the actual times series of environments gave results qualitatively similar to our data (Figure 5b). By contrast, a memoryless tolerance curve reversed the predictions, with the lowest extinction risk predicted under negative or null autocorrelation. The reason for this reversal is that, without environmental memory (and neglecting rapid evolution), environmental autocorrelation does not help phenotypic tracking of the environmental optimum, and therefore has no direct effect on the distribution of population growth rates (dashed lines in Figure 3b-d), but only influences how growth rates are integrated into the population dynamics. Because all things being equal, more autocorrelated growth rates cause larger variance in population size (as confirmed in simulations without memory, p < 10 -3 ), they are expected to increase extinction risk in the absence of environmental memory [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF][START_REF] Heino | Extinction risk under coloured environmental noise[END_REF] . This effect is totally buffered by environmental memory: with higher mean and lower variance in growth rates, populations experiencing positively autocorrelated environments are less likely to fall below a critical threshold for extinction.

Discussion

Random environmental fluctuations largely contribute to natural baseline rates of extinction [START_REF] Ovaskainen | Stochastic models of population extinction[END_REF][START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Lewontin | On population growth in a randomly varying environment[END_REF][START_REF] Lawson | Environmental variation and population responses to global change[END_REF][START_REF] Drake | Population effects of increased climate variation[END_REF][START_REF] García-Carreras | Are Changes in the Mean or Variability of Climate Signals More Important for Long-Term Stochastic Growth Rate?[END_REF] , which are currently aggravated by trends such as global climate change. Furthermore anthropogenic activities, beyond their effects on mean environments, are altering the patterns of these random fluctuations, including their temporal autocorrelation [START_REF] Boer | Changes in interannual variability and decadal potential predictability under global warming[END_REF][START_REF] Wigley | Anthropogenic influence on the autocorrelation structure of hemispheric-mean temperatures[END_REF] . There is thus is a pressing need to empirically characterize and quantify drivers of extinction risk in an autocorrelated stochastic environment. Using a highly replicated experiment with many independent environmental time series, we found that environmental autocorrelation can have a strong impact on population dynamics and extinction risk. In other words, extinction risk depends not only on the mean and variance of environmental fluctuations, but also on the order in which environments are encountered, within and across generations. This implies that the color of environmental noise, as increasingly documented in the ecological literature [START_REF] Vasseur | The color of environmental noise[END_REF][START_REF] Ripa | Noise colour and the risk of population extinctions[END_REF][START_REF] Heino | Extinction risk under coloured environmental noise[END_REF][START_REF] Ruokolainen | Ecological and evolutionary dynamics under coloured environmental variation[END_REF][START_REF] Pike | The effect of autocorrelation in environmental variability on the persistence of populations: An experimental test[END_REF][START_REF] Petchey | Environmental colour affects the dynamics of single-species populations[END_REF][START_REF] Laakso | Environmental noise and population dynamics of the ciliated protozoa Tetrahymena thermophila in aquatic microcosms[END_REF] , is likely to be an important driver of population dynamics. Importantly, the impacts of environmental autocorrelation on population dynamics could be predicted accurately in our experiment using the simple and general concept of the environmental tolerance curve, which summarizes all influences of the environment with a few parameters. This validates the utility of this tool as a way to predict population dynamics in a changing environment, by connecting physiology to ecology [START_REF] Deutsch | Impacts of climate warming on terrestrial ectotherms across latitude[END_REF] . However, we only reached accurate predictions when an influence of past environment was accounted for in these tolerance curves.

The large influence of environmental memory on population dynamics in our experiment implies that the phenotype of an individual cell is partly affected by the environment experienced by its recent ancestor(s), which can be described as trans-generational phenotypic plasticity, or environmentdependent parental effect, a form on non-genetic inheritance [START_REF] Calosi | Thermal tolerance, acclimatory capacity and vulnerability to global climate change[END_REF] . There is abundant evidence in microbes that prior exposure to stress can affect later stress responses [START_REF] Leroi | Temperature acclimation and competitive fitness: An experimental test of the beneficial acclimation assumption[END_REF] , and such acclimation experiments generally span over several cellular generations. In our experiment, genetic evolution could in principle also contribute, because some lines are highly polymorphic, but the time scale of ~3 generations between transfers is not consistent with large evolutionary change between transfers. In addition, clones isolated from mixed BC populations displayed similar tolerance curves as the whole population (Extended Data 5), confirming that (trans-generational) phenotypic plasticity is the major mechanism responsible for this short-term environmental memory.

Our experiment included elements of realism that are uncommon in the laboratory, such as continuously distributed fluctuations of an environmental variable, and a large number of independent realization of the same stochastic process. This represents a large step towards bridging the gap between theory and natural populations on quantitative rather than qualitative grounds, an important goal for population ecology. Note that in our batch culture setting, the environment changed every 3 or 4 days, i.e. every 2 to 5 generations (depending on salinity) for our model species D. salina, while salinity remained constant during the interval between transfers. This means that the autocorrelation we report is not a value per generation or per unit time as in theoretical models, but across transfers. This limitation is common to any experiment that uses organisms with continuous overlapping generations in batch culture, since changing the environment continuously over time in a stochastic way is essentially impossible. Such discretization does not hamper our conclusions about the influence of environmental correlation on population growth and extinction risk in our experiment, but it should be kept in mind when comparing quantitatively the effect of environmental autocorrelation in our experiment versus in wild populations of higher organisms, where censuses and environmental measurements are made at regular intervals (typically a year 2 ).

Microbes are inherently prone to transgenerational effects because they lack soma-germline differentiation, and therefore transmit their cytoplasm content and regulation factors to daughter cells upon reproduction. However, trans-generational effects mediated by non-genetic inheritance are increasingly documented in multicellular eukaryotes [START_REF] Bonduriansky | Extended heredity : a new understanding of inheritance and evolution[END_REF] , and have been shown recently to affect demographic parameters in a fluctuating environment in an animal [START_REF] Proulx | How differing modes of transgenerational inheritance affect population viability in fluctuating environments[END_REF] . This suggests that our results are also relevant to natural populations of larger organisms with conservation concerns, and that plasticity and non-genetic inheritance should generally be considered when investigating extinction risk in a randomly changing environment.

Methods

Dunaliella strains

We used two closely related strains of the facultative sexual microalgae Dunaliella salina,

CCAP

Experimental design

Each ancestral genetic background was exposed to three constant salinities (with 5 replicate lines per salinity), and 156 fluctuating salinity times series. The latter consisted of 39 independent time series (with the first one replicated 3 times) for each of four autocorrelation treatments: ρ = -0.5 (blue noise), ρ = 0 (white noise), ρ = 0.5 and ρ = 0.9 (red noise).

At each transfer, a 15% aliquot of each population was replicated in fresh medium with controlled salinity. The same dilution rate (d = 15%) was used in all treatments, and set so as to compensate for the intrinsic growth rate across variable environments for the slowest growing strain (strain B). However, growth rates were overestimated in our preliminary measurements, so that the dilution rate that we used caused the loss of most B lines before the end of the experiment.

Constant treatments were low ([NaCl] = 0.8M), medium (2.4M), and high (3.2M) salinities. The stochastic treatments involved continuously ranging, (near) normally distributed salinities, with same mean 2.4M and variance 1 among treatments, but temporal autocorrelation set by the treatment. We generated fluctuating salinity time series as follows:

-Generate a first-order autoregressive (AR1) process of length 1000000, mean 2.4, variance 1 and autocorrelation set by the treatment. An AR1 process is a stationary Gaussian process, in which the value at a given time step is a linear combination of the value at the previous time step plus an independent normally distributed noise.

-Truncate this theoretical process so as to keep only reachable salinity values. This involved removing salinities below 0 and above [NaCl] max = 4.8M, as well as salinities that made the transition from the previous time point impossible because of the constraints imposed by the 15% dilution upon transfer.

-Estimate the (reduced) variance of this truncated process -Increase input variance (step += 0.05) and repeat all previous steps, until the reduced variance

V = 1 ± 0.01
-Generate 39 times series of length 37 time steps for each autocorrelation treatment, using as input variance the value from the iterative search described above. The mean, variance, and temporal autocorrelation of these times series were checked to conform to their setpoint values. The distribution of salinities in all time series for the different autocorrelation treatments did not deviate from the expected Gaussian (Extended Data 4).

Transfers and growth

Experimental lines were cultured in 96 deep-well plates with artificial saline water + 2% Guillard's F/2 marine water enrichment solution (Sigma; G0154-500ml). Lines were distributed in 18 Plates were covered with plastic lids, sealed with Parafilm and placed at a fixed position in a growth chamber, with temperature set at 24°C, and light at 200 μmol.m -2 .s -1 for 12:12h LD cycles.

Population density measures

After each transfer, cell density before dilution was measured using 3 complementary methods. We

Statistical analysis

We analyzed population dynamics using an ad hoc state-space model. As typical in state-space models [START_REF] Dennis | Estimating density dependence , process noise, and observation error[END_REF] , our model comprised two main components: (i) an underlying, unobserved process describing the true stochastic dynamics of the populations size; (ii) and an observation model, which describes how this process translates into measurements, with errors that are independent after conditioning by the underlying process. We wrote an explicit likelihood function in C++ and optimized it in R (version 3.5.2), using the TMB package [START_REF] Kristensen | Automatic Differentiation and Laplace Approximation[END_REF] . R and C++ codes are available from the Dryad digital repository [START_REF] Rescan | Code from: Phenotypic memory drives population growth and extinction risk in a noisy environment[END_REF] .

Observation model

In a preliminary experiment, we determined the error distributions of cytometer counts, absorbance and fluorescence error. We made replicated (x3) measured of population sizes in serial dilutions of a Dunaliella culture, with concentration ranging from 30 to 10 6 cells.mL -1 , and fitted generalized linear regression between measurements and the known relative population sizes, expressed as a ratio of the maximal population size. We compared AIC from normal, lognormal, poisson and negative binomial models. Absorbance and fluorescence error were respectively normally and log-normally distributed, while cytometer measures followed a negative binomial distribution, with mean the actual population size. We thus used as the observation model for cytometer counts , fluorescence and absorbance , 5 where denotes a negative binomial, N is the true population density (governed by the unobserved process), and the are parameters that control the error variances (or overdispersion for cytometer counts) of these measurements.

Process model

We assumed that population density followed continuous logistic growth:

6
where N(t) g,ρ,i is the population density for the population from strain g, environmental treatment ρ, and time serie i at time t, K g the carrying capacity of genotype g and r(t) g,i,ρ its intrinsic growth rate in the environment experienced at time t. In the constant treatment analysis, K g was also allowed to vary between salinity treatments (K g,S ).

Given that r(t) g,ρ,i is constant over the time interval between two successive transfer, population size at transfer T + 1 can be implemented from population size at transfer T by integrating equation 6( 3 or 4 days), following a 15% dilution rate (d):

7

A convenient property of state-space models is that observation errors are independent conditional on the process [START_REF] Valpine | Fitting population models incorporating process noise and observation error[END_REF] (and independent of the process). Therefore, the overall log-likelihood in our model was simply the sum of the log likelihood of the process (equation 7, 8, 9, 10) and that of all observations (equation 5), producing an explicit formula amenable to numerical optimization.

Density dependence analysis using constant treatments

As our model was unable to estimate density dependence directly from the time series of fluctuating salinities, we first analyzed the population dynamics from constant time series in an independent model for this purpose. We considered that both the intrinsic growth rate r g,S and the carrying capacity K g,S were constant across replicates and transfers for one strain g in one salinity S. The probability distribution of population size at transfer T + 1, for strain g and salinity treatment S is then 8 where is the continuous logistic function given by equation 8 and is the process error. For each genotype, likelihood ratio test between a model where was constant across salinity and a model where varies with salinity was performed. We found that differed significantly but slightly between salinity treatments (Extended Data 2 d), so we used a constant for each genotype in further analysis.

Fluctuating growth analysis

Environmental fluctuations caused variation in the intrinsic rate of increase r in all our stochastic lines.

The distribution of r results from the interaction between moments of the stochastic environment (mean, variance and autocorrelation) and the growth response of each genotype to the environment.

The distribution of r may be characterized by its mean, variance, but also skewness [START_REF] Chevin | Stochastic Evolutionary Demography under a Fluctuating Optimum Phenotype[END_REF] and potentially by the autocorrelation of r across transfers. We estimated the parameters of this distribution for each autocorrelation treatment x genotype, within the state space model described by equations 5 and 8. We fitted either an autocorrelated normal (equation 9) or a reverse gamma distribution (equation 10) for , 9 or 10 where and are the maximal growth rate, the shape and the scale parameter for the gamma distribution for strain g in the environment. Mean , standard deviation and skewness were simultaneously estimated (estimation and standard error) from , and using the ADREPORT procedure of the TMB package.

Population densities were estimated as random variables with a residual component of variation, and were hence only very weakly affected by the assumed shape of the growth rate distribution (r² > 0.99 between the normal and the reverse gamma estimates, see Extended Data 6). This allowed us to visualize the goodness of fit for the reverse gamma and the autocorrelated normal. Using the population size estimates, we computed the realized intrinsic growth rate, corrected for density dependence:

11 and plotted them against the predicted distribution fitted by the normal and reverse gamma model (Extended Data 1).

To assess the extent to which population fluctuations were driven by responses to fluctuating salinity, we then used salinity as an environmental covariate. We explicitly expressed the growth rate as a function of both the current ( ) and the previous salinity . We modeled a bivariate tolerance curve using a 2 nd order polynomial, with parameters depending on the strain g,

12
leading to equation 2 in the main text.

Population size at T + 1 then followed a logistic growth (equation 7) with growth rate determined by current and previous salinity following equation 10, and we assumed a constant Gaussian distribution of the residuals generated by the process error, with standard deviation β:

13 Similar analysis was performed using only the salinity of growth ( ) in order to fit a tolerance curve without environmental memory (equation 1 in the main text):

14

Moments of population density and intrinsic growth rates.

We used the population density outputs of the gamma model (equation 10) to estimate the mean, variance, skewness and autocorrelation of log population density of non-extinct lines during the stationary phase (day > 40), for each autocorrelation treatment genotype combination (39 independent lines over up to 26 transfersreplicates 2 and 3 of series 1 were excluded from this analysis). Confidence intervals for these estimates were obtained by bootstrap (1000 simulations), where population densities were sampled for each simulation from a normal distribution with mean and standard deviation given by their corresponding estimate and standard error. We then performed mixed linear regression between these computed moments and the environmental autocorrelation, with ancestral genotype as a random factor, using the R package lme4 51 . To include the uncertainty in the estimation of population size moments, we performed 1000 generalized linear regression with population moment sampled from a normal distribution with parameters their estimate and error, and corrected for multiple imputation using the R package MICE [START_REF] Van Buuren | Multivariate imputation by chained equations in R[END_REF] .

Similar regressions were performed on 1000 simulations of the growth rate moments, sampled from a normal distribution with parameters the estimate and standard error obtained from the gamma model.

Mean, variance, skewness estimates, and their associated errors, were directly fitted in the gamma model. Autocorrelation of the intrinsic rate of increase was computed using the realized population sizes estimated in the gamma model, with similar bootstrap and sampling as used for the population log density analysis.

Fitted moments of the intrinsic rate of increase were compared to analytical predictions based on the tolerance curves, with or without memory. Subsequent salinities in our time series closely follow a binormal distribution (except for a weak truncation caused by our experimental design, see Extended Data 4), with mean µ = 2.4M and variance σ 2 = 1 at both steps, and correlation given by the environmental autocorrelation ρ. Combining this with the tolerance curve with memory leads to explicit formulas in the Mathematica 11 notebook available from the Dryad repository [START_REF] Rescan | Code from: Phenotypic memory drives population growth and extinction risk in a noisy environment[END_REF] . In particular, we found that the mean growth rate increases linearly with the environmental autocorrelation ρ(E), while the variance Var(r) is a quadratic function of ρ(E). Analytic formula for the skewness and the variance are more complex, but show notably that environmental memory allows for negative autocorrelation of the growth rate ρ(r), even when the environment is positively autocorrelated, a pattern that we found in our data but that could not be achieved without a memory effect.

Without memory (equation 14), the moments of the distribution of growth rates take a simpler form, summarized in Supplementary Table 2 (details in Mathematica 11 notebook in Dryad digital repository [START_REF] Rescan | Code from: Phenotypic memory drives population growth and extinction risk in a noisy environment[END_REF] ). In the absence of environmental memory, the mean, variance and skewness of the intrinsic rate of increase do not change with environmental autocorrelation ρ. Skewness is always negative for humped-shaped tolerance curves (c univ < 0). The effect of the environmental autocorrelation ρ(E) on the autocorrelation of the intrinsic rate of increase ρ(r) depends on the distance between the salinity optimum and the mean environment. When the phenotypic optimum corresponds to the mean environment, ρ(r) = ρ²(E), while ρ(r) tends rapidly towards ρ(E) as the mismatch between environmental mean and salinity optimum increases.

Survival analysis

We determined the extinction time for each population in the stochastic environment. Because OD, fluorescence and cytometry measures were not precise enough to discriminate between population sizes below 1000 cells.mL -1 , we considered that a population was extinct from the time it fell below an estimated 1000 cells.mL -1 (corresponding to a small number of cytometer counts, and a negligible density compared with the carrying capacity close to 10 6 ) and remained under this threshold for at least 5 transfers. After this delay, observed growth events that occurred only in 8 populations were due to contamination.

The survival responses in each treatment were then assessed by plotting the Kaplan-Meier survival curves (non-parametric model for the ratio of non-extinct populations through time), and analyzing differences in survival by log rank tests (ggsurvplot R package [START_REF] Kassambara | survminer: Drawing Survival Curves using 'ggplot2[END_REF] ). Independent analyses were performed for each genotype, and average extinction rates and autocorrelation effects were found similar in all genotypes except B. To illustrate these effects, we plotted the Kaplan-Meier survival curves and median survival time (time by which 50% of the population were extinct) for the joined genotypes in each autocorrelation treatment (Figure 5).

Simulations

Tolerance curves fitted with (equation 12) and without (equation 14) memory were used to simulate population dynamics under our fluctuating salinity treatments. Process noise and observation errors were not simulated. Population size variance was analyzed using the same bootstrap method as used for the real dynamics, except that here the true simulated population sizes were known and not sampled from their estimation distribution. Simulations did not include demographic stochasticity, so simulated lines were considered extinct after population densities were found below 1000 cells.mL -1 .

Intracellular glycerol dosage

A background culture of genotype C was used to measure variation in intracellular glycerol concentration in response to salinity changes. The ancestral population was cultivated in 2.4M for 10 days until it reached the late exponential phase. It was then split and acclimated for 7 days in 50mL flasks at low (0.5M) and high (3.5M NaCl) salinities. Each of these flasks was then split again and transferred to 0.5M, 2M and 3.5M NaCl in 50mL flasks, and sampled after 16, 45, 72 minutes, 2h15, 4h15 and 1 day for glycerol dosage.

Intracellular glycerol concentration was estimated from the difference between total and extra-cellular (after cell filtration) glycerol concentrations. 800µL Free Glycerol Reagent (Sigma Aldrich) was mixed to 200µL culture (with or without cells), and we measured optical density at 540nm after 5 minutes incubation at 37°C. Glycerol concentration was then interpolated from a (linear) standard curve. Dunaliella densities were also estimated by flow cytometry, and we made two independent replicates of all measures. Intracellular glycerol concentration per cell was computed as ,

where [] denotes a density/concentration. We computed the first-order standard errors of the intracellular glycerol concentration using the delta method, accounting for the (independent) error of the cytometer, the total and the extracellular glycerol measures (Figure 4). We then performed Wald test to assess whether different glycerol levels were found under different treatments (Supplementary Table 1). or without (c) environmental memory, using the same salinity time series as in the experiment.

Figures legends

  plates, such that each combination of ancestral genotype x environmental treatment was equally represented in all plates. Population positions in each plate were then randomized across the 60 central wells (external wells were not used because they experience larger evaporation). A liquid handling robot (Biomek NX from Beckman) was used to serially transfer these lines twice per week (alternating three-day and four-day cycles) into fresh medium. The salinity in this new medium was set independently in each well by adjusting the volumes of a hyposaline ([NaCl] min = 0M) and a hypersaline ([NaCl] max = 4.8M) solution, following 3 where [NaCl] ρ,i,T and [NaCl] ρ,i,T -1 are the setpoint salinities at transfers T and T -1 (respectively) for time series i of autocorrelation treatment ρ, V tot = 800 µL is the total volume of culture, and d = 15% is the dilution rate as defined above.

4 where

 4 measured optical density (absorbance at 680nm) and fluorescence (excitation at 390-80nm, emission at 685-40 nm) in a 200 µL sample of the each population using a BMG ClarioStar spectrophotometer.Cells from the same sample were then counted by flow cytometry (Guava EasyCyte HT). Dunaliella cells were isolated from debris using forward scatter (FSC), side scatter (SSC), red (695/50 nm) and yellow (583/26 nm) fluorescence emissions (excitation 488 nm, See Extended Data 3a., b. and c. for details). Cytometer performance was checked using a standardized fluorescent bead reagent (EasyCheck Kit, Merck Millipore) before each measure. Because of time limitation, only half of the plates went into the flow cytometer after each transfer, while the remaining plates were read in the next transfer, such that each plate was passed through the flow cytometer once per week.To calibrate among our different measures of population density, we first regressed the fluorescence and OD over the population density as estimated from cytometer counts (see Extended Data 3 d. and e.), excluding the data points where spectrometer measures were too low and exhibited high uncertainty (fluorescence < 400 and absorbance < 0.01). We obtained: OD blank and Fluo blank are the fluorescence and absorbance in wells known to be empty (control wells or extinct populations) within each given plate.

Figure 1 :

 1 Figure 1: Salinity, growth rates and population dynamics for two time series with different

Figure 2 .

 2 Figure 2. Population dynamics in a stochastic environment. Colors correspond to the temporal

Figure 3 .

 3 Figure 3. Population growth rate in response to environmental autocorrelation. a. Tolerance

Figure 4 :

 4 Figure 4: Intracellular glycerol concentration in genotype C after transfer to a new salinity. The

Figure 5 :

 5 Figure 5: Population extinctions across environmental autocorrelations. a. Survival curves

  19/12 (A) and CCAP 19/15 (C), and one more distant strain CCAP 19/18 (B). The lines were not axenic nor clonal when received from Culture Collection of Algae and Protozoan (Glasgow).We exposed them to increasing salinities (up to 4.8M) during several weeks to eliminate putative non-Dunaliella species of eukaryotes and finally isolated them by cell sorting (BD FACS Aria IIu, up to 10 5 cells). We then initiated 1074 lines from six ancestral backgrounds, with either low or high genetic diversity. Single strain lines (A, B and C) were used for the low genetic diversity treatment, and high genetic diversity lines were initiated by a 50% mix of two strains (AB, BC and AC).

Acknowledgements

This work was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement STG-678140-FluctEvol). We thank the MRI-IGMM platform for access to the cell sorter, the GenSeq platform (LaBEX CEMEB, Montpellier) for sequencing, and O. Cotto, C.A Klausmeier, R. Lande, and J. Tufto, and three anonymous reviewers for helpful feedback.

Data availability statement

The population dynamics and glycerol acquisition data that support the findings of this study are available from the Dryad digital repository [START_REF] Rescan | Data from: Phenotypic memory drives population growth and extinction risk in a noisy environment[END_REF] .

Code availability statement

R, C++ and Mathematica codes are available from the Dryad digital repository [START_REF] Rescan | Code from: Phenotypic memory drives population growth and extinction risk in a noisy environment[END_REF] .

Author Contribution

L-M and M.R designed the experiment. D.G and M.R conducted the population dynamics experiment.

L.M and M.R performed the statistical analysis. D.G. and M.R extracted and amplified barcode genes of all lines and D.G, E.O-A and M.R analyzed the bioinformatics data.

Competing interests

The authors declare no competing interests.