Transgenerational response to early spring warming in Daphnia

Temperature and photoperiod regulate key fitness traits in plants and animals. However, with temperature increase due to global warming, temperature cue thresholds are experienced at shorter photoperiods, disrupting the optimal seasonal timing of physiological, developmental and reproductive events in many species. Understanding the mechanisms of adaptation to the asynchrony between temperature and photoperiod is key to inform our understanding of how species will respond to global warming. Here, we studied the transgenerational mechanisms of responses of the cyclical parthenogen Daphnia magna to different photoperiod lengths co-occurring with warm temperature thereby assessing the impact of earlier spring warming on its fitness. Daphnia uses temperature and photoperiod cues to time dormancy, and to switch between sexual and asexual reproduction. Daphnia life cycle offers the opportunity to measure the relative contribution of plastic and genetic responses to environmental change across generations and over evolutionary time. We use transgenerational common garden experiments on three populations ‘resurrected’ from a biological archive experiencing temperature increase over five decades. Our results suggest that response to early spring warming evolved underpinned by a complex interaction between plastic and genetic mechanisms while a positive maternal contribution at matching environments between parental and offspring generation was also observed.

role of adaptive transgenerational plasticity is key to predict the consequences of parental effects on population dynamics, and to inform our understanding of how species will respond to rapid environmental change 18 . Presently, the prevalence and strength of TGP in natural systems remains controversial 13 .
Temperature and photoperiod are two pivotal cues that regulate key fitness traits in plants and animals (e.g. flowering time and dormancy). However, whereas climate is changing, photoperiod is stable across latitudes 19,20 . The asynchronous change of temperature and photoperiod cues has disrupted the optimal seasonal timing of physiological, developmental and reproductive events in many species [21][22][23] . Therefore, many organisms have experienced a decline in fitness when preparing to develop, reproduce, hibernate, enter dormancy, or migrate at seasonally inappropriate times 22 . Parthenogenetic zooplankton species (e.g. daphniids, rotifers) use both temperature and photoperiod cues to time life history events such as male production and dormancy 21,24 . Generally, short photoperiod and low temperatures induce a switch from the production of clonal (parthenogenetic) females to the parthenogenetic generation of males and to initiating sexual reproduction, which culminates in dormant embryos 25 . Conversely, long photoperiod and high temperatures are cues for the termination of dormancy in spring 21,22 . Dormant embryos are early stage embryos that arrest their development to overwinter; dormant embryos at the water-sediment interface hatch during the following growing season contributing to the genetic diversity of the local population [26][27][28] . A proportion of these embryos becomes buried in the sediment missing the opportunity to hatch thereby creating biological archives, which represent an excellent resource to investigate population dynamics in natural systems over evolutionary time [29][30][31] . Because of the changes in water temperature due to global warming, the temperature cue thresholds are experienced by some parthenogenetic zooplankers at shorter photoperiods altering population-level investment in dormancy and male formation 32,33 . However, the impact of persistent asynchrony between temperature and photoperiod across generations is poorly understood. Understanding the transgenerational impact of this asynchrony is helpful to predict how rapidly these species, and the community which they sustain, may track climatic change.
Here, we study the impact of the asynchrony between temperature and photoperiod on the fitness of the freshwater crustacean Daphnia magna. Daphnia are keystone species central to the food-web of lentic freshwater ecosystems worldwide 34,35 . We study within-generation plastic and genetic responses, as well as the cross-generational responses of fitness-linked life history traits to environmental conditions that mimic early spring warming (warm temperature and short photoperiod) as compared to a typical spring environment (warm temperature and long photoperiod).
Studies that investigate mechanisms of response to environmental change face two main limitations. If genetically diverse species are studied, the same genotype or a diversity of genotypes cannot be replicated across environmental exposures and transgenerational responses to environmental change are confounded by genetic variation. Conversely, if clonal species are used and genotypes are replicated across experiments, genetic diversity is not representative of sexually reproducing species. Daphnia's life cycle overcomes these limitations to reveal the relative contribution of plastic and genetic responses to environmental change across generations and over extended time periods 36,37 . Daphnia alternates sexual recombination with asexual (clonal) reproduction. Sexual recombination results in early stage embryos that can be 'resurrected' and propagated indefinitely as clonal lineages in the laboratory 38 . These properties and a short life cycle provide the advantages of clonal species while retaining the natural genetic diversity. Taking full advantage of these properties, we measured transgenerational fitness changes in three populations of dormant embryos previously resurrected from the biological archive of Lake Ring and spanning 50 years 39 . We quantified WGP and TGP as changes in trait values of fitness-linked life history traits between environments and across generations. We quantified genetic responses to photoperiod as mean trait differences among populations and as genetic divergence at 15 candidate genes, previously linked to environmentally driven local adaptation 40,41 . We also performed a gene-traits association analysis between the life history traits measured in the parental generation and the 15 candidate genes in an effort to identify the genetic basis of fitness-linked life history traits responding to photoperiod and temperature cues. Our study provides important insights into the transgenerational mechanisms of response to early spring warming in a freshwater keystone species.

