When the woolly rhinoceroses roamed East Asia: a review of isotopic paleoecology of the genus Coelodonta from the Tibetan Plateau to northern Eurasia

1 Introduction

As the final representative of the woolly rhinoceroses, Coelodonta antiquitatis stands as an iconic species thriving on the mammoth steppe during the Late Pleistocene, with its distribution extending widely across northern Eurasia (Kahlke and Lacombat, 2008; Kahlke, 2014). However, East Asia has long been considered the cradle and evolutional center of the genus Coelodonta.

Deng et al. (2011) documented the earliest known ancestor of the woolly rhinoceroses, Coelodonta thibetana, originating from the Zanda Basin in the southwestern Tibetan Plateau at 5.08 Ma (Deng et al., 2011; Wang X. et al., 2013). This species then evolved into Coelodonta nihowanensis, which mainly thrived in northern China during the Early Pleistocene (Deng, 2002, 2008; Jin et al., 2021). C. nihowanensis in the Linxia Basin was considered to be the earliest woolly rhinoceros (2.5–2.2 Ma) prior to the discovery of the C. thibetana (Qiu et al., 2004; Deng, 2006; Zan et al., 2016; Deng et al., 2023). C. nihowanensis then dispersed out the Tibetan region, and it was commonly found in multiple sites in northern China, including the Nihewan Basin and Jinyuan Cave during the Early Pleistocene (e.g. Tong and Wang, 2014; Chen et al., 2017; Dong et al., 2021; Chen et al., 2023). C. nihowanensis gradually evolved into Coelodonta tologoijensis and was discovered in Russia and Mongolia during the early Middle Pleistocene. Its occurrence in Europe is still in debate (Kahlke and Lacombat, 2008; Uzunidis et al., 2022). Coelodonta then terminated in C. antiquitatis and its distribution extended over approximately 190° of longitude and 40° of latitude (Kahlke, 1999; Stefaniak et al., 2023). Through the evolutionary history of the woolly rhinoceroses, the genus Coelodonta was restricted in East Asia for more than 4 million years. Therefore, delving into the paleoecological shifts within Coelodonta in East Asia becomes essential for a deeper comprehension of the interplay between its evolutionary dynamics and the paleoenvironmental changes.

The series of morphological traits of Coelodonta, including special nasal horn, a low position of the head, and complex tooth morphology, were all considered indicators of grazing in cold and dry environment (Fortelius, 1983; Kahlke and Lacombat, 2008; Deng et al., 2011). As the most emblematic member of Coelodonta, the paleoecology of C. antiquitatis was intensively explored through various methods, and it was evidenced as a cold-adapted grazer (e.g. Rivals et al., 2010; Pushkina et al., 2020). The stable isotope analysis was applied on C. thibetana and C. nihowanensis in the Zanda and Nihewan basins. C. thibetana in the Zanda Basin fed on C₃ vegetation in an environment dominated by C₃ plants (Wang Y. et al., 2013), and C. nihowanensis in the Nihewan Basin displayed different foraging ecologies at different sites (Xu et al., 2021; Xu et al., 2023). The presence of C. nihowanensis in the Linxia Basin suggests a connection between the mid-Pliocene Zanda Basin and the Early Pleistocene Nihewan Basin (Deng et al., 2023). This provides an opportunity to explore the paleoecological changes of Coelodonta as they migrated out of the Tibetan Plateau. The previous study explored the paleoecology of Cenozoic rhinoceroses in the Linxia Basin; however, only one C. nihowanensis specimen was analyzed and identified as a grazer in the Longdan fauna (Biasatti et al., 2018). It is evident that more specimens are required to investigate the dietary ecology of C. nihowanensis in the Linxia Basin. Furthermore, all isotopic data from East Asia should be discussed from the perspective of the Coelodonta evolution.

In this study, we performed stable carbon and oxygen isotope analyses on C. nihowanensis and other ungulates from the Longdan fauna to investigate the paleoecology of C. nihowanensis in the Linxia Basin during the Early Pleistocene. Additionally, we conducted an extensive review of all published isotopic data pertaining to Coelodonta and coexisting mammals, aiming to provide a deeper understanding of the paleoecological shifts experienced by Coelodonta in East Asia.

2 Carbon and oxygen isotopes in tooth enamel

Higher plants are generally divided into C₃ (including virtually all trees, most shrubs and herbs, and cool-season grasses), C₄ (warm season grasses), and CAM (succulents) plants according to different photosynthesis pathways (Troughton and Card, 1975; O’Leary, 1981; Farquhar et al., 1989). Modern C₃ plants have an average δ¹³C value of −25‰, while C₄ plants have an average δ¹³C value of −13‰ (Bender, 1968; Smith and Epsten, 1971; Kohn and Cerling, 2002). Additionally, water use efficiency, the canopy effect, natural physiology, and abiotic factors can impact the δ¹³C values of plant tissues (van der Merwe and Medina, 1991; Quade et al., 1995; Peuke et al., 2006; Kohn, 2010). When plants are directly consumed by mammals, isotopic fractionation occurs during the process of tissue formation. The enrichment factors between plant consumption and tooth enamel bioapatite are roughly 14‰ for herbivores, although they may be slightly influenced by variations in their body sizes and digestive strategies (Ambrose and Norr, 1993; Cerling and Harris, 1999; Kohn and Cerling, 2002; Tejada-Lara et al., 2018; Cerling et al., 2021). A tooth enamel with a δ¹³C value of less than −9‰ generally indicates a pure C₃ diet, while a δ¹³C value greater than −2‰ would demonstrate a pure C₄ diet (Cerling et al., 1997). Nevertheless, in water-stressed environments, a conservative enamel-δ¹³C value cut-off for a pure C₃ diet in fossil herbivores should likely be −7‰ or even −6.5‰ (Cerling et al., 1997; Wang et al., 2008a). Therefore, the δ¹³C values of mammals are valuable for reconstructing their paleodiet and paleoecology. Nonetheless, recent researches have uncovered that the relationships between herbivore assemblages and the vegetation landscapes are complicated (Negash and Barr, 2023; Sokolowski et al., 2023), underscoring the importance to avoid simplistic interpretations of paleolandscapes based solely on the δ¹³C values of fossil herbivores.

