Introduction
During the Pliocene, the rise of the Isthmus of Panama connected North and South American lands, triggering the major episode of trans-isthmian mammalian exchanges, an event called the great American biotic interchange (Webb Reference Webb2006; Woodburne Reference Woodburne2010). At that time, the camelids (Mammalia: Camelidae) entered the southern continent and, by the Pleistocene, this group was widely distributed and diversified (Scherer Reference Scherer2013). South American camelids comprise the subfamily Camelinae, tribe Lamini, which diverged from the tribe Camelini during the Miocene, about 16 Ma, in North America (Lynch et al. Reference Lynch, Sánchez-Villagra and Balcarcel2020). South American camelid fossils have been recovered from fossiliferous deposits from Venezuela to southern Chile (see Rocha-dos-Santos et al. Reference Rocha-dos-Santos, Avilla and Scherer2017). According to the fossil record of southern Brazil, three species of camelid Lamini inhabited this region in the recent past: Hemiauchenia paradoxa Gervais and Ameghino Reference Gervais and Ameghino1880, Lama guanicoe Müller, Reference Müller1776, and Vicugna vicugna Molina, Reference Molina1782 (Scherer Reference Scherer2009, Reference Scherer2013; Rocha-dos-Santos et al. Reference Rocha-dos-Santos, Avilla and Scherer2017). While H. paradoxa is currently extinct, L. guanicoe and V. vicugna still live in some regions of South America—from the northern slope of the Andes to the island of Tierra del Fuego in the south (Baldi et al. Reference Baldi, Acebes, Cuéllar, Funes, Hoces, Puig and Franklin2016; Acebes et al. Reference Acebes, Wheeler, Baldo, Tuppia, Lichtenstein, Hoces and Franklin2018; see also Fig. 1), and are popularly known as guanaco and vicuña, respectively. Given this scenario, Lopes et al. (Reference Lopes, Pereira, Kerber and Dillenburg2020) classified the two last taxa as pseudo-extinct. Following the same authors, the climatic changes of the Pleistocene/Holocene transition guided such extinctions in the southern Brazilian Pampa.
Map showing the current distribution of the two endemic wild camelids of South America (Lama guanicoe and Vicugna vicugna; dashed area) (Baldi et al. Reference Baldi, Acebes, Cuéllar, Funes, Hoces, Puig and Franklin2016; Acebes et al. Reference Acebes, Wheeler, Baldo, Tuppia, Lichtenstein, Hoces and Franklin2018) and the fossil collection sites (L. guanicoe, V. vicugna, and Hemiauchenia paradoxa) in the Rio Grande do Sul, southern Brazil.