Methods
Experimental design. We investigate the impact of temperature and photoperiod cues on three populations of D. magna separated in time and previously revived from a biological archive of Lake Ring, a shallow mixed lake (maximum depth is 5 m) in Jutland, Denmark (55°57′51.83″N, 9°35′46.87″E) 42 . According to historical records 43 and paleolimnological analyses 30,36,38 , the lake experienced a steadily, even if modest (∼1 °C), increase in temperature over the last five decades 36 . Further, changes in water chemistry and transparency occurred over time because of eutrophication and agricultural land use. These changes are described in detail elsewhere 39 and summarized in Supporting Information Data 1. The populations are hereafter referred to as: P1: population resurrected from layers of sediment dated >1999; P2: population resurrected from layers of sediment dated 1975-1985; P3: population resurrected from layers of sediment dated 1960-1970. Daphnia reproduces parthenogenetically in the laboratory, allowing experiments across generations. Several parthenogenetic broods are produced in each generation (up to eight) allowing the same genotype to be used in parallel experiments. Between nine and eleven genotypes per population (Table S1; N = 30) were used in a common garden experiment over two parthenogenetic generations (Fig. 1). The sample size per population was chosen based on a previous study that determined the threshold sample size required to assess genetic diversity in Daphnia populations 30 . According to this study, 10-15 genotypes capture the genetic diversity of natural populations both in space and time.
Prior to the experiment, clonal lines established from the thirty genotypes were maintained for two generations in common garden conditions (20 °C, 12:12 h light: dark regime; fed ad libitum with 0.8 mg Carbon/L of Chlorella vulgaris CCAP strain no. 211/11B) to reduce potentially confounding maternal and grandmaternal www.nature.com/scientificreports www.nature.com/scientificreports/ effects (Fig. 1). Individual juveniles of 24 h from the second or following broods were randomly selected to establish experimental generations and assigned to experimental conditions. Different broods from the same generation were used to ensure developmental synchrony among clonal lineages in the experiment. Experimental clonal lineages were kept in individual jars of 100 ml, filled with filtered borehole water, which was refreshed every second day, and fed daily 0.8 mg Carbon/L of Chlorella vulgaris.
After two generations in common garden conditions, 24 h old individual juveniles from the second or following broods were randomly assigned to the experimental conditions for two generations (G1, G2): short photoperiod (SP; 10:14 h light: dark regime) and long photoperiod (LP; 14:10 h light: dark regime) (Fig. 1). The temperature for both photoperiods was 20 °C, so that changes in photoperiod could be studied in a typical spring environment (warm temperature and long photoperiod), acting as a control, and in conditions that mimicked early spring warming (warm temperature and short photoperiod).
The following fitness-linked life history traits were measured in G1 and G2: size at maturity (the distance between the head and the base of the tail spine), age at maturity (first time parthenogenetic eggs are released in the brood pouch), fecundity (number of juveniles across 8 broods), interval between broods (time interval between each pair of broods averaged across eight broods), proportion of male offspring across 8 broods, and mortality. For size at maturity, all animals were measured after releasing their first brood into the brood pouch using imageJ software (https://imagej.nih.gov/ij/index.html).
The experimental design is not full factorial because of the number of genotypes and conditions tested in parallel. However, because genotypes were fixed across experimental conditions and generations, we were able to control for confounding factors e.g. genetic changes occurring from one generation to the next and genetic variation among experimental exposures, enabling WGP and TGP effects to be studied in isolation.
Life history traits: genetic and plastic responses to photoperiod and temperature cues. We assessed evolutionary mechanisms (genetic changes, WGP, TGP, and their interactions) of response to early spring warming on life history traits using linear mixed models (LMMs) in the "lme4" package in R v.3.3.3 44 . We quantified the effects of Generation, Population, Treatment and their interactions on individual life history traits. Genotype was fit as a random effect nested within population. Prior to the analyses, all variables were tested for normality. Except for male proportion (a binomial variable), all life history trait measurements followed a Gaussian distribution. Because the populations separated in time originate from the same genetic pool and because genetic drift is negligible 30 , a significant population term indicates genetic differences among populations. Differences in mean trait values between photoperiods, after controlling for maternal effects, are the expression of WGP, an environmental effect. Differences in mean trait values between generations are the expression of TGP. If the effect of the treatment (photoperiod) differed significantly among populations (genetic effect), we would have evidence of a Pop (population) × Pht (photoperiod) interaction. Similarly, if the effect of the treatment (photoperiod) differed significantly among generations, we would have evidence of Gen (generation) × Pht (photoperiod) interaction. If population means varied by generation, we would have evidence of a Gen (generation) × Pop (population) interaction. We also measured the three-way interaction term (Gen × Pop × Pht). All evidence of interactions or main effects were assessed via Type II analysis of deviance tables using the Anova www.nature.com/scientificreports www.nature.com/scientificreports/ function in the "car" package (R v.3.3.3) 45 . We visualized the main effects of population (P1, P2 and P3), treatment (SP and LP) and generation (G1 and G2) plus their interaction terms on individual life history traits through reaction norms, which describe the pattern of phenotypic expression of each population across treatments and generations 46 .
We performed a principal component analysis (PCA) to quantify the principal modes of variation and covariation among life history traits within generation using the 'prcomp' function in the R "stats" package. The PCA plots were obtained using the "ggbiplot" (R v.3.3.3 44 ). Prior to the PCA analysis, the life history traits variables were log transformed.
Mortality rates per population were calculated with a survival model fit via the "psm" function in the rms R package V.3.3.3 44 . A separate model was fitted to each treatment and generation, in which the day of mortality and the mortality event were combined as the dependent variables (e.g. censoring) and population was treated as fixed effect. All mortality curves were plotted using the "survplot" function form "rms" package in R v.3.3.3 44 .
Candidate gene analysis: genetic responses to temperature and photoperiod cues. The genome of the 30 genotypes used in this study has been sequenced and will be published elsewhere. Here, we use the sequenced genome to identify SNP polymorphisms at 15 candidate genes previously associated with environmentally driven local adaptation, including male formation and the induction of sexual reproduction 40,41 by mapping the sequences of the candidate genes against the reference genome of D. magna 2.4 (NBCI: LRGB00000000.1). The panel of SNPs at the 15 candidate genes is shown in Table S2. Methods used for mapping and variant calling are in Supporting Data 1.
To assess whether the candidate genes were under selection in the populations studied here, we used measures of genetic differentiation (F ST ) and neutrality tests. Previous results identified a small yet significant proportion of neutral genetic divergence among the D. magna populations studied here (1%) 30 . Larger divergence among populations at the 15 candidate genes would identify genes under divergent selection, whereas a smaller divergence at the candidate loci would suggest balancing selection. In the former, different polymorphic sites at the same genes are found in different populations, resulting in higher divergence than at neutral loci. Conversely, balancing selection reduces population differentiation, because of moderate to intermediate frequencies of polymorphic sites 47 .
We quantified pairwise population differentiation as F ST , using MSA Analyser 48 , and the partitioning of genetic variance within and among populations (AMOVA) using Arlequin 49 . The two hierarchical levels used in the analysis are (i) among populations and (ii) within populations. Statistically significant values (P < 0.001) were calculated with permutation tests (10,000 permutations). Further, we performed neutrality tests on the candidate genes to quantify departure from neutrality using DnaSP 50 . Typically, these tests are affected by demography. In our study, the three populations are snapshots of the same genetic pool at different time points across five decades. Therefore, demographic confounding factors are negligible.
To assess whether the candidate genes were associated with life history traits showing responses to photoperiod, we performed a gene-trait association analysis between the candidate gene polymorphisms and life history traits measured in G1, using Mixed Linear Models in TASSEL 51 with the following criteria: minor allele frequency (MAF) > 0.05 and Hardy-Weinberg Equilibrium (HWE) > 0.001, applying a false discovery rate of FDR = 0.05 to correct for multiple testing. We investigated significant association between gene polymorphisms at the 15 candidate genes and the life history traits in long and short photoperiod, separately. We also quantified association between the gene polymorphisms and the plastic change in life history traits as delta values between long and short photoperiod within G1.

Results
Life history traits: genetic and plastic responses to photoperiod and temperature cues. We investigated the impact of photoperiod and temperature cues on fitness-linked life history traits in D. magna. Short photoperiod at 20 °C mimicked early spring warming, whereas long photoperiod at 20 °C mimicked a typical spring environment, acting as reference.
We quantified the effects of Generation, Population, Treatment and their interactions on individual life history traits using an analysis of variance. The three-way interaction term (Generation × Population × Photoperiod) was significant for the proportion of male offspring ( Table 1). The effect of photoperiod varied significantly between generations in four of the five life history traits (Table 1; Gen × Pht). The effect of photoperiod did not vary significantly by population, except for the proportion of male offspring (  Fig. S2).
Short photoperiod caused a decrease in size and in the time elapsing between parthenogenetic broods in the first experimental generation (Fig. 2; G1). In the second generation (G2), short photoperiod induced increase in fecundity, smaller size at maturation, and longer time elapsing between parthenogenetic broods ( Fig. 2; G2). Male proportion varied with generation and population (Fig. S1), showing both genetic and plastic responses (Fig. 2).
Candidate gene analysis: genetic responses to temperature and photoperiod cues. The genetic divergence among populations was higher at neutral loci than at the candidate genes in presence of negligible drift ( Table 2, Table S3). The neutrality test on the candidate genes identified nine out of fifteen genes as putatively under selection; at least two independent tests supported departure from neutrality in these genes ( Table 3). The candidate genes showed significant association with fecundity and male offspring proportion both in short and long photoperiod (Table S4; SP, LP). Fecundity showed significant association with the histone deacetylase complex subunit SAP18 (Table S4). The proportion of male offspring showed significant association with Amino-oxidase, SAP18, Phosphotyrosine (PTB), and Aldo/keto reductase (Table S4). In addition, the candidate genes showed significant association with plastic changes in all traits (Table S4; Delta). Plastic changes in age at maturity were significantly associated with the rhodopsin g-protein coupled receptor; plastic changes in fecundity were significantly associated with the Histone deacetylase complex subunit SAP18; variation in male offspring proportion was significantly associated with Amino-oxidase, SAP18, Phosphotyrosine (PTB), and Aldo/keto reductase; plastic change in size at maturity was significantly associated with Pyridine nucleotide-disulphide oxidoreductase, PTB, SAP18, and an RNA binding protein (Table S3).

Discussion
Investigating both plasticity and genetic adaptation across generations enables to establish a link between evolutionary potential (either additive genetic variation or genotype by environment interactions) and TGP [52][53][54] . We used transgenerational common garden experiments on resurrected D. magna populations separated in time to study the mechanisms of response -genetic, plastic (WGP, TGP) or a combination thereof -to earlier spring warming. Resurrection ecology provided us with the unique advantage of following evolution in situ across decades by studying genetic differences among populations that evolved from an initial genetic pool in response to environmental selection 55,56 . The use of genetically distinct clonal lineages across generations enabled us to perform the concurrent analysis of plastic and genetic mechanisms of adaptation to environmental change within and across generations minimizing confounding factors.
We observed significant difference in the plastic response to spring warming between generations, made evident from a significant generation term and a significant interaction term Generation × Photoperiod in the ANOVA. Trade-offs among life history traits clearly differed between generations. In G1 the time elapsing between parthenogenetic reproductions showed trade-offs with age and size at maturity, and fecundity showed trade-offs with the proportion of male offspring. In G2, size and age at maturation showed trade-offs with all other life history traits; moreover, the investment in male offspring came at a cost of fecundity and time elapsing between parthenogenetic reproductions.
Theory predicts that highly plastic traits show strong maternal effect variance and little to no genetic variance, because highly plastic traits can be influenced by the environment, including parental environment, and because additive genetic variance may be masked by high environmental variation 57 . For example, elytron length in seed beetles, a highly plastic trait, is strongly influenced by parental effects and has no detectable genetic variance 58 ; body and offspring size in stickleback is strongly influenced by the parental environment but does not show high genetic variance 59 . Our results suggest a positive maternal effect on earlier spring warming when the offspring environment perfectly matches the maternal one. This effect is evident from fecundity, which is higher in G2 than G1, and from size at maturity that does not change significantly in G2 but significantly decreases in G1 in short photoperiod. However, a full factorial design measuring the anticipatory maternal effect in matching and not matching environments is needed to confirm positive maternal effect in matching environments 13 . In addition, the photoperiod cue experienced by the second generation during fertilization may have influenced the phenotypic response of the offspring generation, result of developmental plasticity. The timing at which the www.nature.com/scientificreports www.nature.com/scientificreports/ second generation directly experience the parental environment has been shown to influence offspring response in matching environments 60,61 . In particular, exposure during fertilization and early juvenile development may strongly affect offspring response to parental environments 2 . Follow up experiments, in which the parental generation is exposed for different lengths of the life cycle to early spring warming, is needed to exclude whether exposure during fertilization influenced offspring response to early spring warming.
We observe a strong role of WGP in response to earlier spring warming, made evident from significant photoperiod effect on the majority of life history traits. The effect of photoperiod resulted in smaller size at maturation, delayed maturation and a decrease in the time lapse between parthenogenetic reproductions in G1, whereas it resulted in higher fecundity, and an increase in the time between parthenogenetic reproductions in G2.
We observe genetic responses to early spring warming quantified as genetic differences among populations in fecundity. Moreover, we observe a lower than neutral genetic divergence at the candidate loci, which were significantly departing from neutrality expectations. The lake from which the Daphnia populations were sampled has a documented increase in average ambient temperature and recurrence of heat waves over time 36 . Previous studies on these populations showed evidence for temporal evolution of the critical thermal maximum (CT max ), the upper temperature at which animals lose motor function 62 , in presence of warming as a single stress 63 . Here, we observe significant evolutionary differences among populations to early spring warming cues (G1). These evolutionary  Table 1.
www.nature.com/scientificreports www.nature.com/scientificreports/ differences indicate that the most recent population evolved to be less fecund than the historical populations. This may be the result of maladaptation or of evolutionary responses to multiple stressors. It is possible that multiple stressors documented in Lake Ring (e.g changes in water transparency and lake chemistry 39 ) have influenced adaptive responses of the population over time 64 . Indeed, in a previous study on the populations of Lake Ring, higher CT max in the most recent population was observed in presence of warming alone, whereas it was no longer observed in presence of multiple stressors 63 .
In the presence of stable neutral genetic diversity over time 30 , the candidate genes show a departure from neutrality and a lower divergence among populations than neutral markers. These patterns are expected when balancing selection affects the frequency of gene polymorphism across populations 47 Table 1: age at maturity, size at maturity, fecundity, male proportion and average interval between broods.

Among populations Within populations
Neutral µsat 1.08* 98.92* Candidate genes 0 100* www.nature.com/scientificreports www.nature.com/scientificreports/ Rhodopsin has been associated with dormancy in Daphnia 41 ; SAP18 has been associated with embryonic development and environmental stress response in insects 65,66 ; PTB (epidermal growth factor receptor kinase) has been linked to sexual differentiation in insects and crustaceans 67,68 ; MF, a catalyser for the methyl farnesoate and a putative juvenile hormone in daphniids, and AKR (aldo-keto reductase) have been linked to the propensity to form male offspring 69,70 . In the MF pathway, the AKR family catalyses the conversion of farnesal to farnesol, and mutations in the AKR gene have been shown to reduce MF production in favour of the juvenile hormone, resulting in higher male production 70,71 . The function of these genes and their significant association with life history traits in the population studied here suggest their potential role in the response to photoperiod and temperature cues. However, it is possible that other environmental stressors, which have been documented in Lake Ring, may have contributed to the observed response at the candidate genes and/or that additional genes not investigated here are contributing to modulate the life history traits analysed. Indeed, all traits showed significant association with multiple candidate genes suggesting epistasis among genes underpinning the life history traits. Furthermore, single genes showed association with multiple traits, suggesting extensive pleiotropy. It is expected that multiple genes regulate complex fitness-linked life history traits. In a follow up experiment, we generated further evidence supporting this expectation (Supporting Data 1). We used three genotypes that showed divergent patterns in male offspring proportion between photoperiods (Fig. S1) and a reference genotype, which never produces male broods under all tested experimental conditions to date (P-IT, Institute of Ecosystem study, CNR Verbania, Italy) 72 . We focused on the propensity to form males because some of the candidate genes studied here have been previously associated with this trait 40,69,70 . Fixed polymorphisms among the strains with divergent propensity to form male offspring would suggest that the candidate genes are underpinning this trait. Lack of evidence for fixed polymorphisms at the candidate genes would suggest that other genes underpin male offspring formation. We found that the divergent patterns in propensity to form male offspring in the four genotypes did not correspond to fixed polymorphisms at the 15 candidate genes, confirming pleiotropy and the need to investigate genome-wide patterns to identify the genes underpinning complex fitness traits (Supporting Data 1).
Overall, our results point to an adaptive potential for a keystone zooplankter Daphnia to evolve in response to early spring warming mediated by a complex interaction between plastic and genetic mechanisms. The results also suggest a positive maternal effect in presence of matching environments between parental and offspring generations. A positive maternal effect at a life history trait indicates accelerated rates of microevolution that can facilitate adaptation 73 . However, a full factorial design is needed to assess whether non-matching environments provide different responses.
We identified extensive pleiotropy in candidate genes underpinning life history traits responsive to early spring warming. Pleiotropy is commonly observed in complex environments with multifarious selection pressures acting on multiple aspects of the phenotypes, resulting in trade-offs among competing functions [74][75][76] . This finding supports previous experimental results showing that populations in complex environments can overcome fitness cost exhibiting synergistic pleiotropy 74 .  Table 3. Neutrality tests. Statistics testing for departure from neutrality of 15 candidate genes previously associated with environmental driven local adaptation and partial sex determination.