Dietary breadth in kangaroos facilitated resilience to Quaternary climatic variations

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.

Methods
Modern samples

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).

Fossil samples

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.

Data acquisition

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 <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).

Analyses

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 >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<*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*<0.05.