Phenotypic variation in functional traits underpins responses to environmental gradients, influencing thermoregulation, energy balance, and long-term persistence under climate extremes. Climate change is altering these gradients globally, yet in species that have disappeared from much of their range, adaptive phenotypes may have also been lost, potentially limiting population viability following reintroduction to different climates. The greater bilby (Macrotis lagotis) is a threatened Australian marsupial that has undergone severe range contraction since European settlement and has been the focus of subsequent conservation translocations. To test hypotheses of climate-associated morphological divergence, we quantified spatial variation across its historical distribution and in reintroduced populations in traits that are key determinants of heat exchange and thermal buffering: skull size (as a proxy for body size), ear length, and fur morphology (dorsal and ventral hair length, depth, and width). Data were collected from museum specimens held across nine natural history collections in Australia, the United Kingdom, and the United States, and from live-captured individuals at two conservation reserves in South Australia (Arid Recovery Reserve and Venus Bay Conservation Park). Museum specimen records are accompanied by specimen metadata and CHELSA v2.1 climate variables extracted for each locality (Karger et al., 2017, 2021). The dataset comprises seven CSV files and one R Markdown code file. It is suitable for reuse in studies of ecogeographic variation, intraspecific trait–climate relationships, thermoregulation, and conservation translocation outcomes in mammals. There are no legal or ethical restrictions on reuse; no data on human subjects are included.

Museum specimens

We measured skull dimensions and ear lengths of bilby specimens housed in the American Museum of Natural History, Australian Museum, Australian National Wildlife Collection, British Natural History Museum, Melbourne Museum, Museum and Art Gallery of the Northern Territory, Queensland Museum, South Australian Museum, and Western Australian Museum. Cranial measurements were taken using digital callipers (resolution = 0.01 mm) by the same observer following methods outlined by Aplin et al. (2010). For this investigation, we measured skull length (occipitonasal length; greatest length of the skull) and skull width (zygomatic width; maximum width across cranium) as proxies for body size (Aplin et al., 2010; Travouillon, 2016; Umbrello, 2018). From this dataset, we identified 106 bilby specimens with location data allowing climatic assessment, of which 98 specimens were adults from naturally occurring populations (Figure 1). Where possible, we measured ear length from the crown to the tip of the ear using digital callipers to assess whether extremity size showed evidence of divergence consistent with Allen’s rule under differing climatic pressures (Allen, 1877).

We measured fur characteristics of preserved bilby specimens at the Australian Museum, South Australian Museum, and Western Australian Museum. To account for variance in fur properties across the body, we took measurements from six locations on each specimen (Figure 1): three dorsal locations (neck, rump, and outer thigh) and three ventral locations (sternum, belly, inner thigh). All measurements included in the analysis were taken by the same observer to ensure methodological consistency. To minimise potential measurement bias, specimens were measured prior to consultation of specimen metadata, and so the observer was blind to locality and associated climatic information during data collection. Fur length was measured by combing the hair in the opposing direction to its natural orientation and measuring from the skin surface to the tip of the hair. Depth was measured when the fur was combed downwards and recorded as the perpendicular distance between the skin and the outer surface of the fur layer. Measurements were taken by placing the end of the digital callipers against the skin and measuring to the target position to the nearest 0.01 mm.

Fur width was measured by taking opportunistic samples of hairs dislodged during the combing process of museum specimens. At each fur region per individual, three hairs were measured at the midpoint under an Olympus SZX16 microscope using Olympus cellSens Entry v.3.2 imaging software. These three measurements were then averaged to get a measurement for that region. Guard hairs, the longer and sparser component of the coat layer, could easily be distinguished under the microscope and were not measured due to high intra-region variability on the same animal (Appendix S1). All fur measures were taken before reading specimen metadata to help avoid measurement biases. For each specimen, we also measured pes length (length of foot measured from tip of longest toe, excluding nail, to base of heel) using digital callipers (resolution = 0.01 mm) in order to account for differences in body size in skins where the head was no longer intact. Body mass of 45 live bilbies measured independently of this investigation was strongly correlated with both pes length (r = 0.83) and head length (r = 0.92), and pes length was itself strongly correlated with head length (r = 0.82), supporting the use of these metrics as proxies for overall body size within populations (see Appendix S2 for details).

Translocated populations

