Statistical analyses
None of the control populations went extinct, so we only analysed the effect of morph enrichment on extinction in populations subjected to gradual warming after four generations, using the ‘glmer’ function from the ‘lme4’ package in R (Bates et al., 2015; R Core Team, 2021). The extinction was modelled as a binomial response variable (coded 0 for extinction and 1 for survival) with population type (FT vs. ST) entered as a fixed factor and replicate within population as random effect.
To analyse survival of individuals within both GTI and control populations, the effect of treatment (Control vs. GTI populations), population type (F vs. S) and generation (coded as continuous variable) were entered as fixed effects, and population and replicate within population as categorical random effects. A two-variable vector containing the numbers of survived and dead individuals (counted as difference between initial number of juveniles and number of adults at the end of the reproductive period) was a response variable. In order to investigate the effect of founder morph effect on survival of individuals that experienced increased temperature, GTI treatment-specific model with population type and generation as fixed effects was fitted. We similarly analysed sex-specific survival at the adult stage (immature stages cannot be sexed) in GTI populations, with a vector of the numbers of adults which survived and died by the end of reproductive period as a response variable. Finally, we analysed changes in male morph ratio (coded as a two-variable vector containing number of fighters and the number of scramblers) at emergence from nymphs, using GLMM with treatment and generation as fixed factors. The model was run only for F populations data because an evqiuvalent model for S population, as well as a full model with a tree-way interaction, could not be estimated due to high correlation between fixed effects in models including three-way interaction, and quasi compleate separation (a situation when the predictor is associated with nearly uniform value of response variable) in S population. The data were overdispersed in models testing for survival of individuals, sex-specific survival, and changes in male morph ratio. To account for this, beta-binomial GLMMs were conducted using the ‘glmmTMB’ function in ‘glmmTMB’ package (Brooks et al., 2017). The goodness of fit of the models were evaluated based on Akaike information criterion (AIC) by comparing AIC scores between model with the highest order interactions and simplified models, but this only lead to reduction of sex-specific survival model, for which the best fitted model (ΔAIC = 3.2) included two-way interaction between sex and generation, and population type as fixed effects. In all other cases, full models for all response variables had the lowest AIC scores. All statistics were run in R version 4.0.5 (R Core Team, 2021).