Data analysis
The effects of treatments were analyzed for the following dependent variables: (1) body size after summer treatment; (2) sexual development time (days elapsed after the beginning of the winter treatment until the appearance of first gonads); (3) sexual fitness (total number of eggs produced by females and the maximum number of testes seen on a male); (4) asexual fitness (total number of buds produced by an individual from the start to the end of the experiment, i.e. during both the summer and winter treatments) and (5) survival rate. I used Generalized Linear Mixed Models implemented in the glmmTMB package (v1.1.4) in R version 4.2.2 (Brooks et al., 2017; R Core Team, 2020). Gaussian errors were assumed for body size and binomial errors for survival rate. For sexual development time, sexual and asexual fitness (count data), four different error structures were considered: Gaussian, Poisson, negative binomial with linear parametrization (nbinom1 ) and negative binomial with quadratic parametrization (nbinom2 ) (Brooks et al., 2017). The model best fitting to the data was selected with Akaike Information Criteria corrected for small sample size (AICc). After selecting the right distribution, I fitted a model with experimental treatment as a fixed effect and strain ID and batch ID as random effects. This model was compared to a random effects null model via Likelihood Ratio Tests (LRTs). Post-hoc treatments-versus-control comparisons with Dunnett’s correction were obtained using theemmeans R package (Lenth, 2022). For the model analyzing body size after summer treatment only one predictor was included (summer treatment) and I used LRT to evaluate the effect of treatment on size. Finally, I generated model diagnostics plot using the DHARMa R package (v0.4.5; (Hartig, 2022)) to check that model residuals did not show outstanding deviations from model assumptions. Since the model diagnostics for sexual development time showed substantial deviations from the expectation for both males and females, I also analyzed these variables with non-parametric Kruskal-Wallis tests, followed by pairwise Wilcoxon tests with Holm’s correction.