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.