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., Reference Moreno, Denton, Moreno, Lowell, Putnam and Kaplan2015). This record has been used to study the effects of interhemispheric paleoclimatic processes (Denton et al., Reference Denton, Anderson, Toggweiler, Edwards, Schaefer and Putnam2010), 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., Reference Moreno, Videla, Valero-Garcés, Alloway and Heusser2018). 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., Reference Moreno, Denton, Moreno, Lowell, Putnam and Kaplan2015; Moreno, Reference Moreno2020; Supplementary Figure 1).
In northwestern Chilean Patagonia, three late Pleistocene sites have been formally excavated: Monte Verde, Los Notros, and Pilauco (Pino et al., Reference Pino, Martel-Cea, Vega, Fritte, Soto-Bollmann, Pino and Astorga2020). Pilauco exhibits the highest taxonomic diversity (Gomphotheriidae, Equidae, Camelidae, Cervidae, Mephitidae, Xenarthra indet., Myocastoridae, and Cricetidae; Pino et al., Reference Pino, Chávez–Hoffmeister, Navarro-Harris and Labarca2013; González et al., Reference González, Prevosti and Pino2010; González et al., Reference González, Labarca, Chavez-Hoffmeister and Pino2014; Canales-Brellenthin, Reference Canales-Brellenthin2020; Labarca, Reference Labarca, Pino and Astorga2020; 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., Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018).
Some studies undertaken in the area have shown the importance of integrating different proxies beyond the traditional morphological analyses (Sánchez et al., Reference Sánchez, Prado and Alberdi2004; Aguilera, Reference Aguilera2010; Domingo et al., Reference Domingo, Prado and Alberdi2012; González-Guarda et al., Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí2017; Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018). For example, González-Guarda et al., (Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018) 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, Reference Moreno2020, 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., Reference Pino, Abarzúa, Astorga, Martel-Cea, Cossio-Montecinos, Navarro and Lira2019).
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 (δ13C and δ15N) and tooth enamel bioapatite (δ13Cenamel) 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, Reference Meier and Merino2007). 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., Reference Rawlence, Wood, Bocherens and Rogers2016), especially in the absence of modern baselines (Szpak et al., Reference Szpak, White, Longstaffe, Millaire and Sánchez2013; Tejada et al., Reference Tejada, Flynn, Antoine, Pacheco, Salas-Gismondi and Cerling2020). 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., Reference Pino, Martel-Cea, Vega, Fritte, Soto-Bollmann, Pino and Astorga2020), 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, Reference Labarca, Pino and Astorga2020). 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., Reference Pino, Chávez–Hoffmeister, Navarro-Harris and Labarca2013; Pino et al., Reference Pino, Martel-Cea, Vega, Fritte, Soto-Bollmann, Pino and Astorga2020) (PB = Pilauco Bajo). Black silhouettes indicate the taxa of fauna from Pilauco.
According to Pino et al. (Reference Pino, Chávez–Hoffmeister, Navarro-Harris and Labarca2013), 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., Reference Quiroz, Mella, Moreno, Duhart, Carrasco and Miralles2020). 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., Reference Pino, Martel-Cea, Vega, Fritte, Soto-Bollmann, Pino and Astorga2020). 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., Reference Pino, Martel-Cea, Vega, Fritte, Soto-Bollmann, Pino and Astorga2020).
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., Reference Lira, Labarca, Fritte, Oyarzo, Pino, Pino and Astorga2020). 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 δ13C value of plants (Farquhar et al., Reference Farquhar, Ehleringer and Hubick1989). In turn, plant values are fixed in the tissues of herbivores following metabolic processes that involve fractionation in δ13C (Koch, Reference Koch, Michener and Lajtha2007). 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 δ13C values of tooth enamel bioapatite (δ13Cenamel) track the δ13C values of consumed plants (δ13Cdiet), offset by ~14‰ due to fractionation associated with carbonate equilibria and metabolic processes (Cerling and Harris, Reference Cerling and Harris1999). Because carbon isotope values vary with the photosynthetic pathways of plants, and C4 plants, which have much higher carbon isotope values, tend to be more prevalent in open habitats, ranges of δ13C values can be estimated for herbivorous mammals in different habitats. Following the classification by Domingo et al., (Reference Domingo, Prado and Alberdi2012), the ranges of δ13Cbioapatite values for herbivorous mammals can be estimated for pure C3 feeders in different habitats (closed-canopy, −20.5‰ to −14.5‰; mesic/woodland, −14.5‰ to −9.5‰; wooded C3 grassland to open, arid C3 grassland, −9.5‰ to −6.5‰) and pure C4 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., Reference Tejada-Lara, MacFadden, Bermudez, Rojas, Salas-Gismondi and Flynn2018). With these considerations, more precision has been achieved recently regarding habitat ranges for some extinct South American mammals (e.g., Domingo et al., Reference Domingo, Tomassini, Montalvo, Sanz-Pérez and Alberdi2020; Asevedo et al., Reference Asevedo, Ranzi, Kalliola, Pärssinen, Ruokolainen, Cozzuol and Rodrigues do Nascimento2021).
As in tooth bioapatite, bone collagen δ13C values also reflect the diets of consumers, but with an offset of ~5‰ (e.g., Koch, Reference Koch, Michener and Lajtha2007). Thus, a collagen δ13C value of ~ −22.5‰ is a minimum estimate for individuals that consumed closed-canopy plants only (e.g., Hofman-Kamińska et al., Reference Hofman-Kamińska, Bocherens, Borowik, Drucker and Kowalczyk2018). In herbivorous mammals, observed δ15N values in collagen are around 3‰ higher than in consumed plants (Koch, Reference Koch, Michener and Lajtha2007). The δ15N value in plants depends on many factors, including the substrate δ15N, degree of soil development, availability of nutrients, mycorrhizal associations, and soil acidity (Stevens et al., Reference Stevens, Lister and Hedges2006). Globally, δ15N values in ecosystems are climatically controlled (Amundson et al., Reference Amundson, Austin, Schuur, Yoo, Matzek, Kendall, Uebersax, Brenner and Baisden2003) in such a way that generally, lower δ15N values are found in cold and/or moist areas (Fox-Dobbs et al., Reference Fox-Dobbs, Leonard and Koch2008). Warmer, drier, and/or more saline soils raise δ15N 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., Reference Hedges, Stevens and Richards2004). In temperate and semi-arid ecosystems, δ15N values in vegetation vary from 3‰ to 6‰ (Evans and Ehleringer, Reference Evans and Ehleringer1994). In areas where precipitation is greater than 1000 mm/year, δ15N values can range from −2‰ to 0‰ in plants that do not fix nitrogen (Heaton, Reference Heaton1987). δ15N 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., Reference Schwarcz, Dupras and Fairgrieve1999).
Although physiological and body size differences between Pilauco megafauna and P. puda must be considered for any final interpretation, several studies have shown that 15N fractionation is relatively constant in terrestrial mammal herbivores (e.g., Murphy et al., Reference Murphy and Bowman2006; Kuitems et al., Reference Kuitems, van Kolfschoten and van der Plicht2015). 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., Reference Bocherens, Drucker and Madelaine2014). In particular, the turnover rate for bone collagen is relatively slow, thus 13C and 15N values of bone collagen reflect the average isotopic values of dietary protein for several years of the animal's lifespan (Hedges et al., Reference Hedges, Stevens, Koch and Leng2006).
Analysis of microfossils from dental calculus
Plant microremains in dental calculus can provide direct information on feeding habits (Cordova and Avery, Reference Cordova and Avery2017) and a long-term dietary signal (Weyrich et al., Reference Weyrich, Duchene, Soubrier, Arriola, Llamas, Breen and Morris2017). 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., Reference Power, Salazar-García, Wittig, Freiberg and Henry2015). 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., Reference Weyrich, Duchene, Soubrier, Arriola, Llamas, Breen and Morris2017). Considering that older individuals present more microremains (Power et al., Reference Power, Salazar-García, Wittig, Freiberg and Henry2015), 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., Reference Jouy-Avantin, Debenath, Moigne and Moné2003) that can contain a variety of macroscopic and microscopic remains (Martínez and Yagueddú, Reference 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, Reference Reinhard, Bryant and Schiffer1992). Moreover, a variety of macroscopic and microscopic remains allows researchers to obtain interrelated data sets that enable the reconstruction of diets (Piperno, Reference Piperno2006), behavior (Horrocks et al., Reference Horrocks, Irwin, McGlone, Nichol and Williams2003), and paleoenvironments (Jouy-Avantin et al., Reference Jouy-Avantin, Debenath, Moigne and Moné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, Reference Hick1969). All sampled specimens showed typical subadult coat coloration (Hershkovitz, Reference Hershkovitz1982) and permanent teeth had been erupted. Also, because different tissues can show different isotopic values (Carleton et al., Reference Carleton, Kelly, Anderson-Sprecher and del Rio2008), 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., Reference Hedges, Stevens, Koch and Leng2006). To ensure that femur isotopic values were not significantly different from other tissues (Bocherens et al., Reference Bocherens, Cotte, Bonini, Straccia, Scian, Soibelzon and Prevosti2017), dentin from the dental root was measured (Supplementary Table 2). This tissue represents the isotopic composition of the dentin formation period only (Hick, Reference Hick1969). Because studies of bone and dentin have shown that the value of diet-to-tissue trophic discrimination for δ13C is between 3‰ and 5‰ (Drucker et al., Reference Drucker, Bocherens, Bridault and Billiou2003), and between 2‰ and 5‰ for δ15N (Hedges et al., Reference Hedges, Stevens, Koch and Leng2006), 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., Reference Pavez-Fox, Pino and Corti2015). 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., Reference Iacumin, Nikolaev and Ramigni2000; Coltrain et al., Reference Coltrain, Harris, Cerling, Ehleringer, Dearing, Ward and Allen2004) 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 (δ15N and δ13Ccollagen) 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., Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí2017; 2018), new N. platensis samples were analyzed (Table 1). In addition, stable isotope analyses (δ15N and δ13C from root collagen and δ13C from enamel bioapatite) were performed in N. platensis specimens from Los Notros (40°34′S; Lira et al., Reference Lira, Labarca, Fritte, Oyarzo, Pino, Pino and Astorga2020) and Curaco de Vélez (42°26′S; this study; Table 1), which are geographically close and contemporaneous sites.
Table 1. (A) δ13C 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 δ13C (‰, VPDB) and δ15N (‰, 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., Reference Mothé, dos Santos Avilla, Asevedo, Borges-Silva, Rosas, Labarca-Encina and Souberlich2017).
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 N2 air (AIR) scales were done using B2159 sorghum (δ13C VPDB = −13.78 ± 0.17‰ and δ15N AIR = +1.58 ± 0.15‰) calibrated by Elemental Microanalysis to IAEA-CH-6 for δ13C (accepted δ13C VPDB = −10.449 ± 0.033‰) and IAEA-N-1 (accepted δ15NAIR = +0.4 ± 0.2‰). A laboratory soil standard was used as quality control (accepted values δ13C VPDB = −22.69 ± 0.04‰ and δ15N AIR = 7.81 ± 0.07‰), and deviation from accepted values was <0.3‰ and 0.5‰ for δ13C VPDB and δ15N AIR, respectively. Typical standard deviation for replicate samples is δ13C = 0.2‰, and δ15N 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 (Reference Longin1971) and modified by Bocherens et al. (Reference Bocherens, Fizet, Mariotti, Billiou, Bellon, Borel and Simone1991). 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 δ13C and δ15N measures. We used the Ali-j1, Caf-j1, and blk standards, and we used the 3-point linear normalization method (see Werner and Brand, Reference Werner and Brand2001). The VPDB and AIR scales were used as a reference for the δ13C and δ15N 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., (Reference Koch, Tuross and Fogel1997) that were modified in Tornero et al. (Reference Tornero, Bălăşescu, Ughetto-Monfrin, Voinea and Balasse2013). Samples were treated for 4 h in 0.1 M acetic acid [CH3COOH] (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 [H3PO4] at 70°C in individual vessels and purified in an automated cryogenic distillation system. δ13C 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 δ13C; CECC expected values −20.78‰ for δ13C). The mean analytical precision of RC-1 was +0.01‰ for δ13C values and +0.01‰ for CECC.
In this study, a δ13CatmCO2 value of −6.5‰ was used because it is an accepted value for late Pleistocene studies (Tipple et al., Reference Tipple, Meyers and Pagani2010). Therefore, stable isotope data of P. puda and modern vegetation were corrected because the modern composition of 13CatmCO2 has a value of −8‰ (Marino and McElroy, Reference Marino and McElroy1991). Following previous studies (e.g., Koch, Reference Koch, Michener and Lajtha2007; Metcalfe et al., Reference Metcalfe, Longstaffe and Hodgins2013), an average of the diet-to-tissue trophic discrimination was applied at ~3‰ for δ15N and 5‰ for δ13C values in collagen samples. According to diet-to-tissue trophic discrimination studies (e.g., Koch, Reference Koch, Michener and Lajtha2007), 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., Reference Tejada-Lara, MacFadden, Bermudez, Rojas, Salas-Gismondi and Flynn2018), 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 (δ13C = 12‰) value of P. puda (bm: 9.6 kg), it was possible to increase the confidence of the results obtained from the comparisons between δ13Cbioapatite 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., Reference Tejada-Lara, MacFadden, Bermudez, Rojas, Salas-Gismondi and Flynn2018). 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, Reference Cerling and Harris1999; Domingo et al., Reference Domingo, Tomassini, Montalvo, Sanz-Pérez and Alberdi2020). However, according to the estimate of gomphothere body mass (6,000 kg), the study by Asevedo et al. (Reference Asevedo, Ranzi, Kalliola, Pärssinen, Ruokolainen, Cozzuol and Rodrigues do Nascimento2021) 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, Reference DeNiro1985; Ambrose and Norr, Reference Ambrose, Norr, Lambert and Grupe1993; Van Klinken, Reference van Klinken1999), and values within this range in fossil material may be considered as corresponding to unaltered collagen (Clementz et al., Reference Clementz, Fox-Dobbs, Wheatley, Koch and Doak2009). Moreover, it is expected that isotopically well-preserved collagen would exhibit C and N percentages higher than 13% and 4.8%, respectively (Ambrose, Reference Ambrose1990).
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. (Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018). 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. (Reference Wesolowski, de Souza, Reinhard and Ceccantini2007). 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 (Reference Maher1981) as modified by Wesolowski et al. (Reference Wesolowski, de Souza, Reinhard and Ceccantini2010). Coprolite samples were processed according to the methodology described by Katz et al. (Reference Katz, Cabanes, Weiner, Maeir, Boaretto and Shahack-Gross2010). 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), δ13C values ranged from −25.2‰ to −22.2‰, while δ15N values ranged from 5‰ to −1.6‰ (Table 2). Mean δ13C and δ15N values were −23.8‰ ± 1‰ and 3.2‰ ± 2‰, respectively. When applying δ13C diet-to-tissue trophic discrimination of ~5‰, the mean ECP of δ13CECP value was −30.3 ± 0.9‰. When applying δ15N diet-to-tissue trophic discrimination of ~3‰, the mean ECP of δ15NECP was 0.2 ± 2‰. Regarding δ13Cbioapatite 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. δ13C (‰, VPDB) and δ15N (‰, 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 δ13C; 2–5‰ for δ15N). We used the maximum cut-off value (δ13C = −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., Reference Drucker, Bocherens, Bridault and Billiou2003, Reference Drucker, Bridault, Hobson, Szuma and Bocherens2008; Hofman-Kamińska et al., Reference Hofman-Kamińska, Bocherens, Borowik, Drucker and Kowalczyk2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ. All data normalized based on preindustrial atmospheric conditions (δ13CatmCO2 = −6.5‰).
Table 3. δ13C 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 (δ13CatmCO2 = −6.5‰).
When comparing the mean of δ13C and δ15N values between the femur and the dentin of P. puda, a similarity is observed (δ13Cfemur= -23.8‰; δ13Cdentine= -24‰; δ15Nfemur = 3.2‰; δ15Ndentine = 3.6‰; Tables 3, 4). However, looking at δ15N 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. δ13C (‰, VPDB) and δ15N (‰, 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 δ13C; 2‰–5‰ for δ15N). We used the maximum cut-off value (δ13C = −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., Reference Drucker, Bocherens, Bridault and Billiou2003, Reference Drucker, Bridault, Hobson, Szuma and Bocherens2008; Hofman-Kamińska et al., Reference Hofman-Kamińska, Bocherens, Borowik, Drucker and Kowalczyk2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ. All data normalized based on preindustrial atmospheric conditions (δ13CatmCO2 = −6.5‰).
The mean δ13C value of estimated consumed plants (ECP) (+12‰) was −28.4‰. For modern vegetation, the ranges of δ13C 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 δ15N 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 δ13C and δ15N values for each area. The difference among plant δ13C 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 δ15N 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 δ13C 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 δ13C (‰, VPDB) and δ15N (‰, AIR) values in Pilauco Megafauna, modern vegetation, Pudu puda, gomphotheres from the surroundings of Pilauco (González-Guarda et al. Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018), megafauna from the Southern Patagonia region of Chile (Prevosti and Martin, Reference Prevosti and Martin2013; equids, xenarthrans, camelids, ratites) and megafauna from the Pampean region (Bocherens et al., Reference Bocherens, Cotte, Bonini, Scian, Straccia, Soibelzon and Prevosti2016; xenarthrans, equids, cervids, macrauchenids, toxodons, rodents). Vegetation δ13C (‰, VPDB), δ15N (‰, 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) δ13C (‰, VPDB), (B) δ15N (‰, AIR) values. Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ. All data normalized based on preindustrial atmospheric conditions (δ13CatmCO2 = −6.5‰).
Values of stable isotopes in fossil megafauna
The δ13C 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. δ15N 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 (δ13C= −21.3‰; δ15N = 5.7‰; Table 6). Figure 4 shows mean δ13C and δ15N values for each taxon from Pilauco. Results fit well with previous analyses in Pilauco using similar procedures (Aguilera, Reference Aguilera2010, n = 3; González-Guarda et al., Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí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 δ13C values of N. platensis, E. andium, cf. H. paradoxa and Xenarthra indet. (p < 0.001). For δ15N values, significant differences are documented in all taxa (p = 0.009).
Figure 4. Average bone collagen values for δ13C (‰, VPDB), indicating a more closed to open environment, and δ15N (‰, 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 δ13C (‰, VPDB) and δ15N (‰, 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 δ13C; 2‰–5‰ for δ15N). We used the maximum cutoff value (δ13C = −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., Reference Drucker, Bocherens, Bridault and Billiou2003, Reference Drucker, Bridault, Hobson, Szuma and Bocherens2008; Hofman-Kamińska et al., Reference Hofman-Kamińska, Bocherens, Borowik, Drucker and Kowalczyk2018). 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 δ13C (‰, VPDB) values, (B) Collagen δ15N (‰, AIR) values. Stable isotope values of this study were analyzed together with values from Aguilera (Reference Aguilera2010) and González-Guarda et al. (Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí2017). The δ13C (‰, VPDB) and δ15N (‰, 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 δ13C; 2‰–5‰ for δ15N). We used the maximum cutoff value (δ13C = −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., Reference Drucker, Bocherens, Bridault and Billiou2003, Reference Drucker, Bridault, Hobson, Szuma and Bocherens2008; Hofman-Kamińska et al., Reference Hofman-Kamińska, Bocherens, Borowik, Drucker and Kowalczyk2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ.
Once the δ13C diet-to-tissue trophic discrimination of 5‰ was applied, the average value of the estimated consumed plants (ECP = δ13CECP) 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 δ15N diet-to-tissue trophic discrimination of 3‰, mean ECP δ15NECP 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 δ13Cbioapatite values in samples of N. platensis from Pilauco, Curaco de Vélez, and Los Notros, the mean value was −13.8 ± 0.5‰. The δ13Ccollagen and δ15Ncollagen values from Los Notros were −22.2‰ and 6.6‰, respectively, and the δ13Ccollagen and δ15Ncollagen 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 C3-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, Reference Korstanje and O'Connor2002). 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 (C3 plant), Bambusoideae (C3 plant), and Panicoideae (C4 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 C4 plants that are adapted to temperate environments, can tolerate high tree cover (e.g., Imperata condensate; Zuloaga et al., Reference Zuloaga, Belgrano and Zanotti2019) and wide latitudinal and altitudinal distribution ranges (e.g., Paspalum sp.; Rodríguez et al., Reference Rodríguez, Marticorena, Alarcón, Baeza, Cavieres, Finot and Fuentes2018). Diatoms were present in all samples and in variable proportions that were higher in sediment than coprolites. Spherulites found in samples are CaCO3 crystals that form during the digestive process (Canti, Reference Canti1998) in the intestines of certain animals (Brochier et al., Reference Brochier, Villa, Giacomarra and Tagliacozzo1992).
DISCUSSION
Feeding behavior and habitat of extinct species
Assuming a cutoff value of −27.5‰ for δ13Ccollagen values as delineating a plant-diet based on open areas versus forested areas (or δ13Cbioapatite (diet) = −27.2‰; Tejada et al., Reference Tejada, Flynn, Antoine, Pacheco, Salas-Gismondi and Cerling2020), and considering the average values of δ13C 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., Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018). Particularly interesting are the δ13Cbioapatite values (MHMOPI/628a = −14‰ and MHMOPI/627a = −13‰) in the Pilauco gomphotheres, which indicate more closed conditions. However, although slightly different from δ13Ccollagen values, δ13Cbioapatite 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 (δ13C values for woodland-mesic C3 grassland, −14.5‰ to −9.5‰; Domingo et al., Reference Domingo, Prado and Alberdi2012).
Observed differences in the δ15N 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 (δ15Ndentine = 8.5‰; González-Guarda et al., Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí2017; δ15Nbone = 9.2‰; this study), some E. andium samples with δ15Nbone values (8.7‰) differ from modern vegetation. These values could indicate foraging in grassland areas (Bocherens, Reference Bocherens2003), but as stated above, most of the δ13C values obtained suggest feeding in more wooded areas. What factors could have caused high δ15N 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 δ15Ndentine values (from 14.2‰ to 1.3‰; González-Guarda et al., Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018), 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., Reference Metcalfe, Longstaffe and Hodgins2013). 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; δ15Ndentine = 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; δ15Ndentine = 2.3 ‰), 200 km from Pilauco (for a more comprehensive chronological context for the study area, see Supplementary Table 8).
Again, the δ15Ndentine value (Supplementary Table 5) of a new gomphothere (Los Notros) from the Pilauco area does not overlap with δ15N values of modern vegetation (Supplementary Table 4), while the gomphothere sampled in the Curaco de Vélez has a low δ15N value similar to those observed in gomphotheres from the northwestern Patagonia Chilean (Table 8). Therefore, regardless of the specific factors that could affect δ15N 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., Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018 The δ13C (‰, VPDB) and δ15N (‰, 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 δ13C; 2‰–5‰ for δ15N). We used the maximum cutoff value (δ13C = −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., Reference Drucker, Bocherens, Bridault and Billiou2003, Reference Drucker, Bridault, Hobson, Szuma and Bocherens2008; Hofman-Kamińska et al., Reference Hofman-Kamińska, Bocherens, Borowik, Drucker and Kowalczyk2018). Number of samples: n; maximum: Max; minimum: Min; mean values: Mean; standard deviation: σ.
The absence of overlap of δ15N values between Pilauco megafauna and modern vegetation could indicate a diet based on plants found in grassland areas. However, δ13C values indicate a diet from a more wooded environment. When comparing the δ13C 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 δ13Ccollagen 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 δ13Cbioapatite 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, Reference Moreno2020). 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 δ13Cbioapatite 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., Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018).
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, Reference Elton2008). 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, Reference Veit1994). 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., Reference Rivals, Semprebon and Lister2019). 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., Reference Asevedo, Winck, Mothé and Avilla2012). Likewise, H. paradoxa is characterized by grazing behavior, inhabiting more arid habitats (Menégaz and Ortiz Jaureguizar, Reference Menegaz, Ortiz-Jaureguizar, Artiodactilos, Alberdi, Leone y and Tonni1995). However, studies of Hemiauchenia from other areas of America have shown a browser and mixed-feeder diet (Feranec, Reference Feranec2003; Semprebon and Rivals, Reference Semprebon and Rivals2010), which is also demonstrated in specimens sampled from Pilauco. In addition, Equus andium occupied different environments in South America; from mixed C3–C4 grassland (Ecuador; Domingo et al., Reference Domingo, Prado and Alberdi2012) to woodland mesic C3 grassland and closed-canopy forest (Chile; González-Guarda et al., Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí2017). The Xenarthra teeth are morphologically simple and do not provide unambiguous information about their diet (Bargo and Vizcaíno, Reference Bargo and Vizcaíno2008). However, skull and jaw morphology differences among these taxa yielded evidence for possible feeding mechanisms and therefore food composition (Vizcaíno et al., Reference Vizcaíno, Cassini, Fernicola and Bargo2011). Recently, an isotopic study on bulk collagen showed that xenarthrans were exclusively herbivorous (Bocherens et al., Reference Bocherens, Cotte, Bonini, Straccia, Scian, Soibelzon and Prevosti2017). Nevertheless, previous studies have indicated that this group of animals may have had insectivorous (Genise and Fariña, Reference Genise and Fariña2012) and carnivorous representatives (Fariña and Blanco, Reference Fariña and Blanco1996).
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., Reference González, Labarca, Chavez-Hoffmeister and Pino2014), 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 δ13CECPcollagen value for the dataset of the Pilauco megafauna (−27.5 ± 0.8‰) matches the maximum δ13C 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., Reference Drucker, Bocherens, Bridault and Billiou2003, Reference Drucker, Bridault, Hobson, Szuma and Bocherens2008; Hofman-Kamińska et al., Reference Hofman-Kamińska, Bocherens, Borowik, Drucker and Kowalczyk2018). This fits well with δ13C values measured in several tooth enamel bioapatite samples from megafauna (mean δ13Cbioapatite value = −13.5‰ in ECP) from northwestern Chilean Patagonia (González-Guarda et al., Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí2017).
Comparison between δ13CECP values from Pilauco megafauna and δ13C 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 δ13C 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 δ13C values between both types of plants have been reported elsewhere (e.g., Fox-Dobbs et al., Reference Fox-Dobbs, Leonard and Koch2008).
ECP nitrogen values in Pilauco megafauna are not in agreement with mean δ15N 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 δ15N = 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 (δ15N from 3‰ to 6 ‰; Evans and Ehleringer, Reference Evans and Ehleringer1994), 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, δ15N data from modern vegetation and fauna compels us to cautiously apply environmental ranges established based on δ15N.
Comparison between δ13CECPcollagen values from Pilauco megafauna (mean = −27.5 ± 0.8‰) and P. puda (mean = δ13CECPcollagen = −28.8 ± 1‰) indicate very similar results. This similarity is also found in P. puda δ13CECPbioapatite 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 δ13Cmodern vegetation value = −33 ± 2.6‰).
Interestingly, δ15NECP values of P. puda (mean = 0.2 ± 1.8‰) and sampled gomphotheres from Pilauco (δ15NECP = 3.7 ± 0.9‰) show important differences, but δ15NECP 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., Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018). Mean δ13CECPcollagen 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 δ15N 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 δ13C from the Pilauco megafauna derived from collagen and bioapatite indicate that the Pilauco ecosystem was forested, while δ15N 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 δ15N 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 δ15N 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., Reference Tejada-Lara, MacFadden, Bermudez, Rojas, Salas-Gismondi and Flynn2018). As this study and that of González-Guarda et al. (Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018), 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 δ15N 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 δ15N 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 δ15N 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 (δ15N = −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 δ15N values could increase the uncertainty of paleodietary reconstructions based solely on δ15Ncollagen. However, if a particular study is supported by a robust isotopic baseline on current samples, δ15Ncollagen 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.
Acknowledgments
Isotopic samples were prepared at the Biomolecular laboratory of IPHES and measured at the Institute of Environmental Science and Technology (ICTA) with technical assistance by Dr. Pau Comes and scientific supervision from Dr. A. Rosell. Dr. C. Tornero acknowledges the Beatriu de Pinós Post-doctoral fellowship (BP-MSCA Cofound code 2016-00346 from the AGAUR, Government of Catalonia, Spain). EG-G thanks the Postdoctoral Fondecyt 3200806. IRP is the beneficiary of a predoctoral fellowship (2021 FI_B1 00223) funded by AGAUR and the Fons Social Europeu (FSE). Special thanks to Pedro Aburto (Centro de Rehabilitación de Fauna Silvestre de la Universidad Austral de Chile), Karen Moreno (Laboratorio de Paleontología, UACh), and Servicio Agrícola Ganadero de Chile. This research was supported by the Spanish Ministry of Science and Innovation through the “María de Maeztu” excellence accreditation (CEX2019-000945-M) and “CERCA Programme/Generalitat de Catalunya”.
Supplementary Material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2022.6
