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.