Long-Term Isotope Evidence on the Diet and Habitat Breadth of Pleistocene to Holocene Caprines in Thailand: Implications for the Extirpation and Conservation of Himalayan Gorals

Three taxa within the subfamily Caprinae (Himalayan goral Naemorhedus goral, Chinese goral Naemorhedus griseus, and Sumatran serow Capricornis sumatraensis) live in the mountainous upland forests of Southeast Asia, where they are considered as vulnerable or near threatened species. Co-occurrences between these two recognized genera have been documented from some Pleistocene fossil sites in Thailand, suggesting more widely overlapping distribution in the past than today. However, diet and habitat preferences of these Pleistocene and present-day coexisting species have rarely been investigated so far. For the past three decades, stable carbon and oxygen isotope analyses become more commonplace in ecological investigations, allowing us to explore the diets and habitats of ancient and extant animals as well as to reconstruct environmental conditions in the past. We reconstructed diets and habitats of these taxa from five fossil sites in Thailand during the past 400,000 years (from the Middle Pleistocene to the Early Holocene) and from some modern wildlife using the isotopic analysis of carbonate in tooth enamel, in order to test species co-occurrence patterns during the Pleistocene and to examine possible changes of their niche breadths over evolutionary time. Our carbon isotope analysis revealed remarkably different ecological patterns between Naemorhedus and Capricornis. The Pleistocene Sumatran serow has been a greater generalist than both the Himalayan and Chinese gorals that fed on pure C4 or mixed C3 and C4 plants restricted to an open landscape habitat and than its extant population that occupies a closed-canopy forest. This suggests that the habitat contraction of the modern wildlife is likely due to the Holocene climate change and the human impacts on Thai ecosystems. In addition to the loss or reduction of grasslands after the latest Pleistocene when rainforests became dominant and besides the human hunting and predation pressure, the high interspecific competition likely contributed to the extirpation of Himalayan gorals in Thailand. Developing a strategic plan for the future biodiversity conservation, a long-term historical isotope approach allowed us to predict the contrasting habitat suitability, a lowland grassland, for these two threatened goral species as testified by their ecological persistence during the Pleistocene.

Three taxa within the subfamily Caprinae (Himalayan goral Naemorhedus goral, Chinese goral Naemorhedus griseus, and Sumatran serow Capricornis sumatraensis) live in the mountainous upland forests of Southeast Asia, where they are considered as vulnerable or near threatened species. Co-occurrences between these two recognized genera have been documented from some Pleistocene fossil sites in Thailand, suggesting more widely overlapping distribution in the past than today. However, diet and habitat preferences of these Pleistocene and present-day coexisting species have rarely been investigated so far. For the past three decades, stable carbon and oxygen isotope analyses become more commonplace in ecological investigations, allowing us to explore the diets and habitats of ancient and extant animals as well as to reconstruct environmental conditions in the past. We reconstructed diets and habitats of these taxa from five fossil sites in Thailand during the past 400,000 years (from the Middle Pleistocene to the Early Holocene) and from some modern wildlife using the isotopic analysis of carbonate in tooth enamel, in order to test species co-occurrence patterns during the Pleistocene and to examine possible changes of their niche breadths over evolutionary time. Our carbon isotope analysis revealed remarkably different ecological patterns between Naemorhedus and Capricornis. The Pleistocene Sumatran serow has been a greater generalist than both the Himalayan and Chinese gorals that fed on pure C 4 or mixed C 3 and C 4 plants restricted to an open landscape habitat and than its extant population that occupies a closed-canopy forest. This suggests that the habitat contraction of the modern wildlife is likely due to the Holocene climate change and the

