Gene expression in diapausing rotifer eggs in response to divergent environmental predictability regimes

In unpredictable environments in which reliable cues for predicting environmental variation are lacking, a diversifying bet-hedging strategy for diapause exit is expected to evolve, whereby only a portion of diapausing forms will resume development at the first occurrence of suitable conditions. This study focused on diapause termination in the rotifer Brachionus plicatilis s.s., addressing the transcriptional profile of diapausing eggs from environments differing in the level of predictability and the relationship of such profiles with hatching patterns. RNA-Seq analyses revealed significant differences in gene expression between diapausing eggs produced in the laboratory under combinations of two contrasting selective regimes of environmental fluctuation (predictable vs unpredictable) and two different diapause conditions (passing or not passing through forced diapause). The results showed that the selective regime was more important than the diapause condition in driving differences in the transcriptome profile. Most of the differentially expressed genes were upregulated in the predictable regime and mostly associated with molecular functions involved in embryo morphological development and hatching readiness. This was in concordance with observations of earlier, higher, and more synchronous hatching in diapausing eggs produced under the predictable regime.

The variability that characterizes natural environments confronts organisms with the vital challenge of surviving during temporarily adverse periods. To cope with such adversity a widespread strategy among metazoans is diapause 1 . As a kind of dormancy, diapause permits escape in time through the temporal suspension of metabolic activity, arrested development and increased stress resistance 2,3 . In contrast to quiescence (another kind of dormancy controlled exogenously and implying an immediate response to a change in environmental conditions 4,5 ), diapause depends on endogenous control for its initiation and termination 6 . Thus, diapausing organisms enter the passive state even under optimal environmental conditions that would otherwise promote normal metabolism and development, typically preceding the onset of adverse conditions 7 . In diapause, metabolic arrest remains for an internally scheduled period of time, the so-called refractory period 8 , even if suitable conditions resume. Only after the refractory period is completed do specific environmental cues disrupt diapause to restore active development and reproduction 9,10 . However, the decision to leave diapause remains a challenge if the cues anticipating future environmental conditions are unreliable and lead with some probability to zero-fitness events (death and reproductive failure in the case of annual organisms) 11 . To reduce the risk associated with such unpredictability organisms can employ a bet-hedging strategy for waking the diapausing forms whereby only a portion of these resume development at the first occurrence of cues for suitable conditions while the rest wait for later opportunities for reactivation [12][13][14] . The adequacy of diversifying bet-hedging strategies in leaving diapause (a single genotype producing a variety of phenotypes with different diapause durations 15 ) to cope with unpredictable environments has received increasing support in the form of empirical studies [16][17][18] . Access to high-throughput next-generation sequencing (NGS) technologies can motivate studies pursuing better understanding of the molecular mechanisms behind the optimal duration of diapause.
Although it may be not realistic to speak of discrete steps in diapause, a convenient simplification is to divide it into the overlapping ecophysiological phases of initiation, maintenance and termination 19 . Initiation may be driven by differential gene expression leading to altered feeding regimes, changes in the synthesis and activity of cell proteins, and modification of behaviour and/or morphology. Subsequent to initiation, organisms enter the

