DOI: 10.5061/dryad.t4b8gthxv
title: Skull shape of a widely-distributed, endangered marsupial reveals little evidence of local adaptation between fragmented populations

description: 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.
description: <i>Data
 collection (3D acquisition)</i> We
 reconstructed in silico 101 crania of adult individuals of
 <i>Dasyurus hallucatus </i>– 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)<i>.
 </i>Adult status was determined through incisor wear (Oakwood, 2000)
 and P3 eruption (Woolley <i>et al</i>., 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. <i>Template creation </i> 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 <i>et al</i>., 2013; Dumont <i>et al</i>., 2016; Goswami <i>et al</i>., 2019; Gunz & Mitteroecker, 2013; Segall <i>et al</i>., 2016; Watanabe <i>et al</i>., 2019; Weisbecker <i>et al</i>., 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). <i>Landmarking and sliding</i> Each 3D model was landmarked in Viewbox version 4.0 (dHAL software, Kifissia, Greece; http://www.dhal.com; Polychronis <i>et al</i>., 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. <i>Analysis</i> Raw coordinate data were analysed in R version 3.6.1 (R Core Team, 2019) with the ‘geomorph’<i> </i>(version 3.1.2) (Adams & Otárola‐Castillo, 2013) and the ‘Morpho’<i> </i>(version 2.7) (Schlager <i>et al</i>., 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 <i>shape</i> (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. <i>Morphological differences between populations</i> 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 <i>et al</i>., 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 <i>et al</i>., 2016; How <i>et al</i>., 2009; Woolley <i>et al</i>., 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. <i>Sexual Dimorphism and Allometry</i> 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.  <i>Association of shape variation with geography, climate and size </i> 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 <i>varpart</i> function in the <i>vegan</i> package for R version 2.5-6 (Oksanen <i>et al</i>., 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 <i>et al</i>., 2015; Sansalone <i>et al</i>., 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.