Late Pleistocene paleoecology and phylogeography of woolly rhinoceroses

The woolly rhinoceros (Coelodonta antiquitatis) was a cold-adapted herbivore, widely distributed from western Europe to north-east Siberia during the Late Pleistocene. Previous studies have associated the extinction of the species ~14,000 calendar years before present to climatic and vegetational changes, suggesting the later survival of populations in north-east Siberia may have related to the later persistence of open vegetation in the region. Here, we analyzed carbon (d13C) and nitrogen (d15N) stable isotopes and mitochondrial DNA sequences to elucidate the evolutionary ecology of the species. Our dataset comprised 286 woolly rhinoceros isotopic records, including 192 unpublished records, from across the species range, dating from >58,600 to 12,135 14C years before present (equivalent to 14,040 calendar years ago). Crucially, we present the first 71 isotopic records available to date of the 15,000 years preceding woolly rhinoceros extinction. The data revealed ecological flexibility and geographic variation in woolly rhinoceros stable isotope compositions across time. In north-east Siberia, we detected stability in d15N through time, which could reflect long-term environmental stability, and may have enabled the later survival of the species in the region. To further investigate the paleoecology of woolly rhinoceroses, we compared their isotopic compositions with other contemporary herbivores. Our findings suggested isotopic similarities between woolly rhinoceros and both musk ox (Ovibos moschatus) and saiga (Saiga tatarica), albeit at varying points in time, and possible niche partitioning between woolly rhinoceros and both horse (Equus spp.) and woolly mammoth (Mammuthus primigenius). To provide phylogeographic context to the isotopic data, we compiled and analyzed the 61 published mitochondrial control region sequences. The genetic data showed a lack of geographic structuring; we found three haplogroups with overlapping distributions, all of which showed a signal of expansion during the Last Glacial Maximum. Furthermore, our genetic findings support the notion that environmental stability in Siberia influenced the paleoecology of woolly rhinoceroses in the region. Our study highlights the utility of combining stable isotopic records with ancient DNA to advance our knowledge of the evolutionary ecology of past populations and extinct species. © 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Rey-Iglesia), elinelorenzen@ r Ltd. This is an open access article