Results
Diapausing egg hatching experiment. Both predictability regime and diapause condition affected diapausing egg hatching fraction. The effects of both factors and their interaction were significant (Table 1). Hatching fractions were significantly lower under the unpredictable regime for both diapause conditions (Non Forced Diapause (NFD) and Forced Diapause (FD)) ( Fig. 1A,B). Timing of hatching was also significantly affected by selective regime (survival log-rank test; χ 2 = 236, p-value < 0.001) and diapause condition (χ 2 = 28.1, p-value < 0.001). We observed a delayed hatching pattern in populations from the unpredictable regime, with a median hatching time (H 50 ) averaging 9.66 ± 1.20 days and 8.00 ± 0.57 days under NFD and FD conditions, respectively. Under the predictable regime, H 50 values averaged 4.66 ± 0.33 days and 2.66 ± 0.66 days for NFD and FD conditions, respectively. No matter the selective regime, FD condition resulted in significantly lower values for the interquartile range of hatching time (IQR, here taken as a measure of hatching synchronicity) Table 1. Summary of the generalized linear mixed-effect model (GLMM) on diapausing egg hatching fraction, and the linear mixed-effect model (LMM) for the interquartile range (IQR) of hatching times. Note that the significance of the "population" effect for the diapausing egg fraction was tested using a likelihood ratio test (LRT; see "Material and methods").  Table 2. An average of 36,896 ± 388 expressed genes were identified for each combination of predictability regime and diapause condition (U-NFD, U-FD, P-NFD, P-FD). A quality control plot of RPKM density distribution (the squared coefficient of variation (CV 2 ) vs log 10 FPKM) for each combination of predictability regime and diapause condition is shown in Supplementary Figure S1. Overall, the variability in gene expression levels in diapausing eggs produced under the unpredictable regime was always higher than in those produced under the predictable regime, whatever the diapause condition tested.
Variation partitioning analysis of gene expression. Differentially expressed genes (DEGs) associated with the "regime" and "diapause condition" factors, and their interaction, were assessed by (1) the log 2 -fold change (log 2 -FC) values between the levels of each factor and their associated p-values, and (2) by the proportion of expression variance explained by each factor (see "Material and methods"). Our analysis also allowed testing the random effect of population within each selective regime. A total of 18,598 genes were scored and assigned a variance fraction and a differential expression p-value.
The relative amounts of variance in gene expression explained by the predictors of our model are summarized in Fig. 2. Overall, we observed that the "regime" factor accounted for the highest proportion of the variation in gene expression in B. plicatilis s.s. diapausing eggs, with a median percentage close to 40%. In contrast, the  www.nature.com/scientificreports/ "diapause condition" factor and its interaction with "regime" explained approximately 10% of variance in gene expression. The random effect of population within regime was negligible (Fig. 2). Consistent with analysis of variance partitioning, post hoc corrections of the associated p-values of differential gene expression showed that 76% of the whole set of expressed genes retrieved in our analysis were not differentially expressed between selective regimes. No gene was found to be differentially expressed between diapause conditions after p-value correction for multiple comparisons. Since both approaches strongly suggest that diapause condition was not as important in driving differences in the transcriptional profile of diapausing eggs than the selective regime, hereafter we focused only on the latter factor. Variation in gene expression levels among replicates. The estimated biological variation in gene expression levels among replicates follows a very similar trend along the range of average expression values in the two regimes, with much higher variance among genes expressed at very low levels (Supplementary Figure S2). Only among those low-expression genes, expression level is more variable among replicates in the unpredictable compared with the predictable regime. When comparing whole expression profiles between samples, a test of homogeneity of dispersion based on permutations fails to detect any significant difference between the two regimes (p-value = 0.883; Supplementary Figure S3).

Gene function assignment and GO enrichment analysis of DEGs.
To find specific functions related to diapause maintenance or to embryonic development (diapause termination) we inspected GO terms significantly enriched with genes differentially expressed between predictability regimes following the ordination provided by the variation partitioning analysis. Notably, both approaches used for gene ordination yielded similar results (Supplementary Figure S4). A total of 9691 DEGs had GO annotations. In our inspection we included Table 2. Statistics for the filtering and mapping of RNA-Seq reads corresponding to combinations of selective regimes and diapause conditions. P predictable regime, U unpredictable regime, NFD non-forced diapause, FD forced diapause. www.nature.com/scientificreports/ only the most specific GO terms assignable to a gene and skipped those that were too vague or general. Supplementary Tables S1, S2 and S3 show the GO specific terms that were significantly associated with the "regime" factor according to all three algorithms used in the enrichment analysis for biological processes (GO:0008150), molecular functions (GO:0003674), and cellular components (GO:0005575). First, we made a tentative assignment of GO terms to diapause phases on the basis of previous knowledge. Proper references and rationales for GO term assignment to different phases of diapause are shown in Table 3. Next, we used volcano plots for GO specific terms associated with the "regime" factor to better demonstrate the role of DEGs in diapause maintenance (Fig. 3) and termination (Fig. 4). The volcano plots indicate the size of the biological effect (fold change) versus the statistical significance of the result (p-value) in association with the unpredictable regime. Positive log 2 -FC values indicate upregulated genes in the unpredictable regime with respect to the predictable one, thus negative values are indicative of under-expression. These plots were intended to help identify subsets of genes involved in different functions within the entire set of DEGs, and to focus further research on the most promising. A rough estimate of the level of activity of a molecular pathway is the number of genes differentially expressed between the selective regimes. By this criterion, two pathways stand out among the a priori functions assigned to the maintenance of diapause (Fig. 3): G protein-coupled receptor signalling (GPCR) and signal transduction. There are 206 genes differentially expressed in the GPCR pathway (55 with an adjusted p-value lower or equal to 0.1) and 544 genes involved in signal transduction. Of course, not all molecular pathways are equally complex, and a better approach involves comparing the fraction of genes differentially expressed in a given regime. Thus, for instance, the nucleotide-excision repair term is among the most significantly enriched terms with 13 genes annotated. Notably, all 13 were overexpressed in the diapausing eggs from populations evolved in the unpredictable regime. For other functions with a putative role during diapause, including those involved in protection and isolation from external stimuli, the change in expression was not in the same direction for all genes that were regulated by the selective regime; however, we observed that higher numbers of genes were in general expressed in diapausing eggs from populations evolved under the predictable regime.
Functions associated with embryo development included pathways related to (1) cellular association with the extracellular matrix and cell motility, (2) cilium assembly and motility, (3) ion transport, and (4) protein catabolism. We observed that most genes in the pathways were underexpressed in the unpredictable regime ( Fig. 4). Overall, for the set of functions depicted in Fig. 4, the proportion of genes that showed increased expression in the predictable regime averaged 77.8 ± 3.9%.

