Phylogeny and comparative analyses
We built a supermatrix phylogeny of the Amphibolurine Agamidae based on 2 mitochondrial (ND2 and ND4) and 3 nuclear (BDNF, RAG-1 and BACH1) genes, built around a multi-locus nuclear gene backbone taken from the Zheng & Wiens supermatrix dataset (Pyron et al. 2013; Zheng & Wiens 2016). Full details of the supermatrix assembly, alignment and phylogenetic analysis are given in Supplementary Information (Supplementary methods, Table S6, Figure S6). We used a subset of 1300 post-burnin trees (subsampled using logcombiner (Bouckaert et al.2019) and pruned of all non-focal taxa (Phytools R package; Revell 2012) in subsequent phylogenetic comparative analyses.
We tested whether variation in the concentration of carotenoids and pteridines was associated with environmental gradients of habitat productivity (indirect compensation) or indices of sexual selection. The response variables in these models were: 1) total carotenoids; 2) total pteridines; and 3) the ratio of carotenoids to pteridines. The predictor variables were environmental PC1 and PC2, sexual size dimorphism and sexual dichromatism (which are uncorrelated, r2 = 0.05, Estimate = -0.003 – 0.009). Given that information on sexual selection indices only exists at the level of species rather than the individual, we also ran species-level models (27 species). We calculated total carotenoids, total pteridines and the ratio of carotenoids to pteridines based on average pigment concentration per species. We used these measures as the response variables and the two indices of sexual selection as predictors.
We next tested for associations between the concentrations of specific carotenoid and pteridine pigments (direct compensation). The variables in these models were the concentrations of: 1) dietary yellow-orange carotenoids (lutein/zeaxanthin, 3’-dehydrolutein, β-carotene); 2) red ketocarotenoids (astaxanthin, canthaxanthin); 3) yellow pteridines (xanthopterin); 4) red pteridines (drosopterin); and 5) colourless pteridines (pterin, 6-biopterin, isoxanthopterin, pterine-6-carboxylic acid).
Lastly, we tested whether the concentration of pigments present in skin tissue was associated with its colour. The response variables were luminance (the brightness of the colour), saturation (the intensity of the colour) or hue. For luminance and saturation, we ran two models with the following predictors: 1) total carotenoids; and 2) total pteridines. For hue we ran four models with concentrations of pigment subcategories as predictors: 1) dietary carotenoids (yellow-orange); 2) xanthopterin (yellow); 3) ketocarotenoids (red); and 4) drosopterin (red). Using these same subcategories of pigments, we tested for concentration differences between tissues with a yellow to red component (150 tissues, including browns) and those without (36 tissues that were black, grey, white or cream; total 186 samples). Lastly, we tested for associations between colour (luminance, saturation and hue) and environmental gradients of habitat productivity (PC1 and PC2).
All models were run as phylogenetically controlled mixed models in the R package MCMCglmm (Hadfield 2010). We sampled 1300 phylogenies from the posterior distribution of possible phylogenies generated in the Bayesian phylogenetic analyses. The trees employed had 28 tips, which corresponded to the 27 species sampled and two tips from the two populations of Ctenophorus pictus . For all models we used phylogeny as a random factor to control for phylogenetic relatedness between species. Given that we had several individuals per species and all individuals had more than one tissue sampled, we also included as random effects the individual and species ID (except in species-level models). We followed Ross et al. (2013) and sampled a tree at iteration t, and ran the MCMC mixed model for 1500 iterations, saving the last sample. This process was repeated for 1300 iterations (one per tree), and the first 300 runs were discarded as burn-in. Inverse Wishart priors (weakly informative) were used for the covariances and we used parameter expanded priors for the random effects. We ensured that all effective sample sizes were above 1000 and visually assessed convergence in the models using the command plot(model). We used custom code to extract a statistic that quantifies the percentage of variance explained by the fixed factors in our models (equivalent to r2). The graphs presented were generated using ggplot and the predicted fit lines were obtained from simplified mixed models (same as described above but only including significant variables). All pigment concentration variables were loge transformed to facilitate convergence, for variables with concentrations of zero we added 0.1 to all samples to avoid infinite values. We present 95% confidence bounds from the posterior distribution of the estimate based on phylogenetic mixed models run on 1000 phylogenies, where cases in which the upper and lower confidence bounds do not overlap zero indicate a significant effect.