Introduction
The Late Pleistocene was characterized by large-scale climatic and environmental change (Hubberten et al., 2004). Around 33 thousand years before present (ka BP), a progressive cooling started, which led to the cold and dry environment of the Last Glacial Maximum (LGM;~28.6e20.5 ka BP) . This was followed by an increase in temperature, which peaked at the start of the Holocene 11.7 ka BP. These pronounced climatic perturbations led to shifts in the geographic distribution of species and in the composition of entire ecosystems.
During the Late Pleistocene, large mammal species, including cave bear (Ursus spelaeus), cave lion (Panthera spelaea), and woolly rhinoceros (Coelodonta antiquitatis) went extinct in northern Eurasia (Pacher and Stuart, 2009;Lister 2011, 2012). Other Pleistocene megafauna, including giant deer (Megaloceros giganteus) and woolly mammoth (Mammuthus primigenius), experienced strong reductions in their distributions during the Late Pleistocene, but only disappeared from the faunal record later, during the Holocene (Stuart et al., 2002;Lister and Stuart, 2019). These different extinction patterns among species have been attributed to the interplay of several drivers, including climatic, anthropogenic, and/ or genetic factors (Cooper et al., 2015;Pe cnerov a et al., 2017).
In the northern hemisphere, Late Pleistocene glacial phases were characterized by the mammoth steppe ecosystem, which comprised a mosaic of steppe-tundra and shrub vegetation, with high nutrient soils that enabled the growth of plants that sustained grazing species (Guthrie 1982(Guthrie , 2001. The mammoth steppe extended from the northern Iberian Peninsula to Canada, and was characterized by a diverse community of large herbivores. The woolly rhinoceros was one of the iconic inhabitants of this steppe-tundra ecosystem. The species was adapted to cold environments; it was covered with thick and long hair, as documented by mummified remains, but is believed to not have been well adapted to snowfall, due to its massive body on short legs (Boeskorov et al., 2011). Teeth structure and mesowear analysis indicate a grazing diet (Stuart and Lister, 2012;Rivals and Lister, 2016;Stefaniak et al., 2020), which is supported by grass remains recovered between the teeth of some specimens (Guthrie, 1990). Microwear analysis of North Sea woolly rhinoceroses suggests periodical inclusion of woody components in their diet (van Geel et al., 2019). Pollen analysis of stomach contents has primarily identified grasses and sagebrushes (e.g. Boeskorov et al., 2011), and genetic analysis of stomach and gut content supports a diet of primarily grasses, with a contribution of forbs (Willerslev et al., 2014).
Prior to their extinction~14 thousand calendar years before present (cal kyr BP), woolly rhinoceroses ranged from western Europe to north-east Siberia (Fig. 1). Specimens have not been recovered from certain areas in Europe and north-central Siberia, suggesting the species was absent from these regions (Stuart and Lister, 2012). Despite their presence in the fossil record in far northeastern Siberia, adjacent to the Bering Land Bridge, woolly rhinoceroses did not colonise North America; it has been suggested that the environmental conditions of the Bering Land Bridge and the paleoecology of eastern Beringia (present-day Alaska and Yukon territory) were not adequate for the species (Boeskorov, 2011;Stuart and Lister, 2012).
Despite a progressive eastward contraction in the species' fossil record starting~35 cal kyr BP (Stuart and Lister, 2012), fossil evidence indicates woolly rhinoceroses were still present in large parts of Eurasia until~16 cal kyr BP, suggesting an almost synchronous or at least rapid extinction across their range (Lorenzen et al., 2011). Demographic modelling based on genomic data has documented a population increase~30 cal kyr BP, followed by demographic stability until close to woolly rhinoceros extinction (Lord et al., 2020). From 14.6 cal kyr BP, the range of the species had contracted so that woolly rhinoceroses were found only in the Ural Mountains, southern Siberia, and NE Siberia. The most recent fossil records are~14 cal kyr BP from NE Siberia; the youngest fossil has been dated to 12,135 14 C years BP or 14,040 cal yr BP (Boeskorov, 2011;Stuart and Lister, 2012).
Climate and vegetational changes have in combination been invoked as the main causative agents of woolly rhinoceros extinction (Lorenzen et al., 2011;Stuart and Lister, 2012;Lord et al., 2020). Environmental shifts included increased ground moisture due to precipitation and snowfall in the Late Glacial Interstadial (14.7e12.9 cal kyr BP) (Sher, 1997), and the spread of shrubs and trees (Stuart and Lister, 2012). These new environmental conditions would have restricted areas of firm ground and open vegetation that are believed to have been the most suitable habitat for woolly rhinoceroses (Stuart and Lister, 2012). The later survival of the woolly rhinoceros in NE Siberia may be linked to the late persistence of open vegetation in the region, compared with the rest of Eurasia (Stuart and Lister, 2012). NE Siberia also played a key role for other Late Pleistocene megafauna species, and the region was the last mainland area of survival of both steppe bison (Bison priscus) and woolly mammoth prior to their extinctions (Vartanyan et al., 1993;Kirillova et al., 2015a).
In herbivores, d 13 C and d 15 N values reflect the isotopic compositions of plants in their diet (DeNiro, 1985;Bocherens, 2003), the composition of which is in turn influenced by climatic and environmental factors (Hartman, 2011;Bonafini et al., 2013). Carbon isotopic compositions in herbivore bone and tooth collagen reflect the photosynthetic pathways and environmental parameters of plants at the base of the food web (Bocherens, 2003). Nitrogen derives from the plant protein the animal digests; dietary choices may lead to differences in d 15 N, even at the same trophic level (Bocherens, 2003). Furthermore, environmental factors such as water stress, salinity or grazing pressure, are known to affect plant d 15 N values, and changes in vegetation and soil nutrient cycling will therefore be recorded in bone and tooth collagen.
In this study, we analyzed d 13 C and d 15 N retrieved from woolly rhinoceros remains spanning 45 kyr of evolutionary history. We present 192 new isotopic records from across the species range, dating from >58,600 to 12,135 14 C years before present (equivalent to 14,040 calendar years ago). These include 71 isotopic records from the 15 kyr preceding woolly rhinoceros extinction, which represent the first available stable isotope records from this period.
We combined the 192 new records with published d 13 C and d 15 N data of 94 woolly rhinoceroses, totalling 286 samples, to (i) investigate temporal and spatial changes in woolly rhinoceros diet and ecology; (ii) specifically assess changes in the millennia preceding extinction; and (iii) elucidate the later survival of the species in NE Siberia.
To understand the broader ecological context of our findings, we compared the woolly rhinoceros stable isotope records with comparative data from five other contemporaneous megaherbivore species (horse (Equus spp.), musk ox, reindeer (Rangifer tarandus), saiga, woolly mammoth).
To elucidate the phylogeographic context of the palaeoecological data, we analyzed a comprehensive data set of 61 available genetic sequences retrieved from woolly rhinoceros fossils sampled across time and space. The data comprised mitochondrial DNA (mtDNA) control region sequences, a non-coding area with a high substitution rate that is ideal for inferring phylogeographic patterns. The data were analyzed to (i) investigate the distribution of genetic diversity across time and space, which has not previously been done; and (ii) associate patterns of mtDNA diversity and differentiation with isotopic variation.

Radiocarbon and stable isotope data
We present 192 new d 13 C and d 15 N measurements from woolly rhinoceros bone and tooth specimens (Supplementary Table 1), determined by isotope ratio mass spectrometry (CF-IRMS). These measurements were generated during 14 C radiocarbon dating of fossils for two previous studies (Lorenzen et al., 2011;Stuart and Lister, 2012). These two studies only reported 14 C ages, and not the corresponding d 13 C and d 15 N measurements. For the present study, we accessed and compiled the latter data. Although the woolly rhinoceros fossil records presented in (Lorenzen et al., 2011;Stuart and Lister, 2012) are not temporally or spatially exhaustive, and some geographic regions are underrepresented, they include a large number of specimens collected across time and space ( Fig. 1). We calibrated the 14 C dates retrieved from the previous publications using CALIB 8.2 (Stuiver et al., 2021) and the IntCal20 (Reimer et al., 2020) calibration curve. Unless otherwise stated, dates are reported as calibrated years before present (cal yr BP). Dates reported as infinite correspond to samples where the 14 C content was too small to produce a finite radiocarbon age. The samples ranged in age from infinite to 14,040 cal yr BP. The spatial distribution of the material spans northern Eurasia, with records from Spain to NE Siberia (Fig. 1).
Radiocarbon dating and isotopic measurements of samples were conducted at Oxford Radiocarbon Accelerator Unit (University of Oxford; OxA codes) and Aarhus AMS Center (University of Aarhus; AAR codes). At Oxford Radiocarbon Accelerator Unit, d 13 C and d 15 N were determined following the combustion/recycling protocol described in Brock et al., (2010). At Aarhus AMS Center, d 13 C and d 15 N were determined following the protocols described in Olsen et al. (2010)). Reliability of the stable isotopic values was established by measuring the atomic C:N ratio. All our samples fell within the accepted range of 2.9e3.6 (DeNiro, 1985; Ambrose, 1990).
To increase the number and range of d 13 C and d 15 N records in our analysis, we compiled a range-wide database of available d 13 C and d 15 N stable isotope records for woolly rhinoceros (Supplementary Table 1). We used Google Scholar to perform a LGM (24,600e17,000 14 C years BP/28,660e20,520 cal years BP), and post-LGM (<17,000 14 C years BP/< 20,520 cal years BP). (B) Chronology of the 14 C radiocarbon dated samples used in this study for each geographic region. Plotted dates represent uncal 14 C years BP, extended lines represent infinite dates. Symbol colour and shape represent region and data type, respectively. Sample size (n ¼ all samples, n ¼ 14 C dated samples) is indicated for each geographic region. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) literature search with the terms "woolly rhinoceros" and "stable isotopes", and "Coelodonta antiquitatis" and "stable isotopes". We recovered data for 94 specimens, derived from ten published studies and spanning the species' range (Bocherens et al. , 1996(Bocherens et al. , 1997(Bocherens et al. , 2005(Bocherens et al. , 2011(Bocherens et al. , 1996(Bocherens et al. , 2005Higham et al., 2006;Jacobi et al. 2006Jacobi et al. , 2007Jacobi et al. , 2009Jacobi et al. , 2007Kirillova et al., 2015b;Kuc et al., 2012;Yates et al., 2017). All previously published data represented samples dated >28,660 cal yr BP. Our final dataset represented 286 samples (Fig. 1).
To contextualize the woolly rhinoceros isotopic data with other contemporary large herbivores, we performed a literature search and compiled d 13 C and d 15 N data for five other mammoth steppe species from Eurasia: horse, musk ox, reindeer, saiga, and woolly mammoth (Supplementary Table 2). Paleoecological inference based on d 13 C and d 15 N data is highly influenced by regional environmental characteristics, and large isotopic differences have been observed for d 13 C and d 15 N between the North American and Eurasian regions of the mammoth steppe (Szpak et al., 2010). As woolly rhinoceroses were only present in Eurasia, we omitted North American records for the other species.

