Influence of late Quaternary climate on the biogeography of Neotropical aquatic species as reflected by non-marine ostracodes

We evaluated how ranges of four endemic and non-endemic aquatic ostracode species changed in response to long-term (glacial–interglacial cycles) and abrupt climate fluctuations during the last 155 kyr in the northern Neotropical region. We employed two complementary approaches, fossil records and species distribution models (SDMs). Fossil assemblages were obtained from sediment cores PI-1, PI-2, PI-6 and Petén-Itzá 22-VIII-99 from the Petén Itzá Scientific Drilling Project, Lake Petén Itzá, Guatemala. To obtain a spatially resolved pattern of (past) species distribution, a downscaling cascade is employed. SDMs were reconstructed for the last interglacial (∼120 ka), the last glacial maximum (∼22 ka) and the middle Holocene (∼6 ka). During glacial and interglacial cycles and marine isotope stages (MISs), modelled paleo-distributions and paleo-records show the nearly continuous presence of endemic and non-endemic species in the region, suggesting negligible effects of long-term climate variations on aquatic niche stability. During periods of abrupt ecological disruption such as Heinrich Stadial 1 (HS1), endemic species were resilient, remaining within their current areas of distribution. Non-endemic species, however, proved to be more sensitive. Modelled paleo-distributions suggest that the geographic range of non-endemic species changed, moving southward into Central America. Due to the uncertainties involved in the downscaling from the global numerical to the highly resolved regional geospatial statistical modelling, results can be seen as a benchmark for future studies using similar approaches. Given relatively moderate temperature decreases in Lake Petén Itzá waters (∼5 ^∘C) and the persistence of some aquatic ecosystems even during periods of severe drying in HS1, our data suggest (1) the existence of micro-refugia and/or (2) continuous interaction between central metapopulations and surrounding populations, enabling aquatic taxa to survive climate fluctuations in the northern Neotropical region.

1 Introduction

Climate changes are quasi-cyclical natural processes that continuously influence ecosystem dynamics and shape biological diversity worldwide (Blois et al., 2013; Parmesan and Yohe, 2003; Yasuhara et al., 2009, 2017). On inland ecosystems, late Quaternary climate fluctuations such as glacial–interglacial cycles are recognized as the main drivers responsible for past species extinctions (Martínez-Meyer et al., 2004; Nogués-Bravo et al., 2008), speciation events (Peterson and Nyári, 2008; Solomon et al., 2008), delimitation of refugia (Hugall et al., 2002; Peterson et al., 2004) and development of migration pathways (Ruegg et al., 2006; Waltari and Guralnick, 2009) for both plants and animals.

In the northern Neotropics, which include southern Mexico, Central America and the Antilles, late Quaternary climate inferences based on climatic simulations with global climate models (GCMs; Hijmans et al., 2005) and reconstructions from marine and lacustrine sedimentary sequences (Hodell et al., 2008; Pérez et al., 2011; Escobar et al., 2012) have revealed climate fluctuations related to temperature and precipitation, especially during transitions between glacial and interglacial episodes and during climate pulses such as the last glacial maximum (LGM) and Heinrich stadials (HSs; Correa-Metrio et al., 2012b). In the Neotropics, controls of climate fluctuations are related to orbital forcing and internal component variations, such as the position (north–south) of the intertropical convergence zone (ITCZ), strength of Atlantic meridional overturning circulation (AMOC) and changes in Caribbean surface water temperature (Cohuo et al., 2018). Alterations in these features have produced temperature decreases in a range of 3–5 ^∘C, although some estimations suggest decreases up to 10 ^∘C relative to present and large reductions in precipitation, particularly during HS, when most lakes in the region dried completely (Cohuo et al., 2018). Correa-Metrio et al. (2014) found evidence for rapid climate change in terrestrial environments during HSs, which was associated with major ecological and biological shifts (Loarie et al., 2009; Burrows et al., 2011; Sandel et al., 2011). Correa-Metrio et al. (2012a, b, 2014) found that plant survival in the northern Neotropical region during HSs required migrations to refugia. The climatically driven pace and magnitude of changes in aquatic environments can, however, vary considerably relative to effects in terrestrial environments (Sandel et al., 2011; Litsios et al., 2012; Bonetti and Wiens, 2014). It therefore remains uncertain how aquatic species responded to past climate alterations.

To evaluate past biogeographic dynamics of northern Neotropical inland aquatic species, we used freshwater ostracodes (bivalved microcrustaceans) as a model group (Mesquita-Joanes et al., 2012) and two complementary approaches: (1) fossil records (Dawson et al., 2011; McGuire and Davis, 2013) and (2) species distribution models (SDMs; Elith and Leathwick, 2009; Nogués-Bravo, 2009; Veloz et al., 2012; Maguire et al., 2015).

