Introduction
During the Pleistocene epoch, a set of large mammals with a body mass of over 44 kg, the megafauna, populated a large portion of our planet. Most of these animals became extinct in the late Pleistocene and Early Holocene, with some exceptions in Africa and Asia, where a diversity of large mammals still persists (Koch and Barnosky, Reference Koch and Barnosky2006; Stuart, Reference Stuart2015; Svenning et al., Reference Svenning, Lemoine, Bergman, Buitenwerf, Le Roux, Lundgren, Mungi and Pedersen2024). South America showed the highest rate of extinction, where about 83% (approximately 50 genera) of large animals became extinct, as compared with only 25% in Africa (Lima-Ribeiro and Diniz-Filho, Reference Lima-Ribeiro and Diniz-Filho2013; Faith, Reference Faith2014). This extinction event primarily affected large mammals, but also significantly impacted other large vertebrates including crocodiles, giant lizards, and giant flightless birds in regions like Australasia (Allentoft et al., Reference Allentoft, Heller, Oskam, Lorenzen, Hale, Gilbert, Jacomb, Holdaway and Bunce2014; Miller et al., Reference Miller, Magee, Smith, Spooner, Baynes, Lehman and Fogel2016; Saltré et al., Reference Saltré, Rodríguez-Rey, Brook, Johnson, Turney, Alroy and Cooper2016), Madagascar (Goodman and Jungers, Reference Goodman and Jungers2014), and Southeast Asia (Meijer et al., Reference Meijer, Louys and O'Connor2019). This distinguishes it from mass extinctions, which typically have broader effects across a wide range of species, including plants.
Several megafauna extinction triggers have been pointed out (Prates and Perez, Reference Prates and Perez2021), such as the role of climate change (Grayson, Reference Grayson, Martin and Klein1984), human hunting (Martin, Reference Martin, Martin and Wright1967; Burney and Flannery, Reference Burney and Flannery2005), hyperdisease (MacPhee, Reference MacPhee, Goodman and Patterson1997; Rothschild and Laub, Reference Rothschild and Laub2006), and extraterrestrial impact (Firestone et al., Reference Firestone, West, Kennett, Becker, Bunch, Revay and Schultz2007). On the South American continent, the issue of megafauna extinction has generated many controversies, because this loss was synchronous with intense climate change and occurred almost at the same time as the initial human influx over the continent. However, specialists in the research field suggest abandoning debates involving a climate versus human hunting dichotomy as evidence pointing to climate‒human mixed factors, varying across different regions, emerges (Barnosky and Lindsey, Reference Barnosky and Lindsey2010; Lima-Ribeiro et al., Reference Lima-Ribeiro, Nogués-Bravo, Terribile, Batra and Diniz-Filho2013; Metcalf et al., Reference Metcalf, Turney, Barnett, Martin, Bray, Vilstrup and Orlando2016; Villavicencio, Reference Villavicencio2016; Araújo et al., Reference Araújo, Machado, Mothé and Avilla2021; Prates and Perez, Reference Prates and Perez2021). Thus, simultaneously with climate-induced changes in vegetation, impacts from other causes, such as human hunting, are recognized to have enhanced the risk of extinction under the synergistic effect of multiple causes (Cione et al., Reference Cione, Tonni and Soibelzon2003, 2009; Grayson, Reference Grayson2007).
Most of the South American megafauna were associated with open environments that were influenced by Quaternary climate change (Bond et al., Reference Bond, Cerdeño, López, Alberdi, Leone and Tonni1995; Tonni and Scillato Yané, Reference Tonni and Scillato Yané1997; Cione et al., Reference Cione, Tonni, Soibelzon and Haynes2009). Despite their short duration, interglacial periods, characterized by increased temperatures and humidity, led to a shrinking of open vegetation areas and a reduction in biomass of the mammals adapted to this environment (Cione et al., Reference Cione, Tonni, Soibelzon and Haynes2009). On the other hand, during longer glacial periods, open-habitat-adapted mammal populations recovered. This alternation between open and closed vegetation areas versus high and low mammalian biomass is called the Broken Zig-Zag hypothesis. According to the hypothesis, periodic severe climatic changes during the Pleistocene caused dramatic modifications in the distribution and biomass of the biota in South America, and the expression zig-zag is used to stress the periodicity of the phenomenon (Cione et al., Reference Cione, Tonni and Soibelzon2003). Based on this hypothesis, climate change during the late Pleistocene and Early Holocene is suggested to be the primary driver of changes in vegetation that affected food availability, reduced suitable habitats for megafauna, and caused a collapse in the geographic distribution of species, contributing to a reduction in their populations and an increased risk of extinction (Lopes et al., Reference Lopes, Pereira, Kerber and Dillenburg2020; Araújo et al., Reference Araújo, Machado, Mothé and Avilla2021).
Following the conceptual model of habitat restriction, in this study, we aimed to assess climate change effects on open vegetation habitat and, consequently, on megafauna extinction in South America by using a species distribution model (SDM), fossil records, and paleoclimatic projections (Nogués-Bravo, Reference Nogués-Bravo2009; Varela et al., Reference Varela, Lobo and Hortal2011), using taxa identified in the literature as being restricted to open habitats. In contrast to many studies that focus SDMs on one species or a group of related species, here, we intended to model the suitable habitat variation to which open-area megafauna were adapted. In addition, we evaluated the effects of climatic variables (temperature and precipitation) on the distribution of suitable habitats across regions to test whether these effects were homogeneous throughout South America. The hypothesis raised here is that if climate change increased the risk of extinction, it would have caused an important habitat restriction, specifically, considerable loss of open habitats. On the other hand, if climatic changes were not the only predictor, then the climatic conditions suitable for these taxa should not have totally disappeared and should have persisted throughout the Holocene.
Methods
We used an SDM to determine the potential distribution of suitable habitat for extinct South American megafauna taxa, guided by paleoecological studies indicating their restriction to open vegetation areas. These open areas were systems generally constrained by water availability, with a herbaceous layer dominated by grass species, and may have been accompanied by a variable component of woody vegetation (Bilenca and Miñarro, Reference Bilenca and Miñarro2004). We used this approach because species restricted to these open areas would exhibit strong responses to Pleistocene climatic oscillations (Peçanha et al., Reference Peçanha, Althoff, Galiano, Quintela, Maestri, Gonçalves and Freitas2017).
SDMs work as a combination of environmental variables and occurrence records of a species or taxon group to identify the environmental conditions where populations may be maintained (Dalapicolla, Reference Dalapicolla2016). The general objective of SDMs is to obtain an environmental suitability map that defines which locations are more or less suitable for the survival of species, focusing on their ecological requirements, which is called a “modeled geographic distribution” or a “predictive map” (Elith and Leathwick, Reference Elith and Leathwick2009). SDMs are very effective in showing potential habitat suitability for species survival across the landscape; therefore, they have become a popular tool to determine the potential geographic distribution of now-extinct taxa (Nogués-Bravo, Reference Nogués-Bravo2009). By using a multispecies model to assess open vegetation habitat variation throughout the late Pleistocene and Holocene, we intended to minimize dispersal and inter- and intraspecific interactions that tend to generate over- or underestimated distribution models (Araújo et al., Reference Araújo, Machado, Mothé and Avilla2021). Multispecies models consider the ecological relationships, interactions, and shared environmental preferences among multiple species. We used this approach because it describes shared environmental responses (Ovaskainen and Soininen, Reference Ovaskainen and Soininen2011; Pollock et al., Reference Pollock, Tingley, Morris, Golding, O'Hara, Parris, Vesk and McCarthy2014), which is important information in the face of the climatic variation felt throughout the Pleistocene and to understand its role in extinctions.
The predictive maps were generated using the Maximum Entropy (MaxEnt) algorithm, v.3.3.3k (Phillips et al., Reference Phillips, Anderson and Schapire2006). MaxEnt models express a probability distribution where each grid cell has a predicted suitability of conditions for the species or group of species, obtained from a set of environmental variables and georeferenced occurrence locations (Elith et al., Reference Elith, Philips, Hastie, Dudik, Chee and Yates2011). MaxEnt uses only occurrence data and background points of the species because absence data are often rare and unreliable, especially when working with extinct species (Phillips et al., Reference Phillips, Anderson and Schapire2006; Araújo et al., Reference Araújo, Machado, Mothé and Avilla2021). Grayson (Reference Grayson1981) argued that the paleozoological record yields asymmetrical data concerning the presence/absence of taxa. He argued that data asymmetry arises because occurrence data only offer evidence regarding the presence of taxa in the area where their remains are recovered. However, he highlights that the absence of remains does not necessarily indicate the absence of a species. Information on the absence of taxa is often lacking for several reasons, such as the possibility that the remains of a taxon were not preserved, they were not collected due to inadequate sampling strategies, or their failure to be recovered was due to inappropriate collection techniques. Because of these and other shortcomings, Grayson (Reference Grayson1981) recommended that analyses that depend only upon the taxa recorded as present within a fauna are to be preferred. Since then, other authors have provided additional arguments and examples (Lyman, Reference Lyman1995; Ervynck, Reference Ervynck1999; Etnier, Reference Etnier2002). Additionally, Lyman (Reference Lyman2008) estimated the magnitude of data asymmetry and emphasized the importance of using larger samples, as they offer a more accurate database with less asymmetry.
The highest value of the area under the curve (AUC) was used to assess the model performance. The receiver operator characteristic (ROC) statistic was used to assess model accuracy, with 10 repetitions of 10,000 maximum interactions, 10% of the mean replications were randomized as test data, and the remainder were randomized to model training in each replication. The “maximum sensitivity plus specificity threshold” was used as the cutoff value or tolerance threshold and identified the suitable areas (SAs). This is a common approach when using MaxEnt, since it maximizes the cases where the model erroneously assigns unsuitable habitat (true negative) and misses suitable habitat (false positive) (Ferraz et al., Reference Ferraz, de Siqueira, Alexandrino, da Luz and do Couto2012; Jorge et al., Reference Jorge, Galetti, Ribeiro and Ferraz2013; Liu et al., Reference Liu, White and Newell2013; Vale et al., Reference Vale, Tarroso and Brito2014). SDMs were performed in RStudio software (version 1.1.453) using the “dismo” package (Hijmans et al., Reference Hijmans, Phillips, Leathwick and Elith2021).
The taxa were selected considering literature references focusing on grazing behavior, diet of C4 plants (grass), and occurrence in open areas, based on dental morphology and studies involving carbon isotopes. Our database consisted of four taxa: Equus neogeus (grazer with a diet consisting exclusively of C4 grasses) (Dantas et al., Reference Dantas, Cherkinsky, Bocherens, Drefahl, Bernardes and França2017; Pansani et al., Reference Pansani, Muniz, Cherkinsky, Pacheco and Dantas2019), Glyptodon clavipes (grazer adapted to open environments) (Pérez et al., Reference Pérez, Scillato-Yané and Vizcaíno2000; Zurita and Miño-Boilini, Reference Zurita and Miño-Boilini2012), Panochthus spp. (grazer, mesoherbivore specialist, with a diet consisting of 96% C4 plants and skull morphological features adapted to open environments) (Zurita et al., Reference Zurita, Scillato-Yané and Carlini2005; Xavier et al., Reference Xavier, Dantas and Silva-Santana2018; Silva et al., Reference Silva, Leal, Cherkinsky and Dantas2019; Cuadrelli et al., Reference Cuadrelli, Zamorano, Barasoain, Anaya and Zurita2023), and Pampatherium humboldtii (grazer, with skeletal, dental, and masticatory musculature adaptations to herbivory) (De Luliis et al., Reference De Luliis, Bargo and Vizcaíno2001).
We obtained georeferenced information for fossil occurrences from bibliographic surveys of fossil findings in South America, from Prof. Dr. Mario Dantas's (Federal University of Bahia) personal data (Dantas et al., Reference Dantas, Pansani, Asevedo, Araújo, França, Aragão, Santos, Cravo, Waldherr and Ximenes2024), and from the Paleobiology Database (https://paleobiodb.org; accessed February 7, 2024) (Supplementary Tables 1 and 2). We were able to collect a total of 292 georeferenced records (Fig. 1). The data represent the occurrences of fossils from the selected taxa that existed during the last glacial maximum (LGM). Most of the compilated data correspond to stratigraphically dated specimens, while a few correspond to specimens with absolute dating. Therefore, the dating of the majority of occurrences span broad time intervals, but we selected all those demonstrating an occurrence during the LGM and assumed they were representative of that time, which is a method that has been applied to many other previously conducted modeling studies (Varela and Fariña, Reference Varela and Fariña2016; Villavicencio et al., Reference Villavicencio, Corcoran and Marquet2019; Araújo et al., Reference Araújo, Machado, Mothé and Avilla2021). Furthermore, when examining studies with absolute dating conducted in South America, out of 62 dates across more than 40 publications (data compiled by Dantas et al., Reference Dantas, Pansani, Asevedo, Araújo, França, Aragão, Santos, Cravo, Waldherr and Ximenes2024), we observed that the majority, about 71%, are over 17.0 ka (Heinrich stadial), indicating a greater presence of these animals during the LGM and earlier times.
Geographic distribution of the selected megafauna taxa representative of last glacial maximum (ca. 21 ka) open vegetation areas in South America. Abbreviations: FG, French Guiana; SU, Suriname; GU, Guyana; VE, Venezuela; CO, Colombia; EC, Ecuador; PE, Peru; BR, Brazil; BO, Bolivia; PY, Paraguay; CL, Chile; UY, Uruguay; AR, Argentina. The human used for scale has a height of approximately 1.8 m. Map background made with Natural Earth. Free vector and raster map data @ naturalearthdata.com.
We combined the selected occurrence records with bioclimatic data to project the potential geographic distribution of the selected taxa during the LGM; we then used the obtained model to predict their potential distribution during the Heinrich stadial (HS: 17.0–14.7 ka), Bølling-Allerød (BA: 14.7–12.9 ka), Younger Dryas stadial (YDS: 12.9–11.7 ka), Early Holocene, Greenlandian (EH: 11.7–8.3 ka), Mid Holocene, Northgrippian (MH: 8.3–4.2 ka), and Late Holocene, Meghalayan (LH: 4.2–0.3 ka) (Fordham et al., Reference Fordham, Saltré, Haythorne, Wigley, Otto-Bliesner, Chan and Brook2017).
The bioclimatic variables were obtained from www.paleoclim.org, for the periods mentioned above, which provided 19 variables from several general circulation models (Brown et al., Reference Brown, Hill, Dolan, Carnaval and Haywood2018). The bioclimatic variables were downloaded at a spatial resolution of 2.5 arc-minutes. To reduce the use of highly correlated variables, we ran Pearson's correlation within the “vegan package” (Oksanen et al., Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn and Minchin2022) for RStudio. Highly correlated variables (r > 0.8) were eliminated to avoid redundant climate variables (Nascimento et al., Reference Nascimento, Agne, Batalha-Filho and Araujo2021). The variables used were annual mean temperature (BIO1), mean diurnal range (BIO2), isothermality (BIO3), annual temperature range (BIO7), annual precipitation (BIO12), precipitation seasonality (BIO15), precipitation of driest quarter (BIO17), precipitation of warmest quarter (BIO18), and precipitation of coldest quarter (BIO19).
To test whether the effects of temperature and precipitation on the SA distributions were homogeneous throughout South America, we adopted two strategies. First, we grouped regions based on the disjunct distribution of SAs and the geographic distinction between them. Second, we used the generalized linear mixed model (GLMM) to find possible associations between climatic variables and SA distributions. A predictor in the GLMM contains random effects typically represented by the grouping variable to allow for the estimation of variance in the response variable within and among these groups (Harrison et al., Reference Harrison, Donaldson, Correa-Cano, Evans, Fisher, Goodwin, Robinson, Hodgson and Inger2018). As temperature and precipitation are not independent of location in South America, the regions identified in the first step were used for random effects in our analysis. We used the temperature and precipitation averages from each region and period as fixed predictors and the SA distribution from respective periods as the response variable. The GLMM was performed in RStudio using the “lme4” package (Bates et al., Reference Bates, Mächler, Bolker and Walker2015). To check the correct distribution, dispersion, and outliers from our model, we used the “DHARMa” package (Florian, Reference Florian2022).
Results
We obtained an AUC value of >0.905, a very good model performance, which suggested that the obtained model had a good measure of separability and was influenced by environmental variables. Annual precipitation was the variable that most influenced the model (62.9%), followed by the annual mean temperature (14.6%) (Table 1). The highest values of the potential geographic distribution of the selected taxa correspond to ca. 400–500 mm of annual precipitation and between 5 and 15°C of annual mean temperature (Fig. 2).
Relative contributions (%) of the environmental variables to the MaxEnt model for suitable areas of megafauna taxa restricted to open vegetation in South America
The curves show the association between two environmental variables and the presence probability of suitable areas for four megafauna taxa restricted to open vegetation in South America, during the last glacial maximum (ca. 21 ka).
During the LGM, the potential geographic distribution for the selected taxa predicted that open habitats were concentrated mainly in the Andean region, Brazilian Northeast, and the Chaco–Pampas region (Figs. 3 and 4). During all the other periods, the SAs remained in the same regions where they were in the LGM, with some smooth displacements and visible contraction in the Chaco–Pampas region. Significant SA shrinkage in the Brazilian Northeast was observed, with evidence of almost complete disappearance during the HS period (Fig. 4, Table 2).
Suitable areas (SA) calculated according to maximum test sensitivity plus specificity threshold during the last glacial maximum (LGM), Heinrich stadial (HS), Bølling-Allerød (BA), Younger Dryas stadial (YDS), Early Holocene (EH), Mid Holocene (MH), and Late Holocene (LH) and percentage of area change (contraction/expansion) for the predicted periods compared to the LGM
Potential geographic distribution of the selected taxa during the last glacial maximum (LGM) (ca. 21 ka), Heinrich stadial (HS) (17.0–14.7 ka), Bølling-Allerød (BA) (14.7–12.9 ka), Younger Dryas stadial (YDS) (12.9–11.7 ka), Early Holocene, Greenlandian (EH) (11.7–8.3 ka), Mid Holocene, Northgrippian (MH) (8.3-4.2 ka), and Late Holocene, Meghalayan (LH) (4.2–0.3 ka). The area under the curve (AUC) indicates model performance.
Suitable areas for open vegetation during the last glacial maximum (LGM) (ca. 21 ka), Heinrich stadial (HS) (17.0–14.7 ka), Bølling-Allerød (BA) (14.7–12.9 ka), Younger Dryas stadial (YDS) (12.9–11.7 ka), Early Holocene (EH) (Greenlandian, 11.7–8.3 ka), Mid Holocene (MH) (Northgrippian, 8.3–4.2 ka), and Late Holocene (LH) (Meghalayan, 4.2–0.3 ka), considering all values superior to 0.3863 of probability (using maximum test sensitivity plus specificity threshold). The predicted map is divided into three major regions: Brazilian Northeast, Chaco–Pampas, and Andean.
Throughout the modeled periods, alternating moments of expansion and contraction of SAs for the megafauna were noticeable (Table 2). The HS and YDS were characterized by a severe decrease in SAs in the Brazilian Northeast, while in the Andean and the Chaco–Pampas regions, there was an increase in suitable habitat. From the EH onwards, SAs decreased progressively in all regions. Regarding the more recent modeled period (LH), all regions suffered climatic changes that reduced suitable habitat for megafauna: relatively small habitat shrinkage in the Andean and Chaco–Pampas regions (−14.5% and −11.41%, respectively) and significant habitat contraction in the Brazilian Northeast (−75.5%). Therefore, it is clear that climate condition fluctuations were different across different regions of South America, resulting in distinct changes in SAs for open-area megafauna during the late Pleistocene and Holocene (Supplementary Figs. 1 and 2).
Precipitation was negatively correlated with the distribution of SAs for open-area megafauna (Table 3, Fig. 5). This association was not observed for temperature and the interaction between precipitation and temperature (Table 3). The effect of precipitation was responsible for drastically reducing the SA in the Brazilian Northeast, where the previous SA distribution was already lower, and the precipitation variation was broader than that of the Andean and Chaco–Pampas regions (Fig. 5).
Effects of precipitation and temperature and the interaction between them on the distribution of suitable areas for open-area megafauna in South America. The results are from the generalized linear mixed model in which the linear predictor contains as random effects the three different regions: the Andean, Brazilian Northeast, and Chaco–Pampas
Significance: a 0.001; b 0.01.
Effects of precipitation on the distribution of suitable areas for open-area megafauna in the three different regions of South America (Andean, Brazilian Northeast, and Chaco–Pampas) during the last glacial maximum (LGM) (ca. 21 ka), Heinrich stadial (HS) (17.0–14.7 ka), Bølling-Allerød (BA) (14.7–12.9 ka), Younger Dryas stadial (YDS) (12.9–11.7 ka), Early Holocene, Greenlandian (EH) (11.7–8.3 ka), Mid Holocene, Northgrippian (MH) (8.3–4.2 ka), and Late Holocene, Meghalayan (LH) (4.2–0.3 ka).
Discussion
Our results demonstrate that South American climate changes, especially precipitation variations, during the late Pleistocene and Holocene, led to alternating periods of decrease and increase in megafauna-suitable habitat, supporting the Broken Zig-Zag hypothesis (Cione et al., Reference Cione, Tonni and Soibelzon2003). However, the amplitude of these fluctuations varied among different regions: in the Andean and Chaco–Pampas regions, the oscillation was more constrained, while in the Brazilian Northeast, it was broader. The Brazilian Northeast, with an already smaller SA, suffered more intense oscillations in precipitation over the past 21 ka, strongly impacting the extent of adequate habitat for the analyzed megafauna taxa. Yet, despite fluctuations, climatic conditions suitable for open vegetation areas persisted into the Holocene. This suggests that climate change was not the sole predictor of megafauna extinction in South America, but it may have been a more important determinant in certain regions, specifically in the Brazilian Northeast.
An important outcome of our study is that it reveals that we cannot consider the extinction of megafauna in South America to be a uniform or homogeneous event across the continent. Climatic fluctuations throughout the late Pleistocene and Holocene differed in various regions of the continent, leading to distinct variations in SA, which suggests that the extinction of South American megafauna occurred differently in different regions, and its causes may have been diverse as well.
Precipitation, not temperature, had the greatest influence on the SAs for the occurrence of megafauna in open vegetation areas. The influence of precipitation was expected since phenological cycles in fields, pastures, and other open areas are strongly seasonal and governed by rainfall. As it is responsible for water availability in the system, precipitation is generally considered to be the main factor that controls the structure and dynamics of the ecosystem and drives the biological processes (Notaro et al, Reference Notaro, Liu, Gallimore, Williams, Gutzler and Collins2010; Chen et al., Reference Chen, Wang and Fu2020). In addition, in the tropics, changes in precipitation are often more significant and rapid than changes in temperature (Bonebrake and Mastrandrea, Reference Bonebrake and Mastrandrea2010). This is due to large-scale atmospheric processes, such as the El Niño–Southern Oscillation and the North Tropical Atlantic mode, which can cause dramatic changes in precipitation patterns in South America (Marengo et al., Reference Marengo, Alves, Alvala, Cunha, Brito and Maraes2018; Ding et al., Reference Ding, Nnamchi, Yu, Li, Sun, Li and Tseng2023).
The more dramatic oscillation in tropical precipitation may help us to better understand the drastic fluctuations in SAs in the Brazilian Northeast during Pleistocene–Holocene times. The HS period, in which the climatically favorable areas for open habitats almost disappeared, corresponds to the period in which precipitation almost doubled (from 780 to 1450 mm). In the following period, the BA, a combination of lower precipitation and higher temperatures supported the recovery of open areas and, therefore, of SAs for large mammals specifically adapted to this type of environment. In the YDS, SAs shrunk again due to increased precipitation (greater than 1100 mm). Both the HS and YDS periods were characterized by a Northern Hemisphere decrease in temperature, possibly due to inputs of fresh and cold water into the North Atlantic, changing the density and patterns of the thermohaline circulation (Heinrich, Reference Heinrich1988). As a result, atmospheric circulation patterns may also have been affected, as colder waters in the North Atlantic would have created a high-pressure center in the region, forcing the intertropical convergence zone to descend to latitudes that cover the Brazilian Northeast, translating to increased humidity and rainfall.
The influence of North Atlantic sea surface temperatures on precipitation levels in the Brazilian Northeast, which persists to the present day (Mao et al., Reference Mao, Zou, Alves, Macau, Taschetto, Santoso and Kurths2022), may explain the significant increase in precipitation in the region during the HS and YDS and the resulting drastic reduction of SAs. Thus, the drastic reduction of SAs in the Brazilian Northeast leads us to think that in this particular region, climate change could indeed have deeply influenced the disappearance of the open-habitat megamammals over 44 kg. However, it cannot be the only driver, as there is fossil evidence in the region indicating that open-area megafauna survived the HS and YDS, including Toxodon platensis dated to 11,400–11,750 cal yr BP (10,050 ± 30 14C yr BP) (Dantas et al., Reference Dantas, Dutra, Cherkinsky, Fortier, Kamino, Cozzuol, Ribeiro and Vieira2013). Additionally, Faria et al. (Reference Faria, Carvalho, Júnior, Ximenes, Facincani, Machado, Machado and Pohlmann Bulsing2023) inferred very young ages for several megafaunal species (ranging from 3493 to 8189 cal yr BP) recovered from Itapipoca and Miranda (Ceará and Mato Grosso do Sul states, Brazilian Northeast).
Although there is no information about the dating of the fossil taxa used in the present research from the Brazilian Northeast, the dating of other megamammal fossils, such as Catonyx cuvieri, Eremotherium laurillardi, Nothrotherium maquinense, Notiomastodon platensis, T. platensis, and Xenorhinotherium bahiense (data compiled by Dantas et al., Reference Dantas, Pansani, Asevedo, Araújo, França, Aragão, Santos, Cravo, Waldherr and Ximenes2024), can support the chronology appointed by SA reductions shown here. Of all 36 recorded dates for these species in the Brazilian Northeast, 70% were dated to before the HS, and the remaining 30% were dated to between the EH and BA. These results suggest the need for more studies and fossil dating in order to shed more light on the reasons that led to the extinction of megamammals in the Brazilian Northeast.
In the Chaco–Pampas region, despite the reduction of SAs for open-area megafauna, partially due to the evident continental platform reduction, over 2 million km2 of SAs were still modeled for the EH, a reduction of only 11% when compared with the LGM. Additionally, a significant area of suitable habitat was also modeled for the EH in the Andean region. Therefore, our results suggest that climatic changes, which generated dynamic environmental variations between open and closed areas between the late Pleistocene and Holocene, were not able to explain the megafauna extinction in the Andean and Chaco–Pampas regions.
Distribution models based on single species (E. neogeus, N. platensis, and Hippidion spp.) also suggest a reduction in SAs during this time period (Villavicencio et al., Reference Villavicencio, Corcoran and Marquet2019; Araújo et al., Reference Araújo, Machado, Mothé and Avilla2021) but not their total disappearance. Although a reduction of about 61% of potential area of distribution has been suggested for E. neogeus, for instance, 1,043,449 km2 of modeled suitable area remained during the Holocene (Villavicencio et al., Reference Villavicencio, Corcoran and Marquet2019).
Researchers widely agree that, in conjunction with climate change, the detrimental effects of factors such as human hunting exacerbate the risk of extinction through the combined effect of multiple causes (Cione et al., Reference Cione, Tonni and Soibelzon2003, 2009; Grayson, Reference Grayson2007), and this could explain the megafauna extinction in the Andean and Chaco–Pampas regions (Barnosky and Lindsey, Reference Barnosky and Lindsey2010; Prates and Perez, Reference Prates and Perez2021). For instance, human predatory behavior has been suggested as a main factor driving the megafaunal decline in South America through the demonstration of a strong spatial and temporal relationship between fishtail projectile point technology, which is related to large mammal hunting, and the density and distribution of large mammal species, as well as with the distribution and fluctuation in human demography during this critical period (Prates and Perez, Reference Prates and Perez2021). An alternative theory investigated in North America, as suggested by Faith (Reference Faith2011), proposed that late-glacial climate change, particularly increases in atmospheric CO2, shifted herbivore ecosystem dynamics from a nutrient-accelerating mode to a nutrient-decelerating mode at the end of the Pleistocene, leading to reduced forage production, reduced rates of nutrient cycling through the food web, and ultimately reduced megafaunal population densities.
Our paper contributes to the discussion on the extinction of megafauna in South America by showing that climate change, particularly precipitation fluctuations, had an impact on the extent of SAs for large mammals restricted to open areas during the late Pleistocene and Holocene. Therefore, as precipitation fluctuations were more intense in the Brazilian Northeast, which had a smaller suitable habitat to begin with, the climate change factor was more significant in this region. In contrast, the Andean and Chaco–Pampas regions experienced less variability in precipitation and, thus, had a larger SA that persisted until the Holocene, suggesting that climate change was a less critical factor in explaining the extinction of megafauna in these regions.
Conclusions
The considerable amount of area that remains viable in the LH for open habitats and for the megafauna adapted to it in the Andean and Chaco–Pampas regions leads us to believe that habitat reduction induced by climatic changes was not the main factor that contributed to the extinction of these animals. On the other hand, in the Brazilian Northeast, where the SA was already smaller, climatic changes, more specifically precipitation patterns shifts, were more drastic and may have been decisive for the disappearance of these species.
Alternating intervals of expansion and contraction of SAs for megafauna persistence, mainly in response to lower and higher precipitation, have been observed over the last 21 ka in all regions of South America. However, the amplitude of this zig-zag oscillation was more significant in the Brazilian Northeast. In the Andean and Chaco–Pampas regions, greater precipitation stability resulted in greater stability in habitat suitability; therefore, for these regions, other factors must have been predominant for megafauna extinction. We conclude that in the Andean and Chaco–Pampas regions, climate change was not solely responsible for the disappearance of megafauna, but in the Brazilian Northeast, it may have been decisive.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2024.33.
Acknowledgments
A special thanks to Prof. Dr. Mario Dantas, from the Federal University of Bahia (UFBA), for the provided georeferenced information of fossil occurrences. We are grateful to the anonymous reviewers and editors for their valuable suggestions during the drafting of this manuscript.
