Skull shape of a widely distributed, endangered marsupial reveals little evidence of local adaptation between fragmented populations

Abstract The biogeographic 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, behavior, and phenotypic differentiation, allowing us to ask whether populations can be distinguished and whether patterns of variation indicate adaptability to changing environmental conditions. We analyzed 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 component analysis, Procrustes ANOVA, and variation partitioning analyses. The populations harbor 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 genetic‐based conservation studies.

One of the challenges of current conservation efforts is the determination of within-species diversity. Determining population units for management plans ensures the preservation of evolutionary potential in endangered species (Crandall, Bininda-Emonds, Mace, & Wayne, 2000;Moritz, 1994) and the functioning of ecosystems (Des Roches et al., 2018). Population units are largely determined using molecular data (Allendorf, 2017)-for example, researchers of endangered species of squirrels (Finnegan, Edwards, & Rochford, 2008), jaguars (Wultsch et al., 2016), and wolves (Hindrikson et al., 2017) have all relied on genetics to identify their population diversity for conservation purposes. This genetic management has established links between diversity metrics and population fitness; however, it does not assess the phenotypic variation within a species and therefore does not discriminate this aspect of the organismic diversity within a population (Wanninger, 2015). This results in the potential for serious disjuncts between phenotypic intraspecific variation and genotype variability (Boyko et al., 2010;Le Rouzic & Carlborg, 2008;Vogt et al., 2008).
Understanding the phenotypic diversity of fragmented populations can provide valuable information to conservation management.
In particular-in analogy to the interpretation of genetic distancesmorphological differences between populations may indicate local adaptation (Colangelo et al., 2012;Meloro, 2011;Meloro, Guidarelli, Colangelo, Ciucci, & Loy, 2017). Current conservation studies of endangered taxa rarely use morphological data to determine phenotypic differentiation below the species level (Dierickx, Shultz, Sato, Hiraoka, & Edwards, 2015;Wilting et al., 2015); and most quantitatively rigorous assessment of phenotypic differentiation remains the domain of taxonomic studies (Celik et al., 2019;Meloro et al., 2017;Nicolosi & Loy, 2019;Senczuk et al., 2018;Sveegaard et al., 2015). Therefore, quantifying morphological variation within a species represents a largely untapped potential for understanding the phenotypic variation between taxonomic units and testing of hypotheses of adaptation and relatedness within a species. In addition, the morphological diagnosis of populations provides a valuable tool for management, for example, in assessing whether population units may be too morphologically divergent to be crossbred in outbreeding conservation efforts. It can also inform predictions of morphological change during future species' fragmentation events, which is a common consequence of human activity (Bennett & Saunders, 2010;Haddad et al., 2015;Hansen et al., 2013).
In this study, we use geometric morphometric analyses to provide the first population-level study of variation in the morphology of the skull of a mammal of particular conservation concern.
We focus on the endangered northern quoll (Dasyurus hallucatus: Gould, 1842), a small carnivorous marsupial (usually weighing between 300 and 900 g) (Oakwood, 1997) with well-understood morphometrics can provide a comprehensive, statistically rigorous phenomic contribution to genetic-based conservation studies.

