Study organisms
A total of 147 individuals, including 51 Females, 51 males and 45 juveniles representing six species occupying different habitat types were sampled (see Additional file 1: Table S1). We distinguish male agamas from females visually, based on the bulging of the hemipenes in males [22]. Lizards were caught by hand or noose in different localities in South Africa. Agama atra samples (N = 41) were captured in the Muizenberg mountains (34° 05′ S, 18° 26′ E) and the Grootwinterhoek reserve (33° 09′ S, 19° 05′ E) and other parts of the Western Cape in March 2008 and January 2011. Agama aculeata distanti (N = 36) were sampled in Kruger National Park (23° 58′ S, 31° 31′ E) and Welgevonden Reserve (24° 12′ S, 27° 54′ E), Limpopo province, in November 2011 and March 2017 respectively. Both A. anchietae (N = 10) and A. aculeata aculeata (N = 10) were sampled in Tswalu game reserve (27° 17′ S, 22° 23′ E), Northern Cape, in January 2010, with the exception of three A. anchietae from Gobabis (22° 26′ S, 18° 57′ E) and Swakopmund (22° 15′ S, 15° 4′ E), Namibia, and one A. a. aculeata from Zwartskraal farm (33°10′S, 22°34′E), Western Cape. Agama armata (N = 11) were collected at Alicedale Farms (22° 38′ S, 30° 08′ E) and Greater Kuduland Safaris (22° 32′ S, 30° 40′ E) in January 2010 and February 2017. Lastly, Acanthocercus atricollis (N = 39) were caught in the suburban area of Mtunzini (28° 57′ S, 31° 44′ E) and Zululand Nurseries, Eshowe (28° 52′ S, 31° 28′ E), KwaZulu-Natal, in February 2017. The GPS coordinates of each lizard were recorded upon capture. Lizards were marked with a non-toxic marker (to avoid recapture), placed in cloth bags and then transferred back to the field station where they were stomach flushed. Morphological and bite force measurements were also taken. Once all the data were collected, we released the animals at the exact location where they were found.
Morphometrics
Morphological variables were measured for each individual using digital callipers (Mitutoyo; precision 0.01 mm) according to [11] for the following morphological traits (Fig. 2): snout-vent length (SVL); head width (HW) taken at the widest point of the skull; head length (HL) taken from the tip of the snout to the end of the parietal bone; head depth or height (HH) at the tallest part of the head, posterior to the orbital region; lower jaw length (LJL), taken from the snout tip to the end of the retroarticular process; jaw out-lever taken from the posterior end of the quadrate to snout tip (QT), and the distance from the back of the jugal to the tip of the snout (CT). Based on the latter three measurements, two other morphological variables were calculated: the first or closing in-lever of the jaw being the difference between QT and CT; and the second or opening in-lever, being the subtraction of QT from LJL [9, 17]. A longer in-lever for jaw closing provides a higher mechanical advantage and subsequently increases bite force for a given head size [10].
Bite force
Bite forces were measured in vivo following the method of Herrel et al. (2001) [6] using an isometric Kistler force transducer (type 9203, 500 N, Kistler Inc. Winterthur, Switzerland), connected to a Kistler charge amplifier (type 5995A) with all measurements made accurate to 0.1 N. A pair of metal bite plates was placed between the jaws of the lizard which typically results in prolonged and repetitive biting. If needed, sides of the jaw were gently tapped to provoke the lizards to bite the plates. Agamids have solid acrodont teeth implanted onto the jaw and are unlikely to suffer any damage from measuring bite force with metal plates. No audible breaking of teeth was present (contra [39]) and inspection of the teeth showed no damage. Our experience with these and other lizards is not consistent with the findings of Lappin & Jones (2014) [39] and suggest that bites on metal plates provide more accurate measures of maximal bite force. The distance between the plates and the point of application of the bite force were standardised across all animals. Bite force was recorded five times for each animal. The maximum value was then retained as the maximal bite force and used in further analyses. Although air temperatures, humidity and other environmental conditions could not be controlled in this study due to the absence of facilities in the field, we ensured that the lizards were tested at the temperatures at which they are active in the field.
Diet
Data on stomach contents for 67 individuals from four species of agamas (Agama atra, A. aculeata aculeata, A. armata and Acanthocercus atricollis) were extracted from a previously published study (Additional file 1: Table S1) [24]. Stomach contents were classified to the lowest taxonomic level possible. The food items were then blotted dry, measured and weighed using an electronic microbalance (AE100-S, Mettler Toledo GmBH, Zurich, Switzerland; ± 0.1 mg) [11].
For each prey group, we calculated the index of relative importance (IRI) [40]. This compound index indicates the importance of particular prey group and provide a balanced view based on combination of unique individual properties (numbers, mass and occurrence in diet) [17, 24]:
$${\text{IRI}} = \left( {\% {\text{N}} + \% {\text{V}}} \right) \times \% {\text{Oc}}$$
where %N is the percentage of numeric abundance, counted from the number of heads of the prey items, %V is the proportion of mass of that prey group to total prey mass and %Oc is frequency of occurrence of a certain prey group.
Statistical analyses
All morphological and performance variables were logarithmically transformed (log10) before analyses to fulfil the assumptions of normality and homoscedascity. To explore which head variables best explain variation in bite force, stepwise multiple regression analyses were conducted. Pearson correlations were further used to explore relationships between head morphology and bite force.
We grouped the species examined into three habitat groups: rock-dwelling, ground-dwelling, and arboreal. We should point out, however, that these ecological habitat groups may only apply to the populations sampled in our study and are based on our observations in the field. Following the classification of species, we tested whether habitat groups differ in size (SVL) using a univariate analysis of variance (ANOVA). If habitat groups were significantly different in size, multivariate analyses of covariance (MANCOVA) were performed to test for differences in head morphology with SVL as a covariate. We did not test for potential differences between sexes due to the low sample size for each sex per habitat group. Subsequent Tukey’s honest significant difference (HSD) post hoc tests were conducted to test for differences between pairs of habitat groups.
We then investigated whether species assigned to different habitat groups differed in absolute bite force using an ANOVA. A subsequent ANCOVA with the most significant explanatory variables from stepwise regression as covariates was performed to examine whether differences remained when correcting for head dimensions. If so then this would suggest variation in the underlying muscle architecture.
To explore the multivariate association between morphology and diet, we used two-block partial least-squares regressions (PLS) using the two.b.pls function of the geomorph package [41]. Snout vent length and all head variables were computed as the first block of variables while the IRI of the different prey groups were combined in the second block of variables. We first performed the PLS with absolute variables and then repeated the analysis with relative size, in which we used corrected morphological variables using the residuals of each trait following a regression on SVL. We additionally ran a Pearson correlation between absolute bite force and the IRI of all prey and reran the same analysis with size-corrected (residual) bite force. Finally, we tested whether habitat groups differed in their diet composition using a MANOVA and Tukey’s HSD post hoc tests on the IRI of all prey groups.
Using the same statistical analyses above, we tested for ontogenetic differences in the habitat groups. To determine whether differences found between habitat groups are linked to ecological divergence or clade membership, we conducted phylogenetic ANOVAs among the groups using a trimmed phylogeny of Leaché et al. (2014) [21].
All statistical analyses were performed using R v.3.6.2 [42].