Analysis
Raw coordinate data were analysed in R version 3.6.1 (R Core Team, 2019) with the ‘geomorph’ (version 3.1.2) (Adams & Otárola‐Castillo, 2013) and the ‘Morpho’ (version 2.7) (Schlager et al., 2019) packages. A Generalized Procrustes Analysis (GPA) was performed on the raw landmarks to translate, rotate and scale specimens to the same size. This allowed us to extract the size component as Centroid Size (Rohlf & Slice, 1990) and to analyse shape (form minus size) (Kendall, 1989). This GPA step was used for all specimens as well as subsets (e.g., if only specimens of known sex or mainland-only specimens were considered for corresponding analyses). We included specimens from Groote Eylandt as a separate population for all our analyses, but did not include the specimens from four other island populations for population analyses. This is because of the low sample size for these island populations, as well as the uncertain genetic history (Hohnen et al., 2016; How et al., 2009; Woolley et al., 2015) and possible genetic erosion (Cardoso et al., 2009). We did, however, test whether there were differences in shape variation between island and mainland populations, which would occur if divergent selection on the different islands shaped each population differently.
Morphological differences between populations
In order to explore the shape variation of our dataset, we conducted a Principal Component Analysis (PCA) on the landmark coordinates. This method reduces the large dimensionality of the dataset – due to the large amount of variables (i.e., landmark coordinates) – by tracing orthogonal axes along the main variance-covariance of the data, with the by-product being that the first axis (i.e., principal component) represent most of the shape variation. If population is a determining factor of shape variation in the sample, one of the main Principal Components (PCs) is therefore expected to separate specimens according to populations.
We also assessed for shape differences between populations with Procrustes ANOVAs with permutation tests, and then performed permutation-based pairwise comparisons between the shape and centroid size least squares means of each population (Collyer et al., 2015). Heat plot visualizations of mean comparisons between populations were used to understand where the main shape differences were located. We performed all heat plot visualizations with the landvR package (Guillerme & Weisbecker, 2019).
Sexual Dimorphism and Allometry
To assess the degree to which shape variation in the sample was determined by sex differences, we computed the interaction term of Size and Sex on shape to evaluate if both sexes had common allometric relationships. In addition, we corrected for allometric shape differences between sexes by extracting the residuals of allometry of each specimen and adding them up to the consensus shape obtained from the GPA. This allowed us to make heat plots of shape-only sexual dimorphism.
In order to evaluate the influence of size on shape (allometry) in our dataset, we performed a Procrustes ANOVA to quantify the amount of shape variation influenced by the centroid size. We then plotted the Centroid Size vs the projected regression score of shape on size (Drake & Klingenberg, 2008). We used a Homogeneity of Slopes Test to examine if there was a common allometric pattern in all mainland populations plus the Groote Eylandt population.
Association of shape variation with geography, climate and size
To assess the factors influencing cranial shape on the four mainland populations, such as geographical distance among individuals or bioclimatic data, we performed a variation partitioning analysis. For this, we used the varpart function in the vegan R package version 2.5-6 (Oksanen et al., 2018). Latitude and longitude coordinates of each locality corresponding to each specimen were transformed using a principal coordinates of neighbourhood matrix (PCNM) (Borcard & Legendre, 2002) to avoid spatial autocorrelation. The PCNM method presents several advantages. It produces orthogonal (linearly independent) spatial variables over a much wider range of spatial scales (Pandolfi et al., 2015; Sansalone et al., 2015). Bioclimatic data for all specimens (Annual Mean Temperature, BIO1, and Annual Precipitation, BIO12, of each location) were obtained from WORLDCLIM (v. 2.0) (www.worldclim.org/bioclim) and were included in the final model as a hypothesized influence of a climatic factor on shape. Finally, we performed a redundancy analysis (RDA) on the full model, which includes the three factors (Size, Spatial data and Climatic data) that we hypothesized explained the variation on cranial shape in our study system.