INTRODUCTION
Equidae and Proboscidea are the only lineages of South American megamammals with Holarctic origins that went extinct at the Pleistocene/Holocene transition (Avilla and Mothé, Reference Avilla and Mothé2013; Machado et al., Reference Machado, Grillo, Scott and Avilla2018). The South American diversity of equids and proboscideans was low, mainly represented by Equus neogeus (Lund, Reference Lund1840) and Notiomastodon platensis (Ameghino, Reference Ameghino and Torcelli1888), respectively (Mothé et al., Reference Mothé, Avilla, Asevedo, Borges-Silva, Rosas, Labarca-Encina and Souberlich2017a; Machado and Avilla, Reference Machado and Avilla2019). Both species were highly polymorphic, suggesting great genetic variability within the lineages. In addition, E. neogeus and N. platensis were widely distributed across the South American continent, with records from the high altitudes of the Andean Cordillera to the coastal and inland plains (Mothé et al., Reference Mothé, Avilla, Cozzuol and Winck2012, Reference Mothé, Avilla, Asevedo, Borges-Silva, Rosas, Labarca-Encina and Souberlich2017a; Machado et al., Reference Machado, Grillo, Scott and Avilla2018; Machado and Avilla, Reference Machado and Avilla2019).
The Proboscidea's arrival to South America is possibly related to the first migratory pulse of the Great American Biotic Interchange (GABI 1), just after the uplift of the Panamanian Isthmus ca. 3.1–2.8 million years ago (Woodburne, Reference Woodburne2010). Although Proboscidea may have arrived earlier (Campbell et al., Reference Campbell, Frailey and Romero Pittman2000), they are one of the most common megamammals in Quaternary sites, especially those that date to the middle to late Pleistocene (Mothé et al., Reference Mothé, Avilla, Cozzuol and Winck2012, Reference Mothé, Avilla, Asevedo, Borges-Silva, Rosas, Labarca-Encina and Souberlich2017a). Paleoecological inferences suggest that N. platensis had an opportunistic/generalist feeding habit, consuming grasses, branches, leaves, and fruits (Mothé et al., Reference Mothé, Ferretti and Avilla2017b; González-Guarda et al., Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018), which may have influenced its wide geographic distribution in South America.
After an extensive review of South American fossil horses, Machado et al. (Reference Machado, Grillo, Scott and Avilla2018) and Machado and Avilla (Reference Machado and Avilla2019) considered E. neogeus as the only valid species for Equus in this continent (Barrón-Ortiz et al., Reference Barrón-Ortiz, Avilla, Jass, Bravo-Cuevas, Machado and Mothé2019), reinforcing the importance of updated taxonomy for paleoecological and biogeographic studies (de Oliveira et al., Reference de Oliveira, Araújo, Rotti, Mothé, Rivals and Avilla2020). The origin of Equus possibly occurred during the early Pleistocene of North America, with subsequent dispersion to South America during GABI 3 and GABI 4 (between 0.8 Ma and 0.125 Ma; Woodburne, Reference Woodburne2010; MacFadden, Reference MacFadden2013; Rook et al., Reference Rook, Bernor, Avilla, Cirilli, Flynn, Jukar, Sanders, Scott and Wang2019). Based on cranial and dental morphology, the Pleistocene Equus is considered a restricted grazer, and thus an indicator of open and arid environments (MacFadden, Reference MacFadden2005; Semprebon et al., Reference Semprebon, Rivals, Solounias and Hulbert2016). Furthermore, stable isotopes of equids from different localities from the Pleistocene of South America revealed distinct carbon and hydrogen values for E. neogeus, indicating variations in C3 vegetation (possibly grass), C4 grass, and a mix of C3 and C4 vegetation (Sánchez et al., Reference Sánchez, Prado and Alberdi2004).
The combination of these paleoecological features suggests that N. platensis and E. neogeus should have been less vulnerable to negative pressures of natural selection during the Pleistocene/Holocene transition. Indeed, Mann et al. (Reference Mann, Groves, Gaglioti and Shapiro2019) proposed that several characteristics of the megafauna (e.g., low reproductive rates, large size, and body mass) could have been advantageous for their survival under the unstable environmental and climatic conditions that prevailed after the last glacial maximum.
Species distribution models (SDMs) have become a popular tool to determine the potential geographic distribution of now-extinct species (Nogués-Bravo, Reference Nogués-Bravo2009) and Quaternary megafaunal extinction dynamics (Lima-Ribeiro and Diniz-Filho, Reference Lima-Ribeiro and Diniz-Filho2012). In general, species distribution models require occurrence records of the species under study, as well as bioclimatic variables, and employ analytical models to relate species occurrences to bioclimatic predictors. However, SDMs vary from the Hutchinsonian niche concept in that they are based only on abiotic data and occurrence records, and disregard dispersal and inter- and intraspecific interactions. Nonetheless, the result—a potential geographic distribution (PGD) map—can identify areas with environmental conditions similar to the localities of known species occurrences and show the variation in habitat suitability for species survival across the landscape (Guisan et al., Reference Guisan, Edwards and Hastie2002; Soberón and Peterson, Reference Soberón and Peterson2005; Franklin, Reference Franklin2009; Lima-Ribeiro and Diniz-Filho, Reference Lima-Ribeiro and Diniz-Filho2013).
Thus, considering the great resilience of N. platensis and E. neogeus, and the hypothesis that climate changes affected their population dynamics and geographic distribution during the LGM, we applied an SDM to help to understand their ecological niches and to recognize the influence of changing climatic conditions in South America on their extinction during the late Quaternary. We note, however, that the climate change hypothesis does not exclude a possible overkill scenario nor a combination of climate change and human action as the driver of the megafauna extinctions (Cione et al., Reference Cione, Tonni and Soilbenzon2003; Koch and Barnosky, Reference Koch and Barnosky2006; Sandom et al., Reference Sandom, Faurby, Sandel and Svenning2014).
METHODS
We used an SDM to determine the potential geographic distribution for N. platensis and E. neogeus during the LGM (Ray and Adams, Reference Ray and Adams2001) and to estimate their suitable areas (SA) during the middle Holocene (6 ka; Melo and Marengo, Reference Melo and Marengo2008). First, the geographic distributions of N. platensis and E. neogeus were reviewed, considering only specimens identified at the species level by the most recent and reliable taxonomic reviews (Mothé et al., Reference Mothé, Avilla, Asevedo, Borges-Silva, Rosas, Labarca-Encina and Souberlich2017a; Machado and Avilla Reference Machado and Avilla2019). We obtained georeferenced information for the fossil occurrences from previous studies that were referenced to the World Geodetic System 1984 (WGS 84) (Supplemental Table 1). Occurrence records were georeferenced to the WGS 84 and the Geographic Coordinate System in decimal degrees. Our database comprises records dated between 25‒15 ka (LGM interval; Ray and Adams, Reference Ray and Adams2001); records outside this temporal interval were not included.
We combined the selected occurrence records with bioclimatic data to project the potential geographic distribution of the two species during the LGM and middle Holocene (Franklin, Reference Franklin2009; Nogués-Bravo, Reference Nogués-Bravo2009). The predictive maps were generated in Maxent v.3.3.3k (Phillips et al., Reference Phillips, Anderson and Schapire2006), which is based on an algorithm that estimates the species’ distribution through maximum entropy modeling to predict the SA according to the potential niche of the species (Coelho et al., Reference Coelho, Carvalho and Gomide2016). The Maxent program uses only occurrence data and background points of the species (Anderson et al., Reference Anderson, Lew and Townsend Peterson2003) because absence data are often rare and unreliable, especially for extinct species (Amaro and Morais, Reference Amaro and Morais2014).
Ten thousand background points were sampled and 25% of the original occurrence dataset points were randomly selected to test the predictions and to evaluate the performance of the models (see Barve et al., Reference Barve, Barve, Jiménez-Valverde, Lira-Noriega, Maher, Townsend Peterson, Soberón and Villalobos2011; Owens et al., Reference Owens, Campbell, Dornak, Saupe, Barve, Soberón, Ingenloff, Lira-Noriega, Hensz, Myers and Peterson2013). Additionally, five replicates of the models were used to generate the final model. Runs were conducted with a logistic output format, resulting in a map of habitat suitability of the species ranging from 0 to 1 per grid cell. The parameterization was adjusted according to specific characteristics of the species and the data used in the modeling following the recommendations of Merow et al. (Reference Merow, Smith and Silander2013). The models were evaluated by producing the receiver operating characteristic curve and calculating the area under the curve using the maximum training sensitivity plus specificity as the threshold value (Basher et al., Reference Basher, Bowden and Costello2014). Cross-validation was performed to evaluate the model accuracy.
We extracted layers of bioclimatic variables from the WorldClim Project version 1.4 (Hijmans et al., Reference Hijmans, Cameron, Parra, Jones and Jarvis2005) for 21 ka (LGM) and 6 ka (middle Holocene), which provides 19 variables from several general circulation models. In this study, the LGM layers were used as the model predictors and the middle Holocene layers as the model projections. The bioclimatic variables were downloaded at a spatial resolution of 2.5 arc-minutes, and the general circulation model used was the Model for Interdisciplinary Research on Climate-Earth System Model (MIROC-ESM; Watanabe, Reference Watanabe, Hajima, Sudo, Nagashima, Takemura, Okajima and Nozawa2011). Only temperature and precipitation variables were included, resulting in a high correlation between the bioclimatic layers. Because it is recommended against using layers with highly correlated variables, the multicollinearity of climatic data was examined using the Pearson correlation coefficient, and the least correlated layers were then selected (Guisan and Thuiller, Reference Guisan and Thuiller2005; Giannini et al., Reference Giannini, Siqueira, Acosta, Barreto, Saraiva and Alves-dos-Santos2012). The variables used for E. neogeus were mean diurnal range, mean temperature of warmest quarter, annual precipitation, and precipitation seasonality. For N. platensis, the variables used were diurnal range mean temperature of warmest quarter, precipitation of wettest month, and precipitation of warmest quarter.We estimated the continental surface area of South America in each generated map using ImageJ software (Abràmoff et al., Reference Abràmoff, Magelhães and Ram2004), considering sea level variations from the LGM to the middle Holocene and their influence on the South American continental-plain surface area (Stevaux, Reference Stevaux2000). In addition, the suitable areas (SA) on each map were estimated with the Analyze Particles tool at a spatial scale of 1000 km. This method helps understand the changes in maximum suitability areas during the LGM and middle Holocene because it considers the differences in surface area of the continental region of South America in the LGM and middle Holocene due to eustatic movement.
RESULTS
The review of distribution data identified 198 localities for E. neogeus and 142 for N. platensis, all within the LGM interval. Equus neogeus was recorded from northern Venezuela, Colombia, Ecuador, Peru, Bolivia, Chile, Brazil, Uruguay, and the Pampean region of Argentina, whereas N. platensis was recorded across all South American countries except for Guiana, French Guiana, and Suriname (Fig. 1). The co-occurrence of the two species during the late Pleistocene in some localities resulted in similar patterns of geographic distribution, including some overlapping areas (Fig. 1). No SA were predicted below 42°S during the LGM, lower than the Huincul Fault zone (Figs. 2, 3). Data provided by Maxent indicate that predictor variables, average precipitation, and precipitation of wettest quarter increasingly contributed to the probability of occurrence of E. neogeus and N. platensis (see Supplemental Data).
During the LGM, the potential geographic distribution for E. neogeus and N. platensis (Figs. 2, 3) predicted that most of the SA were in the Pampean region, southern and northeastern Brazil, Uruguay, southern Bolivia, northern Paraguay, northern and south-central Chile, Ecuador, Colombia, and Venezuela. The estimated potential distributions also showed suitability for both species in the Pampean region, southern and northeastern Brazil, and some small warmer areas (highly SA) throughout the Andes. However, during the middle Holocene, SA remained in the same regions where they were in the LGM, but contracted, with a southern shift in suitable habitats post-LGM.
DISCUSSION
The ecological niches of E. neogeus and N. platensis are inferred from palynological and paleodietary data, SDMs (such as this study), and comparisons with climatic models. The two species studied here are open and dry environment dwellers based on the superimposition of the patterns of each SA and the South American biomes reconstructions for the LGM (Ray and Adams, Reference Ray and Adams2001; Behling, Reference Behling2002; Cook and Vizy, Reference Cook and Vizy2006). During the LGM, E. neogeus often occurred in open and arid areas that were dominated by tropical and temperate semi-desert biomes (Ray and Adams, Reference Ray and Adams2001; Behling, Reference Behling2002; Cook and Vizy, Reference Cook and Vizy2006). Sánchez et al. (Reference Sánchez, Prado and Alberdi2006) suggested that late Pleistocene Equus populations might have restricted their use of habitats due to their exclusively grazing diet, which is a more selective feeding habit. Modeling of the potential distribution of N. platensis during the LGM suggested that this species occurred in areas with lower temperature and high aridity, such as tropical and temperate semi-deserts.
The Pampean region of Argentina during the LGM was predominantly a temperate semi-desert with a small area of temperate desert, characterized by low and dense vegetation, dry climate, low precipitation, and large daily temperature fluctuations (Iriondo and Garcia, Reference Iriondo and Garcia1993; Ray and Adams, Reference Ray and Adams2001). During the LGM, southern Brazil had a drier and cooler climate than today. It was a period marked by the contraction of tropical and araucaria forests and the great expansion of grasslands (Behling, Reference Behling2002), and, as Ray and Adams (Reference Ray and Adams2001) demonstrated, tropical semi-desert and tropical grasslands were the dominant landscapes in northeastern Brazil. The tropical semi-desert Caatinga biome is an open landscape with 2–10% vegetation cover of scattered grasses and small shrubs. Tropical grasslands are characterized by >20% vegetation cover, mainly grasses, with <5% woody plants (Ray and Adams, Reference Ray and Adams2001). Based on the palynological record, Mayle et al. (Reference Mayle and Power2008) classified the vegetation cover of the region as a mix of Caatinga, xerophytic plants, and grass field. The xerophytic vegetation described in palynological studies suggests the predominance of an arid climate during the LGM. During the Pleistocene/Holocene transition in northern Chile, González-Guarda et al. (Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018) suggested the predominance of a woodland-mesic environment C3 grassland with a significant presence of trees. For the south-central Chile region, González-Guarda et al. (Reference González-Guarda, Petermann-Pichincura, Tornero, Domingo, Agustí, Pino and Abarzúa2018) suggested a temperate environment, with a cold and humid climate, dominated by C3 plants and with a significant woodland cover.
The middle Holocene projection revealed a noticeable reduction in SA for both species, approximately 19% and 16% for E. neogeus and N. platensis, respectively. The fact that E. neogeus was more commonly found in wide-open habitats than N. platensis might explain the larger reduction in its SA because those were the areas more affected by climate change. Additionally, although N. platensis is more frequently found in open areas, it also occurred in dense forests (Campbell et al., Reference Campbell, Frailey and Romero Pittman2000; Codrea and Diaconu, Reference Codrea and Diaconu2007; Mothé et al., Reference Mothé, Avilla, Asevedo, Borges-Silva, Rosas, Labarca-Encina and Souberlich2017a). Thus, N. platensis might have occupied both habitats. Despite the clear reduction in the potential SA for both species, there was a slight increase in the SA in the southernmost part of South America during the middle Holocene. This change might be related to the climatic/environmental variations that occurred at the Pleistocene/Holocene transition, when both species likely searched for colder environments.
The middle Holocene projections showed that the SA were associated with a climate that was more humid than in the LGM (and slightly drier than today), except for northeastern Brazil, where humidity was estimated to be 10% higher during the LGM (Valdes, Reference Valdes2000; Melo and Marengo, Reference Melo and Marengo2008; Silva Dias et al., Reference Silva Dias, Turcq, Silva Dias, Braconnot, Jorgetti, Vimeux, Sylvestre and Khodri2009; Prado et al., Reference Prado, Wainer, Chiessi, Ledru and Turcq2013). Palynological data suggested that the Argentine Pampean region climate was more humid during the middle Holocene (Prieto, Reference Prieto1996). In southern Brazil, the climatic model generated by Silva Dias et al. (Reference Silva Dias, Turcq, Silva Dias, Braconnot, Jorgetti, Vimeux, Sylvestre and Khodri2009) indicated the middle Holocene climate was slightly more humid with vegetation cover of mixed humid forest/tropical forest. In their climate models, Melo and Marengo (Reference Melo and Marengo2008) suggested the occurrence of a warming signal during almost the entire year. According to Behling (Reference Behling1995), the palynological record indicated a reduced forest cover, suggestive of a drier climate. In northeastern Brazil, palynological data indicated that gallery forests were more developed during the middle Holocene (de Oliveira, Reference de Oliveira, Barreto and Suguio1999), which also suggests a wetter climate (Melo and Marengo, Reference Melo and Marengo2008; Silva Dias et al., Reference Silva Dias, Turcq, Silva Dias, Braconnot, Jorgetti, Vimeux, Sylvestre and Khodri2009; Prado et al., Reference Prado, Wainer, Chiessi, Ledru and Turcq2013). The vegetation models by Silva Dias et al. (Reference Silva Dias, Turcq, Silva Dias, Braconnot, Jorgetti, Vimeux, Sylvestre and Khodri2009) indicated that the vegetation of northeastern Brazil was characterized mainly by dry tropical forests/savannas, and a small area covered by seasonal tropical forests.
According to the fossil record and our SDM/PGD (species distribution/potential geographic distribution) models, E. neogeus and N. platensis had large areas of suitable habitat in arid environments with open vegetation during the LGM (Roig-Juñent et al., Reference Roig-Juñent, Domínguez, Flores and Mattoni2006). In both models, these species did not predominantly occupy dense forest environments such as the Amazon. Also, the SA for E. neogeus and N. platensis were more continuous and showed maximum suitability in the regions surrounding the La Plata River, the Pampean region in Argentina, and southern Brazil (Ramos et al., Reference Ramos, Riccardi and Rolleri2004). Moreover, no SA were predicted below 42°S during the LGM. The southern limits of the PGD for both species reach the Huincul Fault zone at the northern limit of Patagonia, which corresponds to the transition from temperate to cold climates. Notiomastodon platensis and E. neogeus seem to have been more associated with tropical and temperate climates, avoiding the colder and more arid environments of southern South America. Climate change caused the remaining glaciers to retreat to the Andes, expanding available areas with favorable weather conditions, such as Patagonia, and enabling both E. neogeus and N. platensis to extend their geographic distribution. However, this process coincided with the timing of their extinction and prevented the two species from extending their distribution in search of cooler and arid habitats (Villavicencio et al., Reference Villavicencio, Corcoran and Marquet2019).
Our results also support the argument that megamammals occupied exposed continental shelf areas during the LGM. The reduction in SA due to the rise in sea level possibly influenced the distribution of suitable areas for N. platensis and E. neogeus at the end of the Pleistocene and the beginning of the Holocene.
The SDMs of both species suggested that climatic/environmental variations played a key role in the decreases in their SA during the Pleistocene/Holocene transition, but would not have been enough to drive them to extinction (Hanski, Reference Hanski2005). Although these decreases appear to be small, the negative effect of any habitat loss is higher for herbivore mammals with large body mass (body mass estimation for E. neogeus is 250 kg, Fariña et al., Reference Fariña, Czerwonogora and Di Giacomo2014; and for N. platensis is 4,500 kg, Larramendi, Reference Larramendi2015), because these losses tend to reduce resource availability (Lino et al., Reference Lino, Fonseca, Rojas, Fischer and Pereira2019), consequently leading to extinction. Additionally, considering the largest SA, the climatic and environmental/biome changes may only have affected the population densities of E. neogeus and N. platensis, whereas the qualitative changes in the biomes may have had a greater effect on these species. The climate in the larger SA changed from cold and dry during the LGM to warmer and more humid during the middle Holocene.
Although the geographic distributions of E. neogeus and N. platensis were similar, their potential distributions were different. In addition to climate change, there was a rapid and considerable modification in the vegetation of the SA, with a rise in dense forests favored by the increases in humidity and temperature during the LGM, and a reduction in open plain habitats, as shown in the survey of past vegetation in this study. Because it had a grazer feeding habit, E. neogeus may have been strongly affected by these vegetation changes due to the reduction of pastures (Sánchez et al., Reference Sánchez, Prado and Alberdi2004; MacFadden, Reference MacFadden2005; Semprebon et al., Reference Semprebon, Rivals, Solounias and Hulbert2016). Sánchez et al. (Reference Sánchez, Prado and Alberdi2006) suggested that late Pleistocene populations of E. neogeus might have adapted to a more selective diet, which restricted its occurrence to certain habitats. This change in diet may explain the greater post-LGM reduction in SA for E. neogeus. Conversely, N. platensis had an opportunistic/generalist diet, feeding on grasses (Mothé et al., Reference Mothé, Avilla, Asevedo, Borges-Silva, Rosas, Labarca-Encina and Souberlich2017a), branches, and other types of vegetation, and thus more access to a wide variety of food resources, which may explain the highly subtle reduction in its SA.
Thus, the changes in climatic conditions modified the habitats of E. neogeus and N. platensis and played a significant role in their extinction. Although climatic changes were relevant to both megamammals’ extinction, we cannot discard a more complex scenario wherein other causes (i.e., human impact, pathologies, etc.), working in synergy, could also explain megafauna extinction in South America.
CONCLUSIONS
Species distribution models revealed similar distribution patterns for E. neogeus and N. platensis during the last glacial maximum, but their ecological niches differed mainly in their diet composition, with E. neogeus being a grazer and N. platensis being an opportunistic/generalist mixed feeder. During the LGM, both species occupied mainly arid and semi-arid open habitats in South America. The total areas suitable for survival decreased after the LGM, suggesting that climatic/environmental changes would have affected E. neogeus and N. platensis populations during the late Pleistocene. Subsequently, the climate in the suitable areas for both species changed from drier and colder in the LGM to more humid and warmer in the middle Holocene. Thus, climatic and environmental changes in South America were an underlying factor that drove these megamammals to extinction.
Supplementary Material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2021.24
Acknowledgments
The authors are grateful to all curators of the paleontological collections where the fossil specimens are housed for allowing the access that supported this study. The authors are also thankful to Claudia Russo and the Programa de Pós-graduação em Biodiversidade e Biologia Evolutiva (PPGBBE/UFRJ) for support with the manuscript review. We deeply thank senior editor Derek Booth and the two reviewers Laura Jiménez and Jeanette Pirlo for their comments and suggestions. We are especially grateful to associate editor Jeff Pigati for the perseverance and precious comments during the review of this manuscript.
Financial Support
This work was funded by Conselho Nacional de Desenvolvimento Científico e Tecnológico (140453/2012-01, 201081/2014-8, PDJ; 153536/2016-0, DM; 248772/2013-9, LSA), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES/PNPD/2016; DM), UNIRIO IC (TA) and Fundação Carlos Chagas de Amparo a Pesquisa do Estado do Rio de Janeiro (E-26/100.246/2014, 202.375/2018 and 202.376/2018, DM; 204036- E_25/2014-Jovem Cientista do Nosso Estado, LSA). HM is thankful to the Department of Earth Sciences, University of Oregon, for financial support, and to the National Science Foundation (NSF) for financial support as a Graduate Employee provided through the project Collaborative Research: ABI Innovation: FuTRES, an Ontology-Based Functional Trait Resource for Paleo- and Neo-biologists (1759821).