The δ¹⁸O values in the bioapatite of the mammal bones and teeth are primarily determined by the source of water consumed and their foraging ecology (Bryant and Froelich, 1995; Pederzani and Britton, 2019). For large herbivores relying heavily on drinking water, the δ¹⁸O values derived from tooth enamel directly reflect the δ¹⁸O values they consume, which can vary due to factors, such as local temperature, precipitation, and humidity (Dutton et al., 2005; Pederzani and Britton, 2019). In general, tree leaves exhibit higher δ¹⁸O values than grasses due to increased transpiration of water (Koch et al., 1989; Quade et al., 1995). Consequently, herbivores relying on such foliage typically demonstrate higher δ¹⁸O values than those consuming grasses.

Mammalian teeth form incrementally from the top to the bottom of the tooth crown during the period of tooth growth (Hillson, 2005). Once tooth enamel is fully mineralized, it remains unaltered throughout the lifespan of mammals, thereby preserving the isotopic sequence formed during its development (Suga, 1979). Consequently, bulk sampling of the entire height of tooth enamel can be used to reconstruct the average diet and ecology during the tooth formation period, while serial sampling of the tooth enamel along the direction of enamel growth can offer detailed insights into the dietary and ecological changes of mammals (Kohn et al., 1998; Balasse, 2002; Blumenthal et al., 2014, Blumenthal et al., 2019; Li et al., 2020). However, enamel mineralization progresses in a non-linear manner, following distinctive trajectories within the enamel layers that do not consistently align with incremental growth lines (Suga, 1979; Smith et al., 2003). Therefore, conventional serial sampling method often yields composite samples from multiple enamel layers, which may potentially dampen seasonal signals (Sharp and Cerling, 1998; Higgins and MacFadden, 2004; Green et al., 2018).

3 Materials and methods

3.1 Geological settings and sample selection

The Linxia Basin, located on the northeastern edge of the Tibetan Plateau, is well known for its thick and continuous Late Cenozoic deposits and its exceptionally rich and diverse mammalian fossils (Deng et al., 2013). The Longdan, the first Early Pleistocene site in China discovered from loess, belongs to the Wucheng Formation, characterized by yellowish hard loess with layers of calcite concretions (Qiu et al., 2002; Qiu et al., 2004). Magnetostratigraphic dating had determined the age of the Longdan fauna to be 2.5–2.2 Ma (Zan et al., 2016). This fauna has yielded abundant fossils of 33 mammalian species, including Equus eisenmannae, Proboscidipparion sinense, C. nihowanensis, Leptobos brevicornis, and a diverse range of carnivores (Qiu et al., 2004). The mammalian fossils in the Longdan fauna exhibit slightly more primitive features than those of the Early Pleistocene Nihewan fauna (Qiu et al., 2004). The Nihewan fauna represents another significant Early Pleistocene mammalian assemblage in northern China, preserving abundant C. nihowanensis fossils, and its geological age are 2.4–1.8 Ma or 2.2–1.7 Ma (Qiu, 2000; Liu et al., 2012; Farjand et al., 2023a; Farjand et al., 2023b). E. eisenmannae in the Longdan fauna was one of the earliest and most primitive Equus in Eurasia, while C. nihowanensis in this fauna was initially regarded as the most primitive woolly rhinoceros before the discovery of the C. thibetana (Qiu et al., 2004; Deng, 2006).

3.2 Materials

27 fossil tooth samples were collected from the Longdan fauna, including the woolly rhinoceros C. nihowanensis (n=7), the deer Nipponicervus longdanensis (n=4), Bovinae (L. brevicornis and Hemibos gracilis) (n=7), and the gazelle Gazella cf. blacki (n=9). These specimens were housed in the Institute of Vertebrate Paleontology and Paleoanthropology and the Hezheng Paleozoological Museum in China. However, the precise stratigraphic levels of these tooth specimens were not recorded. Two modern goats were collected as modern references. Among them, 8 fossil samples and 2 modern goats were serially sampled. In total, 103 specimens were collected.

3.3 Stable isotope analysis

First, the teeth were cleaned by removing a very thin (~0.1 mm) layer of outer enamel using a diamond-tipped dental burr. For bulk samples, enamel powder was collected from the entire range of tooth heights to obtain an average value. Well-preserved teeth with complete tooth crown height were chosen for serial sampling. Each serial sample was obtained by drilling with a 1–2 mm-wide groove mostly perpendicular to the tooth growth axis, penetrating through the enamel layer’s thickness at intervals of approximately 1–2 mm. These sampling methods were designed to encompass the full crown height of each tooth piece. Approximately 12–15 mg of enamel powder was drilled and collected in a 1.5-ml centrifuge tube for each sample.

Second, all the enamel powder underwent pretreatment with 2.5% NaOCl followed by 1 mol/L acetic acid. Subsequently, 3–5 mg of powdered enamel was analyzed at the Laboratory for Stable Isotope Geochemistry, Institute of Geology and Geophysics, Chinese Academy of Sciences using a MAT253 IRMS coupled with a Gasbench unit. The international isotope standards of IAEA CO-8 and NBS-19 were interspersed with every run to establish the standard curve for isotopic calibration. Additionally, the in-house standard GBW04405 (limestone) and sample replicates were analyzed to monitor measurement uncertainty. Carbon and oxygen isotope values were expressed as δ¹³C and δ¹⁸O values relative to VPDB. The analytical precision, based on long-term repeated measurements of the standards and sample replicates, was less than ± 0.15‰ for δ¹³C values and ± 0.20‰ for δ¹⁸O values.

Assessment of fossil tooth preservation in the Longdan fauna, as established in previous studies (Biasatti et al., 2010; Ma et al., 2023), suggests that the specimens in this study from the same fauna and sediments are likely well preserved.

3.4 Data analyses