To assess for potential divergence in translocated populations, we also measured the body size, ears, and fur of live bilbies from two South Australian reserves: the xeric Arid Recovery Reserve and mesic Venus Bay Conservation Park. The Arid Recovery reserve is in a region classified as desert, with a mean annual rainfall of 151 mm. The mean number of days surpassing 40 ℃ is 26.6, and 11.1 days below 0 ℃ annually. Venus Bay has a temperate climate, receiving 425 mm of annual rainfall, and experiences fewer climate extremes than Arid Recovery, with just 1.8 days surpassing 40 ℃ and 0.1 days below 0 ℃ annually. These sites, located ~380 km from each other in different bioclimatic regions, make an ideal comparison to assess phenotypic divergence, having been translocated from the same captive population into different climatic regions.

Bilbies from both populations were reintroduced from a captive-bred population at Monarto Safari Park in South Australia, with the breeding program starting in 1994 from arid zone Northern Territory and Western Australian populations (Appendix S3). Bilbies were reintroduced to Arid Recovery in 2000 and occupy the entire 123 km² fenced reserve. Bilbies were originally translocated to Venus Bay between 2001 and 2005 and were later supplemented with individuals from Thistle Island, South Australia, another temperate site with original stock from Monarto established in 1998 and within 230 km of Venus Bay. These sites were selected as they represent two of the longest-running translocation sites (~25 years) for the species and have not received reinforcement from climatically distinct source populations, allowing potential trait divergence to be assessed without confounding source effects.

Bilbies were caught at both sites using a combination of cage trapping, netting from a vehicle, and burrow trapping (McGregor & Moseby, 2014). The same measurements of head, ear, pes length, fur length, and fur depth were taken immediately upon capture, with the bilby remaining within a catch bag, where only the focal measurement location was unobscured to minimize stress to the animal. Fur samples to assess hair width were taken by shaving a small section of fur from the target area, and fur width was measured using the same method previously described. Cage traps, baited with a peanut butter and oat mixture, were checked before sunrise, and bilbies were released at the capture site. We measured a total of 14 bilbies at Arid Recovery and five bilbies from Venus Bay. Data were collected from both sites during the Austral autumn, May 2024 and April 2025, to reduce any seasonal impacts on fur.

Climate data

We aimed to characterise spatial variation in body size and fur properties to test hypotheses of heat dissipation, heat conservation, and productivity across the full geographic extent of the bilby collections. In arid Australia, weather station records were initially sparse and only became more spatially representative over the course of the twentieth century (Jeffrey et al., 2001), resulting in substantial spatial artefacts in interpolated climate surfaces, particularly in remote regions. To minimise these artefacts and better capture spatial patterns in climate, we extracted data for each specimen location from CHELSA v2.1 climate rasters (Karger et al., 2017, 2021), representing the thirty-year average between 1981–2010 at a native spatial resolution of 30 arc-seconds. Although this approach may introduce some error for older specimens, we judged it preferable to relying on historical interpolations with known spatial biases in sparsely sampled regions, such as inland Australia. To test hypotheses of heat dissipation, heat conservation, and productivity as drivers of variation in fur and skull traits, we used four BIOCLIM variables extracted from CHELSA v2.1 as predictors: maximum temperature of the hottest month (BIO5), minimum temperature of the coolest month (BIO6), annual precipitation (BIO12), and precipitation seasonality (BIO15). For clarity, these are hereafter referred to as MaxTemp, MinTemp, Rainfall, and Seasonality, respectively.

Data preparation

Data preparation and statistical analyses were performed in R v4.4.1 (R Core Team, 2025). The fur dataset contained 41.5% missing measurements, primarily concentrated in ventral traits. Missing values were imputed using predictive mean matching implemented in the mice package (Van Buuren & Groothuis-Oudshoorn, 2011), with 50 imputations run for 50 iterations. To ensure inference was robust to imputation, all analyses were repeated across all iterations, model stability was assessed, and complete-case sensitivity analyses were conducted (see Appendix S4 for full details).

We retained only adult individuals from natural populations for model fitting. After cleaning, we assessed skull size in 98 individuals, ear length in 54 individuals, and fur parameters in 75 individuals from museum collections (Appendix S5). For fur measurements, we averaged length, depth, and width between dorsal and ventral sides, leaving six fur traits as response variables in the final dataset. We averaged climatic predictors over a 10 km buffer around each specimen location to better represent local-scale environmental conditions and account for varying levels of spatial confidence. All predictor variables were scaled prior to modelling to allow for direct comparison of effect sizes.

