Dietary ecological traits of extinct mammalian herbivores from the last glacial termination at the Pilauco Site, Chile

Abstract

Stable isotopes are a powerful tool for reconstructing the past. However, environmental factors not previously considered can lead to misinterpretations. Our study presents a novel analysis of the feeding behavior of the megafauna that inhabited the Pilauco ecosystem in south-central Chile during the last glacial termination. We analyzed a suite of modern plant and animal samples from closed-canopy forests to establish an isotopic baseline with which to compare stable isotope results from fossil megafauna. Using the modern samples as a reference, the δ¹³C results from the Pilauco megafauna indicate feeding behaviors in forested areas. These results were then calibrated with dental calculus samples and coprolites, which suggest the coexistence of graze- and grass-dominated mixed-feeder diets. The δ¹⁵N values found in Pilauco megafauna are not consistent with modern reference data sets or with the low δ¹⁵N values of extinct proboscideans from other contemporaneous and nearby sites. Probably, the δ¹⁵N values of the Pilauco ecosystem were not primarily affected by climate, but rather by disturbance factors (e.g., grazing effect). Our results indicate that the Pilauco megafauna fed mainly on arboreal vegetation; however, non-isotopic proxies indicate that they were also eating open vegetation (e.g., herbs and grasses).

INTRODUCTION

One of the major comprehensive paleoenvironmental records (e.g., pollen, alkenones, etc.) of the late Pleistocene of South America comes from northwestern Chilean Patagonia (38–42°S, 74–71°W; Moreno et al., 2015). This record has been used to study the effects of interhemispheric paleoclimatic processes (Denton et al., 2010), such as changes in landscape composition that took place during the last glacial maximum and the last glacial termination. Records from this region reveal high abundances of Nothofagus trees in the lowlands between ca. 25 and 11.5 cal ka BP, and trees interspersed in a matrix of herbs and shrubs that are currently commonly found above the Andean tree line, between ca. 25 and 18 cal ka BP. This suggests the existence of a variably open parkland during the last glacial maximum (Moreno et al., 2018). Later, ca. 18 cal ka BP, deglaciation in the northwestern Chilean Patagonia led to a gradual increase in the mean annual temperature, which ultimately resulted in a dense canopy forest as the dominant vegetation formation, which persists at present (Moreno et al., 2015; Moreno, 2020; Supplementary Figure 1).

In northwestern Chilean Patagonia, three late Pleistocene sites have been formally excavated: Monte Verde, Los Notros, and Pilauco (Pino et al., 2020). Pilauco exhibits the highest taxonomic diversity (Gomphotheriidae, Equidae, Camelidae, Cervidae, Mephitidae, Xenarthra indet., Myocastoridae, and Cricetidae; Pino et al., 2013; González et al., 2010; González et al., 2014; Canales-Brellenthin, 2020; Labarca, 2020; Supplementary Figures 1, 2). However, in the region, there are ~20 localities with remains of extinct proboscideans (i.e., Notiomastodon platensis) that range in age from ca. 32 to 12 cal ka BP (González-Guarda et al., 2018).

Some studies undertaken in the area have shown the importance of integrating different proxies beyond the traditional morphological analyses (Sánchez et al., 2004; Aguilera, 2010; Domingo et al., 2012; González-Guarda et al., 2017; 2018). For example, González-Guarda et al., (2018) used a multiproxy approach to provide direct evidence of a trend towards a browsing diet in Notiomastodon platensis, a proboscidean that was linked to a closed forest environment during the last glacial termination. This result fits well with general trends described in recent environmental palynological reconstructions of the forest expansion of Nothofagus and dominance of the cold-tolerant hygrophilous conifers Fitzroya cupressoides, Pilgerodendron uviferum, Podocarpus nubigena, and Saxegothaea conspicua between ca. 15.4 and ca. 12.8 cal ka BP (Moreno, 2020, references therein). However, evidence from Pilauco indicates a more open environment between 16 and 12.8 cal ka BP as attested to by 65% of pollen from the site being of non-arboreal and aquatic types, mainly Poaceae, Asteraceae, Cyperaceae, and Sagittaria (Pino et al., 2019).

Because there is a contrast between the local (more open) and regional (more closed) environment, we could expect new feeding behaviors among the region's herbivorous mammals. Consequently, Pilauco could also be expected to have more than one habitat. Therefore, previous dietary interpretations of megafauna recorded at Pilauco need to be reevaluated. This reevaluation is also based on two aspects: 1) megafauna preserved at Pilauco are characterized by different morpho-functional adaptations; and 2) thus far, most dietary interpretations from Pilauco come from indirect evidence such as stable isotopes, which only indicate a dietary trend but not the diet itself (e.g., grasses, shrubs, or trees). These reasons reinforce the ambiguities in the types of diets inferred for herbivorous mammals around Pilauco, which hinders the development of community and ecosystem-level modeling.

To address this problem, new fossil remains from the Pilauco megafauna were analyzed using a variety of techniques. First, collagen (δ¹³C and δ¹⁵N) and tooth enamel bioapatite (δ¹³C_enamel) were examined using stable isotopes. Second, to determine any exceptions regarding the dietary pattern of gomphotheres previously studied at 16 sites in the region, we incorporated new stable isotope data from Los Notros (40°S) and Curaco de Vélez (42°S), two localities in northwestern Chilean Patagonia. Additionally, we used two new proxies to detect direct dietary evidence: phytolith analysis of coprolites (Pilauco) and microfossil analysis of dental calculus (Pilauco, Los Notros, and Curaco de Vélez). To aid in the interpretation of isotopic datasets from fossils, we also sampled and analyzed modern native plant species and modern specimens from the pudu deer (Pudu puda) that were sourced from the closed-canopy forest to provide a comparative modern baseline for this ecosystem type. The pudu deer was used because it is a strong indicator for a closed-canopy forest (Meier and Merino, 2007). This builds robustness into our research because even though isotopic ranges are an excellent guide for interpretations of past diets and climates, if they are not applied with caution there may be significant errors in their interpretation (Rawlence et al., 2016), especially in the absence of modern baselines (Szpak et al., 2013; Tejada et al., 2020). Our research aims to investigate the dietary patterns of the Pilauco megafauna using a multiproxy approach.

Chronological, geological, stratigraphic, and paleontological setting of Pilauco site

Pilauco is an archaeo-paleontological site located in Osorno city, at the Central Depression of south-central Chile (40°34′S–73°70′W; Fig. 1). The site dates from 17.3 to 4.3 cal ka BP (Pino et al., 2020), although megafauna remains are most prevalent in layers (PB-7 and PB-8) dating from 16.4 cal ka BP to ca. 12.9 cal ka BP (Fig. 2). Fossil mammals found at the site include Mephitidae, Myocastoridae, Gomphotheriidae, Equidae, Cricetidae, Cervidae, Camelidae, and Xenarthra indet. (Labarca, 2020). Gomphotheriidae, Equidae, Camelidae and Xenarthra indet. are the taxa that we used for analysis in this study. Based on morphology, the diet of these taxa has been classified as browsers and grazers. However, the fossil record has cast doubt on many preconceptions about their diets.

