DOI: 10.5061/dryad.4j0zpc8nf
title: Dietary breadth in kangaroos facilitated resilience to Quaternary climatic variations

description: Identifying what drove the late Pleistocene megafaunal extinctions on the
 continents remains one of the most contested topics in historical science.
 This is especially so in Australia, which lost 90% of its large species by
 40,000 years ago, more than half of them kangaroos. Determining causation
 has been obstructed by a poor understanding of their ecology. Using Dental
 Microwear Texture Analysis, we show that most members of Australia’s
 richest Pleistocene kangaroo assemblage had diets that were much more
 generalized than their craniodental anatomy implies. Mixed feeding across
 most kangaroos pinpoints dietary flexibility as a key behavioral
 adaptation to climate-driven fluctuations in vegetation structure,
 dispelling the likelihood that late Pleistocene climatic variation was the
 sole or primary driver of their disappearance.
description: <em>Modern samples</em>
   Modern specimens used in this study
 are housed in the Australian Museum, Sydney (prefix AM M); Museum and Art
 Gallery of the Northern Territory, Alice Springs and Darwin (CAM, U);
 Museum Victoria, Melbourne (MV C, DTC); Flinders University Research
 Collection, Adelaide (FUR); Queensland Museum, Brisbane (QM A, J, JM),
 South Australian Museum, Adelaide (SAMA M); Western Australian Museum,
 Perth (WAM M); American Museum of Natural History, New York (AMNH); and
 Papua New Guinea National Museum and Art Gallery, Port Moresby (PNG MR).
 Sixteen extant species were sampled to capture a broad dietary spectrum,
 with an additional five species of *Dendrolagus *from New Guinea analyzed
 as a single unit, because no single species had adequate available samples
 (Table S1, Data S1).  
 <em>Fossil samples</em>
   Paleontological specimens are housed
 in the South Australian Museum, Adelaide (SAMA P, FU). The sample
 originated from excavations of the Main Fossil Chamber deposit of Victoria
 Fossil Cave, Naracoorte, South Australia, which were led by Rod Wells
 (Flinders University) through the 1970s–1990s. The sequence consists of at
 least eight superposed infill sedimentary units (*19*, *42*). A flowstone
 capping the sequence provides a minimum age of ~213 ka (*22*). Using
 existing stratigraphic designations (*42*), unit 8 and upper unit 7 (depth
 bin 7E, see below) are dated to ~220 ka and ~226 ka, respectively
 (*21*).   Here we assess the
 diets of the 14 best-represented VFC macropodid species (Table S1),
 excluding potoroines, because DMTA for modern, largely fungivorous
 potoroines has yet to be investigated. Insufficient samples could be
 attained for *Lagorchestes leporides *and *Procoptodon goliah*. Recently,
 *Protemnodon brehus *and *Prot. roechus *have been relegated to *nomina
 dubia*, and effectively replaced by *Prot. mamkurra *and* Prot. viator
 *(*12*). Unfortunately, the cheek dentition of these two species cannot be
 distinguished, so samples are here referred to as *Prot. mamkurra *due to
 its greater abundance in the VFC deposit based on specimens that can be
 identified to species level (I. Kerr, pers. comm. 21/3/2024). However,
 there remains a possibility that *Prot. viator *is present in the
 sample.   <em>Data
 acquisition</em>  
 Specimens were cleaned and cast using standard procedures
 (*43*), which have been shown to have high fidelity in replicating
 microwear surfaces (*44*). Casts were scanned using a Sensofar Plμ NEOX
 confocal microscope “Bruce” at Flinders University, (100X ELWD objective,
 neural aperture 0.80, blue light 460nm, spatial sampling 0.17µm, step
 height &lt;4nm). Four scans were taken of each surface, (2x2, overlap
 15%) stitched to a total scan size of 242 × 181 µm2.
   To minimize sampling errors, the ‘Soft
 Filter-edited’ data processing template (*45*) was used on all scans in
 SensoMAP 7.1.2.7288 (Digital Surf). Some samples utilized here have been
 previously published upon in methodological papers (*45*, *46*), though
 scanning parameters were identical and all data were recollected here to
 ensure comparability. Photosimulations of representative scans can be
 found in figure S1, and PDFs of all scanned surfaces can be found at the
 Dryad page for this article PDFs_of_Specimens_Scanned.zip.   Intraspecific factors were scored for each specimen (*46*, see table S2). In addition, we also considered the factor *islands*, i.e., whether or not a specimen originated from an island (*47*). Due to ease of scanning, we priorirized larger *tooth facets* (1 and 6 for upper molars, 4 and 9 for lower molars), and *macrowear stages* 2–3, but not exclusively, which allowed effective modeling of interspecific differences where low numbers of specimens were available (*46*). It should also be noted that (*46*) indicates that facet numbering followed (*48*), when in fact it used that of (*49*).   Data were collected using six Scale-Sensitive Fractal Analysis variables (*Asfc, HAsfc9, HAsfc81, Smc, Lsfc, epLsar and*, NepLsar; *43*), 25 International Organization for Standardization (ISO) variables (*S10z, Sa, Sal, Sda, Sdax, Sddx, Sdq, Sdv, Sdvx, Shh, Shhq, Shhx, Shv, Sk, Sku, Smc, Smrk1, Spc, Spd, Spk, Ssk, Svd, Svk, Vmc, *and *Vvv; *(*50*), and four variables relating to Motif, Isotropy and Furrow (*Madf, Metf, Medf, and* Pc;51).   <em>Analyses</em>   Our sample is analyzed in two ways: 1) As a single time-averaged sample, which is ideal for characterizing diets, and modern and VFC species analyzed in the same way, with species present in both samples kept discreet in analysis. 2) by depth (in seven consecutive units; Table S5) to track intraspecific variation through time for the four most abundant VFC macropodids, as well as across subfamilies, and independent of taxonomy to assess broad temporal trends.   All data were transformed for normality prior to analysis using the *BestNormalize *package in R. version 4.3.2 (*52*; Data S1). Many of the 35 variables used measure similar aspects of surface metrology, which was investigated through a correlation matrix (Figure S2). Variables with the highest correlation were systematically removed until a smaller set of 18 variables were retained to best describe the data with minimum collinearity (Fig. S3). The variables *Sku, Smrk1, *and* Ssk *were later removed from analysis having shown no significant differences between species, leaving 15 variables for analysis (Table S3).   Data for final variables were modelled using linear mixed-effect modeling (LME), where factors were included only where they improve the ability of the model to fit the data (*46*). LME also allows sub-sampling (*53*) allowing a larger dataset of *n*=2650 scans to be used to delineate dietary differences. Multiple models were constructed following (*46*), with delineating differences between *species,* being principle and all other factors only included where they improved the ability of the model to fit the data. Models were compared using different measures of model fit, and cross-validation on independent hold-out datasets, in this case using five folds (see table S4; *53*).   LME also allows the inclusion of random factors. In DMTA, this is particularly useful as it accounts for inter-individual sample repetition by using specimen as a random factor (*46*). All remaining variables were also permutated as random effects in the modeling process to see if this better matches the distribution any of the variables. As some factors (e.g., *tooth position* and *tooth facet*) may be better considered as nested factors these were also added into the comparison set (Data S2). LME also requires a reference level for each factor, which acts as a base against which all other levels of each factor are compared (see Table S2). For interspecific comparisons, the mixed feeder,* Thylogale stigmatica*, was chosen as reference level as a generalist mixed-feeding species (*13*).   Comparison used multiple measures of model ‘fit’, (Conditional R2, Marginal R2, Sigma, Performance, AIC, and cross-validated R2 were used; see Table S4 for a glossary of these). Any models with &gt;2 singularities within cross-validated subsets were flagged as suspect and not included for comparison (Table S4; Data S2). To balance remaining models, those which fell in the top ten of the most measures were selected, or where multiple models performed equally at this, models were visually compared to see which best differentiated species. Where these were identical, the model with fewest parameters was chosen (Data S2; Figures S4–6).   ANOVA models were then run to provide statistical support for differences evident in the final models constructed by LME modeling, as well as additional tests considering differences between *species* only. For the LME-based models, non-significant factors were dropped from comparison reiteratively until only significant factors remained, with significance defined as the F statistic *p&lt;*0.05. As ANOVA cannot handle random effects, the effect of *specimen* was dropped from ANOVA comparisons. Because of this, a single scan was used for each specimen, with the most commonly sampled teeth and facets used to determine which scan was used. While removing any subsampling effects, this reduced the sample size to *n*=937 (see Table S1 ‘N. specimens’ for the sample size for each species, and Data S1 for the full dataset). Post-hoc comparisons were undertaken using Tukey’s HSD test of pairwise comparisons (Data S1).   To visualize the dataset across the multiple variables being used, and help attain a simplified consensus, a Principal Components Analysis (PCA) was undertaken using the mean for each of the 15 final variables for each taxon. Subsampled data for modeling was removed prior to PCA visualization as per ANOVA dataset above (Data S1).   A distinct set of analyses considered dietary change through time, restricted to specimens from Pit C in the Main Fossil Chamber of VFC, where the most reliable depth information is available (*42*). Analysis considered stratigraphic units 4–8, with the deepest units 4 and 5 combined due to low sample sizes (*42*). In contrast, the fossil-rich, 0.9-m-thick unit 7 was subdivided by depth into units 7A–7E, respectively (Table S5). Analysis was limited to the four most numerous VFC kangaroos: the macropodines *Macropus giganteus *and *Notamacropus rufogriseus*, and the sthenurines *Procoptodon browneorum *and *Pro. gilli *(Table S5). These four species still did not provide sufficient samples for LME modeling, so analysis was restricted to ANOVA comparisons. Datasets were compiled for each species and tested against stratigraphic unit in ANOVA comparisons for each microwear variable to determine if dietary change across time within a species was detectable. For *Ma. giganteus *and *N. rufogriseus, *this analysis was extended to include modern samples to further investigate dietary change over time. To assess how changes in available foods may impact on fitness at any particular time, relative abundances for these four species within each unit were calculated for comparison by dividing the minimum numbers of individuals for each species by the total minimum numbers of individuals for all macropodid species (excluding potoroines) in each unit. Data were acquired from the South Australian Museum Palaeontology register, with taxonomic identities of all specimens checked.   Dietary analysis through time was then extended to consider all taxa lumped into subfamilies Macropodinae and Sthenurinae to see if broad trends were evident across or between these, using the same methodology. Microwear trends over time were then considered across the entire dataset to see if any taxonomically independent changes in diet indicative of changes in vegetation could be detected.   A final set of analyses was undertaken to compare microwear collected here to an earlier study of macropodid microwear (*6*), which itself contained some data collected previously (*15*). Analysis was undertaken using *Asfc *and* epLsar *values published in (*6*), and limited to taxa that were at a generic level included in (*6, 15*) and here. Data were transformed using *bestNormalize* and analyses were made through ANOVA comparisons (Data S1). These considered firstly if there were broad differences between the studies, following (*45*, *54*) and interspecific comparisons.   All data were collected using SensoMAP 7.1.2.7288 (Digital Surf), and analyzed in R. version 4.3.2 (*52*), using the *bestNormalize, car, factoextra, ggplot2, ggcorrplot, lme4, penxlsx, performance, and* stringr packages. Scripts used can be found at the Dryad archive for this article. In all results and discussion ‘significant’ results indicate *p*&lt;0.05.