We utilized the carbon isotope enrichment factors proposed in the most recent study based on modern animal experiments (Cerling et al., 2021): 13.5‰ for non-coprophagous hindgut fermenter (such as equids, woolly rhinoceroses, and mammoths) and 14.5‰ for ruminant foregut fermenter (including cervids, gazelles, and bovids). For the collagen samples from the Late Pleistocene, we applied 5‰ as the isotope enrichment factor (Cerling and Harris, 1999).

The isotopic data collected from different sites in this study cover a range of regions, necessitating the consideration of factors such as latitude and altitude. Previous studies have identified a positive relationship between elevation and the δ¹³C values of C₃ plants, possibly influenced by air pressure and temperature (Körner et al., 1991; Männel et al., 2007). Conversely, the latitudinal gradient for the δ¹³C values remains less definitive (Krajcarz et al., 2016). The positive correlation between elevation and the δ¹³C values of C₃ plants has been noted in the modern Tibetan region as well (Wang et al., 2008b), which is pertinent to our studied locations. However, the relationship between δ¹³C values of C₄ plants and the elevation gradients is unclear and may exhibit a negative correlation in the modern Tibetan region (Wang et al., 2008b). Additionally, no clear relationship was observed between δ¹⁸O values of precipitation and altitude in this region (Wang et al., 2008b).

Therefore, the relationship between δ¹³C values of C₃ plants and altitude appears to be an important factor in this study when analyzing isotopic values across sites at different elevations. We adopted the altitude gradient as per Männel et al. (2007): δ¹³C-adj-alt= δ¹³C -(0.0011×altitude), where altitude in meters. By applying this correction, we mitigate the impact of altitude, normalizing all data to sea level (0 meter above sea level, asl) for comparability. However, due to the unclear association between δ¹³C values of C₄ plants and elevation, we chose not to adjust data from the regions where C₄ plant components were notably present.

Mann-Whitney U tests were conducted to assess statistical variations in isotopic values across distinct groups. These analyses were conducted in PAST 4.03 (Hammer et al., 2001). Data visualizations were created in EXCEL and RStudio (RStudio Team, 2016) using ggplot2 (Wickham, 2011) and ggpubr (Kassambara, 2020).

4 Results and interpretations

4.1 Bulk carbon and oxygen isotope analyses

The δ¹³C and δ¹⁸O values for each tooth sample (bulk samples and the mean values of each serial sample) are summarized in Table 1 and displayed in Figure 1 (see original data in Supplementary Table S1). In order to capture a complete paleoecological scenario of the Longdan fauna, our previous data of the E. eisenmannae and P. sinense were also included (Ma et al., 2023). Overall, the δ¹³C values of all the specimens in this study ranged from −13.3‰ to −7.9‰, and their δ¹⁸O values ranged between −12.5‰ to −4.6‰. In Figure 1, δ¹³C and δ¹⁸O values of the Longdan equid samples were also fall into these ranges (Ma et al., 2023). This indicated that C₃ vegetation constituted a significant portion of the ungulates’ diet and played a substantial role in the local environment during the Early Pleistocene in the Linxia Basin.

Table 1 Summary of the δ¹³C and δ¹⁸O values of the Longdan mammalian tooth enamel samples.

Figure 1 Stable isotopic data of the bulk samples (blank symbols) and the mean value of each serial sample (solid symbols) of the studied mammals (including data of P. sinense and E. eisenmannae from Ma et al., 2023).

The δ¹³C values of the four ungulates in this study were slightly lower than those of two equid species, which indicated that equids occupied more open environments than other mammals (Table 1; Figure 1). The results of the Mann-Whitney U tests revealed that the δ¹³C values of G. blacki were lower than P. sinense and E. eisenmannae (p-values = 0.015, 0.009), suggesting a possible preference for close environments in G. blacki compared to equids. Overall, four ungulates analyzed in this study exhibited subtle differences in δ¹³C values (Supplementary Table S2), indicating potential overlap in their food resources and habitats. The δ¹³C standard deviation was higher for C. nihowanensis specimens in comparison to other mammals (Table 1), indicating the diverse habitats they inhabited. In a previous study, one C. nihowanensis specimen with a δ¹³C value of −9.9‰ was classified as a grazer (Biasatti et al., 2018), within the range of δ¹³C values observed in our C. nihowanensis specimens (from −13.3‰ to −8.8‰, n=7).

As shown in Table 1, the mean δ¹⁸O values of the G. blacki and N. longdanensis specimens were higher than those of Bovinae and E. eisenmannae specimens, whereas the specimens of C. nihowanensis and P. sinense exhibited intermediate δ¹⁸O values. The results from the Mann-Whitney U tests revealed that the δ¹⁸O values of E. eisenmannae were lower than those of N. longdanensis and C. nihowanensis (refer to Supplementary Table S3), with no discernible differences were observed in other groups. This suggested that E. eisenmannae specimens likely rely more heavily on meteoric water resources (Sponheimer and Lee-Thorp, 1999), while N. longdanensis and C. nihowanensis may access different water resources, potentially obtaining a portion of their water from tree leaves. In general, tree leaves exhibit higher δ¹⁸O values than grasses because trees experience greater transpiration of water through the leaves (Koch et al., 1989; Quade et al., 1995). Meanwhile, high standard deviations in δ¹⁸O values were observed in the specimens of Bovinae, G. blacki, and P. sinense, suggesting potential access to a diverse range of water resources (Table 1; Figure 1). Nevertheless, these varied δ¹⁸O values may be caused by the potential differences in the chronology of these specimens.

A discrepancy was noted in the isotopic results between the averaged serial samples and the bulk samples from the specimens of N. longdanensis, P. sinense, and Bovinae. This difference could be linked to potential sampling bias inherent in the two distinct sampling methods. Despite the bulk samples being designed to encompass the entire height of the tooth, challenges arise when dealing with fragmented or heavily worn teeth. Consequently, the stable isotopic values obtained from partial tooth heights could introduce bias and fail to accurately reflect the averaged isotopic values. Nevertheless, considering the subtle intra-tooth isotopic ranges in these serial samples, we propose that the data disparity probably reflects their varied isotopic ecologies. Furthermore, the potential differences in chronology could be another contributing factor to this data discrepancy, given the absence of precise stratigraphic levels for these tooth specimens.