Discussion
This study focuses on a highly relevant stage of the life cycle of B. plicatilis s.s., the diapausing egg stage, and explores patterns of diapause maintenance and termination. We combined hatching experiments and transcriptomics to evaluate diapausing eggs produced by populations evolving under two divergent regimes of environmental predictability and subjected to conditions either promoting or blocking hatching. Different hatching patterns were observed between selective regimes and diapause conditions and more than 3000 genes with different expression profiles between the studied comparisons were identified. To the best of our knowledge, this is the first study that uses RNA-Seq and a reference genome to relate hatching patterns to transcriptional profiles in diapausing eggs subjected to predictable and unpredictable laboratory environments.
Intermediate hatching fractions and variable diapause duration have been proposed as adaptations to environmental unpredictability in rotifers 11,18,29,33 . Diapausing eggs from populations subjected to a predictable regime showed higher, earlier and more synchronous hatching than those produced by populations under an unpredictable one. Interestingly, these traits have been observed to evolve quickly under the divergent predictability regimes assayed in our experimental evolution approach 18 . Consistent with the initial hypothesis, the timing of diapausing egg hatching was earlier and more synchronous after a forced period of diapause. A striking result was that hatching fractions were slightly lower in FD than under NFD conditions. This may be due to an increased diapausing egg mortality due to the forced diapause period 62,63 , despite the apparent health of all eggs in our assays, a trait that has been shown to correlate with their ability to hatch 64 .
Our results suggest that selective regime is more important in driving differences in the transcriptome profiles of diapausing eggs than the conditions of diapause. Interestingly, while the numbers of up-and downregulated DEGs are similar in either regime, the absolute fold change is on average twice as great in DEGs upregulated in the predictable regime. It is likely that most DEGs overexpressed in the predictable regime are related to the reactivation of diapausing eggs and hatching readiness. This would be in concordance with the results from the hatching experiment in which diapausing eggs produced under the predictable regime showed earlier and higher rates of hatching. Interestingly, after 30 h of exposure to hatching induction conditions embryo movement was observed in most of the diapausing eggs from the predictable regime (personal observation). The incidence of embryonic movement was lower in the diapausing eggs from the unpredictable regime. Correspondingly most genes associated with molecular functions related to morphological development of the embryo (cell-matrix adhesion, cell motility, cilium movement and cilium assembly) and cytoskeleton (motor activity) were upregulated in diapausing eggs obtained from populations subjected to the predictable regime under both diapause conditions.
We also found evidence of upregulation of protein catabolism-related genes (peptidase activity including metallocarboxypeptidases, metalloendopeptidases, serine-type endopeptidases and calcium-dependent cysteinetype endopeptidases) in diapausing eggs from populations subjected to the predictable regime. Protein catabolism may serve as a source of amino acids for development after diapause. These results have points in common with previous studies that explored the expressions of transcripts involved in protein turnover and synthesis of ATP and cytoskeletal elements 36 . Interestingly, research on copepods has also shown that proteins are increasingly Known as a repair mechanism in desiccation tolerant stages 37 . Desiccation/ rehydration cycles and prolonged periods in the dry state are associated with remarkable levels of stress to the embryo genome which can result in mutagenesis of the genetic material, inhibition of transcription and replication and delayed growth and development Response to oxidative stress (10/26) Oxidative stress leads to the accumulation of toxic components called reactive oxygen species (ROS). Therefore, production of antioxidant enzymes is critical for stress resistance during diapause 38 and has been shown in a wide range of organisms such as insects 39 , nematodes 40 and rotifers 34 . The major ROS detoxifying enzymes described in rotifers include peroxidases, thioredoxins, catalases and glutathione-S-transferases 34,36 Oxidoreductase activity (63/147) A number of proteins with oxidoreductase activity have been described as involved in diapause with a role in detoxification processes and degradation of xenobiotics 36 . For example, the cytochrome P450 monooxygenase gene has been shown to be overexpressed in insects during diapause 3,5 . In the case of monogonont rotifers this gene is also related to lipid metabolism (steroidogenesis 41 ) with a potential role in the accumulation of energy reserves for diapause 42 Trehalose biosynthetic process (4/6) Trehalose is related to dormancy in many prokaryotes and eukaryotes, and is known to be synthesized under dehydration conditions to stabilize other proteins and membranes during drying 41,43,44 . Trehalose acts as a water replacement molecule and vitrifying agent serving to stabilize cell membranes 45 . Trehalose is known to accumulate in diapausing cysts of Artemia 45,46 and in diapausing rotifer eggs when they enter anhydrobiosis 47  www.nature.com/scientificreports/ utilized as an energy source by the end of diapause when neutral lipid storage becomes depleted 65,66 . It is known that the diapausing eggs of rotifers store fat reserves 42 and it is likely that harvesting of lipids for the embryo starts at the time of diapausing egg production (diapause initiation); however, no molecular function associated with lipid metabolism had annotated DEGs that related to the predictability regime in our study. This could be because their expression is highly conserved in diapausing eggs independent of environmental predictability regime and conditions during diapause. Indeed, genes related to general lipid metabolism have been suggested to exert an ambivalent role in rotifer diapause 35 . Therefore, during diapause organisms may preferentially catabolize neutral lipids while conserving and synthesizing polar ones, which have an essential role in maintaining cell membrane structure and integrity. Interestingly, we found transcripts for the fatty acid synthase gene complex, known to be involved in lipid accumulation and stress tolerance during diapause in beetles 67 . A hypothesis to be tested is whether the usage of neutral lipids can be arrested during diapause and resumed exactly at the time of reactivation by using the stored lipids that have not been depleted during diapause. This makes sense in rotifers since large numbers of droplets with neutral lipids have been found in females that hatch from diapausing eggs 68 . This finding suggests a role for lipids in rotifer development after hatching and in the improvement of their ability to colonize new environments 42,69,70 .
We found a less probable, delayed, and less synchronous hatching pattern in diapausing eggs under the unpredictable regime than under the predictable regime under both diapause conditions. This hatching pattern is consistent with the general downregulation observed in the diapausing eggs from populations subjected to the unpredictable regime, possibly associated with developmental arrest during diapause 32,69,71 . The maintenance of diapause in a larger fraction of eggs or for a longer time under the unpredictable regime entails a higher Table 3. Gene ontology (GO) term assignment and rationale to the different phases of diapause. Numbers in parentheses denote the number of significant differentially expressed genes (p-value < 0.01) out of the total number of genes annotated with that term.