Ostracodes were selected because they have possessed one of the best fossil records in the region since the late Quaternary (Pérez et al., 2011, 2013) and have demonstrated sensitivity to climatic variation (modern and past) in both terrestrial (Horne et al., 2002) and marine environments (Yasuhara et al., 2008, 2014). Given their intermediate role in trophic chains (Valtierra-Vega and Schmitter-Soto, 2000; Bergmann and Motta, 2005; Cohuo et al., 2016), changes in their abundances and assemblage composition can also reflect changes in primary production and higher trophic levels (Rodriguez-Lazaro and Ruiz-Muñoz, 2012). Paleo-records provide true evidence for the presence of a species in the past at resolutions ranging from decadal to millennial scales, but in the absence of a denser spatial network, this approach is usually limited to the local scale (Maguire and Stigall, 2009; Dawson et al., 2011). Species distribution models are based on the combination of georeferenced species occurrences with environmental information to characterize the range of climate tolerance that a species inhabits (Guisan and Thuiller, 2005; Maguire et al., 2015). By using multiple time periods, species occurrences across different climatic scenarios can be projected to a certain degree (Elith and Leathwick, 2009; Svenning et al., 2011).

The most important limitations and uncertainties of SDMs are the relevant forcing data such as GCMs and the statistical algorithms employed. For instance, simulations of tropical Atlantic climates remain deficient in many climate models due to incomplete characterization of the vertical structure of tropospheric water vapour and humidity. As a consequence, the simulation of temperature and precipitation gradients is afflicted with a high degree of uncertainty in GCMs, especially across regions with irregular and complex topography (Solomon et al., 2010). Statistical algorithms and data parametrization also add another level of uncertainty in the downscaling cascade, including the structure of past surface fields such as topography, vegetation structure and coastline. Moreover, the usage of statistical algorithms for the geospatial mapping also includes uncertainties that are implicitly included in the results (Chen et al., 2010; Neelin et al., 2010).

The combination of paleo-records and SDMs provides a unique opportunity to obtain quantitatively and potentially high-resolution reconstructions of past species dynamics at the local and regional scale during past climate fluctuations in the northern Neotropical region.

In this study, we addressed three overarching questions. (1) Have past climate changes since 155 ka (Hodell et al., 2008; Correa-Metrio et al., 2012a, b, 2014; Cohuo et al., 2018) had profound consequences for aquatic ecosystem stability in the northern Neotropics? (2) Did endemic and non-endemic (widespread) species respond in the same way to climate shifts? (3) Did refugia exist, and if so, what was their spatial distribution?

2 Methods

2.1 Study area and sampling of modern species

Our study area is the northernmost northern Neotropics, an area that extends from southern Mexico to Nicaragua (Fig. 1). We sampled 205 aquatic ecosystems during 2010–2013, including cenotes (sinkholes), lakes, lagoons, crater lakes, maars, permanent and ephemeral ponds, wetlands, and flooded caves. Sampled systems are located at elevations from ∼10 to ∼4000 m a.s.l., and conductivity ranged from 0.1 to 3500 µS cm⁻¹. Most aquatic systems were shallow, with a mean depth <10 m, except for large lakes such as Petén Itzá, Atitlán, Coatepeque, Ilopango, Lachuá, crater and maar lakes, and cenotes, which are mostly >15 m deep. Biological samples were collected at three different sections of the systems: littoral, water column and deepest bottom. At littoral areas, we sampled in between submerged vegetation using a hand net of 250 µm open mesh. The water column was sampled while performing vertical tows and horizontal trawls with a net with a 20 cm wide mouth and 150 µm mesh size. Sediment samples were taken from the deepest part of the systems with an Ekman grab, but only the uppermost centimetres of each grab were used for further analysis. Ostracodes were sorted in the laboratory using a Leica Z4 stereomicroscope, and dissections were carried out in 3 % glycerin. Shells were mounted on micropaleontological slides. Dissected appendages were mounted in Hydromatrix^® mounting media. Taxonomic identification followed Karanovic (2012) and Cohuo et al. (2016). Four ostracode species were selected for this study: Cypria petenensis Ferguson Jr. et al., 1964, Paracythereis opesta (Brehm, 1939), representing taxa endemic to the northern Neotropical region (Cohuo et al., 2016; Fig. 1a, b), Cytheridella ilosvayi Daday, 1905, and Darwinula stevensoni (Brady and Robertson, 1870), which are widely distributed (non-endemic) on the American continent (Fig. 1c, d).

Figure 1 Current ostracode species distributions and predicted distribution based on species niche modelling and two statistical evaluations: true skill statistics (TSS) and area under the receiver operating characteristic curve (AUC). (a) Cypria petenensis, (b) Paracythereis opesta, (c) Cytheridella ilosvayi and (d) Darwinula stevensoni.

2.2 Sediment cores from Lake Petén Itzá and regional paleo-records

Information about fossil occurrences of the target species was obtained from sediment cores retrieved from Lake Petén Itzá (northern Guatemala) by the Petén Itzá Scientific Drilling Project (PISDP). Cores PI-1, PI-2, PI-6 (Mueller et al., 2010) and Petén-Itzá 22-VIII-99 were used. Core chronologies were established independently by radiocarbon dating (Mueller et al., 2010), and for cores PI-1, PI-2 and PI-6, sediments older than 40 ka were dated by identification and correlation of tephra layers (tephrochronology; Kutterolf et al., 2016). The age model proposed by Kutterolf et al. (2016) was used. Correlation of cores was done using lithological markers; stratigraphic boundaries; similarity in magnetic susceptibility patterns; and ash layer correlation such as Congo tephra (53 ka BP), EFT tephra (50 ka BP) and Mixta tephra (∼39 ka BP). Core sampling was done at 20 cm intervals, which is ∼ 100–200 years of temporal resolution (Kutterolf et al., 2016). At sediment transitions and periods of interest such as the LGM and Heinrich stadials, samples were closely spaced at 1 cm, representing ∼ 6–10 years of temporal resolution (Kutterolf et al., 2016). All samples had a volume of 3 cm³ of wet sediments.