Figure 1. Geographic setting of studied samples in southern Chile. Shaded portion of inset map at upper left shows location of larger map. Fossil mammals from Pilauco: Notiomastodon platensis, Equus andium, Xenarthra indet., cf. Hemiauchenia paradoxa. Notiomastodon platensis samples from Los Notros and Curaco de Velez. Sampling areas of modern vegetation: (A) Parque Oncol (Coastal Range); (B) Reserva Biológica Huilo Huilo (Andes Range); (C) Pichirropulli (Central Depression). Modern samples of Pudu puda from: (1) Valdivia; (2, 8) Panguipulli; (3) Río Bueno; (4) Rucatrehua; (5) Máfil; (6) Futrono; (7) Los Lagos; (8) Panguipulli.

Figure 2. Schematic drawing of the stratigraphy of Pilauco site (Pino et al., 2013; Pino et al., 2020) (PB = Pilauco Bajo). Black silhouettes indicate the taxa of fauna from Pilauco.

According to Pino et al. (2013), the basal layers PB-1 to PB-5 at Pilauco correspond with the recently redescribed ignimbrites San Pablo unit (ca. 130–80 cal ka BP; Quiroz et al., 2020). The upper unit at Pilauco includes layers PB-6 to PB-9. The PB-6 bed consists of fluvial pebbles and cobbles representing an alluvial plain that existed between 17.4 and 17.2 cal ka BP. Layers PB-7 and PB-8 are mainly composed of dispersed gravel in a sandy peat matrix. Siliciclastic sediments originated from colluvial processes, while organic matter accumulated in perennial and seasonal wetlands as attested to by pollen and beetle analyses (Pino et al., 2020). PB-7 and PB-8 were deposited in a seasonal wetland following the retreat of the channel of the ancient river Damas (Pino et al., 2020).

The fossil material studied corresponds to two areas: First, where the megafauna layer is found (Fig. 1), there are two paleontological sites (Pilauco and Los Notros) are almost superimposed on each other. The geology of Los Notros is similar to that of the Pilauco, where strata LN-1 and LN-2 are equivalent to layers PB-7 and PB-8 from Pilauco (Lira et al., 2020). Second, the study also includes Curaco de Vélez, a site at the southern limit of the northwestern Chilean Patagonia (42°S) area.

Significance of stable isotope analyses on extinct and extant mammalians

The photosynthetic pathway influences the δ¹³C value of plants (Farquhar et al., 1989). In turn, plant values are fixed in the tissues of herbivores following metabolic processes that involve fractionation in δ¹³C (Koch, 2007). Thus, the bioapatite carbon-isotope composition differs from its source in the diet. This process takes place in ungulate mammals in such a way that δ¹³C values of tooth enamel bioapatite (δ¹³C_enamel) track the δ¹³C values of consumed plants (δ¹³C_diet), offset by ~14‰ due to fractionation associated with carbonate equilibria and metabolic processes (Cerling and Harris, 1999). Because carbon isotope values vary with the photosynthetic pathways of plants, and C₄ plants, which have much higher carbon isotope values, tend to be more prevalent in open habitats, ranges of δ¹³C values can be estimated for herbivorous mammals in different habitats. Following the classification by Domingo et al., (2012), the ranges of δ¹³C_bioapatite values for herbivorous mammals can be estimated for pure C₃ feeders in different habitats (closed-canopy, −20.5‰ to −14.5‰; mesic/woodland, −14.5‰ to −9.5‰; wooded C₃ grassland to open, arid C₃ grassland, −9.5‰ to −6.5‰) and pure C₄ feeders (−1.5‰ to −3.5‰). However, the body mass of each taxon also should be considered because differences have been observed in the standard value of the diet-bioapatite enrichment (Tejada-Lara et al., 2018). With these considerations, more precision has been achieved recently regarding habitat ranges for some extinct South American mammals (e.g., Domingo et al., 2020; Asevedo et al., 2021).

As in tooth bioapatite, bone collagen δ¹³C values also reflect the diets of consumers, but with an offset of ~5‰ (e.g., Koch, 2007). Thus, a collagen δ¹³C value of ~ −22.5‰ is a minimum estimate for individuals that consumed closed-canopy plants only (e.g., Hofman-Kamińska et al., 2018). In herbivorous mammals, observed δ¹⁵N values in collagen are around 3‰ higher than in consumed plants (Koch, 2007). The δ¹⁵N value in plants depends on many factors, including the substrate δ¹⁵N, degree of soil development, availability of nutrients, mycorrhizal associations, and soil acidity (Stevens et al., 2006). Globally, δ¹⁵N values in ecosystems are climatically controlled (Amundson et al., 2003) in such a way that generally, lower δ¹⁵N values are found in cold and/or moist areas (Fox-Dobbs et al., 2008). Warmer, drier, and/or more saline soils raise δ¹⁵N values in plants, which may be the result of metabolic changes in response to water availability and the composition of nitrogen isotopes in individual plants (Hedges et al., 2004). In temperate and semi-arid ecosystems, δ¹⁵N values in vegetation vary from 3‰ to 6‰ (Evans and Ehleringer, 1994). In areas where precipitation is greater than 1000 mm/year, δ¹⁵N values can range from −2‰ to 0‰ in plants that do not fix nitrogen (Heaton, 1987). δ¹⁵N values between −2‰ and 2‰ have also been recorded from atmospheric nitrogen-fixing plants or plants that grow in association with mycorrhizae (Schwarcz et al., 1999).

Although physiological and body size differences between Pilauco megafauna and P. puda must be considered for any final interpretation, several studies have shown that ¹⁵N fractionation is relatively constant in terrestrial mammal herbivores (e.g., Murphy et al., 2006; Kuitems et al., 2015). The isotopic signatures in skeletal tissues are more related to the isotopic variation from consumed plants rather than a change in the isotopic fractionation between herbivores and their diet (Bocherens et al., 2014). In particular, the turnover rate for bone collagen is relatively slow, thus ¹³C and ¹⁵N values of bone collagen reflect the average isotopic values of dietary protein for several years of the animal's lifespan (Hedges et al., 2006).

Analysis of microfossils from dental calculus

Plant microremains in dental calculus can provide direct information on feeding habits (Cordova and Avery, 2017) and a long-term dietary signal (Weyrich et al., 2017). The period involved in formation of dental calculus has not yet been established because the formation processes and their composition can be highly variable between and within individuals (Power et al., 2015). For this reason, it is not possible to determine when, within the lifespan of an animal, a specific plant microremain was ingested (Weyrich et al., 2017). Considering that older individuals present more microremains (Power et al., 2015), it is possible to conclude that dental calculus represents multiple feeding events in the animal's lifespan, assuming there is no replacement or removal of calculus deposits.

Coprolites

Coprolites are fossilized excrements (Jouy-Avantin et al., 2003) that can contain a variety of macroscopic and microscopic remains (Martínez and Yagueddú, 2012). They are one of the most relevant dietary proxies recoverable from sediments because preservation of dietary remains in coprolites is better than in non-fecal deposits (Reinhard and Bryant, 1992). Moreover, a variety of macroscopic and microscopic remains allows researchers to obtain interrelated data sets that enable the reconstruction of diets (Piperno, 2006), behavior (Horrocks et al., 2003), and paleoenvironments (Jouy-Avantin et al., 2003).

MATERIALS AND METHODS

Modern samples