Diapause phase GO term Rationale
Embryo development Cell-matrix adhesion (3/16) The interactions of cells with the extracellular matrix play crucial roles during morphogenesis in developing embryos. Interestingly, Clark et al. 36 have reported the upregulation of genes involved in cell proliferation and adhesions in the transcription profile of non-diapause embryos in their morphological development into neonate rotifers Potassium ion transport (4/56) K + is known to have a role in cell growth through protein synthesis. There is evidence of K + -mediated stimulation of aminoacyl-tRNA (aa-tRNA) synthetase activity in bacteria 54 and insects 55 , with an effect on diapause termination in the latter. Rozema et al. 25 have shown high levels of aa-tRNA synthetases in diapausing compared to asexual eggs in Brachionus Ion transmembrane transport (16/147) Free ions in the cytosol restart the cell machinery by reactivating cellular enzymes 56 . Dumont et al. 57 showed that cyst hatching in some anostraceans increased after the addition of mobile ion carriers. There is also evidence that increased pH induces diapause termination in Artemia salina 26  www.nature.com/scientificreports/ uncertainty concerning the behaviour of the eggs, mirroring the uncertainty of the environment. We tested whether this more unpredictable behaviour in the diapausing eggs translated into more variable gene expression profiles among samples. Only genes expressed at very low level were observed to have more variable expression among samples from the unpredictable regime. The overall biological variation in gene expression levels among samples was of a similar magnitude in both regimes. In general, it is expected that diapausing eggs that have not completed their mandatory diapause period may actively maintain the molecular pathways involved in the processes of repair, protection, and isolation from external stimuli 5,43,72 . Consequently, our specific prediction was that genes belonging to these pathways would be overexpressed in diapausing eggs from populations evolved under the unpredictable regime. Our results partially agreed with this prediction in the case of the expression of nucleotide-excision repair genes to reverse genome damage incurred during diapause. Monogonont rotifers are not unique in exhibiting this ability. Other metazoans that exhibit enhanced DNA repair capacity include tardigrades 73 and bdelloid rotifers 74 . It has been hypothesized that such ability is a consequence of evolutionary adaptation to survive desiccation and exposure to high-energy radiation in ephemerally aquatic habitats 74 .
Functions providing protective mechanisms to the diapausing eggs in rotifers may have also evolved to withstand the physical effects of desiccation 23,34,36 . Since experimental populations of both regimes underwent desiccation-rehydration cycles, it is not surprising that we detected gene expression for some of these functions in their diapausing eggs. Consistent with previous reports on diapausing rotifer eggs, the GO terms associated with protective functions we retrieved in our analysis referred to pathways for trehalose biosynthesis (see also 23,34 ), response to oxidative stress (thioredoxins, peroxidases, and glutathione-S-transferases 35,36 ) and oxidoreductase activity (the cytochrome P450 gene 1 ). However, our results show that most of the genes involved in the three abovementioned pathways were more overexpressed in diapausing eggs from the predictable regime compared with those from the unpredictable one, which remains unclear. The fact that a few genes within these pathways were still overexpressed in the unpredictable regime suggests that diapausing eggs from populations evolved under different selective regimes may be using different sets of genes to protect themselves against desiccation or oxidative stress. This observation may be related to the probability that diapausing eggs from one regime or the other are in a different ecophysiological phase of diapause (maintenance or termination). Naturally, www.nature.com/scientificreports/ both developing embryos and rotifer hatchlings must have mechanisms to protect themselves against oxidative stress 75,76 . It is also conceivable that some genes are annotated as involved in a function that they regulate negatively, a hypothesis deserving further research. While in diapause, rotifer embryos are less sensitive to external stimuli. Therefore, we expected that the molecular pathways responsible for environmental isolation would be more expressed in diapausing eggs that have not completed the refractory period compared with those that have resumed development. In relation to predictability regime, we hypothesized that a higher proportion of diapausing eggs from populations evolved under the unpredictable regime have not passed this period compared with eggs from populations evolved under the predictable regime. Contrary to our hypothesis, most of the genes in the signal transduction and GPCR signalling pathways were underexpressed in the unpredictable regime. Let us here recall the likely ambivalent role of some gene families in relation to environmental information processing. As mentioned above, rotifer diapausing eggs need a light stimulus to re-initiate development. In the absence of light, prostaglandins have been www.nature.com/scientificreports/ shown to cause the same effect 77 . Interestingly, some photoreceptors, such as the opsins, belong to the family of GPCRs, and prostaglandins are known to bind with a GPCR receptor.
In conclusion, our comparative transcriptome analysis, supported by hatching data, provided mechanistic insight into the adaptation of rotifer diapause to habitat unpredictability. Overall, we found evidence that molecular pathways related to diapause maintenance and termination were differentially expressed across rotifer diapausing eggs produced under two laboratory environments: predictable vs. unpredictable. Molecular functions related to embryo development reactivation were in general upregulated in diapausing eggs from predictable environments, consistent with the higher hatching fraction of eggs produced under this selective regime. However, some of the genes related to diapause maintenance were also upregulated in the predictable regime under both diapause conditions assayed. It is possible that expression of some molecular pathways is highly conserved in diapausing eggs independent of the selective regime and conditions during diapause but also that different sets of genes within a pathway are expressed at different phases during diapause 2 . Our study also extends knowledge of the complex molecular and cellular events that take place during diapause. Some of the genes identified here are well known in other anhydrobiotic organisms and resistance forms, but several are new and should be further investigated to determine their involvement in desiccation and stress tolerance. In this sense gene assignment of function to either diapause maintenance or embryonic development processes, although based on previous literature, is still tentative and further research is needed. In fact, several studies have shown that gene families traditionally associated with diapause maintenance, such as some glutathione-S-transferases, aquaporins, ferritins, or trehalose metabolism genes 34 , are highly expressed in diapausing eggs after exposure to 30 h of light 35 , so these genes may also be involved in embryonic development reactivation. Thus, we call for further research to achieve a better understanding of the present results and to further elucidate the genes related to diapause in unpredictable environments. Given future increases in environmental variability due to climate change 78 , understanding the molecular mechanisms underlying differential hatching patterns such as those found in the present and previous studies is essential to understand how organisms cope with environmental unpredictability.