4.2 Serial isotopic analysis

The serial isotopic patterns observed in these mammals indicated notable signals of the ecological and environmental seasonality. All eight fossil teeth examined in this study exhibited stable δ¹³C values, with intra-tooth δ¹³C values ranging from 0.8‰ to 1.6‰ (Supplementary Table S1). This possibly indicated limited seasonal variations in the local vegetation. By comparison, the serial δ¹⁸O profiles of these ungulates displayed interspecific differences (Figure 2). The intra-tooth δ¹⁸O profiles of the large Bovinae (Supplementary Table S1; Figure 2) and equids (Biasatti et al., 2010; Ma et al., 2023) specimens demonstrated stable and monotonous patterns. The seasonal δ¹⁸O values observed in these large obligate drinkers may have been dampened due to their frequent consumption of water from large water bodies, which mitigated seasonal isotopic variations in local rainfall, owning to prolonged water residence times (Roberts et al., 2018; Pederzani and Britton, 2019). Alternatively, the stability and monotony observed in the intra-tooth δ¹⁸O profiles of the equid and bovid specimens might suggest migratory behavior. Consequently, their δ¹⁸O values may reflect a combination of water sources from various locations.

Figure 2 Serial isotopic profiles of the fossil mammals in this study. (Δ¹³C and Δ¹⁸O are the ranges of the intra-tooth δ¹³C and δ¹⁸O values).

The intra-tooth δ¹⁸O ranges of three G. blacki specimens were 3.5‰, 3.5‰, and 4.8‰, exhibiting greater variability compared to those of the Bovinae (1.7‰, 1.3‰, and 1.4‰) and C. nihowanensis (1.7‰, 2.1‰) specimens. The intra-tooth δ¹⁸O ranges of the C. nihowanensis specimens were slightly larger than those of the Bovinae specimens in this study and equid specimens (0.7‰, 1.4‰, 1.5‰, and 2.3‰) in Ma et al. (2023). Regular sinusoidal patterns were observed in the serial δ¹⁸O profiles of the C. nihowanensis specimens (Figure 2). This suggested that C. nihowanensis specimens possibly obtain fewer water resources from the longstanding water bodies, unlike the large Bovinae and equids (Biasatti et al., 2010; Ma et al., 2023). The distinct water resources of the C. nihowanensis specimens, compared to the Bovinae and equids, previously inferred from their higher δ¹⁸O values (Table 1; Supplementary Table S3), were consistent with this finding inferred from the different serial isotopic profiles. Alternatively, C. nihowanensis specimens might be non-migratory and thus could record the δ¹⁸O signals of local environments for more than one year (Figure 2).

Furthermore, an inverse relationship between the intra-tooth δ¹³C and δ¹⁸O values were observed for two G. blacki specimens (LD-29, 31; r =−0.78, −0.92) and one middle-sized Bovinae specimen (LD-33, r =−0.53).

In Supplementary Figure S1, we plotted the serial isotope profiles of two modern goats from the different regions in the Linxia Basin and its vicinity with varying altitudes. The average δ¹³C values of these two samples were −10.1± 2.8‰ (n=6) and −11.9± 0.9‰ (n=5), respectively; and their mean δ¹⁸O values were −5.2± 1.4‰ and −7.1± 1.9‰. The intra-tooth range of δ¹³C values for one goat specimen was higher, and its serial δ¹³C values ranged from −13.4‰ to −6.3‰, clearly indicating seasonal C₄ food resources. The intra-tooth δ¹⁸O ranges of these two goats were 4.5‰ and 3.6‰, slightly higher than those of three fossil G. blacki specimens. Additionally, inverse relationships between δ¹³C and δ¹⁸O values were observed for both goat specimens (r = −0.83, −0.86).

The isotopic results of the modern goats indicated the occurrence of East Asia Summer Monsoon (EASM) and the presence of seasonal C₄ plants. For both specimens, the peaks of δ¹³C values indicated the ingestion of C₄ plants during the summer, and their corresponding δ¹⁸O values were the lowest. In a monsoonal region, the “amount effect” affects seasonal δ¹⁸O values of precipitation, preferentially removing ¹⁸O from vapor, resulting in lower δ¹⁸O values in summer (Dansgaard, 1964; Rozanski et al., 1993). High precipitation and warm temperatures favor C₄ plants, which are naturally more abundant in summer. Therefore, an inverse relationship between δ¹³C and δ¹⁸O values is expected (Higgins and MacFadden, 2004; Wang et al., 2008b; Biasatti et al., 2010; Yang et al., 2012; Zhang et al., 2012; Ciner et al., 2015).

In Figure 2 and Supplementary Figure S1, the serial isotopic profiles of the two modern goats were similar to those of a part of fossil teeth in this study, especially the G. blacki specimens. Through these modern references, this study provides evidence once again that the EASM occurred in the Linxia Basin during the Early Pleistocene (Biasatti et al., 2010). Furthermore, it appears that non-obligate drinkers offer better indicators of seasonal variations in local environments than large obligate drinkers (Roberts et al., 2018; Pederzani and Britton, 2019; Hodgkins et al., 2020).

Based on this and previous studies (Wang and Deng, 2005; Biasatti et al., 2010; Ma et al., 2023), the paleoenvironment in the Linxia Basin, as inferred from mammalian carbon isotopes during the Early Pleistocene was primarily characterized by C₃ plants, with C₄ plants likely representing a minor component that was seasonally consumed by the studied mammals. This finding aligns with the isotopic results of the carbonates from fluvial and lacustrine sediments (Dettman et al., 2003; Hough et al., 2011). The sub-humid climate and the occurrence of EASM are further supported by previous studies, including stable isotope analysis (Biasatti et al., 2010), fossil assemblage analysis (Qiu et al., 2002; Wang, 2005; Xue et al., 2006), and sedimentary analysis (Zan et al., 2016).