Eight specimens of modern P. puda were analyzed (Supplementary Tables 1–3). These specimens were sampled at the Instituto de Anatomía, Facultad de Ciencias Veterinarias, Universidad Austral de Chile (Valdivia, Chile), where they were stored. Specimens had been obtained by the Centro de Rehabilitación de Fauna Silvestre (CEREFAS; Valdivia, Chile) and collected over several years due to domestic dog attacks, and stored for scientific purposes. Only subadults were selected for analyses, avoiding those specimens potentially influenced by breastfeeding in their isotopic signatures. Weaning occurs ~60 days after birth (Hick, 1969). All sampled specimens showed typical subadult coat coloration (Hershkovitz, 1982) and permanent teeth had been erupted. Also, because different tissues can show different isotopic values (Carleton et al., 2008), a single type of bone in P. puda (femur) was selected in all sampled specimens except one (Supplementary Table 1). The femur was selected because it provides dietary information for a long period during the lifespan of the individual (Hedges et al., 2006). To ensure that femur isotopic values were not significantly different from other tissues (Bocherens et al., 2017), dentin from the dental root was measured (Supplementary Table 2). This tissue represents the isotopic composition of the dentin formation period only (Hick, 1969). Because studies of bone and dentin have shown that the value of diet-to-tissue trophic discrimination for δ¹³C is between 3‰ and 5‰ (Drucker et al., 2003), and between 2‰ and 5‰ for δ¹⁵N (Hedges et al., 2006), each of the values that make up these ranges were applied to observe if there were significant differences between mammals and modern vegetation. In addition, the bioapatite of dental enamel of the same eight individuals of P. puda was analyzed (Supplementary Table 3).

Forty-one samples of modern native plant species located in Región de Los Ríos (39°48′30″S, 73°14′30″W) in Chile were sampled for isotopic analyses (Supplementary Table 4) in closed-canopy forest because this is the habitat of P. puda. Leaves were also sampled from up to 1 m above ground level. All plant samples were identified to species. Modern vegetation collected corresponded mainly to plants consumed by P. puda (Pavez-Fox et al., 2015). In February 2016 (summer), the same plant species were collected in three different areas: 1) Coastal Range (Oncol Park, 39°42′S; 715 m asl, Mean Annual Temperature [MAT] = 10°C, Mean Annual Precipitation [MAP] = 2500 mm); 2) Central Depression (Pichirropulli, 40°04′S; 92 m asl, MAT = 12°C, MAP = 1300 mm); and 3) Andes Range (Neltume, 40°01′S; 604 m asl, MAT = 10°C, MAP = 2000 mm). Because previous studies (Iacumin et al., 2000; Coltrain et al., 2004) detected significant isotopic differences between non-nitrogen-fixing and plants that fix atmospheric nitrogen, the sampling design of this study included both types of plants.

Fossil samples

Twenty megafauna samples from the PB-7 layer (~16.4 to 14.0 cal ka BP) at Pilauco were selected for stable isotope analyses from collagen (δ¹⁵N and δ¹³C_collagen) from the tooth root and postcranial bone elements of Notiomastodon platensis (n = 10), Equus andium (n = 7), Xenarthra indet. (n = 2), and cf. Hemiauchenia paradoxa (n = 1) (Supplementary Table 5). To complement previous isotopic data from tooth enamel bioapatite from Pilauco (González-Guarda et al., 2017; 2018), new N. platensis samples were analyzed (Table 1). In addition, stable isotope analyses (δ¹⁵N and δ¹³C from root collagen and δ¹³C from enamel bioapatite) were performed in N. platensis specimens from Los Notros (40°34′S; Lira et al., 2020) and Curaco de Vélez (42°26′S; this study; Table 1), which are geographically close and contemporaneous sites.

Table 1. (A) δ¹³C values in the bioapatite of the molars of the species Notiomastodon platensis. (B) Radiocarbon dating of two specimens of Notiomastodon platensis. Samples were obtained from the collagen in the dental root of Notiomastodon platensis, and their respective values of δ¹³C (‰, VPDB) and δ¹⁵N (‰, AIR) are shown.

Two new radiocarbon dates from sampled gomphotheres from Los Notros and Curaco de Vélez and all radiocarbon data pertaining to megafauna fossils analyzed in previous studies were considered (Table 1). This study considers N. platensis as the only gomphothere species that inhabited south-central Chile at this time (Mothé et al., 2017).

Three molars of N. platensis from Pilauco (sample MHMOPI/14 layer PB-7), one from Los Notros (sample MHMOP/LN/8), and two from Curaco de Vélez (samples CHI1, CH2) were examined to extract dental calculus and recover plant microfossils (Supplementary Table 6). Four coprolite samples from Pilauco and their corresponding control sediments from grid 14AD (layer PB-7) (Supplementary Table 7) were sampled to perform microfossil analysis.

Vegetation sample analysis

Plant samples were dried for 24 hours in a drying oven at 60 °C to ensure elimination of microorganisms. These samples were analyzed at the Isotopes Biosciences Laboratory (ISOFYS) of the Department of Analytical Application and Physical Chemistry, Faculty of Bioengineering of the University of Ghent, Belgium. Samples were first measured using a PDZ Europe Automated Carbon Nitrogen Analyzer-Solids and Liquids (ANCA-GSL) elemental analyzer interconnected with a Sercon 20-20 IRMS with a SysCon electronic system (SERCON, Cheshire, United Kingdom). Normalizations to the Vienna Pee Dee Belemnite (VPDB) and N₂ air (AIR) scales were done using B2159 sorghum (δ¹³C VPDB = −13.78 ± 0.17‰ and δ¹⁵N AIR = +1.58 ± 0.15‰) calibrated by Elemental Microanalysis to IAEA-CH-6 for δ¹³C (accepted δ¹³C VPDB = −10.449 ± 0.033‰) and IAEA-N-1 (accepted δ¹⁵NAIR = +0.4 ± 0.2‰). A laboratory soil standard was used as quality control (accepted values δ¹³C VPDB = −22.69 ± 0.04‰ and δ¹⁵N AIR = 7.81 ± 0.07‰), and deviation from accepted values was <0.3‰ and 0.5‰ for δ¹³C VPDB and δ¹⁵N AIR, respectively. Typical standard deviation for replicate samples is δ¹³C = 0.2‰, and δ¹⁵N AIR = 0.4‰, resulting in a combined uncertainty on VPDB = 0.3‰ and AIR = 0.5‰.

Collagen sample analysis

Collagen extraction was performed at the Biomolecular Laboratory of the Catalan Institute of Human Palaeoecology and Social Evolution (IPHES; Tarragona, Spain) following the original protocol proposed by Longin (1971) and modified by Bocherens et al. (1991). Bone fragments were cleaned mechanically to remove the surface, while bone shards (~300 mg to 350 mg) were demineralized using 1 M HCl, soaked in NaOH (0.125 M) to remove contaminants, then rinsed with distilled water and gelatinized with 0.01 M HCl at 100°C for 17 h. Samples were then filtered, frozen, and freeze-dried at the Institute of Chemical Research in Catalonia (ICIQ). Collagen samples weighing ~300 μg were analyzed in duplicate using a Thermo Flash 1112 elemental analyzer (EA) coupled with a Thermo Delta V Advantage isotope ratio mass spectrometer (IRMS) with a Conflo III interface at the Institute of Environmental Science and Technology (ICTA) at the Autonomous University of Barcelona, Spain. The international laboratory standard IAEA 600 (caffeine) was used as a control. Analytical error was calculated by measuring replicates of IAEA 600 (6 replicates). The average analytical error was <0.15‰ (1σ) calculated separately for each of the isotopic δ¹³C and δ¹⁵N measures. We used the Ali-j1, Caf-j1, and blk standards, and we used the 3-point linear normalization method (see Werner and Brand, 2001). The VPDB and AIR scales were used as a reference for the δ¹³C and δ¹⁵N values.