Material and methods
Experimental design and diapausing egg collection. Diapausing eggs were produced by laboratory populations of B. plicatilis s.s. cultured in an experimental evolution design that tested the effects of hydroperiod length unpredictability on diapause-related traits. Briefly, six genetically diverse laboratory populations were generated by placing together three ovigerous asexual females from each of 30 clonal lines of B. plicatilis s.s. from nine natural populations (a total of 270 clones). Thus, the initial genetic composition of each of the six experimental populations was the same. Clonal lines were founded from the hatchlings of diapausing eggs isolated from the sediments of nine natural ponds in Spain covering a natural range of habitat predictability (see full details in 18 ). Because B. plicatilis s.s. belongs to a cryptic species complex, clones were identified as belonging to B. plicatilis s.s. by genetic analysis of COI based on PCR-RFLP (see 17 ). The six laboratory populations were subjected to two contrasting time-fluctuating laboratory environments, predictable (P) vs unpredictable (U), for eight growth cycles of selection (Fig. 5). These growth cycles were interrupted by periods of habitat unsuitability (dry periods) during which diapausing eggs were collected and then after a fixed period of obligate diapause (28 days) used to re-form each experimental laboratory population. Three laboratory populations were randomly assigned to the predictable selective regime, characterized by growing seasons (hydroperiods) of constant length (28 days). The other three laboratory populations were assigned to the unpredictable regime, characterized by growing seasons whose length varied randomly (between 4 and 53 days with an average of 28 days).
At the end of the last cycle of selection (the eighth growing cycle), diapausing eggs were collected, cleaned and counted. Diapausing eggs obtained from each laboratory population were randomly divided into two groups to be subjected to two different diapause conditions: (1) "non-forced diapause" (NFD) and (2) "forced diapause" (FD) (Fig. 5). For the first condition, a total of 5200 newly produced diapausing eggs per laboratory population were placed in Petri dishes and immediately induced to hatch under standard hatching conditions (constant light, 6 g L −1 salinity and 25 °C 64 ) for 30-36 h. This period of time is considered to be enough to reactivate diapausing embryo development but not to complete hatching (see 35 ; personal observation: no hatchlings were observed when handling the diapausing eggs). Thus, diapausing eggs were exposed to hatching conditions and not forced to remain in diapause (NFD). For the second condition a similar 5200 diapausing eggs per laboratory population were air-dried in darkness at 25 °C and then stored in darkness at 4 °C for 28 days. These conditions typically inhibit hatching 42,79 so diapausing eggs were forced to remain in diapause (FD) during this period of time. After this period of forced diapause, diapausing eggs were placed in Petri dishes and induced to hatch under standard hatching conditions, as for the NFD condition, for 30-36 h. As depicted in Fig. 5 the experiment followed a split-plot design 80 in which the "predictability regime" was the whole-plot factor, the experimental populations were the subplots nested within the whole plots, and "diapause treatment" was the so-called split-plot factor. A total of 12 diapausing egg samples resulting from 2 selective regimes (P vs U) × 3 laboratory populations within regime × 2 diapause conditions (FD vs NFD) were obtained to be further subjected to both hatching and gene expression analyses.
Hatching experiment. Subsamples  www.nature.com/scientificreports/ sample were spun down in 1.5-mL microcentrifuge tubes and any remaining water was removed. Total RNA of each sample was isolated using the TRIzol Plus RNA Purification kit (Ambion, Life Technologies, USA) following the manufacturer's instructions. In the first step of RNA isolation, once TRIzol reagent was added, a pellet pestle motor was used to homogenize the sample and assess the lysis of all diapausing eggs. DNase treatment during RNA purification was performed to reduce the amount of genomic DNA. At the end of the RNA extraction protocol each sample was quantified and the ratios 260/280 and 260/230 were assessed with a NanoDrop spectrophotometer (Thermo Scientific). Since the NFD condition was started one month earlier than the FD condition, the total RNA isolated from NFD samples was kept at − 80 °C. When the 12 samples of RNA were obtained, they were sent frozen to the Servei Central de Suport a la Investigació Experimental (SCSIE) of the University of Valencia (Spain) where RNA-Seq was performed.