We assessed collinearity among climate predictors using both Pearson’s correlation and variance inflation factors (VIFs). Correlations among variables were low to moderate, with the highest correlation observed between MinTemp and Seasonality (r = 0.72) in the fur dataset. As this combination was the only one with a VIF score greater than 4, we did not include any candidate models with both variables in our fur trait model set. VIF scores for all other predictor combinations were between 1 and 4, suggesting moderate correlation but indicating that collinearity was not severe enough to exclude any predictors (Zuur et al., 2010; Dormann et al., 2013).

Statistical analyses

We fitted univariate linear models for skull length and skull width using the lm()** function in R, including combinations of the four climate variables and sex as fixed effects. These response variables were analysed separately because they represent distinct axes of skull morphology, allowing us to distinguish between concurrent changes in skull size (both dimensions responding similarly) and changes in skull shape (change in one measurement or divergent change in both). We ran a similar model set with ear length as a response variable, with the inclusion of pes length as a fixed effect to account for body size. Models were ranked using AICc, and models within 2 ΔAICc units of the top-ranked model were considered competitive (Appendix S6). Where multiple competitive models were identified, we selected the simplest model for parsimony (Burnham & Anderson, 2002; Arnold, 2010). To assess model fit, we simulated and visually inspected scaled residuals using the DHARMa package (Hartig, 2017). To investigate potential overfitting, we performed five-fold cross-validation on all high-ranking models and assessed root mean square error, mean absolute error, and the squared correlation coefficient between predicted and observed values using the Metrics package (Hamner et al., 2018; Jung, 2018; Yates et al., 2023). Marginal means and predicted values were extracted using the emmeans package (Lenth, 2024), and we report Wald t-tests for all linear models.

To measure fur response to climate variables, we fitted multivariate linear models (MANOVAs) using the glmmTMB package for the six fur variables (Brooks et al., 2017). Fur traits were analysed in a multivariate framework as measurements were taken across multiple body regions and represent functionally integrated components of insulation; modelling them jointly allows us to account for covariance among traits and assess coordinated responses to climate. Museum specimens and live captures from translocation sites were analysed separately: museum skins provide a historical baseline across climatic gradients, whereas live bilby measurements reflect contemporary translocated populations that have been moved between climate regions in recent generations, potentially obscuring patterns observed in museum cohorts. For museum specimens, candidate models included the four climate predictors, pes length, and sex, either individually or in combination with each other. Collection year (as years since specimen collection), season (as a factor), or month (fitted with cyclical sine and cosine terms) were also tested in candidate models to assess the possibility of specimen degradation and seasonal moulting.

Models were ranked using multivariate AICc, and those within 2 ΔAICc of the top-ranked model were considered competitive. Once again, where multiple competitive models were identified, the most parsimonious model was selected. Although multiple individuals were sampled from some localities, the majority of sites were represented by just one specimen (median per site = 1; mean = 1.5), making the estimation of a site-level random effect unreliable and producing convergence issues. Climate predictors were expected to capture site-level environmental variation, and inclusion of site as a random effect led to unstable variance estimates. We therefore utilised fixed-effects models and acknowledge potential residual site-level non-independence as a limitation.

We tested whether fur morphology differed between Arid Recovery (n = 14) and Venus Bay (n = 5) populations, again analysing the six fur traits in a multivariate framework. Predictors of interest were translocation site, sex, and pes length, which were included as fixed effects in candidate models and models were selected using the same method described above. To evaluate the stability of the multivariate framework given the limited sample size of the translocated dataset, we implemented a nonparametric bootstrap procedure (500 iterations) to assess model-selection consistency and sensitivity to sampling variation (Davison & Hinkley, 1997; Eck, 2018). We report multivariate test statistics (Wilks’ λ) for terms retained in the top-ranked models. To aid interpretation, we also extracted per-trait coefficients and marginal R² from the best-supported multivariate models. To assess site-level differences in body size and extremity traits, we additionally fitted linear models for head length, ear length, and pes length with translocation site and sex as fixed effects, with ear length models including head length as a covariate to account for body size.