Stable isotope analysis
To investigate temporal variation in the isotopic signature of woolly rhinoceroses and other megaherbivores, we divided the isotopic records into three time bins ( LGM (24,600e17,000 14 C yr BP/28,660e20,520 cal yr BP), and post-LGM (<17,000 14 C yr BP/< 20,520 cal yr BP). The post-LGM cut-off date corresponds to 12,135 14 C yr BP or 14,040 cal yr BP, the age of the youngest woolly rhinoceros specimen in our dataset. The sample size of each time bin was: pre-LGM n ¼ 215, LGM n ¼ 32, post-LGM n ¼ 39 (Table 1).
To accommodate for differences in fossil chronology across Eurasia (such as the later survival of woolly rhinoceros in NE Siberia), and the environmental heterogeneity of Eurasia during the Late Pleistocene, we divided the data into four geographic regions: (i) western and central Europe (n ¼ 97); (ii) European Russia, Urals and Kazakhstan (n ¼ 57); (iii) southern and central Siberia and China (n ¼ 46); and (iv) north-east Siberia (n ¼ 86, Table 1). To simplify, we refer to these regions as W/C Europe, E Europe/Urals (as we only have one sample from Kazakhstan), S/C Siberia, and NE Siberia, respectively.
Stable isotopes were retrieved from both bone and tooth. The isotopic signature in bone material provides information on the average diet of an individual in the years preceding death (Hedges et al., 2007;Hobson and Clark, 1992). In contrast, data from mammalian tooth collagen corresponds to the early years of an individual's life, depending on the chronology of tooth development for the species in question (DeNiro, 1985;Bocherens, 2015).
Furthermore, their d 15 N signature may reflect the suckling 15 Nenriched composition (DeNiro, 1985;Fizet et al., 1995;Bocherens and Drucker, 2013). We have no direct information on the weaning age of woolly rhinoceros, but infer a weaning age of 18 months based on the average age of weaning of a number of extant rhinoceros species (e.g. Owen-Smith, 1974;Osthoff et al., 2008;Gimmel et al., 2018). We do not have information on the age of the individuals used in this study, and assume all samples analyzed belonged to adult individuals.
For the purposes of our analysis, we assumed that the stable isotope composition of bone and tooth collagen are comparable, albeit systematic differences in d 13 C and d 15 N isotopic abundances between bone and dentine have been recognized in some mammals. For instance, lower d 15 N values in bone compared to dentine have been measured in carnivores such as bears and hyenas (Bocherens et al., 1997;Bocherens, 2015), or in herbivores such as reindeer (Fizet et al., 1995). For the species where the isotopic differences between bone and tooth have been directly identified and quantified, d 13 C and d 15 N data are adjusted before comparing bone and tooth isotopic values. However, there has been no systematic study characterising possible d 13 C and d 15 N isotopic differences between bone and tooth in woolly rhinoceros, and our data have thus not been corrected. We are aware of the bias that this might introduce to our data. However, assuming the correction factor for woolly rhinoceros is similar to estimates from other herbivore species, such as reindeer (d 13 C ¼ 0.2‰e0.5‰ and d 15 N ¼ 0.7‰e2‰, lower in bone than in dentine, Fizet et al., 1995), the skeletal element variation in our data set is less than the variation we detect among individuals across time bins and geographic regions. Groups (representing time or space) with low sample size (n < 30) were tested for normality using the Shapiro-Wilkinson test. As the data distribution for some of these groups was not normal, we chose to use non-parametric testing for all analyses comparing isotopic signatures for different time bins and geographic regions. The Kruskal-Wallis test was used to determine statistical differences for each time bin among geographic regions.
To further detect differences between pairs of time bins and/or geographic regions, we used the post-hoc Kruskal Conover test applying the Bonferroni correction to control for Type I errors. Statistical data analysis was conducted in RStudio (Team RStudio, 2018) using the packages nortest (Gross and Ligges, 2015) and Table 1 Summary of d 13 C and d 15 N values of the 286 woolly rhinoceros specimens analyzed in this study. The time bins were defined as: pre-LGM (>24,600 14 C yr BP/> 28,660 cal yr BP), LGM (24,600e17,000 14 C yr BP/28,660e20,520 cal yr BP), and post-LGM (<17,000 14 C yr BP/< 20,520 cal yr BP). Number of specimens in each temporal/spatial bin is indicated in parenthesis.

Time bin (n)
Region ( PMCMR (Pohlert and Pohlert, 2018). Significant differences were reported at the 95% confidence level or p-value of 0.05. Bivariate plots and boxplots were performed in RStudio (Team RStudio, 2018) using the package ggplot2 (Wickham, 2011). Summary statistics were obtained in RStudio. Multispecies comparisons were also performed using SIBER (Jackson et al., 2011), this R package was used to estimate and compare the Standard Ellipse Area (SEA) as a proxy for isotopic niche range.

DNA data compilation and analysis
To contextualize the stable isotope data and evaluate genetic changes across time and space, we compiled all available mtDNA control region sequences of woolly rhinoceros published to date (Lorenzen et al., 2011;Lord et al., 2020). Due to a large amount of missing data for the control region locus, eight of the 14 sequences from Lord et al., (2020) were not included in our analysis (Supplementary Table 3). We did not include the mtDNA sequence published in Willerslev et al. (2009), as the sample was not dated.
The final dataset comprised 61 mtDNA control region sequences; 35 of these samples were also represented in the stable isotope analysis (Supplementary Table 3). Sequences were aligned in Geneious (Kearse et al., 2012), and the length of the final alignment used in our analysis was 541 base-pairs (bp). We generated a haplotype network using the median-joining (Bandelt et al., 1999) algorithm implemented in the program PopART v1.7 (Leigh and Bryant, 2015).
To estimate a phylogeny, we used BEAST v.1.10.4 (Drummond et al., 2012) using the radiocarbon dates of the sequences to estimate a mutation rate and divergence times. This analysis was performed using the same dataset as for the haplotype network. We used a coalescent constant size model, as data were from a single species, and applied a strict clock. The clock rate was set to a normal distribution with an initial value of 6.1 Â 10 À9 substitutions/ site/year, a mean value of 6.1 Â 10 À9 , and a standard deviation of 0.01 as in Lord et al., (2020). As estimated in PartitionFinder (Lanfear et al., 2012), we used the model HKY þ I þ G for all the codons of the control region.
All BEAST analyses were run in two independent MCMC chains of 1 Â 10 7 generations each, sampling trees and model parameters every 1 Â 10 3 generations. Tracer v1.6 (Rambaut et al., 2018) was used to combine and inspect the results of each run and to determine the convergence of each parameter, all of which had ESS values > 200. We identified the Maximum Clade Credibility (MCC) tree in TreeAnnotator v1.8.0, and visualized and graphically edited the MCC tree using Figtree.
The covariation of stable isotopes and haplogroups was also explored. A bivariate plot was performed in RStudio (Team RStudio, 2018) using the package ggplot2 (Wickham, 2011). A cluster dendrogram was generated using the R package cluster (Maechler, 2019). For this, we applied the Euclidean distance method algorithm with Ward's minimum variance criterion to minimize the total within-cluster variance. We set the number of clusters to three, equivalent to the number of haplogroups recovered in the phylogenetic analysis.

Stable isotope analysis: woolly rhinoceros
Summary statistics of woolly rhinoceros stable isotope measurements were compiled in Table 1 and visually explored in Figs. 2e4 d 13 C values ranged from À22.7‰ to À18.1‰, with an average value of À19.9 ± 0.71‰. d 15 N values ranged from þ0.7‰ to þ11.5‰, with an average value of þ5.7 ± 2.04‰. Summary statistics indicated a smaller range of variation for d 13 C than d 15 N. Fig. 2A shows variation in d 13 C through time. We detected an increase in d 13 C starting~34 cal kyr BP, which peaked at the onset of the LGM~27 cal kyr BP, with a subsequent decline. We further explored differences in d 13 C within the three time bins using bivariate plots, and found higher average d 13 C during the LGM than for the other two time bins (Fig. 3). Statistical testing showed significant differences in average d 13 C between pre-LGM and LGM (p < 0.001), and between LGM and post-LGM (p ¼ 0.003).
We investigated the influence of time on d 13 C within each of the four geographic regions (Fig. 4A). For three regions (excluding W/C Europe), we detected significant changes in average d 13 C over time (p-values: E Europe/Urals: p ¼ 0.02; S/C Siberia: p ¼ 0.03; NE Siberia: p ¼ 0.008), with higher values during the LGM. Within each time bin, we identified differences in isotopic values across regions (Fig. 4A, Supplementary Table 4).
For d 15 N, we observed a gradual decline starting at the onset of the LGM~27 cal kyr BP (Fig. 2B). When we split the range-wide dataset into three time bins, statistical testing showed significant differences in d 15 N between pre-LGM and post-LGM (p ¼ 0.002) (Fig. 3). When the data were further investigated across time for each geographic region, we detected significant continuous declines in average d 15 N over time for two regions (p-values: W/C Europe: p ¼ 0.010, E Europe/Urals: p ¼ 0.004), but not for S/C and NE Siberia (Fig. 4B).
For pre-LGM, pairwise comparisons between regions showed significant differences between NE Siberia and both W/C Europe and E Europe/Urals (Fig. 4B, Supplementary Table 4). For LGM, we did not identify any differences in d 15 N between regions. For post-LGM, NE Siberia d 15 N differed significantly from E Europe/Urals.

Stable isotope analysis: other herbivores
To contextualize woolly rhinoceros isotopic variation, we compared the data with stable isotope records from five codistributed herbivore species. To make the data directly comparable across time and space, we restricted the geographic distribution of the data to Eurasia (Supplementary Table 2). Bivariate and boxplots of the data divided by species and by time bin correspond to Figs. 5 and 6. Graphical exploration and statistical testing indicated that pre-LGM d 13 C and d 15 N of rhinoceros differed from all the other species surveyed (Figs. 5 and 6, Supplementary Table 5). Rhinoceros LGM d 13 C and d 15 N were similar to saiga. Post-LGM rhinoceros d 13 C and d 15 N values were similar to musk ox.
To estimate species isotopic ranges, we used the Standard Ellipse Area (SEA) function implemented in SIBER (Supplementary Table 6). For pre-LGM, saiga had the largest isotopic range (SEA ¼ 5.16), followed by rhinoceros (SEA ¼ 4.25). For LGM, we observed a decline in isotopic range for four species, except mammoth (SEA ¼ 2.91) and saiga (SEA ¼ 5.30). For rhinoceros, SEA was reduced to 3.65, still remaining second largest after saiga. For post-LGM, rhinoceros had the largest SEA (SEA ¼ 4.98), with an increase in ellipse area from the LGM (LGM SEA ¼ 3.65). For rhinoceros, post-LGM isotopic range returned to pre-LGM values, albeit with overall lower d 15 N values.

DNA analysis
We identified three haplogroups (A-C) in the mtDNA control region data (Fig. 7, Supplementary Figure 1, and Supplementary Table 3). The alignment contained 66 segregating sites; haplogroups A and B differed from each other at two nucleotide positions, and haplogroups B to C at six nucleotide positions.
Haplogroup A was only found in S/C and NE Siberia; it was present in NE Siberia pre-LGM, in S/C Siberia during LGM, and in Haplogroup C was present across the species range in pre-LGM, extending from central Europe to NE Siberia. During LGM, haplogroup C remained in S/C and NE Siberia. Post-LGM, only woolly rhinoceroses from E Europe/Urals carried this haplogroup.
We also explored the covariation of stable isotopes and haplogroups. We could not identify any relationship between them (Supplementary Figure 2).

Discussion
We analyzed d 13 C and d 15 N stable isotope data from 286 woolly rhinoceros specimens sampled across Eurasia, spanning over 45 kyr of evolutionary history. Our analyses revealed variation across geographic regions; S/C and NE Siberia showed stability in isotopic composition across time, in contrast with other regions, which experienced variation in d 15 N. We found no geographic structuring among the 61 mtDNA sequences analyzed, but did recover three well-differentiated haplogroups with overlapping distributions, all of which showed a signal of expansion during the LGM.

Temporal variation in d 13 C within and among regions
Across samples, we found an average d 13 C value of À19.9‰, indicating woolly rhinoceroses had a diet dominated by C 3 plants (Bocherens, 2003); this is expected in herbivores from cold and temperate climates. We detected significant changes in average d 13 C across time bins; between pre-LGM and LGM, and between LGM and post-LGM (Fig. 3).
Although values remained within the C 3 diet range, we found   1995). Drier environmental conditions during the LGM have been inferred across Eurasia based on paleoenvironmental proxies (Hubberten et al., 2004), and are further supported by bioclimatic models, which show the expansion of tundra and simultaneous reduction of boreal and temperate forests 21 ka BP (Hoogakker et al., 2016).
Post-LGM, average d 13 C declined, reaching values similar to pre-LGM (Fig. 3). Plants from wet environments display lower d 13 C than those from dry environments (Wooller et al., 2007), and our findings of an overall decline in d 13 C after the LGM may reflect an increase in moisture due to precipitation and degrading permafrost during the Late Glacial (Rabanus-Wallace et al., 2017).
We identified significant changes in d 13 C over time for three of the four geographic regions (excluding W/C Europe, Fig. 4A). In the E Europe/Urals and S/C Siberia, we observed a decrease in d 13 C between LGM and post-LGM. In NE Siberia, we detected an increase in d 13 C from pre-LGM to LGM, which remained stable during the post-LGM.
Woolly rhinoceroses ranged across the Eurasian mammoth steppe, a heterogeneous environment encompassing a mosaic of varied local environments (Szpak et al., 2010). This was reflected in our pairwise comparisons between regions, which indicated spatial differences in d 13 C (Supplementary Table 4); for instance, we observed lower average d 13 C in woolly rhinoceroses from NE Siberia compared with the three other regions, in particular compared to E Europe/Urals and S/C Siberia during pre-LGM and LGM (Fig. 4A, Table 1). Lower d 13 C has previously been associated with lower annual temperatures (Szpak et al., 2010), and this may also be the case for NE Siberia. However, this finding is unexpected, as NE Siberia was characterized by consistently arid conditions throughout the Late Pleistocene (Sher et al., 2005), and plants from dry environments have higher d 13 C (Wooller et al., 2007), which should be reflected in the d 13 C composition of herbivores.
Furthermore, the observed regional differences could be driving some of the overall patterns detected in our dataset; for example, the overall higher d 13 C values during the LGM could be caused by the higher d 13 C from E Europe/Urals and S/C Siberia during that period.

Temporal variation in d 15 N within and among regions
Our range-wide dataset showed significant changes in d 15 N over time. A significant, gradual decline in d 15 N was observed across the three time bins in two of the geographic regions (W/C Europe and E Europe/Urals, Fig. 4B). The decline observed from pre-LGM to LGM in S/C Siberia was not significant. A similar pattern of continuous decline in d 15 N from the Late Pleistocene until after the LGM has been observed in European reindeer, attributed to the climatic cooling that culminated during the LGM (Stevens et al., 2008). The post-LGM decline in d 15 N, which is also reported in reindeer, has been linked to an increase in moisture due to precipitation and degrading permafrost (Drucker et al. 2003Stevens et al., 2008). These more recent environmental changes have been attributed a major role in the Late Pleistocene and Holocene megafauna extinctions (Rabanus-Wallace et al., 2017;Drucker et al., 2018).
Between regions, we observed some temporal differences in d 15 N (Fig. 4B) (Bocherens, 2003). Thus, it is not unexpected that we find differences in woolly rhinoceroses d 15 N among regions, which were heterogeneous (Szpak et al., 2010). Moreover, the forms and concentrations of nitrogen that are present in the soil (i.e., mineralized or organic N) will influence the mechanisms by which plants obtain nitrogen for growth. Specifically, the interactions that plants form with symbiotic fungi (mycorrhizae) are significant, as The continuous decline in average d 15 N observed in W/C Europe and E Europe/Urals may reflect changes in the nutrient status of the soil, such as increasing levels of recalcitrant organic nitrogen and decreasing levels of mineralized nitrogen (e.g., NH 4 þ and NO 3 À) that are readily available for plants. If the entire system becomes more dependent on fungi to process organic N and deliver it to plants, the proteolytic capabilities of those fungi impart very low d 15 N values on the plants that receive the mineralized nitrogen (Hobbie and H€ ogberg, 2012). The reduction in the number of herbivores in the landscape, which has been described in previous studies (e.g. Lorenzen et al., 2011;Lister & Stuart, 2019), may also have played a role in the decline in d 15 N values. Herbivores are effective recyclers of nitrogen, by depositing it back into the soil as urine and feces. An increase in the abundance of herbivores would lead to more rapid nitrogen cycling, providing more mineralized nitrogen available for plants (McNaughton, 1979). Thus, there could have been an important feedback loop that was being altered in areas where megafauna was declining.

Temporal stability in d 15 N in NE Siberia
The fossil record, based on published radiocarbon dates, shows the continuous presence of woolly rhinoceroses in NE Siberia from >50 cal kyr BP until their extinction~14 cal kyr BP (Stuart and Lister, 2012). The region was home to the last surviving population of woolly rhinoceros, which survived there~1000 years later than elsewhere.
Later persistence of open vegetation, which may have presented more suitable habitat and food resources, has been suggested to have played a role in the later survival of woolly rhinoceroses in NE Siberia (Stuart and Lister, 2012). This is supported by palaeoecological proxies, which show environmental stability in NE Siberia throughout the Late Pleistocene (Clark et al., 2012;Jørgensen et al., 2012;Willerslev et al., 2014), and is further supported by our findings.
We found no significant changes in d 15 N over time in NE Siberia ( Fig. 4B). This pattern indicates dietary stability, suggesting woolly rhinoceroses fed on resources with similar isotopic compositions through time. Such a scenario would require a stable environment.
Temporal stability in d 15 N in NE Siberia has also been observed in musk ox from Taimyr and in woolly mammoth from NE Siberia (Raghavan et al., 2014;Kuitems et al., 2019), suggesting broad-scale environmental stability, rather than a species-specific pattern.
Environmental stability may have been key to the late survival of woolly rhinoceroses (and other species) in the region. In woolly mammoths, isotopic stability has been reported in both mainland NE Siberia and on Wrangel Island. Both areas were continuously inhabited by the species until their mainland and island extinctions, respectively (Arppe et al., 2019;Kuitems et al., 2019). The survival of mammoths on Wrangel Island well into the Holocene is attributed to dietary stability, and to surviving populations feeding on resources with similar isotopic composition as the Late Pleistocene mammoth population in NE Siberia (Arppe et al., 2019;Kuitems et al., 2019). We suggest a similar scenario in woolly rhinoceroses, where the NE Siberian population retained a stable diet throughout the Late Pleistocene, until its extinction.
Despite being seen as a specialized grazer, woolly rhinoceroses exhibited variable isotopic values that overlap with other large herbivores through time and space (Figs. 5 and 6). This suggests a higher ecological flexibility than expected, and is supported by other studies on woolly rhinoceros diet (Kosintsev et al., 2019;Stefaniak et al., 2020). We found the isotopic composition of woolly rhinoceros overlapped with musk ox and saiga, albeit during different time bins. Woolly rhinoceros pre-LGM average d 13 C and d 15 N differ from the other analyzed herbivores, with intermediate average d 13 C values. During the LGM, woolly rhinoceros average d 13 C and d 15 N were similar to saiga, and post-LGM values were similar to musk ox; all three species are grazers. Saigas feed predominantly on grasses and chenopods (Hopkins et al., 2013). Musk ox are classified as grazers based on the anatomy of their digestive system (Knott et al., 2004), and they feed primarily on grasses, sedges, and willows (Raghavan et al., 2014). The isotopic affinity of woolly rhinoceroses to these two species was not constant across time bins. Standard Ellipse Analysis shows that patterns of variation in isotopic range across time bins is species-specific (Fig. 5). Thus, it is expected that isotopic niches changed across time periods, reflecting species-specific responses to environmental change.
Overall, our results showed that woolly rhinoceros d 13 C values are higher than horse and woolly mammoth, the two other monogastric hindgut fermenters included in our analysis, and we suggest this may reflect ecological niche partitioning or even dietary competition between woolly rhinoceros and woolly mammoth (Figs. 5 and 6). Such dietary competition, and at times  ecological exclusion, has been demonstrated in ecological studies of Asiatic elephants and Indian rhinoceroses (Pradhan et al., 2008), and of African elephants and black rhinoceroses (Landman et al., 2013).

mtDNA shows absence of phylogeographic structuring and LGM lineage expansion
Previous studies have used mtDNA to address the phylogenetic placement of woolly rhinoceros relative to other rhino species (Orlando et al., 2003;Willerslev et al., 2009), to explore the mitogenomics of NE Siberian woolly rhinoceroses (Lord et al., 2020), and to elucidate the demographic history of the species (Lorenzen et al., 2011;Lord et al., 2020). However, a comprehensive phylogeographic study across time and space has so far been lacking. To investigate spatial and temporal patterns of genetic variation across Eurasia, we compiled and analyzed 61 radiocarbon dated DNA sequences from Late Pleistocene woolly rhinoceroses.
We identified three well-differentiated mtDNA haplogroups (A-C) (Fig. 7). Despite the haplotype network analysis separating haplogroups A and B by only two nucleotide differences, our BEAST analysis showed strong posterior support differentiating them (Supplementary Figure 1). Our analysis included six control region sequences from Lord et al. (2020), which in that study, based on entire mitochondrial genomes, were assigned to woolly rhinoceros mitogenome clade 1 (n ¼ 2) and clade 2 (n ¼ 4). Mitogenome clade 1 corresponded to our haplogroup B, and mitogenome clade 2 corresponded to our haplogroup A, further supporting our haplogroup division based on the shorter control region sequences.
The mtDNA control region haplogroups were not geographically structured in time and/or space, and extended from Central Europe to NE Siberia; the only mtDNA sequence included from W/C Europe grouped in haplogroup C (Fig. 7). During the pre-LGM, haplogroup A was present in NE Siberia only, with an LGM expansion into S/C Siberia; the first appearance of haplogroup A in S/C Siberia was dated to 25,870 cal yr BP. A similar pattern was observed in haplogroup B; pre-LGM it was limited to S/C and NE Siberia, and expanded into E Europe/Urals during the LGM. The expansion of haplogroups A and B into new regions occurred at a time where demographic analysis of subsets of the data indicates an increase in population size (Lorenzen et al., 2011;Lord et al., 2020).
The presence of spatially overlapping yet divergent haplogroups likely reflects the allopatric divergence of lineages, with subsequent population expansion leading to the distribution of woolly rhinoceros populations across northern Eurasia. Siberia harbored the highest levels of genetic diversity, and was the only area where all three haplogroups were found. However, DNA analysis of additional specimens from the western part of the species range are needed to confirm this pattern.
For other megaherbivore species, such as bison, high genetic diversity has also been identified in NE Siberia (Kirillova et al., 2015a;Massilani et al., 2016). Genetic data from wolves (Canis lupus) suggest Late Pleistocene populations likely originated and expanded from NE Siberia into Eurasia and North America (Loog et al., 2020). Climatic stability has been suggested as a possible explanation for this pattern of maintained genetic diversity in NE Siberia (Loog et al., 2020); a scenario supported by our study.
Demographic analysis of 55 of the 61 mtDNA sequences analyzed here showed changes in genetic connectivity through time, with increasing fragmentation prior to extinction leading to population isolation (Lorenzen et al., 2011). The phylogeographic pattern revealed in our analysis supports this hypothesis of a decline in connectivity over time; we observed some changes in the geographic distribution of haplogroups over time, although we did not detect a loss of mtDNA haplogroups from pre-LGM to post- LGM. Two of the youngest samples included in our study (dated to 14,390 and 15,340 cal yr BP) were from NE Siberia and represented haplogroups A and B, respectively (Fig. 7). The youngest haplogroup C sample was from E Europe/Urals and is dated to 16,390 cal yr BP.
Previous analysis of 14 mitochondrial genomes from NE Siberia also support lineage continuity from pre-LGM to post-LGM in this region (Lord et al., 2020). However, further data are needed to confirm this pattern, as inferences in (Lord et al., 2020) were based on a small number of individuals, with an overrepresentation of sequences from NE Siberia.

Implications for woolly rhinoceros extinction
We investigated a comprehensive d 13 C and d 15 N dataset of woolly rhinoceroses, and included the first 71 stable isotope records spanning the 15 kyr preceding species extinction. Our study underscores the applicability of combining range-wide stable isotope studies with ancient DNA analysis, to obtain a broader understanding of the evolutionary ecology of past populations. We uncovered ecological flexibility and geographic variation in woolly rhinoceros, and suggest spatial and temporal variation in isotopic composition was driven by environmental and vegetational factors.
Our analysis showed Late Pleistocene stability in d 15 N in NE Siberia, which we suggest reflects long-term environmental stability that may have favoured the later survival of woolly rhinoceros in the region.
Our comparative analysis with contemporaneous herbivores suggested possible niche partitioning of woolly rhinoceroses. We detected temporal variation in the isotopic profile of woolly rhinoceroses compared with other herbivores, which further supports the ecological flexibility of the species.
Ancient DNA analyses showed a lack of geographic structure at the mitochondrial level, with divergent lineages overlapping in time and space, as has been shown with a more spatially and temporally limited dataset from NE Siberia (Lord et al., 2020). We did not detect lineage loss prior to species extinction.
Stable isotope analyses have been used to investigate causes of extinction in other herbivores. For woolly mammoths, isotopic evidence suggested either a population decline (due to human encroachment) allowing other species (horse) to occupy their niche , or a niche change forced by climate-induced environmental change (Drucker et al., 2018). Our data did not support either of these scenarios in woolly rhinoceroses.
Radiocarbon chronologies and genetic data have correlated the extinction of woolly rhinoceroses with climatic and environmental factors (Lorenzen et al., 2011;Stuart and Lister, 2012;Lord et al., 2020). Demographic analysis based on genetic data has indicated an increase in population size~30 cal kyr BP followed by demographic stability until~18.5 cal kyr BP, where a decline in woolly rhinoceros population size started. Our results indicated that environmental stability enabled the late survival of woolly rhinoceroses in NE Siberia, and thus, environmental changes around the time of extinction may have had detrimental effects on the remaining populations. Pollen records for the region have documented a spread of shrub tundra communities and increased precipitation 13.5e12.7 cal kyr BP, compared to the open steppetundra community that was previously found in the region (until 15 cal kyr BP) (Müller et al., 2009). We suggest the interplay between environmental instability and fragmentation/isolation of populations, as presented in Lorenzen et al. (2011), played a major role in the extinction of the woolly rhinoceroses.

Funding
The work was supported by Villum Fonden Young Investigator Programme, grant no. 13151 to EDL.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.