RNA-Seq.
In the Genomic core facility at SCSIE the samples were quantified using a Qubit 2.0 Fluorometer (Life Technologies) and the purity and integrity of the RNA was assessed using electrophoresis in a Bioanalyzer 2100 (Agilent). After this quality control, cDNA libraries were constructed using TruSeq stranded mRNA (Illumina) with an enrichment of mRNA using the standard poly-A based method followed by chemical fragmentation and cDNA synthesis. Libraries were sequenced using 75 nt single-read sequencing on an Illumina NextSeq 500 platform. The quality of the raw reads was assessed using FastQC version 0.11.5 (www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/). Adapter sequences were trimmed and reads < 20 nt were discarded.

Data analysis.
Hatching data. The fraction of hatched eggs for each laboratory population and diapause condition assayed (NFD and FD) was estimated excluding deteriorated eggs 64 . The pure effects and interaction of the fixed-effect factors "regime" and "diapause condition" on the hatching fraction of diapausing eggs were analysed using a generalized linear mixed-effect model (GLMM) for a split-plot design considering binomial errors and logit as link function 81 . Factor "population", nested within "regime", was treated as a random effect whose significance was tested by means of a likelihood ratio test (LRT) comparing a model with all effects versus a model in which the "population" effect was not included. www.nature.com/scientificreports/ the interquartile range (IQR) as a measure of dispersion. Thus, the higher the IQR, the greater the hatching asynchrony. Differences between selective regimes (P vs U) and between diapause conditions (FD vs NFD) in the IQR were assessed using a linear mixed-effect model (LMM). All analyses were performed in R, version 3.2.2 83 . GLMM and LRT were implemented by means of the glmer and anova functions of the "lme4" package 84 . LMM was performed using the lmer function of the same package. For Kaplan-Meyer survival analysis and Cox's model fitting, the Surv and coxph functions from the "survival" package 85 were used.
Transcriptome assembly. Transcriptome assembly was performed using the Tuxedo protocol 86 . RNA-Seq reads were mapped to a draft genome of B. plicatilis s.s. 87 via TopHat v. 2.1.1 88 using the following parameters: -I/-maxintron-length = 15,000; -i/-min-intron-length = 10; -max-multihits = 1; -read-gap-length = 1; -read-realign-editdist = 0; -library-type fr-firststrand. For the rest of parameters default values were used. To assemble the aligned reads into transcripts, Cufflinks software 86 was used with default parameters modified as in TopHat's mapping. Maximum and minimum intron length in TopHat and Cufflinks was assessed using the reference genome of B. plicatilis s.s. The values applied in the abovementioned parameters intended to follow a conservative approach to detect differential gene expression. Finally, RSEM (RNA-Seq by Expectation Maximization 89 ) was used for quantifying both gene and isoform abundances based on the read mappings. Although the abundance of gene transcripts is typically expressed as RPKM 88 , the raw counts of reads mapped to transcripts were also used in subsequent analyses (see below).
Differential expression of genes. Expression data for genes in the dataset were fitted to the split-plot design described above (Fig. 5) using the "variancePartition" package 90 in R, which allowed for the handling of fixed ("regime" and "diapause condition") and random effect ("population") factors. The dream function within this package was used to test for differential expression associated with each factor. Prior to model fitting, raw counts of reads were transformed to log 2 -counts per million mapped reads using the voom function in the "limma" package 91 ]. We considered a gene to be expressed if it had at least five reads, and we only chose for analysis the genes expressed in at least four samples. Differentially expressed genes (DEGs) associated with the "regime" and "diapause condition" factors were assessed separately in two alternative ways. On one hand, DEGs were ordered as usual by taking into account the log 2 -fold change (log 2 -FC) values between the levels of each factor and their associated p-values. To estimate the proportion of true null hypotheses among the whole set of statistical tests performed we used the least slope method 92 implemented in the estim.p0 function of package "cp4p" in R 93 . On the other hand, DEGs were also ordered by the proportion of expression variance explained by each factor. The latter was the preferred criterion for ordination in our analyses since genes with low total variability in expression may still be highly associated with a given factor even if the fold change between the levels of that factor is low. Notwithstanding, the differential expression data based on log 2 -FC and associated p-values were still used to validate the results of the variation partitioning approach.
Comparison of among-replicate dispersion of gene expression levels between selective regimes. We used the "DESeq2" package 94 in R to estimate the dispersion of gene expression levels among biological replicates separately for samples of the two selective regimes. Gene-wise dispersion values are expected to reflect biological variation in true proportions of transcripts among replicates. They are estimated from few replicates by assuming a similar mean-variance relationship among genes 95 . In addition, we computed dissimilarities of gene expression profiles among samples using Poisson distance 96 . We tested for homogeneity of dispersions between the predictable and unpredictable regimes using the betadisper function of the "vegan" package 97 , that is, whether samples from the unpredictable regime were more different in terms of gene expression profiles within regime than samples from the predictable regime.
Functionality assignment of DEGs and Gene Ontology (GO). The gff2fasta tool of CGAT toolkit 98 was used to extract the sequences of the DEGs. Then sequentially, we used TransDecoder (https ://githu b.com/Trans Decod er/Trans Decod er/wiki) to identify the proteins encoded within each transcript, and InterProScan 99 to assign functional annotations to the proteins identified by TransDecoder.
We used the R package "TopGO" 100 to run a gene ontology (GO) enrichment analysis among the genes that were differentially expressed between selective regimes or hatching conditions using a Kolmogorov-Smirnov test and a combination of three algorithms (weight01, lea and elim). GO terms that were significantly associated with the factor in question according to all three algorithms were chosen for further analysis. Because the association of expression with each classification factor was quantified in two different ways (see "Differentially expressed genes" section), the functional enrichment analysis was run separately on the dataset of DEGs ordered by the proportion of expression variance explained by the factor concerned, and the dataset ordered by the p-values of the corresponding tests of differential expression. Both approaches were compared to validate and generalize the results of differential gene expression associated with the study factors.
Finally, we used information from previous studies that have provided very useful genetic resources for the discovery and function assignment of genes associated with diapause and embryonic development in B. plicatilis [34][35][36] , Brachionus calyciflorus 101 , and Brachionus manjavacas 58 as well as in other aquatic and terrestrial invertebrates. This information was used to (1) perform a tentative assignment of GO terms to diapause maintenance and termination, and (2) establish links between the observed hatching phenotypes and their corresponding RNA expression profiles.

Data availability
The sequence data generated during the current study are available in compressed fastq format in the Zenodo repository (https ://zenod o.org/recor d/39536 65#.Xx63n qYp5h E). The whole analysis described here is registered and documented in executable scripts available in the GitHub repository https ://githu b.com/ignas iluca s/Brach ionus .