Ostracode separation methods and counting can be found in Cohuo et al. (2018). We looked at near-continuous ostracode fossil occurrences in the sediments over the last 155 kyr. There was, however, a gap in sediment availability during the period 83–53 ka. We also compiled fossil data for our target species from 19 other studies in the northern Neotropical region to obtain past spatial distributions of the target species (Supplement, Table S1). These studies were restricted to the LGM and middle Holocene.

Shells of the target species were measured and photographed using a Canon PowerShot A640 digital camera attached to a Zeiss Axiostar Plus light microscope. Abundances of the target species in each core were plotted using C2 software version 1.5 (Juggins, 2007).

2.3 Species niche modelling (SNM): modern projections and reconstruction of past distributions

We determined modern macro- and micro-ecological preferences for our target species using our dataset (multivariate approach) and the literature (Pérez et al., 2010). Given the ecological preferences of the species, we used seven environmental variables related to temperature and precipitation that show the lower Pearson correlation coefficient within 19 regional environmental variables (Supplement, Table S2) and are known to have the strongest relationships with ostracode distribution: (1) mean annual temperature, (2) mean diurnal temperature range, (3) isothermality (day-to-night temperature oscillation relative to summer and winter), (4) temperature seasonality, (5) annual temperature range, (6) total annual precipitation and (7) precipitation seasonality, all available from the WorldClim database version 1.4 (Hijmans et al., 2005; http://www.worldclim.org, last access: 20 May 2019). Variables of importance were analysed to identify those with the greatest influence on each ostracode species distribution.

Environmental conditions of the present corresponded to the interpolation of average monthly climate data from weather stations of various locations of the world and major climate databases such as the Global Historical Climatology Network (GHCN) and the Food and Agricultural Organization of the United Nations (FAO). Grids had a spatial resolution of 30 arcsec. Although modern climatic data are generated at very high resolution, one should note that modelling of tropical climate and circulation is still afflicted by a comparatively high degree of uncertainty, especially the realistic simulation of the hydrological cycle and precipitation. In this context, the purpose of the study is also to investigate the extent to which differences in profound background climatic changes during glacial–interglacial periods are responsible for lateral and/or vertical changes in ecological niches of the respective species.

Past species distributions were investigated using climate conditions inferred for three time periods: ∼120 ka (last interglacial), ∼22 ka (LGM) and ∼6 ka (middle Holocene). For environmental data corresponding to ∼120 ka (Otto-Bliesner et al., 2006), grids have a spatial resolution of 30 arcsec, which represents ∼1 km² in the northern Neotropical region. Environmental conditions at ∼22 and ∼6 ka were obtained from downscaled paleo-climatic simulations forced with the coarsely resolved output fields of two global circulation models (GCMs), the MIROC-ESM 2010 (Watanabe et al., 2011) and CCSM4 (Gent et al., 2011).

These GCMs were selected because they yield slightly varying temperature and differences in precipitation fields (Fig. 2). At ∼22 ka, the MIROC-ESM model shows colder and drier conditions in the region than the CCSM4 model (Fig. 2). At ∼6 ka, the CCSM4 model simulates slightly cooler and drier conditions than the MIROC-ESM model (Fig. 2). These differences enable assessment of model uncertainty with respect to global climate simulations.

Figure 2 Estimated mean annual temperature and mean annual precipitation values for ∼120, ∼22 and ∼6 ka and present. Estimates for ∼22 and ∼6 ka were based on general circulation models CCSM4 (grey line) and MIROC-ESM (black line).

The target grids at the lower end of the downscaling cascade have a spatial resolution of 2.5 arcmin, which represents ∼5 km² in the study area. For all periods, grids with global information were trimmed to match the extent of our study area. The SDM toolbox (Brown, 2014), implemented in Arc GIS, was used for this purpose.

The modelling framework was constructed using five presence-based and absence-based algorithms because of true species absences in our database. We used the generalized linear model (GLM; McCullagh and Nelder, 1989), the generalized additive model (GAM; Hastie and Tibshirani, 1990), the generalized boosting model (GBM; Ridgeway, 1999), maximum entropy (MAXENT; Tsuruoka, 2006) and the surface range envelope (SRE; Busby, 1991). The first three algorithms, GLM, GAM and GBM, are regression-based models, which are flexible in handling a variety of data response types (linear and non-linear) and are less susceptible to overfitting than other algorithms such as multivariate adaptive regression splines (MARSs; Guisan et al., 2002; Franklin, 2010). MAXENT is a general-purpose machine-learning method which predicts a species probability occurrence by finding the distribution closest to uniformity (maximum entropy); it requires previous knowledge of the environmental conditions at known occurrence localities (Elith et al., 2011). The SRE algorithm is an envelope-type method that uses the environmental conditions of locations of occurrence data to profile the environments where a species can be found (Araújo and Peterson, 2012). All these modelling techniques are, to a different degree, limited by several numerical factors, such as missing values, outliers, sampling size, overfitting and interaction between predictors. Special attention therefore must be paid to producing reliable models which maximize the agreement of the predicted species occurrences with the observed data (Guisan et al., 2002; Franklin, 2010). In most cases the combination of methods (e.g. GLM and GAM) is recommended to assess the robustness of according results of individual models (Guisan et al., 2002).

For our study, settings for all modelled techniques, such as the number of trees, number of permutations, iteration depths, Bernoulli distribution normalization and node size, follow Georges and Thuiller (2013). Records were split randomly into a training (calibration; 70 %) and a test (validation; 30 %) dataset, with 10 replications for each model type. A total of 50 models (5 algorithms and 10 replications) were generated for each ostracode species and time period. All projections were evaluated using three statistical approaches to reduce uncertainty in species niche models: (1) the true skill statistics (TSS), (2) the area under the receiver operating characteristic curve (AUC) and (3) Cohen's kappa statistics (Thuiller et al., 2009, 2015). For all algorithms, best-fit model runs above critical values (TSS values >0.4, AUC > 0.7 and kappa > 0.4) were used to construct consensus maps for each modelling technique. Final maps were constructed using an ensemble of all techniques. The combination of methods reduces the effect of inter-model variability and uncertainties that arise from using single algorithms (Araújo and New, 2007; Marmion et al., 2009; Thuiller et al., 2009). The final distribution maps thus indicate areas simulated by most modelling techniques. All calculations were done using the “biomod2” v.3.1-64 package (Thuiller et al., 2015), implemented in R v.3.2.1 software (R Development Core Team, 2015).

3 Results

3.1 Northern Neotropical paleo-records, species permanence and displacement

Records of the period corresponding to the last interglacial (130–115 ka) were obtained from core PI-7 (155–83 ka). Abundances of our four target species were generally low, with <60 adult shells g⁻¹, and frequencies (relative abundances) varied considerably (Fig. 3). The endemic C. petenensis was the most frequent species (Fig. 3). Paracythereis opesta and C. ilosvayi, which are bottom-dwelling organisms, were recovered only from sediments deposited in ca. 87–85 ka, where high abundances of C. petenensis were observed (Fig. 3). Darwinula stevensoni showed a sole occurrence at ∼ 155–153 ka.

Figure 3 Fossil record of the period 155–83 ka and species niche modelling results for ∼120 ka (representing last interglacial climate) for four ostracode species: Cypria petenensis, Paracythereis opesta, Cytheridella ilosvayi and Darwinula stevensoni. (a) Fossil record of four ostracode species from core PI-1 in Lake Petén Itzá, and (b) maps from niche modelling, showing the probability of species distributions.

Records of the last glacial and deglacial periods were obtained from Lake Petén Itzá core PI-2 (Fig. 4a) and published data from core PI-6 (Fig. 4b; Pérez et al., 2011). Pérez et al. (2011) found nearly continuous presence of endemic species in core PI-6 during the interval 24–10 ka. Gaps of millennial duration are, however, evident for the periods 24–22 and 13–10.5 ka. The record from PI-2 shows a complementary pattern to that of PI-6 because species presence in PI-2 coincided with species absence in core PI-6. Cypria petenensis in the PI-2 record, for example, shows high abundances at the onset of the LGM (23–21 ka), and P. opesta displays high abundances around 22 and 19 ka (Fig. 4a). Thus, the two records suggest the continuous presence of endemic species in Lake Petén Itzá during both the LGM and deglacial periods.

Figure 4 Fossil record of the period 53–10 and species niche modelling results for the last ∼22 kyr (representing last glacial maximum climate) for four ostracode species: Cypria petenensis, Paracythereis opesta, Cytheridella ilosvayi and Darwinula stevensoni. Fossil ostracode record from Lake Petén Itzá. (a) Core PI-2 for the period 53–14 ka, (b) core PI-6 for the period 24–10 ka (taken from Pérez et al., 2011) and (c) map showing the probability of species distributions based on the CCSM4 climate model.

Non-endemic species show intermittent distributions in both the PI-2 and PI-6 cores (Fig. 4a, b). Darwinula stevensoni was recorded exclusively at ca. 23, 22–20 and 19–18 ka. Similarly, Cytheridella ilosvayi was present in very low abundances during two short episodes at about 20 and 14 ka. We recorded low abundances of both species during the LGM (<250 adult shells g⁻¹), compared to periods immediately before and after, when temperatures are thought to have been warmer. For example, during the deglacial period, abundances were always >250 adult shells g⁻¹.

Fossil records from the middle Holocene were obtained from core Petén-Itzá 22-VIII-99 and 11 regional studies (Fig. 5a). The record from core Petén-Itzá 22-VIII-99, retrieved from 11.5 m water depth, shows that endemic species were present continuously during the last 6.5 kyr (Fig. 5a). Most regional records came from cenotes and lakes on the Yucatán Peninsula (Supplement, Table S1). All fossil records show that endemic species were spatially distributed throughout the current ranges of extant populations (Fig. 5b).

Figure 5 Fossil record of the last 14 kyr and species niche modelling results for the ∼6 ka (representing middle Holocene climate) for four ostracode species: Cypria petenensis, Paracythereis opesta, Cytheridella ilosvayi and Darwinula stevensoni. (a) Ostracode fossil record from core Petén Itzá 22-VIII-99. (b) Map showing the probability of suitable species distribution based on the CCSM4 climate model. Numbers in maps represent regional fossil records. Numbers correspond to those in the Supplement (Table S1).

For non-endemic species, regional fossil records from the middle Holocene revealed their presence in areas ranging from the northern Yucatán Peninsula to northern Guatemala and Belize (Supplement, Table S1). Core Petén-Itzá 22-VIII-99 highlights an almost continuous presence of C. ilosvayi in the lake, characterized by high abundances, except for the period 11–8.5 ka, when the species was absent (Fig. 5a). Darwinula stevensoni was present continuously during the last 9 kyr, but in the lower section of the core, dated to 14–10 ka, the species was absent (Fig. 5a).

3.2 Species niche modelling: distribution hindcasting for time slices ∼120, ∼22 and ∼6 ka

For the 205 aquatic ecosystems sampled, 145 had at least one of the target species present: C. petenensis, P. opesta, C. ilosvayi and D. stevensoni. Forty-nine systems contained C. petenensis, 37 had P. opesta, 79 were inhabited by C. ilosvayi, and 61 contained D. stevensoni. Analysis of variables of importance showed that environmental variables with the greatest influence on species distribution are precipitation seasonality and mean annual temperature (Table 1). For individual species, however, variables received different scores, indicating that each species' optimal climate niche is controlled by a particular combination of variables (Table 1). Diagnostic tests of the reconstructions (TSS, AUC and kappa) show good performance for all algorithms and periods evaluated (Table 1). There were, however, differences in predictive accuracy within species. Modelled distributions of endemic species have the highest evaluation scores (AUC = 0.8, TSS = 0.49 and kappa = 0.45). Non-endemic species models (AUC = 0.75, TSS = 0.46 and kappa = 46) have slightly lower values but also fall within the acceptable range.

Table 1 Ostracode species niche modelling, input data and evaluation scores. Variables of importance (mean of 10 evaluation runs) and evaluation model performances based on true skill statistics (TSS) and area under the receiver operating characteristic curve (AUC). Variable importance scores ≥0.30 are shown in bold.

Abbreviations are as follows: annual mean temperature (BIO 1), mean diurnal range (BIO 2), isothermality (BIO 3), temperature seasonality (BIO 4), temperature annual range (BIO 7), annual precipitation (BIO 12) and precipitation seasonality (BIO 15).

Reconstructions for the period ∼120 ka suggest very broad distributions of endemic taxa, as the climate enabled the species to expand their ranges. Probability values, however, were relatively low (<80 %; Fig. 3b). For the non-endemic species, reconstructions for ∼120 ka show different areas of climatic suitability, with species presence probabilities reaching 60 %. Zones of higher probability (>80 %) are dispersed throughout the region. The most extensive zones of species distribution suitability are located along the Caribbean coast of the Yucatán Peninsula and in northern Guatemala (Fig. 3b).

Inferences for endemic taxa distributions at ∼22 ka, based on the CCSM4 model, suggest that these species remained in the core area but that they may have been displaced somewhat to the northern portion of the Yucatán Peninsula (Fig. 4c). This estimate has probability values of >75 %. The MIROC-ESM model suggests areas of distribution similar to those presented by the CCSM4 model but slightly more restricted areas for C. petenensis and more widespread areas for P. opesta. Probability values were low in this model (<65 %; Supplement, Fig. S1). Models for non-endemic species reveal fragmented and discontinuous distributions (Fig. 4c). At ∼22 ka, corresponding to the LGM, both the CCSM4 and MIROC-ESM models suggest that non-endemics moved northward on the Yucatán Peninsula to the Gulf of Mexico (>65 % probability) and/or were displaced southward to Central America (85 % probability; Fig. 4c; Supplement, Fig. S1).

For ∼6 ka, the CCSM4 model suggests discontinuous areas of distribution on the Yucatán Peninsula (Fig. 5b) for endemic species, whereas the MIROC-ESM shows more continuous distributions, particularly along the eastern portion of the Peninsula (Supplement, Fig. S1). For non-endemic species, the CCSM4 and MIROC-ESM models show very similar patterns. Extensive regions of climatic suitability were identified for C. ilosvayi, but those with higher probability are located along the Caribbean Coast (Fig. 5b; Supplement, Fig. S1). For D. stevensoni, areas of maximum probability are discontinuous. The maximum probability was found in isolated regions such as the southern part of the northern Yucatán Peninsula, Belize and eastern Honduras (Fig. 5b).

4 Discussion

4.1 Congruence between paleo-records and modelled paleo-distributions of freshwater ostracodes in the northern Neotropical region

Our study highlights the fact that accuracy and congruence between paleo-records and modelled paleo-distributions of freshwater ostracodes in the northern Neotropical region were influenced by multiple factors such as the climate model used, modelling algorithm employed, sediment core characteristics and target species.

For instance, distribution models and the modelling cascade were characterized by a high degree of uncertainty with regard to precipitation and temperature estimations of climate models (GCMs). This limited the full estimation of spatial distribution of target species, especially during older periods such as the LIG and LGM, where fossil evidence (spatial and temporal) was scarce.

The simulation of precipitation of GCMs is afflicted with high degrees of uncertainties because the vertical structure of stratospheric water vapour and humidity profile have large biases, especially in the tropics (Gettelman et al., 2010). This implies that GCMs commonly reproduce a large-scale pattern of precipitation with high confidence, but models tend to underestimate the magnitude of precipitation change at the regional or local scale (Stephens et al., 2010). Similarly, GCM temperature estimations in the tropics may display large biases because changes in climate drivers of the continental temperature of the northern Neotropics such as Atlantic sea surface temperature and the Atlantic warm pool are usually underestimated (Liu et al., 2013). Simulations of temperature variations during the LGM, for example, tend to overestimate cooling in tropical regions (Kageyama et al., 2006; Otto-Bliesner et al., 2009).

In our study, reconstructed maps based on MIROC-ESM and CCSM4 models simulate slightly different areas of distribution for the target species. This is associated to differences in precipitation and temperature estimations between models. The most important difference between their respective reconstructions pertains to the extent of suitable areas of distribution of the species, being generally broader in the MIROC-ESM model than in the CCSM4 model.

The scarcity of fossil records also limited the full reconstruction of distribution dynamics of species, especially during the LIG and LGM, because records were obtained only from Lake Petén Itzá and were relatively scarce. The period 24–14 ka was highly informative because the comparisons between cores PI-2 and PI-6, and specifically the compensation effect between them (the presence of species in a core in periods were absences were determined in the other), highlight the fact that gaps in the fossil record may be related to core location in the lake, shell preservation and individual species ecology and not only by species absence. This therefore suggests that short gaps, lasting less than 10 ka, cannot be considered evidence for species absence.

In general, the comparison between species distribution models and paleo-records shows a quite high degree of similarity. This is especially evident for the middle Holocene, as the individual species distribution models (SDMs) output of the target species was compared with the fossil records at the regional scale. In all cases, SDM reconstructions show distributional areas where fossil records were recovered. This congruence may be supported by the agreement between estimations of temperature in climate models and paleo-records.

4.2 Endemic and non-endemic species responses during long-term climatic fluctuations: glacial and interglacial cycles and marine isotope stages

Paleo-climate inferences derived from Lake Petén Itzá sediments suggest that glacial–interglacial cycles in the northern Neotropical region did not have profound consequences with respect to the spatial distribution of isotherms in terrestrial environments (Hodell et al., 2008; Pérez et al., 2011, 2013; Escobar et al., 2012). Most paleo-climate studies in the region based in different proxies such as ostracods, pollen and δ¹⁸O fluid inclusion data from speleothems suggest that temperatures during the last glacial period and start of the deglacial period may have been up to 5 ^∘C lower than today (Correa-Metrio et al., 2012a, b; Arienzo et al., 2015; Cohuo et al., 2018), although some estimations suggest a temperature depression of about 10 ^∘C compared with modern records (Hodell et al., 2012; Grauel et al., 2016). Precipitation was affected more profoundly, but not consistently, during glacial–interglacial cycles and likely fluctuated in response to changes in local atmospheric circulation. For instance, the position of the Hadley cell and ITCZ, together with climate forcing, such as Heinrich stadials, seems to drive precipitation fluctuation locally. During the LGM, for example, humid conditions has been estimated to the region (Bush et al., 2009; Cohuo et al., 2018).

Our results, however, suggest that temperature fluctuations affected aquatic species associations to a higher degree compared to reductions in precipitation (changes in lake water chemistry) because the presence and absence of species and fluctuations in total abundances match periods of temperature change rather than times of lake level shifts.

Endemic and non-endemic species responded similarly to glacial and interglacial cycles and transitions. Fossil records from Lake Petén Itzá sediment cores PI-1, PI-2, PI-6 and Petén-Itzá 22-VIII-99 reveal that endemic species were almost continuously present during the last 155 kyr. Short gaps, lasting less than 10 ka, were not considered evidence for species absence.

Non-endemic species show patterns of expansion and contraction that track temperature fluctuations. Modelled paleo-distributions and paleo-records show that distributions of non-endemic species were widespread during the LIG and fragmented during the middle Holocene, when climates were warmer. During the last glacial, non-endemic species were absent or sporadically present. This may result in response to lower temperatures that characterized the last glacial. Modelled paleo-distributions for the LGM also show that non-endemic species were displaced from their current ranges toward the northern Yucatán Peninsula and/or southward toward Central America, where a warm climate likely persisted. This scenario suggests migrations of regional magnitude, as species were lost from areas such as southern Mexico and northern Guatemala but persisted within their current range of distribution in fragmented populations, such as areas of southeastern Honduras and northeastern Nicaragua.

The presence of endemics and absence of non-endemic species during the LGM reveal a clear ecological signal, which may be associated to the degree of adaptation to ecological niches. Endemic species seem to be highly resilient to long-term natural disturbances, whereas non-endemic ones demonstrated higher sensitivity. There is increasing evidence that biological communities, particularly terrestrial taxa, display strong resilience in the face of natural and human disturbances in the northern Neotropical region. Hurricane impacts, widespread pre-Columbian agricultural activities and decadal-to-centennial climate changes are recognized as the main disrupters of Holocene ecosystem composition and function in the region. Such perturbations, however, did not severely and permanently alter plant associations such as moist forests (Bush and Colinvaux, 1994; Cole et al., 2014) and dry tropical forests (Van Bloem et al., 2006; Holm, 2017), which persisted in the region despite these disturbances. Plant taxa of Panama demonstrated a recovery time of just 350 years after strong deforestation by pre-Columbian agriculture (Bush and Colinvaux, 1994). Similarly, the rainforest in Guatemala recovered from Mayan alterations in a time span of 80–260 years (Mueller et al., 2010). Bird composition has also demonstrated rapid recovery time after hurricane impacts, species compositions affected in Central America and the Caribbean returned to pre-hurricane conditions in time periods ranging from months to years (Will, 1991; Wunderle et al., 1992; Johnson and Winker, 2010).

The continuous presence of both endemic and non-endemic (except during the LGM) ostracode species in the northern Neotropics during glacial–interglacial cycles also reflects the fact that aquatic ecosystem functionality was altered little during the last 155 kyr. High abundance of ostracodes, which belong to intermediate trophic levels, suggests high rates of primary production and ample food sources for higher consumers, especially during the LIG and middle Holocene. During the LGM, the presence of endemics and absence of non-endemics, along with lower total ostracode abundances, suggest moderate alteration of aquatic ecosystem dynamics. Reduced primary production and the loss of poorly adapted species might also be inferred for this period.

Marine isotope stages (MISs), which describe shorter periods of climate variability than glacial–interglacial cycles, were also used to evaluate the distribution dynamics of aquatic species. During MISs, ostracode composition remained relatively constant even across MIS boundaries (Fig. 6). Sediments from Lake Petén Itzá that correspond to warmer periods MIS 3 (57–29 ka) and MIS 1 (14 ka to present) were characterized by abundant fossils. MIS 2 (29–14 ka) shows lower species abundances (total adult and juvenile valves), likely related to persistent cold temperatures. The absence of Cytheridella ilosvayi during most of MIS 2 illustrates the sensitivity of non-endemics to cool temperatures (Fig. 6). Similar to glacial–interglacials in terrestrial environments, during MISs, northern Neotropical endemic species showed high resilience to changes between cold and warm phases, whereas non-endemic species proved to be more sensitive to cold periods, especially the LGM.

Figure 6 Master profile of the fossil ostracode record during marine isotope stages of the last 155 kyr in Lake Petén Itzá. Zone delimited by dashed lines represents a period of data absence. Grey peaks during the period of 24–10 represent results from core PI-6, whereas black peaks during the same period represent results from core PI-2.

4.3 Species responses during abrupt climate shifts, and refugia for aquatic taxa

Sedimentological and fossil records from Lake Petén Itzá suggest that the periods of the strongest climatic fluctuations during the last 155 kyr in the northern Neotropics occurred around 85 ka (Mueller et al., 2010) and Heinrich stadials (Correa-Metrio et al., 2012b; Cohuo et al., 2018). Those episodes were characterized by dramatic decreases in the lake level, suggesting intense aridity in the region. The lowest estimated temperatures for the entire record correspond to HS1.

Correa-Metrio et al. (2013) estimated high climate change velocity in the region during HS1, which produced large changes in terrestrial plant communities. Correa-Metrio et al. (2012b, 2014) estimated that one of the consequences of such ecological instability was the substantial migration of tropical vegetation and development of refugia. The high velocity of climate change inferred for the northern Neotropical region is, however, opposite to trends observed elsewhere in the tropics, which suggests that high biodiversity and endemicity are associated with low climate change velocities and high species resilience (Sandel et al., 2011). It remains uncertain how climate change velocity during periods of abrupt climate change affected aquatic communities in the northern Neotropical region. It is also unclear whether aquatic taxa were as dramatically affected as local terrestrial species during these abrupt episodes or if they simply displayed high resilience.

We used the fossil record of freshwater ostracodes from HSs published in Cohuo et al. (2018) to analyse the HS1 structure in detail (Fig. 7) because that was the period of the coldest temperatures and extreme drought during the last 85 kyr (Mueller et al., 2010; Correa-Metrio et al., 2012a; Cohuo et al., 2018).

Figure 7 Fossil record of two endemic (Paracythereis opesta and Cypria petenensis) and two non-endemic (Cytheridella ilosvayi and Darwinula stevensoni) ostracode species during the period 53–14 ka. Grey horizontal bars represent temporal extent of Heinrich stadials (HS5a–HS1). Modified from Cohuo et al. (2018).

Estimated paleo-temperatures based on δ¹⁸O fluid inclusion data and biologically based (ostracodes and pollen) transfer functions suggest an overall temperature decrease of about 5 ^∘C in comparison with today's temperatures during the HS1. With respect to temperature, environmental conditions probably remained suitable for tropical species (especially endemics) distribution across large areas of the Yucatán Peninsula and in northern Central America. Mueller et al. (2010) estimated that the Lake Petén Itzá water level decreased by ∼50 m during HS1, which would imply that lakes in the region with maximum depths <50  dried completely.

We assume that lakes that held water during HS1 served as “refugia” for aquatic taxa, as temperature apparently did not limit species distributions (Cohuo et al., 2018).

Systems such as cenotes and lakes that are not directly dependent on precipitation to maintain the water level but are instead controlled by large subterranean aquifers (Perry et al., 2002; Schmitter-Soto et al., 2002; Vázquez-Domínguez and Arita, 2010) may serve as “refugia” for aquatic species, enabling native species to remain in the region during periods of low rainfall. To date, it remains uncertain whether lakes and cenotes (approximately 7000 in the Yucatán Peninsula) held water during HS1, and little is known about their spatial distribution. Isolated water bodies (refugia) may explain the high percentage of endemicity and micro-endemicity (species distributed in a single or limited group of lakes) for aquatic taxa on the northern Yucatán Peninsula (Mercado-Salas et al., 2013). Species that inhabited such systems may have remained isolated and adapted to specialized environmental niches.

Deevey et al. (1983) studied sediment cores from Salpetén (z_max = ∼30 m) and Quexil lakes (z_max = ∼ 30 m), Guatemala, and inferred that most lakes, including cenotes, in the northern Neotropics dried out during the deglacial period because of the hydrological sensitivity of the region. They also found that most lake sediment cores from the region bottom out at ∼8 ka, which means that the lakes probably first filled in the early Holocene in response to wetter conditions and a rising sea level, which raised the local water table. The authors therefore suggested that only large lakes in the region, with maximum depths >50 m (e.g. Petén Itzá, Macanché, Atitlán, Coatepeque and Ilopango), held water during the dry deglacial period but possessed water chemistry much different from today, which limited habitats for aquatic species.

This second scenario favours the hypothesis of central populations (meta-populations) in one or more large lakes, which enabled species exchange with surrounding aquatic environments, thereby preventing species losses in small populations by demographic stochasticity. The two scenarios are not mutually exclusive, and it is possible that both account for the success of aquatic tropical taxa through periods of abrupt or prolonged climate fluctuations. Lake Petén Itzá may have played an important role for aquatic species survival and dispersal in the northern Neotropical region because it held water for at least the last 400 kyr (Kutterolf et al., 2016).

Our findings contrast with results from terrestrial environments, which show that HS1 drove plant species to migrate and retreat to a few well-defined micro-refugia (Cavers et al., 2003; Dick et al., 2003; Correa-Metrio et al., 2013). Burrows et al. (2011) demonstrated that the pace of climate shifts in aquatic and terrestrial systems can be very different. They estimated that vegetation responds rapidly to climate change, especially to precipitation and temperature shifts. Indeed, changes in these variables can alter the composition of vegetation abruptly, within a few years. Conversely, in aquatic environments, the velocity of climate change tends to be slower. For instance, given the geomorphology of water systems in the region such as cenotes (small area <1 km² and deep waters >10 m) and large lakes such as Petén Itzá (>120 m deep), dramatic changes in air temperatures are needed to alter the temperature of the water column and thus impact species niche stability. Our study suggests that the velocity of change in aquatic environments remained low in the northern Neotropical region, enabling local species to adapt and specialize to their environments instead of migrating and/or remaining isolated in refugia, as observed in tropical areas elsewhere.

5 Conclusion

Our study integrates species distribution models and paleo-records to reconstruct aquatic species distribution dynamics during the last 155 kyr in the northern Neotropics. Both approaches show strengths and limitations. Species distribution models were afflicted by a degree of uncertainty due to uncertainties of general circulation models MIROC-ESM and CCSM4 simulations related to precipitation and temperature. Although these uncertainties can be considered to be systematic errors, it remains uncertain whether the lower-end simulations based on SDMs generated in this study fully reconstruct suitable areas of distribution of aquatic species, especially because in tropical regions the larger biases in simulated values of precipitation and temperature have been estimated.

The most important limitations of paleo-records relate to the scarcity of fossil evidence spatially and temporally, especially for the older periods evaluated. Low abundances in ostracodes were associated to species ecological preferences, core location and preservation processes. The integration of fossil evidence from two long cores of the Lake Petén Itzá was highly informative, as the full range of the temporal presence and absence of the target species was recovered.

In spite, limitations of both approaches, the comparison of SDM outputs and fossil records, resulted in congruent patterns. For the older periods such as LIG and LGM, temporal agreement between approaches was observed. For the most recent period (middle Holocene), temporal and spatial agreement was observed.

Given the congruence between approaches, our study highlights the following conclusions:

1. Distribution dynamics of endemic and non-endemic species result in similar patterns throughout long-term climatic fluctuations such as glacial–interglacial cycles and marine isotope stages.

2. More divergent patterns can be observed during episodes of profound climatic alterations such as the LGM and HS1.

3. Endemic species are highly resilient and remained in the core area during periods of strong alteration of temperature and precipitation.

4. Non-endemic species are sensitive to decreases in temperature, being displaced to Central America to track climates compatible with their tolerance ranges.

This study represents, to our knowledge, the first insight into the magnitude of ecological alteration of aquatic ecosystems during different past climatic scenarios in the northern Neotropical region. Further studies may therefore consider refining the spatial and temporal resolutions of the analyses and incorporate additional lines of evidence such as molecular data. The understanding of historical species dynamics can help with generating strategies for the protection of the biota which can be highly threatened by the future emergence of non-analogous climates.