The biogeographical distribution of diversity among populations of threatened mammalian species is generally investigated using population genetics. However, intraspecific phenotypic diversity is rarely assessed beyond taxonomy-focused linear measurements or qualitative descriptions. Here, we use a technique widely used in the evolutionary sciences – geometric morphometrics – to characterize shape diversity in the skull of an endangered marsupial, the northern quoll, across its 5,000 km distribution range along Northern Australia. Skull shape is a proxy for feeding, behaviour, and phenotypic differentiation, allowing us to ask if populations can be distinguished and if patterns of variation indicate adaptability to changing environmental conditions. We analysed skull shape in 101 individuals across four mainland populations and several islands. We assessed the contribution of population, size, sex, rainfall, temperature, and geography to skull shape variation using Principal Components Analysis, Procrustes ANOVA, and variation partitioning analyses. The populations harbour similar amounts of broadly overlapping skull shape variation, with relatively low geographic effects. Size predicted skull shape best, coinciding with braincase size variation and differences in zygomatic arches. Size-adjusted differences in populations explained less variation with far smaller effect sizes, relating to changes in the insertion areas of masticatory muscles, as well as the upper muzzle and incisor region. Climatic and geographic variables contributed little. Strikingly, the vast majority of shape variation - 76% - remained unexplained. Our results suggest a uniform intraspecific scope for shape variation, possibly due to allometric constraints or phenotypic plasticity beyond the relatively strong allometric effect. The lack of local adaptation indicates that cross-breeding between populations will not reduce local morphological skull (and probably general musculoskeletal) adaptation because none exists. However, the potential for heritable morphological variation (e.g. specialization to local diets) seems exceedingly limited. We conclude that 3D geometric morphometrics can provide a comprehensive, statistically rigorous phenomic contribution to genetics-based conservation studies.

Data collection (3D acquisition)

We reconstructed in silico 101 crania of adult individuals of Dasyurus hallucatus – including males and females from four mainland populations: Queensland (n = 35), Northern Territory (n = 31), Pilbara (n = 15), and Kimberley (n = 8); and island populations: Groote Eylandt (n = 7) and other small islands (n = 5). Adult status was determined through incisor wear (Oakwood, 2000) and P3 eruption (Woolley et al., 2013). We 3D-scanned most of the specimens from museum collections (Queensland Museum, Australian Museum, Western Australian Museum, Australian National Wildlife Collection and American Museum of Natural History) with a GoMeasure 3D HDI109 blue light surface scanner (LMI Technologies Inc., Vancouver, Canada). Each cranium was placed in 3 different orientations on a motorised rotary table that turned every 45 degrees (8 rotations per round). The 24 resulting 3D images (8 rotations x 3 orientations) were then meshed together with the scanner’s software (FlexScan3D 3.3) to export a complete 3D image of each skull. This file was then treated for post-processing cleaning (so as to not affect the biological shape of the structure), mesh decimation (to facilitate computation) and mesh reformatting (as “.ply” files need to be in binary format for subsequent importations of the mesh into R). Several photographs of each specimen were also taken to help identify landmarks by distinguishing biological structures from 3D artefacts in the landmarking process. Seven fully fleshed specimens from the Groote Eylandt population were CT-scanned at the Centre for Advanced Imaging at The University of Queensland in a micro CT-scanner. In order to obtain the 3D model, segmentation of the DICOM grayscale images provided by the micro CT-scan was performed with Mimics Research version 20.0. All 3D models can be accessed through MorphoSource. The University of Queensland animal ethics committee (permit number SBS/009/16/ARC) and the Northern Territory Parks and Wildlife Commission (permit number 58566) approved the research methods and the collection of the Groote Eylandt specimens.

Template creation

