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.