5 Discussion

5.1 Foraging ecologies and environmental contexts of C. nihowanensis in Early Pleistocene northern China

Studying the foraging ecologies of C. nihowanensis in northern China from west to east can provide valuable insights into the relationship between their ecological behaviors and paleoenvironments during the Early Pleistocene. This study provides a comprehensive insight into the paleoecology of the earliest C. nihowanensis in northwestern China, situated on the edge of the Tibetan Plateau. Within the Linxia Basin, C. nihowanensis specimens were likely a mixed feeder inhabiting diverse habitats, as indicated by their wide distribution of the δ¹³C values and intermediate δ¹⁸O values compared to other coexisted mammals. Meanwhile, recent stable isotopic studies on mammalian faunas in the Nihewan Basin aiming to explore the paleoenvironmental settings of early humans can serve as a strong evidence for the paleoecology of C. nihowanensis in eastern China (Xu et al., 2021; Xu et al., 2023). Therefore, it is a great chance to compile all these data and explore the paleoecological changes of the C. nihowanensis (Figure 3; Supplementary Table S4). Besides, the same pretreatment methods of the stable isotopic analysis for these two faunas in the same laboratory could guarantee the validity of the cross comparison. Considering the potential bias in δ¹⁸O values stemmed from variations in phosphoric acid reaction temperatures (Passey et al., 2007) during the isotopic measurements conducted in different laboratories, the isotopic data from Wang and Deng (2005) as well as Biasatti et al. (2010, 2018) were deliberately not included in Figure 3. It is worth noting, however, that the intra-tooth isotopic ranges from these studies have been collated and are presented in Figure 4 (Supplementary Table S5). The altitude of the Linxia and Nihewan basins during the Early Pleistocene were already situated at their current elevations and therefore, δ¹³C_diet values were corrected with the modern altitudes of 2248 m and ~940 m (Deng et al., 2013; Pei et al., 2019), respectively. δ¹³C_diet-adj-alt were corrected using the formula described in section 3.4. This correction was applied at the Longdan and Shanshenmiaozui sites, where the ungulates’ diets predominantly consisted of almost pure C₃ plants. However, in light of the significant components of the C₄ plants in the diet of these ungulates, this correction was not applied in other three sites during ca. 1.2–1.1 Ma in the Nihewan Basin.

Figure 3 δ¹³C_diet-adj-alt and δ¹⁸O values of the Longdan fauna and Nihewan fauna (~2.2–1.7 Ma: Shanshenmiaozu; ca. 1.2–1.1 Ma: Madigou, Feiliang, and Donggutuo) (Xu et al., 2021, 2023; Ma et al., 2023; this study) (Supplementary Table S4). (A) Scatter plot of δ¹³C_diet-adj-alt and δ¹⁸O values of different taxa. (B) Violin plot and boxplot of δ¹³C_diet-adj-alt values of all the mammalian species. (C) Violin plot and boxplot of δ¹⁸O values of all the mammalian species. (These data consist of bulk samples and the mean values of each serial sample. For panels B and C, the width of the “violin” at a particular value reflects the data density at that point. Within the boxplot, the lower hinge, median line, and upper hinge denote the 25th percentile, median, and the 75th percentile, respectively. The lower whisker extends to the smallest observation at or above the first quartile (Q1) minus 1.5 times the Interquartile Range (IQR), while the upper whisker extends to the largest observation at or below the third quartile (Q3) plus 1.5 times the IQR. The potential outliers in the dataset are marked by separate points).

Figure 4 Intra-tooth δ¹³C (A) and δ¹⁸O ranges (B) of the serial samples in the Linxia Basin (Longdan) and Nihewan Basin (Shanshenmiaozu, ~2.2–1.7 Ma; Madigou, Feiliang, and Donggutuo are ca. 1.2–1.1 Ma) (Biasatti et al., 2010; Xu et al., 2021, 2023; Ma et al., 2023; this study) (Supplementary Table S5). (The boxplot conveys the same meaning as Figure 3, but without annotating the outliers to avoid redundancy with the original data points represented by the jitter scatter plot).

Among the four sites in the Nihewan Basin, Shanshenmiaozui has yielded the richest mammalian fossils. Based on the analysis of the faunal assemblage, the age of the Shanshenmiaozui fauna likely corresponds closely with the classic Nihewan fauna during the early Early Pleistocene, with an approximate age of 2.2–1.7 Ma (Tong et al., 2021). The other three sites–Madigou, Donggutuo, and Feiliang–are classified as Early Paleolithic, with the age of ca. 1.2–1.1 Ma, late Early Pleistocene (Wang et al., 2005; Deng et al., 2007; Pei et al., 2019).

In Figure 3, the δ¹³C_diet-adj-alt values of the Shanshenmiaozui fauna (−25.9± 1.6‰, n=65) were higher than those of the Longdan fauna (−27.1± 1.4‰, n=39) (p-value of Mann-Whitney U test: 0.0005) (Ma et al., 2023; Xu et al., 2023; this study). However, such variation may not be directly linked to the difference in landscape between these two sites (Negash and Barr, 2023). The δ¹⁸O values of the Shanshenmiaozui fauna (−9.8± 1.4‰, n=65) were slightly lower than those of the Longdan fauna (−8.8± 1.7‰, n=39) (p-value of Mann-Whitney U test: 0.001), possibly indicating a wetter and/or more humid environment in the Nihewan Basin during the early Early Pleistocene (Ma et al., 2023; Xu et al., 2023; this study).

The intra-tooth δ¹³C and δ¹⁸O ranges of the single equid specimen in Shanshenmiaozui were higher than those of the Longdan equid specimens (Figure 4) (Biasatti et al., 2010; Ma et al., 2023; Xu et al., 2023; this study). This might suggest a more pronounced climatic seasonality in the Nihewan Basin compared to the Linxia Basin, despite the limited serial sample amount in Shanshenmiaozui. Nonetheless, it seems that C₄ grasses may not have become a significant component of the local ecosystems in the Nihewan Basin during this period.

