Statistical analyses
Before the statistical tests, we evaluated the assumptions of normality and equality of variances using Kolmogorov–Smirnov and Levene tests. To compare critical thermal limits and demographic variables among experimental groups, we included thermal treatment as a factor with six levels (C, V, CV, CC, VC and VV) that describe the thermal experience of flies (direct experience through ontogeny, and indirect thermal experience, through parental thermal exposition; see Figure 1). This factor allowed us to compare phenotypic response between P and F1, and also to perform comparisons across P and F1 groups. To compare critical thermal limits among experimental groups, we employed linear mixed model with trial as random effect (random intercept) and sex and treatment as predictor variables. Also, to test the potential of trade-off between fitness related traits ( \(R_{0}\) and \(T_{g}\))and CTmax the one tailed correlation analysis was conducted.
Because the population density effected significantly \(R_{0}\) and\(T_{g}\) (Supplementary Table S1), we assessed the global response of\(R_{0}\) and \(T_{g}\) to temperature and density. We performed a nonparametric regression analysis using a generalized additive model (GAM) incorporating populational density (D), parental thermal environment (T P), offspring thermal environment (T F1), and thermal treatment (treat ) as predictors. We performed GAM since it does not make any a priori assumptions about the shape of relationships between variables, which is key to our evaluation of the effects of population density. Moreover, the main difference between GAMs and linear models is that linear functions of the variables in GAM are replaced by unknown smooth functions, giving additional flexibility to the modeling process (Wood, 2017) .The complexity of the curve (the number of degrees of freedom) and the smoothing terms were determined by penalized regression splines and generalized cross-validation to avoid overfitting (Wood, 2017). Also, we allowed the shrinkage of the smoothers. This technique allows for an extra penalty to be added in the model, and if the penalty is high enough, it will shrink all smoothing coefficients to zero. Model selection was done using the AIC criterion (ΔAICc < 2; Burnham & Anderson (2002)). To perform pairwise comparisons between experimental groups, we performed a posteriori Tukey test following the linear mixed models or GAMs.
All analyses and visualizations were performed in the R statistical environment (http://www.R-project.org/ ). The datasets generated during the current study will be available in the DRYAD repository