The Pleistocene climate was marked by several glaciation events that favored the development of grasslands in southern Brazil due to the colder and drier conditions (Behling Reference Behling2002). However, from ca. 18 ka onward, the last deglaciation process began (Cruz et al. Reference Cruz, Burns, Karmann, Sharp, Vuille and Ferrari2006; Denton et al. Reference Denton, Anderson, Toggweiler, Edwards, Schaefer and Putnam2010), increasing temperatures and humidity since mid-Holocene times (Behling et al. Reference Behling, Pillar, Orlóci and Bauermann2004, Reference Behling, Pillar and Bauermann2005; Medeanic and Corrêa Reference Medeanic and Corrêa2010). Climate amelioration led to changes in vegetation composition: a progressive decrease of grassland pollen types, following an increase in forest taxa (Cordeiro and Lorscheitter Reference Cordeiro and Lorscheitter1994; Lorscheitter and Dillenburg Reference Lorscheitter and Dillenburg1998; Behling et al. Reference Behling, Pillar, Orlóci and Bauermann2004, Reference Behling, Pillar and Bauermann2005; Macedo et al. Reference Macedo, Cancelli, Bauermann, Bordignon and das Neves2007; Masetto and Lorscheitter Reference Masetto and Lorscheitter2019; Roth et al. Reference Roth, Lorscheitter and Masetto2021). With the loss of grassland habitats, the extant Lamini became restricted to the Andean region, where cold and temperate weather is still predominant (González et al. Reference González, Palma, Zapata and Marín2006; Scherer Reference Scherer2013). Guanacos and vicuñas show ecophysiological adaptations to harsh environments that allowed them to survive and colonize the remaining open habitats of South America after the megafauna extinction of the Pleistocene/Holocene transition (González et al. Reference González, Palma, Zapata and Marín2006).
Pleistocene colder and drier conditions likely supported more heterogeneous vegetation, which allowed the co-occurrence of several herbivorous including closely related taxa (e.g., Yann Reference Yann2014). Forage availability and quality are related to population dynamics, spatial distribution, and feeding habits of South American wild camelids (González et al. Reference González, Palma, Zapata and Marín2006; Shaw et al. Reference Shaw, Galaz and Marquet2012; Rivals et al. Reference Rivals, Rindel and Belardi2013). Thus, assessing the paleoenvironment and the paleodiet of Pleistocene camelids can give us insights into the factors involved in the extinctions and pseudo-extinctions in southern Brazil. In addition, determining the diet of extinct animals contributes to the understanding of the function and structure of the south Brazilian paleocommunity, such as resource partitioning and competition between co-occurrent species (e.g., Feranec Reference Feranec2003).
The ecology of extinct species is typically analyzed through their dental morphologies and tooth wear (i.e., microwear and mesowear) (e.g., Hayek et al. Reference Hayek, Bernor, Solounias and Steigerwald1992; Kaiser et al. Reference Kaiser, Solounias, Fortelius, Bernor and Schrenk2000; Rivals and Semprebon Reference Rivals and Semprebon2011; Saarinen and Lister Reference Saarinen and Lister2016). However, the specialized morphology observed in some herbivores does not always indicate a specialized diet (Bocherens and Drucker Reference Bocherens, Drucker, Elias and Mock2013). An alternative method for studying such species is stable isotope analysis (see Clementz Reference Clementz2012), because it allows for more accurate reconstructions of the actual diet of fossil mammals and, consequently, of the environments they inhabited (Bocherens and Drucker Reference Bocherens, Drucker, Elias and Mock2013). This analysis is based on the premise that the organism's isotopic composition changes in predictable ways, following the biogeochemical processes of the trophic web (Peterson and Fry Reference Peterson and Fry1987). In paleontological studies, the most common isotopes are carbon and oxygen, which can be measured in the inorganic fraction (biological apatite) of skeletal remains (e.g., Yann et al. Reference Yann, DeSantis, Koch and Lundelius2016; Ma et al. Reference Ma, Wang, Jin, Hu and Bocherens2019; Domingo et al. Reference Domingo, Tomassini, Montalvo, Sanz-Pérez and Alberdi2020; Asevedo et al. Reference Asevedo, Ranzi, Kalliola, Pärssinen, Ruokolainen, Cozzuol, do Nascimento, Negri, Souza-Filho, Cherkinsky and Dantas2021). The isotopic ratios of oxygen (δ18O) in inert tissues are mainly related to water ingested by drinking and eating (Luz et al. Reference Luz, Kolodny and Horowitz1984; Kohn Reference Kohn1996). Thus, the oxygen isotope composition of mammalian herbivores is strongly correlated with local meteoric water compositions and average humidity (Kohn Reference Kohn1996). The isotopic ratios of carbon (δ13C), in turn, are derived from the diet, which implies that the δ13C values of herbivores reflect the isotopic ratios of the consumed plants (plus predictable changes—called isotopic enrichment), being considerably different between producers with physiological differences (e.g., plants with C3 and C4 photosynthetic pathways) (Kelly Reference Kelly2000).
Therefore, stable isotope analysis can provide proxies for habitat and diet information, regardless of morphological characters (Clementz Reference Clementz2012). Using this valuable tool, we assessed the diet of South American camelids from southern Brazil with insights into the paleoenvironment of this region.
Materials and Methods
Samples
Teeth and bones of Lamini camelids were obtained from the paleovertebrate collection of the Museu de Ciências Naturais (Secretaria do Meio Ambiente e Infraestrutura–SEMA) (MCN–PV), located in Porto Alegre, RS, Brazil.
Because enamel is less prone to diagenetic alteration (Hedges et al. Reference Hedges, Stevens and Richards2004; Clementz Reference Clementz2012; Bocherens and Drucker Reference Bocherens, Drucker, Elias and Mock2013), sampling focused on this tissue. However, when enamel was not available or to avoid damage to fossil specimens, bone or dentin was sampled (see Table 1).
Sample information (collection site and sampled tissue) and stable isotope results of Lamini camelids: Hemiauchenia paradoxa, Lama guanicoe, and Vicugna vicugna. VPDB, Vienna Pee Dee Belemnite; VSMOW, Vienna Standard Mean Ocean Water.