In the Shanshenmiaozui fauna, C. nihowanensis and bovids exhibited intermediate δ¹³C_diet-adj-alt values, which were lower than those of equids but higher than cervid ones (Figure 3; Supplementary Table S6) (Xu et al., 2023). This indicated that C. nihowanensis and bovid specimens in Shanshenmiaozui likely inhabited mixed grasslands and woodlands. Statistical analyses displayed no differences in the δ¹⁸O values among these groups (Supplementary Table S7). However, δ¹⁸O values were lowest for C. nihowanensis among all the coexisting mammals (Figure 3), possibly indicating that these specimens obtained abundant meteoric water and possibly inhabit areas close to lakes and/or rivers (Xu et al., 2023).

Until the late Early Pleistocene, the mean δ¹³C_diet-adj-alt and δ¹⁸O values of the three sites during ca. 1.2–1.1 Ma were higher than those of the Shanshenmiaozui fauna (p-values of Mann-Whitney U tests for δ¹³C and δ¹⁸O: 0.0005, 0.0001) (Figure 3) (Xu et al., 2021; Xu et al., 2023). This indicated that the C₄ vegetation became significant in the Nihewan Basin during the late Early Pleistocene. The intra-tooth ranges in δ¹³C and δ¹⁸O values of the equid and bovid specimens at these three sites were greater than those of the Longdan and Shanshenmiaozui specimens (Figure 4; Supplementary Table S5) (Biasatti et al., 2010; Xu et al., 2021; Ma et al., 2023; Xu et al., 2023; this study), possibly indicating seasonal ingestion of a substantial amount of C₄ plants in their diets and more pronounced seasonal fluctuations in precipitation. An inverse relationship between δ¹³C and δ¹⁸O values was observed in some serial equid and bovid specimens in Donggutuo and Feiliang sites (Supplementary Table S5) (Xu et al., 2023), which was the signal of the EASM as aforementioned. Therefore, the paleoenvironments in the Nihewan Basin during the late Early Pleistocene were possibly varied and experienced pronounced seasonality under the effect of intensive EASM and the winter monsoon (Xu et al., 2021; Xu et al., 2023).

The δ¹³C_diet-adj-alt values of the C. nihowanensis specimens during ca. 1.2–1.1 Ma were also lower than those of equid and bovid specimens but higher than cervid ones (Figure 3; Supplementary Table S8) (Xu et al., 2021; Xu et al., 2023), indicating that C. nihowanensis specimens fed on mixed C₃/C₄ vegetation, and consumed less C₄ grasses than equids and bovids. Furthermore, C. nihowanensis specimens exhibited the lowest δ¹⁸O values (Figure 3; Supplementary Table S9), and the intra-tooth ranges in δ¹³C and δ¹⁸O values of C. nihowanensis were found to be smaller than those of the equids and bovids (Figure 4; Supplementary Table S5) (Xu et al., 2021; Xu et al., 2023). This possibly indicated that C. nihowanensis occupied relatively close and well-watered environments which experienced fewer seasonal changes. Therefore, C. nihowanensis in the Nihewan Basin during the late Early Pleistocene was likely a mixed feeder and showed a preference for close environments (Xu et al., 2021; Xu et al., 2023).

During the early Early Pleistocene, the paleoenvironments in the Linxia and Nihewan basins were primarily dominated by C₃ vegetation, with the latter likely exhibited a stronger EASM compared to the Linxia Basin. Previous studies found that C₄ signals were more abundant in the Linxia Basin at ~1.2 Ma (Wang and Deng, 2005; Biasatti et al., 2010). The increasing prevalence of C₄ plants and a pronounced EASM in the Linxia and Nihewan basins around 1.2–1.1 Ma potentially aligns with the widespread aridification in northern China (Da et al., 2015; Zhou et al., 2018). Against this backdrop, C. nihowanensis owned a notable dietary flexibility and were able to thrive in a variety of environments across northern China from west to east throughout the Early Pleistocene. Such ecological adaptability of C. nihowanensis is presumed to have been influenced by the intricate terrain of Early Pleistocene East Asia and its diverse environments, which were situated along the migration path of the woolly rhinoceroses from the Tibetan Plateau to northern Eurasia.

5.2 Ecological adaptation of Coelodonta and coexisting mammals

After thoroughly examining the isotopic paleoecology of C. nihowanensis in Early Pleistocene East Asia, it is imperative to delve deeper into the paleoecology of its ancestral and derived species to capture the paleoecological evolution of Coelodonta. We complied all δ¹³C_diet-adj-alt values of the Coelodonta genus and other coexisting mammals in East Asia from the mid-Pliocene to Late Pleistocene in Figure 5 (Supplementary Table S4). During the mid-Pliocene, it was estimated that the elevation of the Zanda Basin was similar to its modern altitude, approximately 4,000 meters asl (Deng et al., 2012). Given that the modern altitudes of these Late Pleistocene sites were around 100–200 meters, it is assumed that the altitude effect could be disregarded in this context.

Figure 5 δ¹³C_diet-adj-alt values of the Coelodonta genus from the mid-Pliocene to Late Pleistocene (Wang and Deng, 2005; Wang Y. et al., 2013; Ma et al., 2017; Biasatti et al., 2018; Ma et al., 2021; Xu et al., 2021; Ma et al., 2023; Xu et al., 2023; this study) (the box plot conveys the same meaning as Figure 3).