Bioapatite sample analysis

Tooth enamel was sampled horizontally in bands. Tooth surfaces were first cleaned with a tungsten abrasive drill bit. Then, samples of enamel were removed by drilling with a diamond bit. One band for each tooth (n = 12 teeth) (fossil and modern samples) were sampled for oxygen and carbon isotope analyses. Powdered enamel samples were chemically treated at the Biomolecular Laboratory of the Institut Català de Paleoecologia Humana i Evolució Social (IPHES). Samples weighed from 3.5 mg to 9.5 mg. Chemical treatment of samples was based on protocols originally proposed by Koch et al., (1997) that were modified in Tornero et al. (2013). Samples were treated for 4 h in 0.1 M acetic acid [CH₃COOH] (0.1 ml solution/0.1 mg of sample), neutralized with distilled water, and freeze-dried. Pretreated powders were analyzed individually on a Thermo Kiel III device interfaced with a MAT Finnagan 253 at the Scientific and Technological centers of the University of Barcelona (CCiTUB), Spain. The samples were reacted in a vacuum with 100% phosphoric acid [H₃PO₄] at 70°C in individual vessels and purified in an automated cryogenic distillation system. δ¹³C values are expressed relative to VPDB. Accuracy and precision of the measurements were checked using two internal laboratory calcium carbonate standards (RC-1 and CECC) normalized to NBS18 and NBS19 international standards. A total of 16 RC-1 and CECC samples were measured (RC-1 expected values +2.83‰ for δ¹³C; CECC expected values −20.78‰ for δ¹³C). The mean analytical precision of RC-1 was +0.01‰ for δ¹³C values and +0.01‰ for CECC.

In this study, a δ¹³C_atmCO2 value of −6.5‰ was used because it is an accepted value for late Pleistocene studies (Tipple et al., 2010). Therefore, stable isotope data of P. puda and modern vegetation were corrected because the modern composition of ¹³C_atmCO2 has a value of −8‰ (Marino and McElroy, 1991). Following previous studies (e.g., Koch, 2007; Metcalfe et al., 2013), an average of the diet-to-tissue trophic discrimination was applied at ~3‰ for δ¹⁵N and 5‰ for δ¹³C values in collagen samples. According to diet-to-tissue trophic discrimination studies (e.g., Koch, 2007), it is possible to estimate an approximate isotopic average value (estimated consumed plants = ECP) for the Pilauco megafauna. This average value was compared with the values of modern vegetation sampled in this study.

Because it has been suggested that body mass (bm) affects the physiological values of carbon enrichment (Tejada-Lara et al., 2018), the equation ɛ* = 2.4 + 0.034 (bm) was applied to obtain the enrichment between bioapatite and the diet of P. puda (ɛ*_{diet-bioapatite}) (Supplementary Table 8). When obtaining the ɛ*_{diet-bioapatite} (δ¹³C = 12‰) value of P. puda (bm: 9.6 kg), it was possible to increase the confidence of the results obtained from the comparisons between δ¹³C_bioapatite values of mammals with different body masses. Thus, it was possible to obtain the ECP value for each mammal by using a specific ɛ*_{diet-bioapatite} value, which is contingent upon body mass, metabolism, and phylogeny (Tejada-Lara et al., 2018). For the analysis of newly obtained bioapatite values in gomphotheres (four samples; Table 1), an enrichment of 14.1‰ (ɛ*_{diet-bioapatite}) was used (Cerling and Harris, 1999; Domingo et al., 2020). However, according to the estimate of gomphothere body mass (6,000 kg), the study by Asevedo et al. (2021) used enrichment of 15‰ (ɛ*_{diet-bioapatite}). Therefore, it was necessary to apply a multiproxy approach to the present study, and to consider the temporal resolution scale of the isotopic proxy to reduce the uncertainty offered by the value of enrichment (ɛ*_{diet-bioapatite}).