Finally, we tested whether insulation in reintroduced populations was converging towards historical measurements for each release site. Site-specific predictions were generated from the best-supported museum model (using MaxTemp as the focal climate predictor). Predictions were evaluated at each site’s MaxTemp (35.1 and 23.9 ℃ at Arid Recovery and Venus Bay, respectively). Observed site means for translocated bilbies were obtained as estimated marginal means from trait-wise linear models using the emmeans package. For each trait and site, Δ was computed on the original scale, and its uncertainty was estimated as the square root of the summed squared standard errors of the predicted and observed means. Two-sided Z-tests were used to assess Δ = 0, with Benjamini–Hochberg FDR adjustments made to p-values across traits and sites.

References

Allen, J. A. (1877). The influence of physical conditions in the genesis of species. Radical Review, 1, 108-140.

Aplin, K. P., Helgen, K. M., & Lunde, D. P. (2010). A review of Peroryctes broadbenti, the giant bandicoot of Papua New Guinea. American Museum Novitates, 2010(3696), 1-41.

Arnold, T. W. (2010). Uninformative parameters and model selection using Akaike's Information Criterion. The Journal of Wildlife Management, 74(6), 1175-1178.

Brooks, M. E., Kristensen, K., Van Benthem, K. J., Magnusson, A., Berg, C. W., Nielsen, A., Skaug, H. J., Machler, M., & Bolker, B. M. (2017). glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal, 9(2), 378-400.

Burnham, K. P., & Anderson, D. R. (2002). Practical use of the information-theoretic approach. In Model Selection and Inference (2nd ed.). New York, NY: Springer.

Davison, A. C., & Hinkley, D. V. (1997). Bootstrap Methods and their Application. Cambridge University Press.

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G., Gruber, B., Lafourcade, B., & Leitão, P. J. (2013). Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography, 36(1), 27-46.

Eck, D. J. (2018). Bootstrapping for multivariate linear regression models. Statistics & Probability Letters, 134, 141-149.

Hamner, B., Frasco, M., & LeDell, E. (2018). Package 'Metrics'. R Foundation for Statistical Computing.

Hartig, F. (2017). Package 'DHARMa'. R Package, 531, 532.

Jeffrey, S. J., Carter, J. O., Moodie, K. B., & Beswick, A. R. (2001). Using spatial interpolation to construct a comprehensive archive of Australian climate data. Environmental Modelling & Software, 16(4), 309-330.

Jung, Y. (2018). Multiple predicting K-fold cross-validation for model selection. Journal of Nonparametric Statistics, 30(1), 197-215.

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., & Kessler, M. (2017). Climatologies at high resolution for the earth's land surface areas. Scientific Data, 4(1), 1-20.

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., & Kessler, M. (2021). Climatologies at high resolution for the earth's land surface areas [Raster dataset]. EnviDat. https://doi.org/10.16904/envidat.228

Lenth, R. (2024). emmeans: estimate marginal means, aka least-squares means: version 1.10.3. https://cran.r-project.org/web/packages/emmeans/index.html

McGregor, H., & Moseby, K. (2014). Improved technique for capturing the greater bilby (Macrotis lagotis) using burrow cage traps. Australian Mammalogy, 36(2), 259-260.

R Core Team. (2025). R: A Language and Environment for Statistical Computing (Version 4.4.1): R Foundation for Statistical Computing. https://www.r-project.org/

Travouillon, K. J. (2016). Investigating dental variation in Perameles nasuta Geoffroy, 1804, with morphological evidence to raise P. nasuta pallescens Thomas, 1923 to species rank. Zootaxa, 4114(4), 351-392.

Umbrello, L. (2018). Evolution and diversification of dasyurid marsupials of the Australian arid zone. (PhD), University of Western Australia.

Van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate imputation by chained equations in R. Journal of Statistical Software, 45, 1-67.

Yates, L. A., Aandahl, Z., Richards, S. A., & Brook, B. W. (2023). Cross validation for model selection: a review with examples from ecology. Ecological Monographs, 93(1), e1557.

Zuur, A. F., Ieno, E. N., & Elphick, C. S. (2010). A protocol for data exploration to avoid common statistical problems. Methods in Ecology and Evolution, 1(1), 3-14.