During the mid-Pliocene, the δ¹³C_diet-adj-alt values of C. thibetana in the Zanda Basin ranged narrowly between −27.7‰ to −26.0‰ (n=22), with their mean δ¹³C_diet-adj-alt (−27.0± 0.4‰) being higher than those of the bovid (−28.6± 1.3‰, n=19), equid (−27.6± 0.5‰, n=30), and cervid (−28.0± 0.1‰, n=3) specimens (Supplementary Table S10; Figure 5) (Wang Y. et al., 2013). Additionally, both carbon and oxygen isotopes of the C. thibetana specimens showed lower standard deviations than those of the bovids and equids (Figure 5; Wang Y. et al., 2013). This suggested that the ancestral C. thibetana primarily fed on C₃ grasses and occupied a confined ecological niche in open environments at the high-altitude Zanda Basin. C. thibetana occupies the most basal position of the Coelodonta lineage and was considered the earliest representative of this genus (Deng et al., 2011). The morphological traits of the C. thibetana, such as the highly elongated skull, large and roughened nasal horn base, and high tooth crowns with rough enamel, are indicative of an adaptation to efficiently graze in cold and arid environments. These were in consistent with its foraging ecology inferred from stable isotopes.

C. thibetana moved down from the Tibetan Plateau toward northeastern Linxia Basin (~2,500 m asl) at the initiation of the Quaternary and evolved as C. nihowanensis. The average δ¹³C_diet-adj-alt value of the Longdan fauna (−26.9± 1.4‰, n=55) (Wang and Deng, 2005; Biasatti et al., 2018; Ma et al., 2023; this study) was higher than that of the Zanda fauna (−27.7± 1.0‰, n=74) (p-value of Mann-Whitney U test: 0.0008) (Wang Y. et al., 2013), possibly indicating a drier/warmer climate in the Early Pleistocene Linxia Basin. The broader range of δ¹³C_diet-adj-alt values observed in the C. nihowanensis specimens in the Linxia Basin (Biasatti et al., 2018; this study) compared to coexisting mammals, in contrast to the narrower range of the δ¹³C_diet-adj-alt values seen in the ancestral C. thibetana in the Zanda fauna, suggested that C. nihowanensis diversified its ecological niche during this period. In the Zanda and Longdan faunas, both woolly rhinoceroses and equids had higher δ¹³C_diet-adj-alt values than those of artiodactyls (Figure 5) (Wang Y. et al., 2013), suggesting a possible ecological niche partitioning between artiodactyls and perissodactyls.

In Section 5.1, a thorough discussion has been conducted on the paleoenvironmental and paleoclimatic changes in the Nihewan Basin during the Early Pleistocene. In the Shanshenmiaozui fauna, equids occupied open grasslands and exhibited the highest δ¹³C_diet-adj-alt values (−24.8± 1.1‰, n=18) among all the mammalian species (Supplementary Table S6) (Xu et al., 2023). Similar δ¹³C_diet-adj-alt values were observed for woolly rhinoceroses (−25.9± 1.2‰, n=21) and bovids (−25.7± 1.7‰, n=13) in this fauna (Supplementary Table S6), indicating similar habitat preferences for both species (Xu et al., 2023). The lowest δ¹³C_diet-adj-alt values (−27.8± 0.8‰, n=13) of the cervid specimens reflected a preference for browsing in a more close environment (Supplementary Table S6) (Xu et al., 2023). Therefore, C. nihowanensis and bovids in Shanshenmiaozui were possibly mixed feeders and likely inhabited relatively open environments in Shanshenmiaozui. These ecological niches within the Shanshenmiaozui fauna inferred from the δ¹³C_diet-adj-alt values were similar to the pattern observed in the Longdan fauna.

In the Feiliang, Madigou, and Donggutuo sites, bovids and equids exhibited much higher δ¹³C_diet-adj-alt values (−19.3± 2.6‰, n=29; −19.6± 1.9‰, n=51) than those of the woolly rhinoceroses (−21.7± 1.6‰, n=34) and cervids (−25.6± 2.0‰, n=13) (refer to Supplementary Table S8) (Xu et al., 2021; Xu et al., 2023). This suggested that the former consumed a greater amount of C₄ grasses and inhabited more open environments compared to the latter during the late Early Pleistocene. C. nihowanensis likely had a mixed diet and favored a relatively close habitat during ca. 1.2–1.1 Ma. This contrasted slightly with the preference for open environments observed in earlier C. nihowanensis specimens from Shanshenmiaozui within the same region. Consequently, it seems that C. nihowanensis in the Nihewan Basin displayed varying ecological preferences during different periods of the Early Pleistocene.

The ecological niche partitioning among these mammalian taxa during the late Early Pleistocene differed from the pattern observed in previous faunas, indicating notable ecological shifts among these mammalian taxa in eastern China. Throughout the mid-Pliocene to late Early Pleistocene, equids consistently favored open environments, reflecting a grazing preference similar to their modern counterparts (Ransom and Kaczensky, 2016). Within the Nihewan Basin, bovids exhibited a highly diverse foraging ecology, gradually expanding into open environments compared to the Zanda and Linxia basins (Xu et al., 2021; Xu et al., 2023). However, the habitats of the C. nihowanensis slightly shifted towards more close environments in the Nihewan Basin (Xu et al., 2021; Xu et al., 2023). This suggested distinct foraging strategies between bovids and woolly rhinoceroses in response to changing paleoenvironments.

C. nihowanensis is estimated to have a larger body size (approximately 1700 kg according to the NOW database¹) compared to the ancestral C. thibetana (~1470 kg¹), and its morphological traits were also considered as a grazing adaptation in cold and dry environments. However, the foraging ecologies of C. nihowanensis were more flexible and even displayed a preference for close environments. Additionally, the limb bones of C. nihowanensis are more slender than those of the derived C. tologoijensis, indicating a better running ability (Deng, 2008). However, the isotopic results of the C. nihowanensis specimens from both Linxia and Nihewan basins likely exhibited a potential of non-migratory behaviors as aforementioned. This highlights the distinction between individual behavioral selection and population dispersal as well as evolutionary history.

C. tologoijensis (~2000 kg¹) has a remarkably larger body size than the primary C. nihowanensis, and its paleoecology was rarely explored. According to the paleoenvironmental records of the sites in Asia which preserved C. tologoijensis fossils, Kahlke and Lacombat (2008) presumed that C. tologoijensis possibly inhabited meadow to forest steppes, similar to those found in northern Mongolia today.