K E Y W O R D S
conservation, Dasyurus hallucatus, geometric morphometrics, intraspecific variation, procrustes ANOVA, shape variation, variation partitioning genetic differentiation among populations (Cardoso et al., 2009;Firestone, Houlden, Sherwin, & Geffen, 2000;Hill & Ward, 2010;Hohnen et al., 2016;How, Spencer, & Schmitt, 2009;Woolley, Krajewski, & Westerman, 2015) but no information on morphological adaptation of the skull. Northern quolls appear to have had a precolonial distribution over 5,000 km across northern Australia (Braithwaite & Griffiths, 1994). They are now separated by major biogeographic breaks into four mainland populations with no apparent gene flow (Hill & Ward, 2010) and several island populations (Woinarski et al., 1999). Northern quolls are also a suitable study system for this investigation because they inhabit a wide range of habitats, ranging from rainforests to deserts (Begg, 1981;Moore et al., 2019;Oakwood, 2002). They are opportunistic foragers of small vertebrates, invertebrates, fruit, and carrion (Dunlop, Rayner, & Doherty, 2017). The species is also expected to evolve quickly because, as a semelparous species, most males die off in their first year after mating (Oakwood, Bradley, & Cockburn, 2001).
This process also has the ability to provide statistical analyses of shape variation patterns alongside visualizations of the major axes of shape variation, thus permitting a finely resolved examination of the drivers of shape variation that is not possible with conventional linear measurements.
We examine several potential factors that might impact on northern quoll skull shape variation. Given the discrete distribution of populations across Northern Australia, we expect shape differences between populations to be a main part of overall skull variation. If this variation relates to heritable adaptation to local environments, for example, through dietary differences between high and low rainfall areas (Dunlop et al., 2017), we expect differences between populations to increase with local environmental differences (such as rainfall and temperature). Alternatively, if local adaptation is limited by a constraint due to the quoll's early timing of birth, most shape variation is expected along a size gradient, rather than a population or climatic gradient, or might be unexplained. Such a result would suggest a one-to-many mapping scenario, possibly under stabilizing selection, where similarly sized individuals can fulfill multiple functions specific to their environment with one shape. It is also possible that most variation relates to the biomechanical use of the cranium due to their generalist feeding habits; this would be apparent from a strong first main axis of variation (principal component 1) which highlights areas involved in mastication-such as the zygomatic arch or incisor region-as varying most.