Setting and Geological Context
The specimens analyzed in the present study are from two geological localities, both situated in southern Brazil: Touro Passo and Santa Vitória Formations (Table 1, Fig. 1).
Touro Passo Formation is located in the western region of Rio Grande do Sul, near the borders of Uruguay and Argentina (Fig. 1). It comprises sedimentary deposits generated on the Serra Geral Plateau of the Paraná Basin (Da-Rosa Reference Da-Rosa, Ribeiro, Bauermann and Scherer2009). Fossils were recovered from the homonymous creek (i.e., Touro Passo Creek, 29°42′08″S, 56°45′16″W), a tributary of the Uruguay River, which exposed late Pleistocene lithologies (conglomeratic facies overlapped by silt and sand levels) along its margins (Kerber et al. Reference Kerber, Pitana, Ribeiro, Hsiou and Oliveira2014). Thermoluminescence dating of sediments indicates a deposition age of 42 up to 6.4 ka (Da-Rosa Reference Da-Rosa, Ribeiro, Bauermann and Scherer2009). Most vertebrate remains accumulated during the same period. According to electron spin resonance dating, fossils from Ponte Velha I locality, Touro Passo Creek, have ages between 42.6 and 10 ka and are therefore from the last glacial maximum (Kerber et al. Reference Kerber, Kinoshita, José, Figueiredo, Oliveira and Baffa2011).
Fossils were also recovered from Hermenegildo Beach (33°40′00″S, 53°15′35″W), typically found ex situ, after erosion of the submerged fossiliferous deposits by storm waves (Buchmann Reference Buchmann2002). Thus, they lack any precise stratigraphic context. However, it is supposed that specimens recovered on the beach are from the same stratigraphic provenance as those found in situ in a near creek (Chuí Creek), Santa Vitória Formation (Soliani Reference Soliani1973), as they have similar faunas. Santa Vitória Formation is located in the southern coastal plain, on the emerged portion of the Pelotas Sedimentary Basin (Fig. 1). Its formation is related to fluvial and eolian deposits from the Middle–Upper Pleistocene (Lopes et al. Reference Lopes, Dillenburg, Savian and Pereira2021). Pleistocene sediments are mainly composed of quartz sediments as well as silty-clay sediments (Buchmann et al. Reference Buchmann, Caron, Lopes, Ugri, de Lima, Ribeiro, Bauermann and Scherer2009). As a result of rework and deposition guided by sea-level fluctuations, vertebrate fossils span a large temporal window. Submerged terrestrial remains are associated with glacial periods, when the continental shelf was exposed, and have ages between 700 and 18 ka, as revealed by electron spin resonance dating (Lopes et al. Reference Lopes, Oliveira, Figueiredo, Kinoshita, Baffa and Buchmann2010).
The current climate of the study areas is warm temperate (Cfb) and subtropical (Cfa), according to Köppen climate zones, without a pronounced dry season. The Touro Passo Formation region's average summer and winter temperatures are 24°C and 14°C, respectively, while mean summer temperatures are below 24°C and winter temperatures vary around 12°C along the coast (Wrege et al. Reference Wrege, Steinmetz, Júnior and Almeida2012). During the Pleistocene glacial periods, however, the climate was drier and temperatures were 5°C–7°C colder than today (Behling Reference Behling2002), probably as a result of more frequent, intense, and persistent Antarctic cold fronts over the Southern Hemisphere (Behling Reference Behling2002; Cruz et al. Reference Cruz, Burns, Karmann, Sharp, Vuille and Ferrari2006).
Both Touro Passo and Santa Vitória Formations are currently located in the Pampa biome, with the landscape dominated by grasslands similar to those found in Uruguay and Argentina (Behling Reference Behling2002; Overbeck et al. Reference Overbeck, Boldrini, do Carmo, Garcia, Moro, Pinto, Trevisan, Zannin, Pillar and Lange2015). Despite the greater contribution of the grasslands, this region is quite heterogeneous concerning its plant physiognomies, with the presence of forest formations, although they do not dominate the landscape (Overbeck et al. Reference Overbeck, Boldrini, do Carmo, Garcia, Moro, Pinto, Trevisan, Zannin, Pillar and Lange2015). Conforming to the climatic characteristics of Rio Grande do Sul (e.g., temperature and rainfall), models predict that the natural vegetation of this region should be composed of forests, rather than grasslands (Buriol et al. Reference Buriol, Estefanel, de Chagas and Eberhardt2007). However, due to the Quaternary's history of changes, vegetation is likely still developing toward forest environments. Both forest and grassland environments of southern Brazil are strongly influenced by human activity for agriculture and livestock and, consequently, only small areas of natural vegetation remain (Behling Reference Behling2002).
Several botanical families are found in the Pampa biome, such as Poaceae, Asteraceae, Fabaceae, Cyperaceae, and Verbenaceae, among others (Boldrini et al. Reference Boldrini, Overbeck, Trevisan, Pillar and Lange2015). The proportion of C3 and C4 grasses varies according to the climate. Toward the south, where the climate becomes temperate, there is an increase in the number of C3 grasses, whereas from east to west, there is a greater proportion of C4 grasses, due to higher temperatures (Overbeck et al. Reference Overbeck, Boldrini, do Carmo, Garcia, Moro, Pinto, Trevisan, Zannin, Pillar and Lange2015). Grasslands were also the main biome in southern Brazil during the last glacial maximum, including in regions where different forest biomes exist today (Behling Reference Behling2002).
Stable Isotope Analysis
Fossil bone, dentin and enamel were sampled using a small drill and then pretreated according to Bocherens et al. (Reference Bocherens, Koch, Mariotti, Geraads and Jaeger1996), with some modifications.
First, to remove any residual organic matter, the powder was soaked in a 2.5% NaOCl solution for 24 h at room temperature and rinsed with distilled water. To remove the exogenous carbonate, a 0.1 M acetic acid solution was added to the samples. After 24 h at room temperature, the samples were rinsed again with distilled water and dried in an oven at 60°C for 48 h. Approximately 20 mg of the resulting powder from each sample of bone, dentin, and enamel was packed into small tubes for further isotopic analysis at the University of Georgia, Center for Applied Isotope Studies (UGA-CAIS, Athens, GA, USA). Carbon and oxygen isotope compositions were measured through a Thermo Gas Bench coupled to a Delta V Isotope Ratio Mass Spectrometer with a GC Pal autosampler (GasBench–IRMS). For this, the solid carbonate was converted to CO2 gas by dissolution with phosphoric acid in a helium atmosphere.
The resulting isotopic ratio (R) of each element (13C/12C and 18O/16O) in the samples was reported in delta notation (δ) of the element in its heavy form (X) in parts per thousand (‰) relative to the international standard Vienna Pee Dee Belemnite (VPDB), as expressed by the following equation: δX = [(RSAMPLE / RSTANDARD) – 1]. Internal laboratory control standards used for calibration were fish (δ13C = −0.64‰; δ18O = −14.90‰) and A1296 (δ13C = 2.56‰; δ18O = −0.60‰). The analytical uncertainty (σ) of the 12 repeated measurements of the standards was <0.02‰ for δ13C and <0.10 for δ18O.
Data Analysis and Interpretation
The diet of the Lamini camelids that inhabited the southern Brazilian Pampa during the Pleistocene was reconstructed based on carbon isotope ratios. Such assessment is possible due to the carbon isotope fractionation during the photosynthetic process. For example, C3 plants fix CO2 through the action of the enzyme ribulose 1,5-biphosphate carboxylase, which discriminates against the heavy isotope, 13C; whereas, in the C4 photosynthetic pathway, the CO2 is initially carboxylated in organic acids by the enzyme phosphoenolpyruvate carboxylase, with fractionation poorly expressed (O'Leary Reference O'Leary1988). Because the biomass productivity of crassulacean acid metabolism plants is low at the community level and they are mostly present in stressful environments like deserts (Lüttge Reference Lüttge2004), it was assumed that they did not represent an important food source, and therefore this specific photosynthetic pathway was not addressed here.
Because carbonate bioapatite is derived from blood bicarbonate (Sullivan and Krueger Reference Sullivan and Krueger1981; Gannes et al. Reference Gannes, de Rio and Koch1998), there is no significant fractionation between bone, dentin, and enamel (Tejada-Lara et al. Reference Tejada-Lara, MacFadden, Bermudez, Rojas, Salas-Gismondi and Flynn2018). Thus, correction factors were not applied to account for differences between tissue types.
The isotopic enrichment (ε*) between carbonate bioapatite of ungulate mammals and their diet is about 14‰ (Cerling and Harris Reference Cerling and Harris1999); however, enrichment is not fixed, because it depends on the digestive physiology of herbivorous species (Tejada-Lara et al. Reference Tejada-Lara, MacFadden, Bermudez, Rojas, Salas-Gismondi and Flynn2018). For a more accurate ε* value, the formula proposed by Tejada-Lara et al. (Reference Tejada-Lara, MacFadden, Bermudez, Rojas, Salas-Gismondi and Flynn2018) for foregut fermenters was applied: ε* = 2.34 + 0.05 (BM), where BM is the natural logarithm of body mass in kilograms. For extant species, the body mass of living populations was adopted. Weight estimations for Lama guanicoe range from 80 to 120 kg (González et al. Reference González, Palma, Zapata and Marín2006 and references therein), therefore a conservative body mass value of 100 was considered here, whereas for a smaller camelid, Vicugna vicugna (weighing between 20 and 50 kg; Ccora et al. Reference Ccora, Condori, Contreras, Curasma, Cordero, Valencia, Mayhua and McGregor2019), a value of 40 was considered. The estimated weight reported in the literature for the extinct Hemiauchenia paradoxa is 1000 kg (Prevosti and Vizcaíno Reference Prevosti and Vizcaíno2006). Thus, after calculation of the natural logarithm, ε* values were estimated to be 13.1‰, 12.5‰, and 14.7‰ for L. guanicoe, V. vicugna, and H. paradoxa, respectively.
As a consequence of differences in carbon isotope fractionation that occurs during photosynthesis, the mean δ13C values of C3 plants is −27.1 ± 2.0‰ and that of C4 plants is −13.1 ± 1.2‰ (O'Leary Reference O'Leary1988). However, these isotopic ratios were determined for modern plants, and therefore carbon isotope ratios must be corrected for the Suess effect (i.e., atmospheric 14C decreases due to the burning of fossil fuels after the Industrial Revolution; Baxter and Walton Reference Baxter and Walton1970). For this, a correction value of 1.3‰ was applied to the values of modern plants, following the calculation proposed by Dombrosky (Reference Dombrosky2020).
The resulting δ13C values of C3 and C4 plants (−25.8 ± 2.0‰ and −11.8 ± 1.2‰, respectively), the estimated ε*, and δ13C values of the consumers were then applied to stable isotope mixing models in R (simmr package; Parnell and Inger Reference Parnell and Inger2016) to estimate the contribution of C3 and C4 plants to the diet of the Lamini camelids. By definition, simmr is a Bayesian mixing model incorporating a Markov chain Monte Carlo method that uses the isotopic data to calculate the probability of contributions from potential sources to a mixture.
Bioapatite δ18O values of herbivores with different behaviors and physiologies track distinct aspects of the environment. Herbivores that supply their water requirements by forage reflect changes in aridity, whereas those that are obligated drinkers trace local meteoric water (Levin et al. Reference Levin, Cerling, Passey, Harris and Ehleringer2006).
Different molecular species (isotopologues) of H2O have distinct vapor pressures that are proportional to their masses (Ehleringer and Rundel Reference Ehleringer, Rundel, Rundel, Ehleringer and Nagy1989). Consequently, water molecules composed of light isotopes are more prone to evaporation, so more positive δ18O in meteoric waters indicate warmer or drier climates, whereas more negative δ18O can represent colder or wetter climates (Bocherens and Drucker Reference Bocherens, Drucker, Elias and Mock2013). Conversely, when herbivores fulfill their water requirements through food, the interpretation of δ18O values must take into account the fractionations that occur in plant organs and tissues. Although there is no fractionation during the passage of water through root membranes, leaf transpiration can result in enrichment in 18O isotopes according to environmental conditions (e.g., atmospheric humidity and air temperature) (Dawson et al. Reference Dawson, Mandelli, Plamboeck, Templer and Tu2002; Barbour Reference Barbour2007; Simonin et al. Reference Simonin, Roddy, Link, Apodaca, Tu, Hu, Dawson and Barbour2013). Hence, differences in δ18O values indicate climate variability within a local population.
All oxygen isotope ratios were converted from VPDB to Vienna Standard Mean Ocean Water (VSMOW) according to the equation: δVSMOW = 1.03086 * δVPDB + 30.86 (Friedman and O'Neil Reference Friedman and O'Neil1977) before analysis.
Differences in carbon and oxygen isotope ratios between H. paradoxa, L. guanicoe, and V. vicugna were assessed through parametric analysis of variance (ANOVA), as the data followed a normal distribution according to a Shapiro-Wilk test. ANOVA was also applied to compare stable isotopic ratios of specimens from distinct collection sites. Unfortunately, this last analysis was only suitable for H. paradoxa, as it had a minimum of three observations per site, unlike L. guanicoe and V. vicugna (see Table 1). Finally, to verify whether the diet of South American camelids changes according to the climate, a linear regression was applied using δ13C and δ18O values.
All analyses were performed in the R statistical environment v. 4.0.3 (R Core Team 2020). The significance level for two-sided tests was α = 0.05.
Results
A total of 17 specimens were analyzed in the present study: 7 Hemiauchenia paradoxa, 6 Lama guanicoe, and 4 Vicugna vicugna (Table 1). δ13CVPDB values ranged from −11.5‰ to−4.7‰, with mean and standard deviation (±SD) of −8.5 ± 1.9‰, while δ18OVSMOW values varied from 25.0‰ to 31.8‰ (mean 29.8 ± 1.8‰). Although the variation in δ18O between individuals was high, the mean values were quite similar among the species analyzed (H. paradoxa, δ18O = 29.6 ± 2.4‰; L. guanicoe, δ18O = 29.9 ± 2.5‰; V. vicugna, δ18O = 30.2 ± 1.0‰), whereas the mean carbon isotope ratios varied broadly. Vicugna vicugna showed higher mean values of δ13C (δ13C = −6.8 ± 0.6‰), compared with H. paradoxa (δ13C = −8.7 ± 1.2‰) and L. guanicoe (δ13C = −9.4 ± 2.5‰). Differences among species were nonsignificant in both δ13C (p-value = 0.10) and δ18O (p = 0.86) values. Isotopic ratios of H. paradoxa also indicated no differences between the collection sites (i.e., Touro Passo and Santa Vitória Formations) in δ13C (p = 0.28) and δ18O values (p = 0.28).
Considering the distinct isotopic enrichment (ε*) of camelids, the δ13C values of each species instead of those of C3 and C4 plants were corrected for graphical representation and are shown in Figure 2. According to the scatter plot, specimens of H. paradoxa and L. guanicoe (except for one individual) are isotopically similar to C3 plants. Specimens of V. vicugna, in contrast, are lying on the boundary between C3 and C4 plants, which points to a mixed diet. Mixing models indicated the same pattern. H. paradoxa from southern Brazil had a diet comprising mostly C3 vegetation (>80% of the diet, credibility interval [CI]: 70%–90%), with a minimal contribution of C4 vegetation (20%, CI: 8%–30%) (Fig. 3A). Similarly, the proportional contribution of C3 plants to the diet of L. guanicoe was higher (~75%, CI: 55%–90%) than that of C4 plants (24%, CI: 8%–45%) (Fig. 3B). Conversely, V. vicugna had a mixed diet; C3 and C4 vegetation comprise about 54% (CI: 40%–70%) and 46% (CI: 30%–60%) of its diet, respectively (Fig. 3C). No significant relationship was found between climate (represented by δ18O values) and diet (δ13C values) of the Lamini camelids (H. paradoxa, r 2 = 0.02, p = 0.78; L. guanicoe, r 2 = 0.05, p = 0.68; V. vicugna, r 2 = 0.004, p = 0.94) (see also Fig. 2).
Scatter plot of δ13C and δ18O values for the specimens of Hemiauchenia paradoxa, Lama guanicoe and Vicugna vicugna from the Pleistocene of southern Brazil. Note that the carbon isotope data are corrected for trophic enrichment. VPDB, Vienna Pee Dee Belemnite; VSMOW, Vienna Standard Mean Ocean Water.

