Phenotypic divergence
Phenotypic divergence was analysed using multivariate and univariate Bayesian linear mixed models. All traits were standardised to have a mean of zero and standard deviation of one to improve comparability and model convergence. All models described below were fit with a Gaussian error distribution in the MCMCglmm package in R . Fixed effects were given weakly informative flat priors, random effects were given default inverse Wishart priors, and we fit full unstructured covariance matrices for the multivariate models. Each model was run for a total of 1,020,000 iterations with a burn-in of 20,000 iterations and thinning of 1000 iterations, which resulted in low autocorrelation. Convergence of models was assessed by examining traceplots to visualise sampling mixing and by assessing effective sample sizes and autocorrelation.
Effects of site and habitat –To investigate the extent of phenotypic variation among sites and habitats, we compared three linear mixed effects models with different random effects structures using the deviance information criterion (DIC) to identify the model with the most support, at ΔDIC = 2 : a model with no random effects (i.e., a “null” model), a model that included site as a random effect, and a model that included habitat type as a random effect. These models were designed to identify whether there was phenotypic variance between all habitats or all sites. We grouped traits in the response variable to investigate phenotypic divergence in functionally correlated traits. TL and gut length were each fit as a univariate response; the four defence traits (plate number, DS1, DS2 and PS) and the four gill raker traits (GRL2, GRL3, GRW, GRN) were fit as multivariate response traits, respectively. This resulted in a total of 12 models (three models with different random effects structures for each of the four responses), all fitted with sex as a fixed factor. Models for gut length, defence traits, and gill raker traits included TL as a covariate. Fitting the interaction between sex and TL did not improve model fits, so we present results from models with sex and TL fit as single term effects. Owing to varying levels of replication for each sex at each site (Table 1) we cannot estimate or exclude sex differences in the patterns of divergence. The multivariate models fit full variance-covariance matrices for random effects, allowing us to estimate phenotypic covariances at both individual (residual) and habitat/site levels . Note that because of the different sample sizes for each group of traits (see below), it was not possible to run a single multivariate model with all traits to measure the full covariance matrix. Fish from sites were considered phenotypically divergent if 95% credible intervals of the posterior distributions of predicted trait values did not overlap.
Association with ecological parameters – We next ran a suite of univariate, linear mixed effects models that investigated the effect of ecological predictors on each phenotype independently. All models were fitted with TL and sex as fixed effects and site as a random effect (except for model on total length, which only had sex as a fixed factor). Because many of our eight ecological variables were highly correlated (Table 1), we ran one model per ecological predictor per phenotype (total = 7 models per phenotypic trait) and compared each to a null model without any ecological variables using DIC (as above). In these models, all predictor variables were standardised to have a mean of 0 and standard deviation of 1 to improve comparability between models. We present the mode and 95% credible intervals of the posterior distribution for the linear coefficients for each ecological predictor, unless ΔDIC to the null model was > 2 because this indicated that the null model was not improved by fitting the ecological variable.