| Template creation
The template consists of 900 landmarks: 101 fixed landmarks, 93 curves (271 semilandmarks), and 18 surfaces (528 semilandmarks; Figure S1 and Table S1). 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 1,000 landmarks, have been demonstrated empirically to successfully capture genuine To ensure the repeatability of landmarking of the manually placed fixed landmarks, three morphologically close specimens were digitized ten times. Measurement error was much lower than interindividual variation, confirming the high repeatability of the template used in this study ( Figure S2).
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 analyzed in R version 3.6.1 (R Core (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 analyze 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.

| 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 (1,000 iterations) implemented in the "geomorph" R package (Adams & Otárola-Castillo, 2013) and then performed permutation-based (1,000 iterations) pairwise comparisons between the shape, size-corrected shape, and centroid size least squares means of each population (Collyer, Sekora, & Adams, 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 . We performed UPGMA clustering analyses based on the morphological data (Figures S5 and S6) 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 ( Figure S5) and among the specimens ( Figure S6). We also estimated the morphological disparity (Procrustes variance) of island and mainland individuals and of each population and performed pairwise comparisons among these groups.

| 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 whether 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 versus the projected regression score of shape on size (Drake & Klingenberg, 2008). We used a homogeneity of slopes test to examine whether there was a common allometric pattern in all mainland populations plus the Groote Eylandt population.

| 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 geographic 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 neighborhood 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, Maiorino, & Sansalone, 2015;Sansalone, Kotsakis, & Piras, 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.world clim.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 1,000 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.

| Morphological differences between populations
The principal component analysis reveals no visually obvious shape differentiation among populations along the two main axes of variation. The first two principal components together account for 35% of the total shape variation (PC1 = 24.58%, PC2 = 11.58%) ( Figure   S3). An example of shape changes along the first principal component axis with allometric shape changes put together with other unrelated shape changes is illustrated in Figure 1.
Despite the low differentiation of populations within the PCA, the Procrustes ANOVA indicates that at least one of the five populations differs significantly in shape with low effect (Table 1)

| Sexual dimorphism and allometry
We first confirmed that known sexual dimorphism in animal weight and skeletal measurements (Oakwood, 1997;Schmitt et al., 1989) are reflected in cranial size ( Figure 3) and shape (Table 1; Figure S4).
Size differences are significant between males and females, island and mainland populations and populations. Pairwise comparisons of size means between populations are shown in Table S2.
Females and males show no significant difference between allometric slopes (Table 1), such that small males and large females overlap on the allometric slope (see Figure 4). Most of the sex-related differences in shape are due to sex differences in size. We found statistically significant but low effects of sex-related differences in allometry-corrected shape (Table 1 and Figure S4). In other words, although there is some non-size-related variation between the sexes, small males and large females are similarly shaped according to their common size. We therefore included individuals of both sexes in our analyses of population differences.
In the full dataset, and among all variables tested, size manifests as the strongest determinant of shape variation in northern quolls, F I G U R E 1 Shape changes associated with the First Principal Component (above), allometry (middle), and precipitation (below). Spheres are the landmarks used in this study. Warmer colours represent higher landmark variation between minimum and maximum estimated configurations based on the three linear predictors. Black vectors show direction and magnitude of shape variation. Red arrows indicate anatomical zones of major landmark variation associated with allometry. Blue arrow indicates anatomical zone of major landmark variation associated with precipitation  Figure 4), meaning that the hypothesis of populations following the same allometric slope is not rejected. Allometry-corrected shape analysis also reveals the shape differences between populations; Procrustes ANOVA performed on the residual shapes of allometry revealed significant differences between populations but with low effect sizes (Table 1)

| Association of shape variation with geography, climate, and size
Significant shape differences were observed along both latitudinal and longitudinal gradients on mainland specimens (for effect sizes and significance levels, refer to Table 2). Size is significantly different only along the latitudinal gradient but not the longitudinal gradient.
Temperature and precipitation have a significant, but small, effect on shape. Size differences are only explained by precipitation and not by temperature. Other biogeophysical and environmental variables tested did not contribute significantly to size or shape variation (Table 2).
We dissected the influence of size, geography and climate (precipitation + temperature) with a variation partitioning analysis ( Figure 5).
The full model including all parameters ([a + b + c + d + e + f + g] on Finally, in accordance with our predictions, size alone [b] contributes mostly to the model by accounting for 17% of the total cranial shape variation (F 1,75 = 17.482, adjusted R 2 = 0.165, p = .001).

| D ISCUSS I ON
We expected that the biogeographic adaptive and genetic divergences between northern quolls should be evident across their 5000-kilometer longitudinal range. Our aim was to improve our F I G U R E 4 Allometry plot, centroid sizes (proxy for body size) versus shape scores obtained from the regression of shape on size (Drake & Klingenberg, 2008 understanding on whether the species represent morphological or functional conservation units, which could be used alongside population genetic approaches in the conservation management of the species. Surprisingly, however, we found little structure in northern quoll shape variation (~76% of shape variation remains unexplained) and no strong evidence that any of the populations have evolved into discrete, possibly locally adapted, morphotypes. In particular, population differences have low effect sizes and explain less variation on shape than size. This appears to be a one-to-many mapping case in which similarly sized individuals from opposite ends of the biogeographic distribution and under distinct climatic conditions with likely different functional constraints are likely to be similar in shape.
It also seems that most variation is evenly distributed within each of the populations, such that more localized western populations appear just as disparate in shape as individuals across the length of the eastern Queensland seaboard.
There is some limited support for the hypothesis that developmental constraints might reduce the adaptability of northern quoll skulls, as size explains most shape variation and the populations seem to differ in body size. Given the low amount of variation (~16%) that size explains, and the broad overlap of populations in size and shape, any such constraint is unlikely to be strong. However, there is an intriguing indication that at least some of the larger-scale evolutionary association between skull shape and size among marsupials may be visible at the within-species level. This is contrary to findings in other marsupials where feeding behaviors clearly influenced craniofacial morphology (Mitchell, Sherratt, Sansalone, et al., 2018;Weisbecker et al., 2019), and might represent one of several ways to shape morphological traits within the species.
Shape appears to have not been affected by population identity, even when size is taken into account. This might suggest a stochastic, possibly heritable, shallow divergence between populations which, however, does not appear to reflect local adaptation. These effects also demonstrate the ability of skull shape to vary independently of size based on genetic factors related to sex or population, again contradicting the developmental constraints hypothesis.
Thus, the population divergences do not appear to coincide with adaptive morphological differentiations. This provides an indication that genetic fitness benefits of outbreeding populations (Cardoso et al., 2009;Kelly & Phillips, 2019) would not risk any adverse effects due to differential local adaptation, although of course this would need to be further investigated based on nonmorphological (behavioral or physiological) traits. However, these assumptions need to be interpreted carefully because the low effect sizes of population identity on shape variation indicate a very low biological impact on cranial shape of individuals (see also Weisbecker et al., 2019).
It is possible that much of the variation in the northern quoll skull derives from a remodeling process based on individual uses of the skull. Perhaps, the best example of this is the sagittal crest, which varies widely in length among northern quoll individuals. It is common for mammals-particularly males-to display larger sagittal crests with age (Flores, Giannini, & Abdala, 2006), but this is a purely behavioral consequence of the pulling action of the temporalis muscle (Washburn, 1947). Similarly, bone remodeling during an individual's lifetime has been suggested for the zygomatic arch of mammals (Abdala & Giannini, 2000;Ravosa, 1991), and is also suspected in wombats and kangaroos Mitchell, Sherratt, Sansalone, et al., 2018;Weisbecker et al., 2019).
It is therefore possible that the emphasis of shape variation on the masticatory apparatus, found along PC1 and the allometric pattern, does not arise from heritable adaptation. Rather, they might derive from highly generalist masticatory behaviors and possibly the "mating bite" (Braithwaite & Griffiths, 1994;Oakwood, 2000) of the larger males. The nonallometric main shape variation found along PC1 (shortening of the muzzle) could be explained by differences in precipitation (proxy for dietary behaviors) (blue arrow in Figure 1).
To a limited extent due to the small effect, this size-corrected main shape variation could also be explained by nonallometric sexual differences ( Figure S4), consistent with the weak tendency of males to have shorter muzzles than females, related to higher bite forces (Wroe & Milne, 2007).
In line with the size component being the strongest determinant of shape, the "island rule" (Flannery, 2002;Foster, 1964;Lomolino, 1985) (Hohnen et al., 2016;Woolley et al., 2015), the Groote Eylandt population is most differentiated from all other mainland populations. This was also revealed in the cluster analyses, supporting the high genetic erosion observed by Cardoso et al. (2009) over 8,000 years of isolation from the mainland (Woinarski et al., 1999). The Groote Eylandt morphology reveals more compacted zygomatic arches combined with longer muzzles, which are changes associated along the allometric relationship and possibly reveal evolutionary shape changes tightly constrained by body size.
The genetic proximity of the Northern Territory and Queensland populations relative to the Pilbara and Kimberley populations (Hohnen et al., 2016;How et al., 2009;Woolley et al., 2015) appears to parallel our clustering analyses on population mean shapes. However, the unsubstantial morphological differentiation and very short branches limit the interpretation of the cluster analysis.
Aside from sexual and allometric shape variation, slight rearrangements of the zygomatic arch and the muzzle appear to provide some diffuse distinction between northern quoll populations.
This variation is consistent with well-known evolutionary adaptations to changes in mastication (Meloro, 2011;Mitchell, Sherratt, Sansalone, et al., 2018;Mitchell & Wroe, 2019;Weisbecker et al., 2019). The variation in muzzle length and zygomatic arch placement among some populations may be adaptations to a particular diet or feeding habit within each population. For instance, drier environments such as the Pilbara or Groote Eylandt show a shortening in the muzzle, when compared to wet environments, such as Northern Territory or Queensland. Shorter faces are known to imply greater bite forces, and thus might indicate the mastication of tougher foods (Wroe & Milne, 2007). However, although precipitation is a main predictor of quoll diets (Fancourt et al., 2015), this variable explains little variation in shape. Whether the shortening of the muzzle is truly an adaptive effect requires further research due to the lack of association between climatic factors and shape, and the abovementioned small effect sizes and extensive overlap in shape between populations. This provides a good starting point for discussion and future investigations to identify whether the differences among populations are or not a decisive factor for the individuals' survival, and rather originate from potentially nonadaptive factors such as genetic drift or individual variation in feeding (Weisbecker et al., 2019) between populations.
An apparent limitation of our variation partitioning analyses is the high levels of unexplained variation. While these are reasonably common in geometric morphometrics, they are also generally not expected (Cardini, Filho, Polly, & Elton, 2010;Cardini, Jansson, & Elton, 2007;Hendges, Bubadué, & Cáceres, 2016;Pandolfi et al., 2015). Where unexplained variation is large, it is commonly presumed that other, unknown variables are responsible for this variation. This may also be the case in our analysis; alternatively, the effect of factors not measured for this study, such as genetic, physiological, or behavioral traits, could contribute to this proportion of unexplained variation.

ACK N OWLED G M ENTS
We thank all museum collections who granted us access to scan We would finally like to thank Ariel Marcy, Jessica Ivory-Church, and Andrew Baker for bringing up landmarking insights.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

O PE N R E S E A RCH BA D G E S
This article has been awarded Open Materials, Open Data Badges.
All materials and data are publicly accessible via the Open Science