C. antiquitatis first occurred during the late Middle Pleistocene and thrived in northern Eurasia during the Late Pleistocene as an emblematic member in the Mammuthus-Coelodonta Faunal Complex, and its paleoecology was intensively explored. C. antiquitatis was graviportal and has a massive body size (~2900 kg¹), as a highly progressive species of Coelodonta (Boeskorov, 2001; Stuart and Lister, 2012). A series of anatomical features observed in C. antiquitatis, including wide lips, low posture of the head, long nasal horn, high tooth crown, long hair, and dense underwool, all indicating an adaptation for grazing in dry and frosty environments (Kahlke and Lacombat, 2008; Boeskorov et al., 2011b; Boeskorov, 2012).

The increasing number of paleoecological studies using multiple methods has revealed that C. antiquitatis from the northern Eurasia was a typical grazer with a highly abrasive diet, including dental microwear and mesowear analyses (Rivals et al., 2010; Rivals and Lister, 2016; Rivals and Álvarez-Lao, 2018; Pushkina et al., 2020), stable isotope analysis (e.g. Bocherens et al., 1991; Bocherens et al., 1997; Pushkina et al., 2014; Schwartz-Narbonne et al., 2019), and the examination of plant remains from gastro-intestinal content (Boeskorov et al., 2011a; Stefaniak et al., 2021) and teeth (van Geel et al., 2019; Stefaniak et al., 2021). Nonetheless, some studies have indicated a more mixed diet, signifying dietary plasticity and seasonality of C. antiquitatis (Rivals et al., 2010; Tiunov and Kirillova, 2010; Rivals and Lister, 2016; Rey-Iglesia et al., 2021).

The stable carbon and nitrogen isotope analyses of the Mammuthus-Coelodonta Faunal Complex in Heilongjiang suggested that C. antiquitatis may have been a grazer in East Asia, while also exhibiting dietary flexibility (as indicated by their intermediate δ¹³C and relatively high δ¹⁵N values, see Supplementary Tables S4, S11, S12; the isotopic pattern of the plants on the mammoth steppe was different from other biomes discussed here) (Ma et al., 2017; Ma et al., 2021). Within this fauna, delineation of ecological niche using carbon and nitrogen isotopic signatures roughly suggested a grazing ecology for mammoths, woolly rhinoceroses, and bovids, while equid specimens appear to have consumed more shrub or tree leaves, and cervid specimens likely ingested lichens, mosses, and leaves (Ma et al., 2017; Ma et al., 2021). This suggests that various taxa exhibited distinct ecological strategies when confronted with a completely distinct paleoenvironments during the Late Pleistocene.

The abundance and wide distribution of C. antiquitatis in Late Pleistocene East Asia have been well-documented (Tang et al., 2003; Li and Dong, 2007; Dong et al., 2014; Chen et al., 2016; Li et al., 2018). C. antiquitatis coexisted with woolly mammoths in northeastern China and the straight-tusk elephant Palaeoloxodon in warmer zones, with the southernmost boundary at 33° latitude (Tong, 2004; Tong et al., 2008; Dong et al., 2014; Wang et al., 2021). Coelodonta-Palaeoloxodon sympatry also occurred in England (Schreve, 2001, 2019). The faunal mixture and the wide distribution of C. antiquitatis may be linked to glacial/interglacial cycles, yet they also suggest the high ecological adaptability of C. antiquitatis, enabling this species to flourish in both cold and temperate environments.

This study conducted a comprehensive review of the paleoecological and evolutionary history of the woolly rhinoceroses from a spatiotemporal perspective. The ancestral C. thibetana and the terminating C. antiquitatis were more likely grazers that inhabited the high-altitude Tibetan Plateau during the mid-Pliocene and the middle to high latitude Eurasia during the Late Pleistocene. In contrast, C. nihowanensis, an intermediate form of the woolly rhinoceroses, exhibited flexible foraging ecologies across varied environments from west to east in northern China. The wide distribution of C. antiquitatis in temperate regions also indicates its high environmental tolerance. As an iconic cold-adapted mammal, Coelodonta adjusted its ecological strategies upon departing the Tibetan Plateau, ultimately dispersing from East Asia and colonizing middle to high latitude Eurasia. Further exploration of the paleoecology of the Middle Pleistocene C. tologoijensis will enhance our understanding of the complete ecological evolution of Coelodonta. While this paper reconstructed the paleoecology of the genus Coelodonta from an isotopic perspective, it is evident that further multidisciplinary investigations including systematic stratigraphy analysis, microwear and mesowear analyses, ecomorphological analysis, and locomotive behavior analysis are required to gain a more comprehensive understanding of the paleoecology of the ancestors of the woolly rhinoceroses.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

JM: Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Writing – original draft, Writing – review & editing. SW: Project administration, Resources, Writing – original draft, Writing – review & editing. TD: Conceptualization, Investigation, Resources, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by State Key Laboratory of Palaeobiology and Stratigraphy (Nanjing Institute of Geology and Palaeontology, CAS) (No. 223119), the National Natural Science Foundation in China (No. 42002005), Award of Youth Innovation Promotion Association (No. 2022070), and National Key Research and Development Program of China (No. 2023YFF0804501, 2022YFF0903800).

Acknowledgments

We appreciate BY Sun, QQ Shi, and Y Chen for identifying the specimens and engaging in useful discussions. We extend our gratitude to L. Zhang from the Hezheng Paleontological Museum for assisting in collecting fossil specimens, and to M.Q. Ma (the first author's father) for gathering modern samples. We thank H.W. Li and H.Y. Rao for their help in the lab. We thank two reviewers (Maciej Tomasz Krajcarz and Naoto Handa) and the Editor Florent Rivals for valuable suggestions, which have highly improved this manuscript. We would like to thank ChatGPT 3.5 for its valuable inputs and assistance for improving the language and checking grammar of this manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2024.1377000/full#supplementary-material

Footnotes

1. ^ The NOW Community. New and Old Worlds Database of Fossil Mammals (NOW). Licensed under CC BY 4.0. https://nowdatabase.org/now/database/.

References