INTRODUCTION
To allow the conservation of endangered mammal species, it is necessary to find a suitable habitat where they could live and prosper. Looking for habitats similar to the ones where they currently live seems a reasonable approach (e.g., Lekagul and McNeely, 1988;Raia et al., 2012;Rovero et al., 2014;Cavada et al., 2019), but many mammal species today are restricted to relic populations occasionally occupying suboptimal niches that do not reflect the habitat preference of the species, and they would possibly have better chances of long-term survival in a different habitat (e.g., Kuemmerle et al., 2012;Bocherens et al., 2015;Lea et al., 2016). Investigating the ecology of past population of mammal species can therefore provide important information about their preferred habitats and diets, and help to define a successful strategy for their long-term conservation (e.g., Bocherens et al., 2015;Jürgensen et al., 2017;Archer et al., 2019).
Southeast Asian caprines are fundamentally represented by four species within two genera, gorals (Naemorhedus goral and Naemorhedus griseus) and serows [Capricornis sumatraensis (possibly including three subspecies: Capricornis sumatraensis thar, Capricornis sumatraensis milneedwardsii, and Capricornis sumatraensis), and Capricornis rubidus] based on the recently genetic evidence (Dou et al., 2016;Bover et al., 2019;Mori et al., 2019). Their taxonomic status is, however, under debate today. All of these Southeast Asian caprine taxa are currently listed in Appendix I (threatened with extinction) of CITES and categorized as near threatened or vulnerable species on the IUCN Red List of Threatened Species (Duckworth and MacKinnon, 2008;Duckworth and Zaw, 2008;Duckworth et al., 2008a,b,c) ( Table 1). The Chinese goral N. griseus and Sumatran serow C. sumatraensis are additionally classified as the Class II-protected species in China (Grubb, 2005).
These caprines usually inhabit a wide range of habitats from high to low montane forests, along steep slopes and/or near cliffs (Lekagul and McNeely, 1988;Grubb, 2005;Groves and Grubb, 2011) (Table 1). Their current geographic ranges mostly overlap in distributions (except for N. goral), especially C. sumatraensis that coexists widely with N. griseus and possibly C. rubidus in some parts of Myanmar, Thailand, and China (Figure 1). Some previous studies have demonstrated a clear coexistence pattern between extant N. griseus and C. sumatraensis in Western Sichuan, China, where the large area of a nature reserve has an elevation ranging from 1,100 and 4,800 m above sea level and is protected (Chen et al., 2009(Chen et al., , 2012. This evidence suggests that these two sympatric caprine species reduced intergeneric competition by occupying different habitats, in terms of canopies, altitudes, and food resources, which may allow their coexistence today. Steep slopes adjacent to the high-altitude area with a high shrub density are preferred by C. sumatraensis, while gentle slopes in the lower altitude area with a lower shrub density are the habitat of N. griseus (Chen et al., 2009(Chen et al., , 2012. To better understand occupancy patterns and niche partitioning within these extant caprine communities, investigating their past ecological history that resulted in the basis for those processes is therefore required. In a similar way, an effective future wildlife conservation also needs to take into consideration a substantial difference between the extant and past ecological patterns of the species. Based on fossil records during the Early to Late Pleistocene, both Naemorhedus and Capricornis are considered to be members of the typical Southeast Asian mammal assemblages, also known as the "Ailuropoda-Stegodon fauna complex." Their fossil remains have been found in large areas of Asia including China (e.g., Colbert and Hooijer, 1953;Kahlke, 1961;Hu and Qi, 1978;Han and Xu, 1985;Rink et al., 2008;Zhu et al., 2015), Vietnam (Long et al., 1996;Bacon et al., 2008aBacon et al., , 2018a, Laos (Bacon et al., 2008b(Bacon et al., , 2011(Bacon et al., , 2012, Thailand (Tougard, 1998(Tougard, , 2001Zeitoun et al., 2010;Filoux et al., 2015;Suraprasit et al., 2016;Wattanapituksakul et al., 2018), Cambodia (Bacon et al., 2018c), Peninsular Malaysia (Ibrahim et al., 2013), and Indonesian islands (Java and Sumatra) (Hooijer, 1958) (Figure 1).
In Thailand, previous studies on several Pleistocene mammal assemblages have documented the occurrence of at least three caprine taxa: N. goral, N. griseus, and C. sumatraensis and have sometimes shown coexistences between these two genera [see Suraprasit et al. (2016): appendices 15-17 for further detailed fauna lists of all caprine-bearing sites in Southeast Asia and China]. For instance, the late Middle Pleistocene of Tham Wiman Nakin (Tougard, 1998(Tougard, , 2001 has yielded both N. goral and C. sumatraensis, whose distribution ranges do not overlap today in the area (Figure 1). The co-occurrence of these three species has also been documented from the late Pleistocene of Tham Lod Rockshelter (Wattanapituksakul et al., 2018) where only two caprine taxa coexist today in the area (Figure 1). This implies that their past biogeographic distributions were latitudinally broader and probably overlapped more widely than today. Moreover, the fossil records have also suggested that the extant Southeast Asian caprines are represented by relict populations, which were more widespread during the Pleistocene but became restricted  Grubb, 2005;Duckworth and MacKinnon, 2008;Fakhar-I-Abbas et al., 2008Thapa et al., 2011Lovari, 1997Buranapim and Sitasuwa, 2014;Grubb, 2005;Duckworth et al., 2008a;Chen et al., 2009Chen et al., , 2012Grubb, 2005Duckworth et al., 2008b,c;Thuc et al., 2012 Duckworth andZaw, 2008;Groves and Grubb, 2011;Castelló, 2016 to the highlands and the high-or even low-altitude mountain ranges during the Holocene interglacial. It is then presumable that these caprines were abundant in tropical lowlands and that some habitat partitioning between these closely related species might have occurred in the past coexistence areas. This fossil evidence emphasizes the need to address the evolutionary issues of niche differentiation and competition between these three caprine taxa through space and time, and to understand better the past ecological factors driving the distribution patterns of their modern populations. Moreover, the understanding of the diet and habitat use of these sympatric species is crucial, not only for improving our knowledge on the refugial hypothesis but also to achieve more effective and efficient conservation plan of the species for the future, as it is currently implemented for other endangered species (e.g., Birks, 2012;Louys, 2012;Archer et al., 2019). Among several ecological approaches, stable isotope tracking on tooth enamel is recognized as a powerful tool for exploring several aspects of the life history of living and fossil animals. This approach has successfully been used to document shifts or restrictions of ecological niches for mammalian species that survived the final Quaternary extinction wave, such as a Madagascar lemur (Crowley et al., 2012), European bison (Bocherens et al., 2015;Hofman-Kamińska et al., 2019), and saiga antelope (Jürgensen et al., 2017). An analysis of stable carbon isotopes ( 13 C/ 12 C) is an efficient method for investigating the diet and habitat use of animals, whether they foraged mostly C 3 plants in a closed canopy, C 4 plants in an open-habitat, or mixed C 3 and C 4 plants in an intermediate area (e.g., DeNiro and Epstein, 1978;MacFadden, 1998;Cerling and Harris, 1999;MacFadden et al., 1999;Cerling et al., 2004;MacFadden and Higgins, 2004;Feranec and MacFadden, 2006;Stacklyn et al., 2017). Measurements of stable oxygen isotope composition ( 18 O/ 16 O) in tooth enamel are also useful to explore the drinking or plant consumption behavior of animals due to isotopic differences in local meteoric water, depending on climatic conditions such as air temperature, humidity, and precipitation amount (e.g., Dansgaard, 1964;Bryant et al., 1996;Fricke and O'Neil, 1996;Kohn, 1996;Higgins and MacFadden, 2004). In some Pleistocene Southeast Asian sites, previous carbon isotope measurements on tooth enamel have revealed variabilities in dietary signatures, either pure C 3 or C 4 plants or even both types of vegetation, for the species of C. sumatraensis [e.g., data from Tham Wiman Nakin (Pushkina et al., 2010), Khok Sung (Suraprasit et al., 2018), Boh Dambang (Bacon et al., 2018c), and Nam Lot (Bacon et al., 2018b)], leading to the question of whether it was fundamentally a grazer or a browser. Moreover, diet and habitat preferences of both species of Naemorhedus (N. griseus and N. goral) have rarely been investigated so far.

Study Areas and Sample Collection
We designed a hierarchical isotope sampling of caprine tooth enamel from five different-aged fossil sites (Table 2), namely, Pha Bong (400 to 200 ka), Khok Sung (217 or 130 ka), Tham Wiman Nakin (>169 ka), Tham Lod Rockshelter (32 to 12 ka), and Ban Rai Rockshelter (10 to 6 ka) (Figure 1), as well as from the modern representatives (see Supplementary Material S1 for more detailed information on geological backgrounds of aforementioned sites and see Supplementary Material S2 for the fauna list). A total of 94 bulk tooth enamel samples were extracted mostly from the second or third molars (fully formed adults) of caprine individuals (Supplementary Material S3). The bulk enamel powder was sampled along the entire length of a tooth in order to obtain average isotopic compositions over the time period of dental growth. The analyzed fossils comprise eight specimens of N. goral from Pha Bong, three specimens of C. sumatraensis from Khok Sung, four specimens of N. goral and nine specimens of C. sumatraensis from Tham Wiman Nakin, eight specimens of N. goral, 19 specimens of N. griseus, and 12 specimens of C. sumatraensis from Tham Lod Rockshelter, and six specimens of C. sumatraensis from Ban Rai Rochshelter. Modern specimens including eight specimens of wild N. goral from the Indian Himalayan region, one specimen of wild N. griseus from Northern Thailand, and 16 specimens of wild C. sumatraensis from Thailand and Sumatra were collected from several natural history museums and institutes including Chulalongkorn University Museum of Natural History (CUMNH, Bangkok), Natural History Museum (THNHM, Pathum Thani), Naturhistorisches Museum Basel (NMB, Basel), and Zoologische Staatssammlung München (ZSM, Munich). In this study, we additionally performed the stable isotope analysis of tooth enamel of other mammals from Tham Lod (12 samples) and Ban Rai (10 samples) rockshelters, used as baseline values for isotopic comparisons, in order to understand the environmental contexts of these sites (Supplementary Material S4).

Taxonomic Identification
To conduct isotope analysis of tooth enamel, we selected almost complete fossil tooth specimens in order to have access to morphological features that are useful for the taxonomic identification. All analyzed fossil teeth were identified based on their morphological and dimensional comparisons with modern samples. The teeth of N. goral differ from N. griseus in having a less rounded outline on P4, more distinct paracone on upper molars, more distally curved crown on M3 in labial view (straight in N. griseus), U-shaped anterior valley on p3, more prominent labial valley on p4, and more developed conids and entoconulids on lower molars (Wattanapituksakul et al., 2018). The teeth of N. goral are similar in size to those of N. griseus but are clearly smaller than those of C. sumatraensis (Suraprasit et al., 2016;Wattanapituksakul et al., 2018). In addition, Capricornis is distinguished from Naemorhedus by several dental characters: more developed cusps and labial styles on upper premolars and molars (especially the mesostyle and metastylar wing on M3) and more distinct lingual stylids on lower premolars and molars (Wattanapituksakul et al., 2018). In the case of modern samples analyzed here, we followed the taxonomic identification labeled by the field collectors and/or curators.

Isotopic Preparation and Analysis
All caprine enamel samples were pretreated and analyzed at the Department of Geosciences, University of Tübingen, Germany (see Supplementary Material S5 for more detailed information on pretreatment and isotope protocols). Stable carbon and oxygen isotopic values are expressed as the following standard delta (δ)-notation: X = [(Rsample/Rstandard)-1] * 1000, where X is referred to δ 13 C and δ 18 O values and R is equivalent to 13 C/ 12 C or 18 O/ 16 O, respectively. The recorded delta values follow the international reference standards, "Vienna PeeDee Belemnite" (VPDB) for the carbon and oxygen. Additionally, δ 18 O values relative to "Vienna Standard Mean Ocean Water" (VSMOW) are given. The carbonate content (CaCO 3 %) was calculated using the ratio between amount of CO 2 released by the reaction, as detected from the peak intensity for mass 44 and the weight of pure carbonate used as a standard, with an analytical error of 0.3%, based on the multiple analysis of reference enamel samples. Before comparing the stable isotope ratios between fossil and modern faunas, we took into account the changes in the carbon composition due to anthropogenic CO 2 emission ( 13 C Suess effect) since the beginning of the industrial revolution in 1850 and ongoing today (Keeling et al., 1979). The atmospheric δ 13 C values in relation to increased CO 2 emissions have decreased approximately from −6.5 to −8.0 over the past few centuries due to fossil fuel combustion and biomass destruction (e.g., Friedli et al., 1986;Marino and McElroy, 1991). To convert δ 13 C values to modern plant isotope compositions due to the Suess effect, the δ 13 C values of modern carbonate bioapatite are therefore adjusted upward by 1.5 for comparing with the Pleistocene enamel carbon isotope compositions (Cerling et al., 1997a;Cerling and Harris, 1999;Passey et al., 2005). However, modern samples analyzed herein do not require the Suess effect δ 13 C correction because these wild mammals died or were collected before or around 1950 (Supplementary Material S3).

Statistical Analysis
We performed statistical analyses of the isotopic data (δ 13 C and δ 18 O) to show a significant difference among our analyzed specimens. Based on the Shapiro-Wilk test, our samples do not correspond to a normal distribution with unequal variances. In the case for which a sample of at least five individuals is available, we used non-parametric tests to examine differences in median δ 13 C and δ 18 O values among our analyzed caprine samples (Kruskal-Wallis test) and between the species and the locality (Mann-Whitney U-test). Significant differences are statistically attained when p-values are equal to or less than 0.05.

RESULTS AND INTERPRETATIONS
To interpret the δ 13 C and δ 18 O results of the caprine bioapatite carbonates, we followed the fundamental principles of stable isotope applications on the ecological study of mammals (e.g., Biasatti et al., 2012, pp. 53, 57;Suraprasit et al., 2018, pp. 26, 27). The carbonate content (CaCO 3 ) of caprine tooth enamel varies from 2.4% to 6.9% (Supplementary Material S3), corresponding to the normal range of modern samples that have not been altered during the diagenetic processes (e.g., Ecker et al., 2013). The results thus do not show any evidence of remarkable diagenetic contamination. The enamel carbonate δ 13 C values of all analyzed caprines exhibit a median of −3.0 and a range between +3.0 and −16.6 , whereas the δ 18 O VPDB results make on average −5.2 for a median and range from +0.1 to −11.5 (Figure 2 and

Naemorhedus goral
The δ 13 C values of Himalayan gorals from Pha Bong (including material previously analyzed by Bocherens et al., 2017) displayed a median of −0.5 and ranged from −2.6 to +3.0 (Figures 2, 3A), indicating the consumption of mostly pure C 4 plants with a limited signal of mixed C 3 and C 4 diets. Compared to the previously isotopic-analyzed samples of Pha Bong Himalayan gorals by Bocherens et al. (2017), our results showed similar tooth enamel δ 13 C values (Mann-Whitney U: NS; Supplementary Material S6), supporting the occurrence of a single species in this locality. The Himalayan gorals from Tham Wiman Nakin had a similar dietary signature (a median of −1.2 and a range between −3.5 and +0.4 ) to those of Pha Bong (Figures 2, 3A) but exhibited a shift of the δ 13 C values toward the C 3 /C 4 mixed signals. A trend of the δ 13 C shifts into the direction of C 3 plants is also observed from the Tham Lod Rockshelter Himalayan gorals, which yielded a median of −3.3 and a range between −5.3 and −0.1 (Figures 2,  3B), indicating their diet selection chiefly on both C 3 and C 4 plants and only sometimes on pure C 4 plants. In contrast, the δ 13 C values of modern Himalayan gorals displaying a median of −12.2 and varying from −15.2 to −6.0 (Figures 2,  3C) fall into the range of pure C 3 and mixed C 3 /C 4 ecosystems, reflecting an adaptation for occupying either more closed habitats or C 3 -dominated grasslands at the Indian Himalayan Region where the conditions at such an extremely high elevation are unfavorable for C 4 plants (Dobremez, 1976;Blasco et al., 1996;Galy et al., 2008;Bhattacharyya et al., 2019).
The Himalayan gorals from Pha Bong showed more positive δ 18 O values, compared to those from other localities (Mann-Whitney U: P < 0.05; Supplementary Material S6), ranging FIGURE 2 | Box plots of enamel δ 13 C and δ 18 O values of three caprine species, Naemorhedus goral (green), Naemorhedus griseus (red), and Capricornis sumatraensis (blue) from all analyzed fossil sites and modern representatives. The abbreviation of "n" indicates the number of analyzed samples. and a range from −6.5 to −1.7 ) and Tham Lod Rockshelter (a median of −6.4 and a range from −11.4 to −1.8 ), as well as from the modern representatives (a median of −6.2 and a range from −10.6 to −3.3 ) (Figures 2, 3), likely a consequence of the consumption of water depleted in 18 O compared to those from Pha Bong. The most variable δ 18 O values recorded from the Tham Lod Rockshelter Himalayan gorals might be explained by more opportunistic feeding on both 18 O enriched and depleted plant water. In general, the wide range of tooth enamel δ 18 O values observed in the Pleistocene and extant Himalayan gorals might reflect the dietary habit of each individual in response to seasonal differences because the modern population has been suggested to feed in the variable frequency on plant species between the summer and winter (Ashraf et al., 2015(Ashraf et al., , 2016. Otherwise, the wide δ 18 O ranges in the tooth enamel of both Pleistocene and extant Himalayan gorals were possibly due to their seasonal migration across different altitudes in response to temperature changes and threats from predators and hunters (Perveen and Khan, 2013). Among the Himalayan gorals, some higher δ 18 O values corresponding to more dominated C 4 signals are observed, probably reflecting a consumption of 18

Naemorhedus griseus
Pleistocene Chinese gorals documented only in Tham Lod Rockshelter yielded a median δ 13 C values of −0.4 with a range between −7.1 and +1.9 (Figures 2, 3B), wider than those of Himalayan gorals from Pha Bong and Tham Wiman Nakin as well as in the same locality (although non-statistical significance tested by Mann-Whitney U; Supplementary Material S6) (Figure 2). This implies that the Chinese gorals had the ability to forage more exclusively on some C 3 plants and/or to occupy more closed habitats than the Himalayan gorals. However, this was possibly due to the sampling bias for which more or less specimens were analyzed among these localities. One modern sample of the wild Chinese goral individual collected from the region of Northern Thailand, displays a δ 13 C value of −3.3 (Figure 3C), suggesting the consumption of mixed C 3 and C 4 plants, similar to some Pleistocene populations in Tham Lod Rockshelter (Figure 2).
The Tham Rod Rochshelter Chinese gorals had a median δ 18 O value of −5.6 , similar to the modern sample (−5.7 ), and exhibited broad δ 18 O values varying from −8.9 to −0.3 (Figures 2, 3B), suggesting that they derived their water from highly variable food or local resources, similar to what is observed in the congeneric species (N. goral) in the same locality. Otherwise, it is conceivable that the Chinese gorals might have ingested water from local resources in different altitudes in relation to their seasonal habitat use or migration, as seen in the modern population (Chen et al., 2009(Chen et al., , 2012.

Capricornis sumatraensis
The δ 13 C values of Khok Sung Sumatran serows including three samples in this study and one previously analyzed material from Suraprasit et al. (2018, Table 1) displayed a median of −13.9 and a range between −15.1 and −11.8 (Figures 2, 3A), indicating diet and habitat preferences in a closed C 3 -dominated environment. The Sumatran serows recovered in Tham Wiman Nakin and Tham Lod Rockshelter showed a median of −11.0 and −12.2 (Mann-Whitney U: NS; Supplementary Material S6) and exhibited considerably wide δ 13 C values ranging from −12.8 to +1.7 and from −14.3 to +1.9 , respectively (Figure 2). The Sumatran serows from both localities inhabited the various types of tropical ecosystems from exclusively closed C 3 to open C 4 vegetation landscapes. The carbon isotope compositions of tooth enamel of these Pleistocene Sumatran serows clearly support a broader foraging niche than those for Naemorhedus (Mann-Whitney U: P < 0.05; Supplementary Material S6). However, the Early Holocene Ban Rai and modern populations showed more limited diets of either pure C 3 or C 3 /C 4 mixed vegetation ( Figure 3C) and more restricted habitat use in either closed environments or areas adjacent to forest/woodland landscapes based on the median δ 13 C values of −14.1 with a narrow range between −16.6 and 12.3 , for the former site and on the median δ 13 C values of −13.3 with a narrow range between −17.8 and −8.8 for the latter (Mann-Whitney U:

NS; Supplementary Material S6).
Highly variable δ 18 O values are observed among most of the Pleistocene and modern Sumatran serows including a median of −5.7 and a range between −8.9 and −3.6 for Tham Wiman Nakin, a median of −6.3 and a range between −8.1 and +0.1 for Tham Lod Rockshelter, and a median of −7.0 and a range between −10.5 and −1.9 for the modern population (Mann-Whitney U: NS; Supplementary Material S6) (Figure 2). Most of the Pleistocene and extant Sumatran serows with the exception of Khok Sung samples thus indicated their independent habits in consuming various leaf water or ingesting meteoric water across variable local resources. The narrower range of the δ 18 O values, with a median of −5.0 and a range between −6.4 and −3.7 (Figure 2), for Khok Sung serows might be explained by fewer isotopic-analyzed samples from this locality. Compared to other localities, C. sumatraensis from Ban Rai Rockshelter likely displayed more limited and negative δ 18 O values with a median of −9.0 and a range between −11.5 and −8.3 (Figure 2) (Mann-Whitney U: P < 0.05; Supplementary Material S6), corresponding to its occupation of a closed forest canopy with a more humid condition. Some C 4 grazers tended to have higher δ 18 O values compared to C 3 browsers, consistent with their dietary intake and habitat use (Figure 3 and Supplementary Material S7).

Evolutionary Niche Differentiation and Adaptation of Naemorhedus and Capricornis Through Time
Analyzing stable carbon isotope composition of fossil caprines from different time periods shows that none of the analyzed Pleistocene N. goral revealed a significant shift in their tooth enamel δ 13 C values over time prior to the Holocene (Figure 4). N. goral therefore had an ecological persistence on the diet of FIGURE 4 | Bulk δ 13 C values of mammals from several Pleistocene to Holocene caprine-bearing sites in Southeast Asia. In the tropics, the cut-off δ 13 C values of −10 between pure C 3 and C 3 /C 4 mixed diets and of +2 for pure C 4 diet (dashed lines) are applied to herbivorous ungulates according to Cerling et al. (1997a), Cerling and Harris (1999), and Biasatti et al. (2012). The cut-off δ 13 C value of lower than −13 is adopted here for the mammals that lived in closed forests according to Cerling et al. (1997aCerling et al. ( , 2004, Kohn et al. (2005), and Biasatti et al. (2012). For the isotope comparison between large herbivores and predators, the fractionation factor of δ 13 C values for carnivores is adjusted by +1.3 (Clementz et al., 2009;Domingo et al., 2013). The faunal isotope baseline includes previously analyzed enamel samples from Pha Bong , Khok Sung (Suraprasit et al., 2018), and Tham Wiman Nakin (Pushkina et al., 2010;Bocherens et al., 2017) in Thailand, Nam Lot (Bacon et al., 2018b) in Laos, and Boh Dambang in Cambodia (Bacon et al., 2018c) (see Supplementary Material S4 for raw data). exclusively C 4 plants but sometimes mixed C 3 /C 4 vegetation and on the habitat in open environments through the Middle to latest Pleistocene. In contrast, stable carbon isotope compositions of modern samples of N. goral showed that their δ 13 C values became lower, indicating a diet with significantly more C 3 food items in more closed habitats or C 3 -dominated grasslands (Figure 4). This isotopic shift corresponds to the extant Indian Himalayan population that is surviving in the extreme condition of highland or mountain climate refugia where the C 3 -grasslands are dominant (Dobremez, 1976;Blasco et al., 1996;Galy et al., 2008;Bhattacharyya et al., 2019) and where the anthropic effects are presumably weaker.
The isotopic analysis of tooth enamel of N. griseus from the Late Pleistocene of Tham Lod Rockshelter indicated its occupation of resources in the same habitats of N. goral (Mann-Whitney U: NS; Supplementary Material S6). Although only one modern isotope sample of N. griseus has been analyzed for this study due to the inaccessibility and rareness of wild individual specimens, the Chinese goral usually inhabits steep landscapes and plateaus in mountainous areas or sometimes subtropical mixed and evergreen-deciduous forests near cliffs based on the field observation (Duckworth et al., 2008a;Chen et al., 2009Chen et al., , 2012Groves and Grubb, 2011;Buranapim and Sitasuwa, 2014). This information recorded in the field indicates an ecological adaptation struggling to live in a more closed habitat for the extant population of N. griseus, despite having mixed C 3 /C 4 plants in diet (as seen in the single analyzed sample from this study), different from what has initially been observed in some Pleistocene inhabitants that mainly occupied an open habitat.
Our stable isotope results also indicate that C. sumatraensis has been well-adapted to variable environmental conditions expanding from closed C 3 to open C 4 landscapes since the Middle Pleistocene and has been a generalist with a greater ecological flexibility in diet and habitat than both Pleistocene and extant Naemorhedus (Figures 2, 4). Although some isotopic data from well-dated localities such as Khok Sung (MIS7 or MIS5) (Suraprasit et al., 2018;Duval et al., 2019;this study) and Nam Lot (MIS5) (Bacon et al., 2018b) as well as from the Early Holocene of Ban Rai Rockshelter have documented the preferred habitat restriction of the Sumatran serows in closed forest environments (Figure 2), this might be explained by the sampling bias for which few specimens are available or analyzed or by the ecological adaptation of C. sumatraensis in response to more homogeneous and closed conditions during interglacial stages. In the heterogeneous environments during the Pleistocene, C. sumatraensis preferably occupied in both closed forest and open grassland landscapes but rarely in an intermediate area between these two canopies. The rareness of C. sumatraensis as being a mixed C 3 /C 4 feeder in the intermediate habitat area was possibly a consequence of its competition avoidance for food resources with other local wildlife such as megaherbivores (e.g., Stegodon) and omnivores (e.g., Sus) (Figure 4). Previous isotope studies have revealed that Stegodon was mainly a C 3 /C 4 -mixed feeder and significantly held an intermediate canopy landscape in some coexisting area (e.g., Khok Sung, Suraprasit et al., 2018) (Figure 4). Moreover, the preferred habitat of suids in the intermediate habitat canopy between closed and open ecosystems remained relatively stable over time (e.g., Pha Bong, Khok Sung, Tham Wiman Nakin, Nam Lot, and Boh Dambang) (Pushkina et al., 2010;Bocherens et al., 2017;Bacon et al., 2018b,c;Suraprasit et al., 2018) (Figure 4 and Supplementary Material S4). The limited δ 13 C range of numerous modern enamel samples suggests that the ecological flexibility of C. sumatraensis has become restricted after the latest Pleistocene (Figures 3C, 4). This rests on the assumption that the present-day wildlife habitat in Thailand is likely constrained by the Holocene climate change and the human impact on the ecosystem restricted its population to more closed canopy habitats such as dense forests in higher altitude areas.

Tooth Morphology and Body Size Influences on the Ecological Adaptation of Caprines
As Naemorhedus and Capricornis did coexist in several Thai fossil sites during the Pleistocene, niche partitioning between them might have occurred on other aspects that are not discernible through the stable carbon isotope analysis.
Regarding the dental morphology, these three caprine taxa possess a distinctly selenodont appearance with similar occlusal patterns. In the evolutionary history of mammals, the development of hypsodont teeth allows them to feed on more abrasive food such as grasses and to therefore expand the range of dietary resources (e.g., Mendoza and Palmqvist, 2008;von Koenigswald, 2011). Both Naemorhedus species have the same degree of hypsodonty as an index of 4.03 (Mendoza and Palmqvist, 2008), while C. sumatraensis has comparatively lower (but still high) crowned teeth according to the hypsodonty index of 3.39 (Kaiser et al., 2013). The dental morphologies of these three caprine species are congruent with our stable carbon isotope results, which indicate their life history adaptations to grazing during the Pleistocene.
Additionally, differences in the average body size between Naemorhedus and C. sumatraensis, as exemplified from the tooth dimensions on fossils (e.g., Suraprasit et al., 2016;Wattanapituksakul et al., 2018) and from the observation of recent populations (e.g., Nowak, 1999;Grubb, 2005) may permit them to feed at different elevations. The body size of N. goral and N. griseus is identical and ranges from approximately 22-35 kg for an adult individual (Nowak, 1999;Grubb, 2005). The Sumatran serows, in general, having a larger body size (about 50-140 kg) (Nowak, 1999), are more capable of foraging the higher level of plant bodies, such as leaves of shrubs or small trees, than the Himalayan and Chinese gorals. In an open habitat landscape around Tham Lod Rockshelter, C. sumatraensis might indeed have had more intense competition with other similar-sized grazers such as an Eld's deer Rucervus eldii (about 95-150 kg; Nowak, 1999) than these two smaller gorals, as suggested by our carbon isotope data (Figure 4 and Supplementary Material S4).
Although the tooth morphology supports a grazing adaptation for both genera, differences in body size between Naemorhedus and C. sumatraensis might have led to niche partitioning and minimized intergeneric competition, which have likely contributed to their prevalent coexistence during the Late Middle Pleistocene (i.e., Tham Wiman Nakin) and during the Last Glacial Maximum (LGM) (i.e., Tham Lod Rockshelter).

Causes of Local Extinction of Himalayan Gorals in Thailand
Today Chinese gorals and Sumatran serows are surviving in the wild despite having been considered as being vulnerable species due to a population decline (Duckworth et al., 2008a,b), while Himalayan gorals went locally extinct from the Thai territory likely at the end of the Pleistocene based on the lack of their fossil remains in the Holocene stratigraphic sequence records, as deduced from the rockshelters of Tham Lod and Ban Rai in Mae Hong Son, Northwestern Thailand. We thus proposed other possible extinction-related causes that are discernible through these two archeological sites.
It is likely that the environments around the Mae Hong Son Province where the three caprine taxa did coexist became more homogeneous and closed after the end of the Pleistocene based on the possible δ 13 C shift of C. sumatraensis and Bubalus arnee into more exclusively C 3 -dietary signals (Figure 4), and the high diversity and abundance of primates (Supplementary Material S2) in the Early Holocene of Ban Rai Rockshelter. Although an open landscape might have existed around the region of Ban Rai Rockshelter based on the C 4 signals of Sambar deer Rusa unicolor (Figure 4), large bovid remains are very scarce in Ban Rai Rockshelter and cervid subfossils are much less abundant than those found in Tham Lod Rockshelter. Supporting this argument, the sequential δ 18 O records of fresh water bivalves from both localities have indicated that substantially increased precipitation and decreased climate variability occurred in Northwestern Thailand after 9,800 BP (Marwick and Gagan, 2011). The climatic and environmental changes related to an increasing development of rainforests starting at the Late Pleistocene/Holocene transition likely contributed to the local extinction of grassland-dwelling taxa as Naemorhedus but other ecological factors also need to be considered.
The major role of anthropic activities was characterized by the accumulation of teeth and bones in Tham Lod Rockshelter where the animal remains have been modified (e.g., traces of fire, stone artifacts, and cut marks) by human hunter-gatherers at that time. Carnivores such as Asiatic black bear Ursus thibetanus and tiger Panthera tigris were major predators in the Tham Lod Rockshelter area but likely did not play a significant role in modifying or accumulating the skeletal elements in the assemblage. According to the observed fire modification of bones, it is assumed that prehistoric humans also hunted these carnivores for food and resources (Wattanapituksakul et al., 2018). Wattanapituksakul et al. (2018) demonstrated that prime adult individuals of Naemorhedus and Capricornis were preferentially selected by the hunter-gatherers based on the mortality profile analysis. The highest number of caprine individuals in the accumulation of the Tham Lod Rochshelter assemblages is observed in N. griseus (represented by the minimum number of individuals: MNI = 40), followed by C. sumatraensis (MNI = 15) and N. goral (MNI = 13). This evidence may provide general information on the abundance of caprine populations in this area (N. griseus > C. sumatraensis > N. goral), reflecting that N. goral was likely more scarce than N. griseus in the area and probably started decresing in its population prior the latest Pleistocene. Moreover, the mortality profile analyzed in Tham Lod caprines implies a hunter-gatherer subsistence strategy with selective hunting of easy preys such as these gorals and serows (Wattanapituksakul et al., 2018). As exemplified by the Tham Lod Rockshelter assemblage, we suggest that the impacts of the human hunting on the Himalayan and Chinese goral community were possibly one of the major factors leading to their population decline in that region during the latest Pleistocene.
Concerning the effects of predation pressure on these gorals, the diversity of carnivores in Southeast Asia appears to have been higher in the closed canopy than the open area throughout the Pleistocene based on numerous fossil records and previous carbon isotope data so far (e.g., Pushkina et al., 2010;Bocherens et al., 2017;Bacon et al., 2018b,c) (Figure 4). Our isotope results analyzed on the Tham Lod carnivores indicate their habitat preferences in closed forest ecosystems far away from the ones where the Himalayan and Chinese gorals lived at that time (Figure 4). Accordingly, the predation probably had a direct effect on which both of the goral species were forced to avoid closed habitats occupied by more diverse predators. Other predators that occupied an open habitat have been poorly known from the Late Pleistocene of Southeast Asia since spotted hyaenas possibly went extinct prior to MIS2 (see discussion in Suraprasit et al., 2015Suraprasit et al., , 2019. These gorals consequently maintained their preferred habitats in the different area as an open landscape where predation pressure was lower or minimal. Regarding the coexistence patterns, other sympatric herbivore species such as large bovids and cervids mainly occupied the open habitat around Tham Lod Rockshelter based on our isotope results (Figure 4). Accordingly, they possibly competed for available food resources with these three caprine taxa. However, this study has demonstrated that the interspecific exploitative competition between sympatric N. goral and N. griseus in Tham Lod Rockshelter was more intense than the intergeneric overlap for Capricornis and Naemorhedus due to the inflexibility of shared food resources and habitats in an open landscape for both goral species. Indeed, N. griseus might have been better adapted to a broader realized niche than N. goral did according to our stable carbon isotope analysis (Figures 2, 3B). The considerable competition for resources between these two closely related species that have a similar dental morphology and body mass and that occupied the same habitat might have led to their reduced population over time and to an extirpation of one species (i.e., N. goral) at the end of the Pleistocene.
Overall, we propose that these four factors: (1) the loss or reduction of grassland habitats after the latest Pleistocene when rainforests became dominant, (2) the negative effects of human hunting, (3) the predation pressure, and (4) the high interspecific competition between both goral species were possibly main causes driving the patterns of local extinction of N. goral in Thailand.

The Ecological Proposal for the Species Conservation of Extant Himalayan and Chinese Gorals
The pace of decline population of extant threatened Himalayan and vulnerable Chinese gorals are accelerating and their distribution ranges are constantly decreasing (Mead, 1989;Duckworth and MacKinnon, 2008) due to the habitat loss caused by several disturbance factors of natural processes (such as forest fire) and human activities (such as deforestation, hunting, and agricultural and commercial use) (Lydekker, 1907;Anwar, 1989;Duckworth and MacKinnon, 2008;Duckworth et al., 2008a;Ashraf et al., 2015;Shakeel et al., 2015). Some studies on the feeding biology of modern Himalayan and Chinese gorals directly targeted at the field have indicated that these two species are both grazers and browsers based on the food availability in their stomach, gizzard, and rumen content (e.g., Anwar and Chapman, 2000;Awasthi et al., 2003;Junaid et al., 2012;Buranapim and Sitasuwa, 2014). However, the most preferred food of both Himalayan and Chinese gorals are grasses, followed by shrubs, trees, twigs, nuts, and herbs (Duckworth and MacKinnon, 2008;Duckworth et al., 2008a;Fakhar-I-Abbas et al., 2008Buranapim and Sitasuwa, 2014). Our stable isotope results of tooth enamel of the Himalayan gorals reflected that this refugee species usually consumes the C 3 vegetation or grasses nowadays on the Himalayan regions where no or few C 4 grasses are able to develop under such an extremely high altitude condition due to an increase of atmospheric CO 2 and humidity levels after the LGM (Galy et al., 2008). As is the case for N. griseus, the current geographic distribution clearly shows the habitat restriction of this species in the high-elevation mountain areas (e.g., the regions of Northwest Thailand and South China). Although the natural or anthropogenic grasslands may exist somewhere in these regions, most of which are covered by dense forests such as deciduous dipterocarp and coniferous forests (Duckworth et al., 2008a;Groves and Grubb, 2011; (Table 1). In subtropical regions where Chinese gorals are living and where vegetation types vary along altitudinal gradients, alpine shrubs and grasslands generally occur at the highest altitudinal zonation of mountainous regions (e.g., Sichuan in China), followed by coniferous and mixed evergreen and deciduous forests in the lower-elevation canopy landscapes (Chen et al., 2009(Chen et al., , 2012 (Figure 5). In contrast, based on paleobotanical evidence, the plant zonation in Southeast Asia during Pleistocene glacials has been suggested to drop about 1,000 m in elevation lower than today, implying the expansion of grasslands in low-altitude mountain ranges and lowland floodplains (e.g., Kershaw et al., 2001;Morley, 2012) (Figure 5). A long-term ecological stasis in the occupation of both N. goral and N. griseus based on our isotope analysis of tooth enamel through several periods of glacial-interglacial cycles therefore provides a dietary consistency on (C 3 or C 4 ) grass components but a contrasting traditional view of their FIGURE 5 | A comparison between Pleistocene and Early Holocene to present-day caprine altitude ranges in relation to major ecosystem distributions in subtropical and tropical regions. For this graphic (as an example), the boundaries of altitudinal vegetation zones are briefly given according to this study for the Pleistocene in Thailand and to Chen et al. (2009Chen et al. ( , 2012  habitat use during the Pleistocene when the anthropic impacts on the ecosystems were lower than today. This study has demonstrated that the populations of both N. goral and N. griseus were C 4 grazers in lower altitude grasslands on mountainous areas or lowland floodplain during the Pleistocene glacials (e.g., Pha Bong, Tham Wiman Nakin, and Tham Lod Rockshelter) ( Figure 5). It is thus presumable that, during the Holocene, the warming climate and intensely anthropic effects have forced these two goral species to retreat upslope in pursuit of their comfort zone with the food availability (grasses). Particularly, the critical anthropic effects still happen on a regular basis upon the modern wildlife goral community with or without intention and may drive these species to the population decline or extinction risk in the future. Besides minimizing the ecological effects caused by human activities, we propose an alternative conservation-concern plan about the sustainable biodiversity management by supplying an optimal open lowland area with the availability of grass resources to the rest of Himalayan and Chinese goral populations.

CONCLUSION
Applying a stable isotope approach on tooth enamel of extirpated and extant caprines provides empirical evidence on the diet and habitat widths between Naemorhedus and Capricornis and allows us to better understand their adaptation to modern environments. It is clear that Sumatran serows have been generalist browsers and grazers, unlike Himalayan and Chinese gorals that had a diet specialist (not a specialized C 3 -plant feeder but C 4 grazer prior to the Holocene). In contrast with the traditional view, our study has revealed that both Himalayan and Chinese gorals likely had a preferred open habitat, such as low-elevation grasslands on mountainous areas or lowland floodplain, with a subsistence strategy in more variable climates over long time duration through the Pleistocene. In addition to human activities, the ecological status of these three taxa from what we have observed today is likely a consequence of the habitat contraction in response to the climatic change during the Holocene interglacial. Environmental changes as turning into more closed forest canopies during the Holocene, human impacts on the ecosystems over a long time, effects of predation pressure, and high interspecific competition between these two same-sized gorals advocate the demise of Himalayan gorals in Thailand and likely force the rest of their populations into the refugial stages. As grasses have been major dietary components of both N. goral and N. griseus throughout long-term historical time, this study finally provides an exemplary case of merging paleontological application with conservation biology to guide the future habitat restoration for these two goral species.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because we do not damage the specimens of vertebrates and do not work on the living animals. All the specimens analyzed for this study have been allowed by the curators and/or owners.

AUTHOR CONTRIBUTIONS
KS, J-JJ, RS, and HB conceived the research project. J-JJ, RS, and YC conducted the field excavations and provided the fossil samples. KS collected the modern samples and was a major investigator for the project. KS and AW identified the fossil material. KS and HB performed and directed the sampling, preparation, and isotopic analysis of the fossil and modern material, and wrote the manuscript. All authors discussed the results and commented on the manuscript.

FUNDING
For this research, the cost of stable isotope analyses is funded by Georg Forster Research Fellowship (Alexander von Humboldt Foundation). All transportation and living expenses during the sample collection process are supported by MESA RU (Department of Geology, Faculty of Science, Chulalongkorn University) and by Grants for Development of New Faculty Staff (Ratchadaphiseksomphot Endowment Fund, Chulalongkorn University). We would like to express our gratitude for funding supports by the DMR (Bangkok) and the CNRS-TRF (Caenozoic Biodiversity) Program for the fossil expeditions in Pha Bong, Khok Sung, and Tham Wiman Nakin and by the Thailand Science Research and Innovation (TSRI) grant RTA6080001 (The Prehistoric Population and Cultural Dynamics in Highland Pang Mapha Project) for the archaeological excavations in Tham Lod Rockshelter and Ban Rai Rockshelter, all of which make the fossil material available for this study.