The template consists of 900 landmarks: 101 fixed landmarks, 93 curves (271 semilandmarks), and 18 surfaces (528 semilandmarks) (Supplementary Figure 1 and Supplementary Table 1). The number of semilandmarks on curves or surfaces was determined by the complexity (inflection points) of the curves or area covered. High density landmark and semilandmark configurations, ranging from 800 to more than 1000 landmarks, have been demonstrated empirically to successfully capture genuine biological signal (Cornette et al., 2013; Dumont et al., 2016; Goswami et al., 2019; Gunz & Mitteroecker, 2013; Segall et al., 2016; Watanabe et al., 2019; Weisbecker et al., 2019)

To ensure the repeatability of landmarking of the manually placed fixed landmarks, three morphologically close specimens were digitised ten times. Measurement error was much lower than inter-individual variation, confirming the high repeatability of the template used in this study (Supplementary Figure 2).

Landmarking and sliding

Each 3D model was landmarked in Viewbox version 4.0 (dHAL software, Kifissia, Greece; http://www.dhal.com; Polychronis et al., 2013). One operator (first author) manually placed the fixed landmarks and curves. Curve semilandmarks were placed equidistantly and then were allowed to slide along their respective curves. Surface semilandmarks were placed following a thin-plate spline interpolation between the template and each specimen, followed by a projection to the surface of the 3D-model and the sliding procedure. Sliding was performed by minimizing bending energy.

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).  Despite its small sample size, we included specimens from Groote Eylandt (n = 7) as a separate population for all our analyses, taking into consideration this caveat in our interpretation of the results. We did not include the specimens from four other island populations for population analyses (total n = 5). We did, however, test whether there were differences in shape variation among all island (including Groote) and mainland individuals, which would occur if divergent selection on the different islands shaped each population differently.

1. Morphological differences between populations

In order to explore the variation of shape in 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 number of variables (i.e., landmark coordinates) – by tracing orthogonal axes along the main variance-covariance axis of the data, with the result being that the first axes (i.e., principal components) represent most of the shape variation. If the identity of a population determines 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 using the permutational procedure (1000 iterations) implemented in the ‘geomorph’ R package (Adams & Otárola‐Castillo, 2013), and then performed permutation-based (1000 iterations) pairwise comparisons between the shape, size-corrected shape and centroid size least squares means of each population (Collyer et al., 2015). We adjusted p-values of pairwise comparisons with the Bonferroni method.  Heat plot visualizations of mean comparisons between populations were used to understand where the main shape differences were located. We executed all heat plot visualizations with the ‘landvR’ package (Guillerme & Weisbecker, 2019). We performed UPGMA clustering analyses based on the morphological data (Supplementary Figures 5 and 6) to contrast with the genetic-based phylogenies of the populations (Hohnen et al., 2016; How et al., 2009; Woolley et al., 2015). We reconstructed the phenetic relationships based on the Euclidean distances among the population mean shapes (Supplementary Figure 5) and among the specimens (Supplementary Figure 6). We also estimated the morphological disparity (Procrustes variance) of island and mainland individuals and of each population and performed pairwise comparisons among these groups.

2. 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 to the consensus shape obtained from the GPA. This allowed us to make heat plots of sexual dimorphism of shape and size-corrected shape.

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. 

3. Association of shape variation with geography, climate and size

We performed a partitioning analysis of cranial shape variation test for factors that may influence cranial shape, such as geographical distance among individuals or bioclimatic variables for the four mainland populations. For this, we used the varpart function in the vegan package for R 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). We retained the PCNM axes showing positive eigenvalues, then we checked for significance for the selected axes and these (n = 10) were the ones included in the analyses. Environmental and biogeophysical data for all specimens (elevation, distance to permanent water, primary productivity and vegetation type; annual mean temperature and annual precipitation) were obtained from the Atlas of Living Australia (www.ala.org.au) and WORLDCLIM (v. 2.0) (www.worldclim.org/bioclim), respectively. Those variables that had a clear effect on size and/or shape variation were included in the final model. Finally, we performed a redundancy analysis (RDA) on the partial and full models with 1000 permutations, which includes the three factors (Size, Spatial data and Climatic data) that we hypothesized to explain the variation on cranial shape in our study system.