description: # Dietary breadth in kangaroos facilitated resilience to Quaternary
 climatic variations
 [https://doi.org/10.5061/dryad.4j0zpc8nf](https://doi.org/10.5061/dryad.4j0zpc8nf) The dataset includes two Excel (.xlsx) spreadsheets to accompany the manuscript, and seven *R* scripts to replicate the results. #### File: S1\_data.xlsx The Data S1 spreadsheet contains all of the data collected and more detailed results of the various analyses included in the manuscript. In all worksheets, significant (*p*&lt;0.05) results are highlighted in yellow. The contents of each worksheet is: | Worksheet | Description | | :--------------------------- | :--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | | Raw data LME | All data used in LME analysis. For explanation of factors included see Table S2, and variables in Table S3. Empty Cells for some intraspecific factors, particularly Depth and Ecoregion are where that information was not available. 'NULL' and 'NA' are not used for empty cells as some R functions would not work if included. | | Raw data ANOVA | All data used in the ANOVA analysis (limited to one scan per individual). Includes ‘VFC depth layer’ indicating unit for depth analysis. Empty Cells for some intraspecific factors, particularly Depth and Ecoregion are where that information was not available. 'NULL' and 'NA' are not used for empty cells as some R functions would not work if included. | | How Normalised | Transformation used for normalization, lambda (where applicable) and Pearson’s P / DF (see Table S4) | | ANOVA summary interspecific | Summary of the interspecific ANOVA models. | | ANOVA PHC interspecific | Tukey’s post-hoc comparisons of the interspecific ANOVA models. Lower left of each intraspecific table is the species-only ANOVA; upper left is intraspecific comparisons where models include intraspecific factors (based on LME models). Lowest table is a count of significant (*p*&lt;0.05) interactions. Post-hoc comparisons for intraspecific factors are on the right of interspecific tables. | | ANOVA summary depth analysis | Summary of the depth analysis ANOVA. | | ANOVA PHC depth | Tukey’s post-hoc comparisons of the depth analysis ANOVA. Lowest table of each column is a count of significant (*p*&lt;0.05) interactions. | #### File: S2\_data.xlsx Data S2 presents the outputs of the model comparisons undertaken in the LME modelling. See Table S4 for definitions of the various model comparison terms. Models with &gt;2 singularities (which were removed from consideration) highlighted in red. The ten best performing model for each parameter used to refine model selection are highlighted in green. Final model (Table S3) selected is in bold. Each worksheet is named by the dependent DMTA variable being modelled (see table S3 for full variable names). ##### File: LME process for macropodid microwear This file is the *R* script used to run the LME analysis. The file is marked down with instructions describing how to replicate the results. Users will need to save the 'Raw Data LME' worksheet from Data S1 in an appropriate directory which they can use as the working directory to carry out the analysis. This analysis is somewhat complicated and takes several hours to process for each variable. A more detailed explanation of the method used can be found in Arman et al. (2019) [https://royalsocietypublishing.org/doi/full/10.1098/rsif.2018.0957](https://royalsocietypublishing.org/doi/full/10.1098/rsif.2018.0957). The *R* script also contains instructions for exporting the models for comparison (replicating Data S2), as well as calculating the final modelled DMTA variables for each taxon, and producing the plots to replicate figures S4-6. ##### File: ANOVA interspecific.R This file is the *R* script used to run the interspecific ANOVA comparisons, including Tukey's post-hoc comparisons. The file is marked down with instructions describing how to replicate the results. Users will need to save the 'Raw Data ANOVA' worksheet from Data S1 in an appropriate directory which they can use as the working directory to carry out the analysis. ##### File: Depth\_Analysis\_individual\_species.R This file is the *R* script used to run the ANOVA comparisons between depth layers for the four most abundant species in VFC: *Macropus giganteus, Notamacropus rufogriseus, Procoptodon browneorum *and *P. gilli.* The file is marked down with instructions describing how to replicate the results. The user will need to filter the 'Raw Data ANOVA' worksheet in S1 Data by the 'species' column to create the datasets for analysis. ##### File: Depth\_Analysis\_subfamily.R This file is the *R* script used to run the ANOVA comparisons between depth layers for the subfamilies Macropodinae and Sthenurinae. The file is marked down with instructions describing how to replicate the results. The user will need to filter the 'Raw Data ANOVA' worksheet in S1 Data by the 'Subfamily' column to create the datasets for analysis. ##### File: Depth\_Analysis\_full\_no\_taxon.R This file is the *R* script used to run the ANOVA comparisons between depth layers independant of taxonomy. The file is marked down with instructions describing how to replicate the results. The user will need to use the 'Raw Data ANOVA' worksheet in S1 Data to undertake the analysis. ##### File: PCA.R This file is the *R* script used to run the Principal Components Analysis features in Figures 3 and S7. The file is marked down with instructions describing how to replicate the results, including calculating the means for each taxon for each variable. The user will need to use the 'Raw Data ANOVA' worksheet in S1 Data to undertake the analysis.