A non-parametric Kruskal-Wallis test was used to compare isotopic values when data were not normally distributed. A one-way ANOVA test was used to compare mean values, and any differences were detected by performing a post-hoc Tukey (HSD) test when there was no statistically significant difference between the variances of the groups analyzed (Levene's test). The significance level was set at p = 0.05.

The carbon to nitrogen (C:N) atomic ratio in fossil bone-collagen samples is the most widely used method to determine the degree of diagenesis. According to previous studies, this ratio ranges between about 2.9 and 3.6 in living mammals (DeNiro, 1985; Ambrose and Norr, 1993; Van Klinken, 1999), and values within this range in fossil material may be considered as corresponding to unaltered collagen (Clementz et al., 2009). Moreover, it is expected that isotopically well-preserved collagen would exhibit C and N percentages higher than 13% and 4.8%, respectively (Ambrose, 1990).

Dental calculus analysis

Three molars were examined to extract dental calculus and recover plant microfossils following the methodology described in González-Guarda et al. (2018). Three slides were prepared for each of the samples using Entellan®. The extraction of microfossils from calculus samples was carried out using the chemical processing method defined by Wesolowski et al. (2007). To estimate the quantities of microfossils in dental calculus, a Lycopodium tablet was added to each sample. Lycopodium spores were counted and recorded to calculate the concentration of microfossils using the method of Maher (1981) as modified by Wesolowski et al. (2010). Coprolite samples were processed according to the methodology described by Katz et al. (2010). The slides both for tooth calculus and coprolites were examined under a polarized light microscope at 200× and 630× magnification.

RESULTS

Isotope sample preservation

Results of stable isotope analyses and collagen quality indicators are reported in Table 1 and Supplementary Tables 1, 2, and 5. Collagen was successfully extracted from all samples analyzed. All samples presented higher %C and %N values than what is recommended as acceptable limits. The atomic C:N ratio ranged from 3.0 to 3.4 in Pilauco fossil samples, and from 3.0 to 3.5 in the modern P. puda specimens.

Stable isotope values in modern samples

In femur samples of P. puda (n = 8), δ¹³C values ranged from −25.2‰ to −22.2‰, while δ¹⁵N values ranged from 5‰ to −1.6‰ (Table 2). Mean δ¹³C and δ¹⁵N values were −23.8‰ ± 1‰ and 3.2‰ ± 2‰, respectively. When applying δ¹³C diet-to-tissue trophic discrimination of ~5‰, the mean ECP of δ¹³C_ECP value was −30.3 ± 0.9‰. When applying δ¹⁵N diet-to-tissue trophic discrimination of ~3‰, the mean ECP of δ¹⁵N_ECP was 0.2 ± 2‰. Regarding δ¹³C_bioapatite values in the enamel of P. puda, the mean value was −16.4 ± 1.1‰ (Table 3). The minimum and maximum values were −18.3‰ and −15.2‰, respectively.

Table 2. Stable isotope analysis of the femur of the deer species Pudu puda. δ¹³C (‰, VPDB) and δ¹⁵N (‰, AIR) values of estimated consumed plants (ECP) were calculated by applying trophic corrections of 5‰ and 3‰, respectively, and ranges of values of diet-to-tissue trophic discrimination are incorporated (3–5‰ for δ¹³C; 2–5‰ for δ¹⁵N). We used the maximum cut-off value (δ¹³C = −27.5 ‰) commonly used to identify specimens of herbivorous mammals that consume a significant number of plants growing in closed-canopy conditions (Drucker et al., 2003, 2008; Hofman-Kamińska et al., 2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ. All data normalized based on preindustrial atmospheric conditions (δ¹³C_atmCO2 = −6.5‰).

Table 3. δ¹³C enamel isotopic results from the deer species Pudu puda. Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ; ECP: estimated consumed plants. All data normalized based on preindustrial atmospheric conditions (δ¹³C_atmCO2 = −6.5‰).

When comparing the mean of δ¹³C and δ¹⁵N values between the femur and the dentin of P. puda, a similarity is observed (δ¹³C_femur= -23.8‰; δ¹³C_dentine= -24‰; δ¹⁵N_femur = 3.2‰; δ¹⁵N_dentine = 3.6‰; Tables 3, 4). However, looking at δ¹⁵N values, a difference of up to 3‰ between the dentin and the femur of the same individual was observed (sample 16139; Supplementary Tables 1, 2).

Table 4. Stable isotope analysis of dentin from the deer species Pudu puda. δ¹³C (‰, VPDB) and δ¹⁵N (‰, AIR) values of estimated consumed plants (ECP) were calculated by applying trophic corrections of 5‰ and 3‰, respectively, and ranges of values of diet-to-tissue trophic discrimination are incorporated (3‰–5‰ for δ¹³C; 2‰–5‰ for δ¹⁵N). We used the maximum cut-off value (δ¹³C = −27.5 ‰) commonly used to identify specimens of herbivorous mammals that consume a significant number of plants growing in closed-canopy conditions (Drucker et al., 2003, 2008; Hofman-Kamińska et al., 2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ. All data normalized based on preindustrial atmospheric conditions (δ¹³C_atmCO2 = −6.5‰).

The mean δ¹³C value of estimated consumed plants (ECP) (+12‰) was −28.4‰. For modern vegetation, the ranges of δ¹³C values for Oncol Park (Coastal Range), Neltume (Andes Range), and Pichirropulli (Central Depression) were −36.6‰ to −26.3‰; −38.7‰ to −28.9‰; and −33.1‰ to −28.4‰, respectively. The ranges of δ¹⁵N values for Oncol Park (Coastal Range); Neltume (Andes Range), and Pichirropulli (Central Depression) were −5‰ to 3.5‰; −9.9‰ to 2.1‰; and −4.9‰ to 3.5‰, respectively (Table 5). Figure 3 shows mean δ¹³C and δ¹⁵N values for each area. The difference among plant δ¹³C values from Oncol Park, Neltume, and Pichirropulli are statistically significant (one-way ANOVA, F = 3.9, p = 0.03) and the post-hoc Tukey test revealed that all three areas differ from each other. Regarding δ¹⁵N values, the only group in which no significant differences were found was the Neltume-Pichirropulli group (one-way ANOVA, F = 3.1, p = 0.06). The δ¹³C values in modern vegetation did not show significant differences between nitrogen fixation and non-fixing plants, except in the sample of a lichen species (Stereocaulon ramulosum) that is nitrogen-fixing plant, and showed the highest negative value of the samples (−26.3‰).

Figure 3. The plot of results of isotopic analyses from various sources discussed herein showing the mean of δ¹³C (‰, VPDB) and δ¹⁵N (‰, AIR) values in Pilauco Megafauna, modern vegetation, Pudu puda, gomphotheres from the surroundings of Pilauco (González-Guarda et al. 2018), megafauna from the Southern Patagonia region of Chile (Prevosti and Martin, 2013; equids, xenarthrans, camelids, ratites) and megafauna from the Pampean region (Bocherens et al., 2016; xenarthrans, equids, cervids, macrauchenids, toxodons, rodents). Vegetation δ¹³C (‰, VPDB), δ¹⁵N (‰, AIR) values were calculated applying trophic corrections of 5‰ and 3‰ respectively. The error bars represent the values ± 1 standard deviation.

Table 5. Summary of stable isotope data from modern plant samples collected in northwestern Chilean Patagonia. (A) δ¹³C (‰, VPDB), (B) δ¹⁵N (‰, AIR) values. Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ. All data normalized based on preindustrial atmospheric conditions (δ¹³C_atmCO2 = −6.5‰).

Values of stable isotopes in fossil megafauna

The δ¹³C values for N. platensis, E. andium, and Xenarthra indet. ranged from −23.4‰ to −22.8‰, −22.2‰ to −20.9‰, and −22‰ to −21.8‰ (VPDB), respectively. δ¹⁵N values for N. platensis, E. andium, and Xenarthra indet. ranged from 4.9‰ to 9.2‰, 3.9‰ to 8.7‰, and 7.2‰ to 8.1‰ (AIR), respectively. Only one sample was analyzed for H. paradoxa (δ¹³C= −21.3‰; δ¹⁵N = 5.7‰; Table 6). Figure 4 shows mean δ¹³C and δ¹⁵N values for each taxon from Pilauco. Results fit well with previous analyses in Pilauco using similar procedures (Aguilera, 2010, n = 3; González-Guarda et al., 2017, n = 17). Table 7 shows the statistical analysis of data from previous studies and this study. The Kruskal-Wallis test shows significant differences in δ¹³C values of N. platensis, E. andium, cf. H. paradoxa and Xenarthra indet. (p < 0.001). For δ¹⁵N values, significant differences are documented in all taxa (p = 0.009).

Figure 4. Average bone collagen values for δ¹³C (‰, VPDB), indicating a more closed to open environment, and δ¹⁵N (‰, AIR), indicating dryer/warmer to wetter/colder climate, for the taxa studied at Pilauco site (Notiomastodon platensis, cf. Hemiauchenia paradoxa, Equus andium, Xenarthra indet.). The error bars represent the values ± 1 standard deviation.

Table 6. Summary of stable isotope data obtained from the Pilauco site megafauna. The δ¹³C (‰, VPDB) and δ¹⁵N (‰, AIR) values of estimated consumed plants (ECP) were calculated by applying trophic corrections of 5‰ and 3‰, respectively, and ranges of values of diet-to-tissue trophic discrimination are incorporated (3‰–5‰ for δ¹³C; 2‰–5‰ for δ¹⁵N). We used the maximum cutoff value (δ¹³C = −27.5 ‰) commonly used to identify specimens of herbivorous mammals that consume a significant number of plants growing in closed-canopy conditions (Drucker et al., 2003, 2008; Hofman-Kamińska et al., 2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ.

Table 7. Summary of stable isotope data and statistical tests results (non-parametric Kruskal-Wallis test) from the Pilauco site megafauna. (A) Collagen δ¹³C (‰, VPDB) values, (B) Collagen δ¹⁵N (‰, AIR) values. Stable isotope values of this study were analyzed together with values from Aguilera (2010) and González-Guarda et al. (2017). The δ¹³C (‰, VPDB) and δ¹⁵N (‰, AIR) values of estimated consumed plants (ECP) were calculated by applying trophic corrections of 5‰ and 3‰, respectively, and ranges of values of diet-to-tissue trophic discrimination are incorporated (3‰–5‰ for δ¹³C; 2‰–5‰ for δ¹⁵N). We used the maximum cutoff value (δ¹³C = −27.5 ‰) commonly used to identify specimens of herbivorous mammals that consume a significant number of plants growing in closed-canopy conditions (Drucker et al., 2003, 2008; Hofman-Kamińska et al., 2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ.

Once the δ¹³C diet-to-tissue trophic discrimination of 5‰ was applied, the average value of the estimated consumed plants (ECP = δ¹³C_ECP) was −28.2 ± 0.5‰ for N. platensis; −26.7 ± 0.6‰ for E. andium; −27.1 ± 0.9‰ for cf. H. paradoxa, and −26.7 ± 0.3‰. for Xenarthra indet. When applying δ¹⁵N diet-to-tissue trophic discrimination of 3‰, mean ECP δ¹⁵N_ECP values were 3.7 ± 0.9‰ for N. platensis, 2.1 ± 1.6‰ for E. andium, 2.3 ± 0.3‰ for cf. H. paradoxa, and 3.5 ± 0.3‰ for Xenarthra indet.

Regarding δ¹³C_bioapatite values in samples of N. platensis from Pilauco, Curaco de Vélez, and Los Notros, the mean value was −13.8 ± 0.5‰. The δ¹³C_collagen and δ¹⁵N_collagen values from Los Notros were −22.2‰ and 6.6‰, respectively, and the δ¹³C_collagen and δ¹⁵N_collagen values from Curaco de Vélez were −22.1‰ and 2.3‰, respectively (Table 1).

Analysis of dental calculus and fossil coprolites

The dental calculus analysis was conducted on the molars (n = 4) of N. platensis. The microfossil count included phytoliths, sponge spicules, fragments of diatoms, and non-silica tissues. The most abundant microfossils were phytoliths, including herbaceous and arboreal morphotypes (total phytoliths for samples: MHMOPI/14 = 123681 mf/g; MHMOP/LN/8 = 787736 mf/g; CHI1 = 34639 mf/g; CHI2 = 399249 mf/g; where mf/g indicates microfossils per gram; Supplementary Table 6). Molars MHMOPI/14, CHI1 and CHI2 presented a higher percentage of arboreal morphotypes. Herbaceous morphotypes predominated in the MHMOP/LN/8 sample (Fig. 5), followed by diatoms and sponge spicules.

Figure 5. Microfossils from dental calculus, coprolites, and sediment samples. (A) Stipa sp. morphotypes; (B) Chusquea sp. ‘dumb-bell’ morphotypes; (C) Panicoideae ‘dumb-bell morphotypes; (D) Chusquea sp. ‘bulliform’ cell; (E, F) Arboreal irregular morphotypes; (G) sponge spicule; (H) diatom; (I) spherulite.

Coprolite analysis shows the coprolites contained mostly herbaceous plants. The record of phytoliths contained in the sediments of the grid 14AD confirms an open landscape in the area of Pilauco, with C₃-type grasslands (e.g., Stipa spp.) and Chusquea spp. shrubland. Phytoliths and diatoms were present in coprolite and control sediment samples (Pi494, Pi505, Pi507), whereas spherulites were only present in coprolite sample C588 (Supplementary Table 7), and were identified as spherulites of Camelidae according to the description of microfossils for this taxon (Korstanje, 2002). Phytoliths in both the control sediments and coprolites corresponded to arboreal and herbaceous morphotypes, although a portion of the phytoliths found in the coprolites was unidentifiable due to either taphonomic processes or digestive degradation. The identified phytoliths in coprolites were mainly herbaceous morphotypes (>95%) belonging to Poaceae, Cyperaceae, and Pterydophyta. Three subfamilies of Poaceae were identified, including Pooideae (C₃ plant), Bambusoideae (C₃ plant), and Panicoideae (C₄ plant). Control samples showed a similar proportion of arboreal/herbaceous morphotypes, and presence of the same subfamilies. Panicoideae was residual, occurring only in samples C585 and PI494 (<1%). Panicoid species in modern Chile are C₄ plants that are adapted to temperate environments, can tolerate high tree cover (e.g., Imperata condensate; Zuloaga et al., 2019) and wide latitudinal and altitudinal distribution ranges (e.g., Paspalum sp.; Rodríguez et al., 2018). Diatoms were present in all samples and in variable proportions that were higher in sediment than coprolites. Spherulites found in samples are CaCO₃ crystals that form during the digestive process (Canti, 1998) in the intestines of certain animals (Brochier et al., 1992).

DISCUSSION

Feeding behavior and habitat of extinct species

Assuming a cutoff value of −27.5‰ for δ¹³C_collagen values as delineating a plant-diet based on open areas versus forested areas (or δ¹³C_{bioapatite (diet)} = −27.2‰; Tejada et al., 2020), and considering the average values of δ¹³C from P. puda (femur= −28.8‰; dentine = –29.0‰; dental enamel= −28.2‰), results suggest that the megafauna species analyzed preferred different environments. Equids, xenarthrans, and camelids may have regularly occupied more open areas, while gomphotheres may have occupied more forested areas in the Pleistocene ecosystem of the Pilauco area (Fig. 4; Table 7). Previous studies of dental microwear and dental calculus confirm this interpretation for gomphotheres in northwestern Chilean Patagonia (González-Guarda et al., 2018). Particularly interesting are the δ¹³C_bioapatite values (MHMOPI/628a = −14‰ and MHMOPI/627a = −13‰) in the Pilauco gomphotheres, which indicate more closed conditions. However, although slightly different from δ¹³C_collagen values, δ¹³C_bioapatite values (−14 ± 0.8‰) in the dental enamel of E. andium also indicating feeding in a more closed environment. This does not translate into the absence of grass consumption during the lifespan of the specimens (δ¹³C values for woodland-mesic C₃ grassland, −14.5‰ to −9.5‰; Domingo et al., 2012).

Observed differences in the δ¹⁵N values between N. platensis (6.7‰ ± 0.9) and E. andium (5‰ ± 1.7) may indicate the consumption of different types of plants (Table 7). However, as in N. platensis samples from Pilauco (δ¹⁵N_dentine = 8.5‰; González-Guarda et al., 2017; δ¹⁵N_bone = 9.2‰; this study), some E. andium samples with δ¹⁵N_bone values (8.7‰) differ from modern vegetation. These values could indicate foraging in grassland areas (Bocherens, 2003), but as stated above, most of the δ¹³C values obtained suggest feeding in more wooded areas. What factors could have caused high δ¹⁵N values? This anomaly is not only observed at Pilauco but is also seen in the gomphotheres of central Chile (31–42°S), which exhibit high variability in δ¹⁵N_dentine values (from 14.2‰ to 1.3‰; González-Guarda et al., 2018), with values as low as in mastodons and as high as in mammoths in some regions of the northern hemisphere (e.g., Metcalfe et al., 2013). Interestingly, these differences are also observed in the two gomphotheres sampled from elsewhere in the region (Table 1), in Los Notros (40° 34′S; 16,021–15,643 cal yr BP, 2 sigma range; δ¹⁵N_dentine = 6.6 ‰), which is located only 60 meters from Pilauco, and in Curaco de Vélez (42°26′S; 13,569–13,396 cal yr BP, 2 sigma range; δ¹⁵N_dentine = 2.3 ‰), 200 km from Pilauco (for a more comprehensive chronological context for the study area, see Supplementary Table 8).

Again, the δ¹⁵N_dentine value (Supplementary Table 5) of a new gomphothere (Los Notros) from the Pilauco area does not overlap with δ¹⁵N values of modern vegetation (Supplementary Table 4), while the gomphothere sampled in the Curaco de Vélez has a low δ¹⁵N value similar to those observed in gomphotheres from the northwestern Patagonia Chilean (Table 8). Therefore, regardless of the specific factors that could affect δ¹⁵N variability (e.g., grazing intensity, animal manure, fire regimes, coprophagy, starvation, etc.), future studies on these values should also consider aspects such as volcanic soils, proximity to the Pacific Ocean, the pronounced orographic climates typical of Chile or a combination of these factors.

Table 8. Statistical summary of molars at gomphotheres that do not come from the Pilauco site (38°–42°S) (from González-Guarda et al., 2018 The δ¹³C (‰, VPDB) and δ¹⁵N (‰, AIR) values of estimated consumed plants (ECP) were calculated by applying trophic corrections of 5‰ and 3‰, respectively, and ranges of values of diet-to-tissue trophic discrimination are incorporated (3‰–5‰ for δ¹³C; 2‰–5‰ for δ¹⁵N). We used the maximum cutoff value (δ¹³C = −27.5 ‰) commonly used to identify specimens of herbivorous mammals that consume a significant number of plants growing in closed-canopy conditions (Drucker et al., 2003, 2008; Hofman-Kamińska et al., 2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ.

The absence of overlap of δ¹⁵N values between Pilauco megafauna and modern vegetation could indicate a diet based on plants found in grassland areas. However, δ¹³C values indicate a diet from a more wooded environment. When comparing the δ¹³C values for the megafauna (ECP) with those of modern vegetation, it is very likely that the megafauna fed in open forest areas, where the canopy effect is less intense than that of a closed forest. Nevertheless, it is not possible to state with certainty that megafauna from Pilauco had a leaf-browsing diet for two reasons: first, there were equids with high δ¹³C_collagen values (e.g., −20.9‰), which could indicate the consumption of herbs because they lived in more open areas. Second, results from two proxies from this study show unequivocal evidence of food consumption in grassland areas. Considering the temporal resolution of the dietary proxies, in the shorter term, coprolites show grazing behavior; and in the longer term, the analysis of microfossils of dental calculus from one gomphothere sample shows a grass-dominated mixed-feeder diet (MHMOP/LN8, Los Notros).

Interestingly, the δ¹³C_bioapatite value of the MHMOP/LN8 (Los Notros) sample shows a value of −14‰ (Table 1), which indicates a wooded environment. This is inconsistent with the dental calculus microfossils from the same specimen (MHMOP/LN/8), but this last scenario may indicate different dietary behaviors in the life history of the MHMOP/LN/8 specimen, perhaps under the influence of an environment characterized by the assemblage Nothofagus-Myrtaceae-Poaceae (ca. 16.3–15.4 cal ka BP; Moreno, 2020). However, microfossil analysis of dental calculus in another contemporary gomphothere (sample MHMOP/PI/14; Pilauco) indicates a leaf-browsing diet. In addition, a new analysis of δ¹³C_bioapatite from sample MHMOP/PI/16 (Pilauco) shows a value of −14‰ (Table 1), consistent with results from dental microwear and stable isotopes analysis from Pilauco (samples MHMOP/PI/627, MHMOP/PI/628; González-Guarda et al., 2018).

These results are compatible with the environment interpreted between ca. 16 cal ka BP to 12.8 cal ka BP being a domain of non-arboreal pollen in Pilauco surroundings (~35% arboreal, 65% non-arboreal), which would characterize its landscape as having been a ‘mosaic habitat,’ defined as a range of different habitat types, scattered across and interspersed within a given area (Elton, 2008). Thus, a more complex feeding behavior is revealed in northwestern Chilean Patagonia. Perhaps the new dietary behavior related to the appearance of grasslands, environments that were not exclusive to the Pilauco area (Supplementary Figure 2). In a closed-canopy context, many of these less-wooded areas can be explained by the low quality of the soil, flooding, and geomorphology, and not necessarily due to a cold event such as the Antarctic Cold Reversal. Particularly, periglacial processes would have disturbed the soils of the Central Depression (40°S–43°S; Veit, 1994). Therefore, for the region to support high animal biomass, megafauna should have needed heterogeneous environments and high plant diversity to avoid some aspects of competitive exclusivity.

Finally, the dietary categories discovered for the late Pleistocene ecosystem of the area can now be understood with greater complexity, and are in line with the increase in the fossil record of mammals with generalist behavior but with specialized morphologies (Rivals et al., 2019). Here, the studied taxa are good examples of this variability. N. platensis is considered a browser (Fox and Fisher, 2004), but previous studies have shown that this species could also exist as a grazer and mixed-feeder elsewhere in South America (e.g., Asevedo et al., 2012). Likewise, H. paradoxa is characterized by grazing behavior, inhabiting more arid habitats (Menégaz and Ortiz Jaureguizar, 1995). However, studies of Hemiauchenia from other areas of America have shown a browser and mixed-feeder diet (Feranec, 2003; Semprebon and Rivals, 2010), which is also demonstrated in specimens sampled from Pilauco. In addition, Equus andium occupied different environments in South America; from mixed C₃–C₄ grassland (Ecuador; Domingo et al., 2012) to woodland mesic C₃ grassland and closed-canopy forest (Chile; González-Guarda et al., 2017). The Xenarthra teeth are morphologically simple and do not provide unambiguous information about their diet (Bargo and Vizcaíno, 2008). However, skull and jaw morphology differences among these taxa yielded evidence for possible feeding mechanisms and therefore food composition (Vizcaíno et al., 2011). Recently, an isotopic study on bulk collagen showed that xenarthrans were exclusively herbivorous (Bocherens et al., 2017). Nevertheless, previous studies have indicated that this group of animals may have had insectivorous (Genise and Fariña, 2012) and carnivorous representatives (Fariña and Blanco, 1996).

Therefore, according to previous results and results from this study, the presence of mammals that were morphologically adapted to a grassland environment while inhabiting a more wooded environment would not be rare, although their presence in these environments must be explained by knowledge of their biogeographic history. It is likely that a dietary diversity characterized the Pleistocene fauna in the northwestern Chilean Patagonia because populations of Nothofagus trees were abundant in the lowlands, and because the trees were interspersed within a matrix of herbs and shrubs during times of cold climate. During warmer episodes, the closed-canopy forest could have covered most of the region. Indeed, the fossil record of the P. puda, considered a strong indicator of closed forest (i.e., cf. Pudu; González et al., 2014), and the presence of herbivorous mammalian grazers, such as camelids and equids in Pilauco, are examples of the dietary complexity/flexibility that this fauna ensemble could have had throughout the last ice ages of the northwestern Chilean Patagonia.

Landscape reconstruction from isotopic dietary patterns

Both collagen and bioapatite carbon isotope data from the Pilauco megafauna point to closed-canopy conditions. The mean δ¹³C_ECPcollagen value for the dataset of the Pilauco megafauna (−27.5 ± 0.8‰) matches the maximum δ¹³C cutoff value (−27.5‰) commonly used to identify herbivorous mammals that consume a significant number of plants growing under closed-canopy conditions (Drucker et al., 2003, 2008; Hofman-Kamińska et al., 2018). This fits well with δ¹³C values measured in several tooth enamel bioapatite samples from megafauna (mean δ¹³C_bioapatite value = −13.5‰ in ECP) from northwestern Chilean Patagonia (González-Guarda et al., 2017).

Comparison between δ¹³C_ECP values from Pilauco megafauna and δ¹³C results from modern vegetation (mean −33 ± 2.6‰) show differences (Fig. 3), although some low results from the megafauna overlapped with the high results from modern vegetation. The absence of differences in δ¹³C values between nitrogen fixation and non-fixing in modern plants indicates it is not possible to detect the consumption of both types of plants by fossil megafauna even though differences in δ¹³C values between both types of plants have been reported elsewhere (e.g., Fox-Dobbs et al., 2008).

ECP nitrogen values in Pilauco megafauna are not in agreement with mean δ¹⁵N values of modern vegetation (Parque Oncol: −1 ± 2.7‰; Reserva Huilo Huilo: −3.6 ± 3.7‰; Pichirropulli: −1.2 ± 2.7‰; Table 5). In contrast, the mean value of modern P. puda is very close (mean δ¹⁵N = 0.2 ± 2.2‰) to that of modern vegetation. The mean ECP nitrogen isotope values in the Pilauco megafauna are closest to those observed in modern temperate and semi-arid ecosystems (δ¹⁵N from 3‰ to 6 ‰; Evans and Ehleringer, 1994), suggesting that the Antarctic Cold Reversal did not significantly alter the Pilauco ecosystem. Nevertheless, some exceptions are observed: when applying ECP, some values of E. andium are very close to 0‰ (cold/humid sites); some values of modern vegetation indicate warm or semi-arid ecosystem; and modern individuals of P. puda included values similar to the megafauna bones. Consequently, δ¹⁵N data from modern vegetation and fauna compels us to cautiously apply environmental ranges established based on δ¹⁵N.

Comparison between δ13C_ECPcollagen values from Pilauco megafauna (mean = −27.5 ± 0.8‰) and P. puda (mean = δ¹³C_ECPcollagen = −28.8 ± 1‰) indicate very similar results. This similarity is also found in P. puda δ¹³C_{ECPbioapatite} values (mean −28.4 ± 1.1‰). However, P. puda results show that this species inhabits more open areas than those where the modern vegetation was collected (mean δ¹³C_{modern vegetation} value = −33 ± 2.6‰).

Interestingly, δ¹⁵N_ECP values of P. puda (mean = 0.2 ± 1.8‰) and sampled gomphotheres from Pilauco (δ¹⁵N_ECP = 3.7 ± 0.9‰) show important differences, but δ¹⁵N_ECP P. puda values are very similar to those retrieved from gomphotheres (n = 19; mean = 0.6 ± 1.8‰; Table 8) from northwestern Chilean Patagonia (38°– 42°S, González-Guarda et al., 2018). Mean δ¹³C_ECPcollagen values of Pilauco gomphotheres (mean values = −28.2 ± 0.5‰; n = 19) and the other gomphotheres samples (mean value = −27.8 ± 0.6‰; n = 19) are very similar. This matching pattern in carbon values but not in nitrogen, also documented between both groups of gomphotheres, may be a strong indication that the influence of a non-climatic factor operating on δ¹⁵N signatures from the fossil megafauna of Pilauco.

CONCLUSIONS

Our multiproxy approach has found new feeding behaviors in late Pleistocene megafauna in Chilean Patagonia, and we interpret that the Pilauco area was an ecosystem that had more than one habitat. Stable isotope values of δ¹³C from the Pilauco megafauna derived from collagen and bioapatite indicate that the Pilauco ecosystem was forested, while δ¹⁵N values only indicate that the soil of the Pilauco ecosystem was relatively altered, probably due to grazing and trampling effects. We developed our interpretation based on data from modern vegetation and fauna. Therefore, we might expect that Pleistocene plants from both open and closed habitats in Pilauco could have had high δ¹⁵N values.

Carbon isotope data from the Pilauco megafauna support previous interpretations that these animals were leaf browsers. Nonetheless, differences between fossil megafauna and both the modern P. puda and plant δ¹⁵N values suggest that the Pilauco megafauna were using more open woodland habitats than the purely closed-canopy environments in which modern P. puda are found. It is likely that Pleistocene megafauna visited the ecosystem of the Pilauco area mainly to consume water and to feed on open forests.

However, coprolites from Pilauco and dental calculus from Los Notros reveal new dietary categories for the area, such as grazer and mixed-feeder. This confirms that the megafauna also consumed food from open environments. The gomphothere sample from Curaco de Vélez reinforces the browsing pattern of the taxon in northwestern Chilean Patagonia. These results indicate that it is crucial to incorporate proxies that show direct evidence of paleodiet. Recent proposals suggesting specific enrichment for each species of mammal have been an advancement to aid in determining a more realistic ranges of habitats (e.g., Tejada-Lara et al., 2018). As this study and that of González-Guarda et al. (2018), have shown, regardless of the value of the enrichment applied (+14‰ or +15‰) in some gomphotheres (i.e., from Los Notros), coherence in the interpretation of the diet is not observed when more proxies are applied to the same specimen (i.e., dental calculus).

Regarding nitrogen signatures, samples from Pilauco specimens show higher δ¹⁵N values than those observed in modern samples representing closed-canopy woodland ecosystems. This difference is also observed between gomphotheres from the Pilauco area and other known gomphothere samples from northwestern Chilean Patagonia with very low δ¹⁵N values, commonly interpreted as evidence of colder and wetter climatic conditions during the last glacial termination. The challenge of future biogeochemical studies will be to determine the factors (perhaps synergistic) driving these large variations in δ¹⁵N values (particularly in gomphotheres) that vary significantly between localities, e.g., from 9.2‰ (40°S, Pilauco) to 1.3‰ (41°S, Monte Verde). This high nitrogen isotope variability is not only found in the paleontological samples, it is also evident in modern P. puda samples (δ¹⁵N = −1.6‰ to 5‰), which is surprising since the modern puda sampled inhabit the same environment and climate. Therefore, evidence of these large differences in δ¹⁵N values could increase the uncertainty of paleodietary reconstructions based solely on δ¹⁵N_collagen. However, if a particular study is supported by a robust isotopic baseline on current samples, δ¹⁵N_collagen values can be a very powerful tool to detect a microhabitat or alterations in a paleoecosystem.

Finally, given the scarcity of formally excavated and studied sites in northwestern Chilean Patagonia and the excellent preservation of collagen, Pilauco will continue to have special relevance for future studies at different levels of the paleoecological organization and those related to methodological and taphonomic biases.