Dietary contribution of C3 and C4 plants to the diet of Hemiauchenia paradoxa (A), Lama guanicoe (B), and Vicugna vicugna (C) according to simmr Bayesian models.

Discussion
Isotopic analysis of bone, dentin, and enamel bioapatite revealed the diet of camelid populations that inhabited southern Brazil during the Pleistocene. Hemiauchenia paradoxa and Lama guanicoe had diets predominantly composed of C3 vegetation, whereas Vicugna vicugna had a mixed diet, as it consumed C3 and C4 plants in very similar proportions.
Studies with modern guanacos and vicuñas have shown that the diet of these species comprises a wide range of plant resources, especially grass species, depending on the seasonal and spatial availability of primary producers (Puig et al. Reference Puig, Videla, Monge and Roig1996, Reference Puig, Videla and Cona1997, Reference Puig, Rosi, Videla and Mendez2011; Torres and Puig Reference Torres and Puig2010; Rivals et al. Reference Rivals, Rindel and Belardi2013; Castellaro et al. Reference Castellaro, Orellana, Escanilla, Bastías, Cerpa and Raggi2020). Hence, based on their modern diet, it is possible to infer that the results found in the present study reflect the consumption of grasses rather than shrubs or trees. In addition, pollen analyses indicate a predominance of grassland environments in southern Brazil during glacial periods (Behling et al. Reference Behling, Pillar and Bauermann2005; Roth et al. Reference Roth, Lorscheitter and Masetto2021), which corroborates the high consumption of grasses by the Lamini camelids.
Similar to modern studies, paleontological and archaeological studies have reported high dietary plasticity, in which populations oscillate between C3 and C4 feeders (MacFadden and Schockey Reference MacFadden and Shockey1997; Feranec Reference Feranec2003; Gil et al. Reference Gil, Ugan, Otaola, Neme, Giardina and Menéndez2016; Yann et al. Reference Yann, DeSantis, Koch and Lundelius2016; González-Guarda et al. Reference González-Guarda, Domingo, Tornero, Pino, Fernández, Sevilla, Villavicencio and Agustí2017; Domingo et al. Reference Domingo, Tomassini, Montalvo, Sanz-Pérez and Alberdi2020; Gayo et al. Reference Gayo, Martens, Stuart-Williams, Fenner, Santoro, Carter and Cameron2020; Szpak and Valenzuela Reference Szpak and Valenzuela2020). Given their opportunistic behavior, the feeding habits of southern Brazilian camelids may have been the result of forage availability, ecological opportunity, and intra- and interspecific competition, as discussed in the following paragraphs.
Variations in camelid diet are mainly linked to climate (Gil et al. Reference Gil, Ugan, Otaola, Neme, Giardina and Menéndez2016). Climatic conditions are responsible for the distribution of grass species and characterize their spatialization at a regional scale (Moreira et al. Reference Moreira, Fontana, Kuplich and Guasselli2019; Preston and Fjellheim Reference Preston and Fjellheim2020). Due to the greater consumption of C3 grasses compared with C4 grasses, the former likely occurred in higher abundance in southern Brazil during the Pleistocene. Grasses that use the C4 photosynthetic pathway require warm growing-season temperatures, and their growth is therefore limited by low temperatures and the occurrence of frost during winter, whereas C3 photosynthesis occurs in cool-season grasses (Boldrini et al. Reference Boldrini, Overbeck, Trevisan, Pillar and Lange2015; Moreira et al. Reference Moreira, Fontana, Kuplich and Guasselli2019). So, the low temperatures of glacial periods may have favored the development of cool-season grasses over tropical grasses. Indeed, grasses that dominate temperate regions are from the subfamily Pooideae (family Poaceae) (Preston and Fjellheim Reference Preston and Fjellheim2020), all of which are C3 photosynthesizers (Longhi-Wagner Reference Longhi-Wagner2012). Pollen records from the coastal plain of northern Rio Grande do Sul state (municipality of Torres: sedimentary profiles dated from 24 ka onward) (Roth et al. Reference Roth, Lorscheitter and Masetto2021) and western Rio Grande do Sul (municipality of São Francisco de Assis: 22 ka to present) (Behling et al. Reference Behling, Pillar and Bauermann2005) corroborate the predominance of Poaceae grasses during the late Quaternary. Nonetheless, taxonomic identification of Poaceae grasses was not provided beyond the family level in both Roth et al. (Reference Roth, Lorscheitter and Masetto2021) and Behling et al. (Reference Behling, Pillar and Bauermann2005), due to pollen morphological similarities among members of this family. The comparable isotopic composition of specimens of Touro Passo (continental fluvial) and Santa Vitória (coastal plain) Formations also indicate similarities in the floristic composition of these two localities, with a predominance of grasslands during glacial periods, as stated by Bauermann et al. (Reference Bauermann, Behling, Macedo, Ribeiro, Bauermann and Scherer2009).
Besides ecological differences, C3 and C4 grasses differ in their nutritional quality. Generally, C4 grasses have lower nitrogen and higher cell wall concentrations than their C3 counterparts (Ehleringer and Monson Reference Ehleringer and Monson1993). The South American camelids are foregut ferments with a great capacity to digest low-quality forages (López et al. Reference López, Maiztegui and Cabrera1998). They are even able to digest C3 and C4 plants more efficiently than hindgut fermenters (Sponheimer et al. Reference Sponheimer, Robinson, Roeder, Hammer, Ayliffer, Passey, Cerling, Dearing and Ehleringer2003). Among them, the vicuña shows the highest degree of specialization for grazing, mainly evidenced by the adaptation of its lower incisors (Menegaz et al. Reference Menegaz, Goin and Jaureguizar1989), which allowed this species to include large proportions of C4 grasses in its diet, avoiding competition with other sympatric camelids. Contemporary vicuñas and guanacos also differ concerning food preferences, despite their generally similar diet (Cassini et al. Reference Cassini, Borgnia, Arzamendia, Benítez, Vilá and Gordon2009). For instance, in western Argentina, where populations of these species occur in sympatry, both guanacos and vicuñas consume several plant species in common; however, the guanacos’ diet comprises a greater diversity of plant species and, consequently, more generalist feeding habits (Cajal Reference Cajal1989). Thus, different foraging behaviors may have favored niche partitioning and the coexistence of species with similar ecologies, as observed in paleocommunities worldwide (e.g., DeSantis et al. Reference DeSantis, Feranec and MacFadden2009; Ma et al. Reference Ma, Wang, Jin, Hu and Bocherens2019).
A factor that contributes to niche partitioning between phylogenetically similar taxa is different habitat use. Direct encounters between living populations of guanacos and vicuñas result in agonistic behavior, which in turn leads to differences in habitat use (Cajal Reference Cajal1989). Unfortunately, such behaviors are almost impossible to assess in nonliving populations.
In summary, the results of this study likely reflect the opportunistic behavior of South American camelids. The higher nutritional quality of C3 grasses, associated with their greater availability in the environment, may explain their high consumption by H. paradoxa and L. guanicoe during the Pleistocene of southern Brazil. Even so, among the guanaco specimens, one individual had an isotopic signature indicating a mixed-feed diet, probably a result of individual variation in its foraging ecology to minimize intraspecific competition (Bolnick et al. Reference Bolnick, Svanbäck, Fordyce, Yang, Davis, Hulsey and Forister2003). Other studies also reported generalist feeding habits in guanaco populations, with some individuals differing from the general population (e.g., Gil et al. Reference Gil, Ugan, Otaola, Neme, Giardina and Menéndez2016), but further efforts are needed to fully understand past ecological dynamics.
The isotopic ratios of oxygen can provide information about paleoclimate, but such results should be interpreted with caution. Despite their strong phylogenetic proximity (Lynch et al. Reference Lynch, Sánchez-Villagra and Balcarcel2020), South American camelids are known to have some differences in their behavioral ecology (see Franklin Reference Franklin1983). For example, vicuñas are obligate water drinkers, whereas guanacos supply their water requirements by forage and are periodic water drinkers (Franklin Reference Franklin1983). The physiology and behavior of the extinct H. paradoxa are difficult to assess; nevertheless, by comparing the oxygen isotope ratios of proboscideans (obligated drinkers) and camelids of this genus, Yann et al. (Reference Yann, DeSantis, Haupt, Romer, Corapi and Ettenson2013) suggested that Hemiauchenia was more sensitive to δ18O variations in plant leaves that they ate. These contrasting behaviors can lead to differences in the oxygen isotope ratios recorded in these species’ tissues (Kohn Reference Kohn1996; Levin et al. Reference Levin, Cerling, Passey, Harris and Ehleringer2006). In addition, C4 feeders tend to have higher δ18O values than C3 feeders, following the isotopic composition of the consumed plants (Kohn Reference Kohn1996). Although South American camelids show different behaviors regarding water drinking and diet, significant differences in δ18O values were not detected between guanacos, vicuñas, and the extinct H. paradoxa. One possible explanation is that both meteoric waters and the ingested plants had similar oxygen isotope compositions. The isotopic composition of plant leaf tissues is strongly influenced by leaf transpiration, which can vary with environmental conditions (Dawson et al. Reference Dawson, Mandelli, Plamboeck, Templer and Tu2002; Barbour Reference Barbour2007; Simonin et al. Reference Simonin, Roddy, Link, Apodaca, Tu, Hu, Dawson and Barbour2013). So, glacial conditions likely provided slightly evaporative environments, leading to a weak fractionation in leaf water compared with source water.
According to oxygen and carbon isotope ratios from specimens analyzed in this study, climate does not influence the diet of the Lamini camelids. It is important to emphasize, however, that all samples were from glacial periods, and there may not have been substantial changes in climate and vegetation to affect the feeding habits of these animals. On the other hand, during interglacial periods of global warming, cool-season grasses likely became restricted to montane areas, such as the South American Andes (Preston and Fjellheim Reference Preston and Fjellheim2020). Likewise, the distribution of the South American camelids followed the distribution of these grasses, as they are specialized grazers (Menegaz et al. Reference Menegaz, Goin and Jaureguizar1989).
Grazing is an important factor in maintaining the Pampa biome in southern Brazil, as stated by Overbeck et al. (Reference Overbeck, Boldrini, do Carmo, Garcia, Moro, Pinto, Trevisan, Zannin, Pillar and Lange2015). Intense grazing pressure favors the development of vegetation composed of grasses and increases the diversity of low-stratum species. In contrast, when the grazing pressure is low or absent, shrubs and tree species develop, reducing the diversity of smaller species. Therefore, the extinction of the Brazilian Pampa megafauna may have led to strong changes in the landscape through reduced grazing pressure. Camelids also play other ecological roles in grasslands, as their excrement and urine increase soil fertility, which contributes to higher plant diversity and forage production (Franklin Reference Franklin1983).
Finally, it is important to recognize the possible limitations associated with stable isotope analysis in fossilized tissues. To date, there are no methods available to assess the preservation of carbon isotope ratios in bioapatite (Domingo et al. Reference Domingo, Tomassini, Montalvo, Sanz-Pérez and Alberdi2020), and to test the preservation of oxygen isotope ratios, it is necessary to analyze both phosphate and carbonate compositions of bioapatite (Prado et al. Reference Prado, Sánchez and Alberdi2011), which requires a greater quantity of samples. However, as noted by Domingo et al. (Reference Domingo, Tomassini, Montalvo, Sanz-Pérez and Alberdi2020), a consistency in isotopic ratios among Pleistocene taxa suggests that they were derived from animal paleobiology rather than diagenesis. In addition, Demény et al. (Reference Demény, Gugora, Kesjár, Lécuyer and Fourel2019) suggested that only δ18O values of enamel should be included in past climate interpretations, due to large variations in OH related to δ18O of bone. Despite this caveat, δ18O values of bones were similar to those of enamel. Thus, we assume that bioapatite from bone samples is also reflecting natural changes in oxygen isotope ratios. Furthermore, all bioapatite δ18O values were obtained under the same analytical conditions, decreasing the chances of possible bias and ensuring that the values are comparable (Demény et al. Reference Demény, Gugora, Kesjár, Lécuyer and Fourel2019).
Conclusions
This was the first study to assess the paleobiology of extinct and pseudo-extinct camelids from southern Brazil. Based on stable isotope analysis, it was possible to reconstruct the diet of animals that inhabited the region during the Pleistocene. Carbon stable isotopes revealed that Hemiauchenia paradoxa (n = 7) fed mainly on C3 grasses. Similar results were found for Lama guanicoe (n = 6), but it showed a broader diet due to one individual that consumed C3 and C4 plants equally (i.e., mixed feeder). Conversely, Vicugna vicugna (n = 4) had a mixed diet. Distinct feeding habits may have contributed to niche partitioning and the coexistence of Lamini camelids. Combined oxygen and carbon isotope data indicated a consistent diet in glacial periods, probably marked by greater availability of C3 grasses compared with C4 grasses, as the former are adapted to cold weather. These results contribute to the understanding of the dynamics of paleocommunities of southern Brazil and the causes of Pleistocene/Holocene extinction. Given their adaptations to grazing, the climate amelioration followed by loss of grasslands likely had a major impact on camelid populations, leading to their extinction. For future studies, we recommend analyzing a larger number of samples to confirm the patterns reported here and also to account for individual variation.
Acknowledgments
We thank the Museu de Ciências Naturais (Secretaria do Meio Ambiente e Infraestrutura–SEMA) for providing fossil samples. This work was supported by the Paleontological Society International Research Program Sepkoski Grant awarded to T.S.C., who also received a scholarship from the Coordination for the Improvement of Higher Education Personnel–Brazil (CAPES) (88882.436025/2019-01). A.M.R. received financial support from the Brazilian National Council for Scientific and Technological Development (CNPq) (306951/2017-7). We thank the reviewers and editors for their contributions, which improved the quality of this article.


