Ice Age effects on genetic divergence of the American crocodile (Crocodylus acutus) in Panama: reconstructing limits of gene flow and environmental ranges: a reply to O’Dea et al.

Abstract

O’Dea et al. (2022) (Pleistocene sea level changes and crocodile population histories on the isthmus of panama: a comment on Avila-Cervantes et al. (2020). Evolution, 76(11), 2778–2783. https://doi.org/10.1111/evo.14610) question our hypothesis that sea-level changes during the past glaciation played a role in restricting gene flow between Pacific and Caribbean Crocodylus acutus in Panama. They argue that an error in sea-level high-stand reconstruction during the last interglacial period (118–130 ka) does not support our hypothesis. Although they are correct in our high-stand reconstruction error, overlooked the point in that we were presenting a model of restricted gene flow across the Panamanian Isthmus during low sea levels. We review the assumptions of gene demographic methods, emphasizing that we were focusing on times of genetic divergence. We expand here why gene flow between these coastal populations could have been restricted during the last glacial maximum (19–26.5 ka) and the 50,000 years preceding it when sea levels were lower than today. O’Dea et al. suggest local climates may have played larger roles than sea levels. We demonstrate that paleoclimate estimates for the past 3.3 Ma in Panama are within the bounds of extant C. acutus. The importance of Ice Age Sea-level dynamics on Neotropical species was likely profound and should be incorporated into evolutionary studies of these taxa.

Crocodylus, Ice Age, relative sea levels, gene flow

We thank O’Dea et al. (2022) for their comments on our publication and for highlighting an error we made using maximum sea levels during the last interglacial (LIG) that occurred about 130,000 years ago (ka) to 118 ka. We take this opportunity to clarify key aspects of our paper, develop a more extensive discussion on the importance of including environmental changes to population demographics, and emphasize that any LIG sea-level change above current levels has no bearing on our results.

Our study analyzed the genetic variation of the American crocodile (Crocodylus acutus) on both sides of the Central American Isthmus. We sampled several populations from the Caribbean and Pacific coasts of Panama to estimate the impact of the formation of the Isthmus of Panama (about 3 million years ago) on the restriction of gene flow and construction of the Panama Canal (about 100 years ago) on the presumed resumption of gene flow between these coastal populations. These populations were hypothesized to be genetically separated at least until the construction of the Panama Canal. Species distribution maps of C. acutus in the Atrato Basin in northern Colombia sometimes reconstruct a continuous range between the Caribbean and Pacific coast (e.g., Briggs-Gonzalez et al., 2017; Rossi et al., 2020). However, these populations are restricted to the Atrato River with a range extending only to the Caribbean coast and no reported populations on the adjacent Pacific coast (Thorbjarnarson et al., 2006). In addition, we inferred the demographic history and calculated the divergence times of the sampled populations.

Demographic models use genetic and genomic data to estimate population size, migration rates, and divergence times between populations (Gutenkust et al., 2010). Divergence times are estimated from generation times, so the models require some assumptions of average generation times in those populations. We used two methods for the demographic models: Approximate Bayesian Computation (ABC) on DIYABC v.2 (Cornuet et al., 2014) and the Diffusion approximation demographic inference (δaδi) (Gutenkunst et al., 2009). It is important to note that the calculated times of restricted gene flow can only be the most recent isolation and, in cases where secondary contact is inferred, the penultimate isolation (Table 1).

Table 1.

Estimated divergence times for populations using DIYABC and δaδi. Only the most supported models from Avila-Cervantes et al. (2021) are presented.

Divergence	Barrier	Population pairs	DIYABC Model 2	δaδi
			Divergence in thousands of years (ka)	Divergence in thousands of years (ka)	Secondary contact in thousands of years (ka)
Caribbean	Sympatric	GAL-BCI	7.82–9.77	5.8–7.3	3 - 3.7
Pacific	Sympatric	LAG-COIB	11.04–13.8	1.1–1.4	1.1 - 1.4
Caribbean—Pacific	Allopatric	GAL-LAG	32.6–40.75	30.9–38.6	20.7 - 25.9
Caribbean—Pacific	Allopatric	BCI-LAG	NA	53.2–66.8	26.7 - 33.4

The results identified several periods of restricted gene flow with the most realistic ranging between 26.7 and 66.8 ka (Table 1). All were about 2 orders of magnitude younger than the formation of the Isthmus of Panama. Our primary conclusion was that the Isthmus was not a barrier at all, and at best semipermeable to these large-bodied amphibious reptiles. In an attempt to explore other factors that may have been responsible for these relatively recent divergence dates between the Pacific and Caribbean populations, we looked for temporal associations between sea-level lows and restricted gene flow. Extreme sea-level lows are expected to have made the Isthmus which was otherwise semipermeable to crocodile gene flow, less permeable.

The ABC and δaδi divergence times coincide within the period of the LGM (26.5–19 ka) and the penultimate glacial maximum (140 ka) (Colleoni et al., 2016). These periods were characterized by extreme rises and falls in eustatic sea levels. These levels varied from +6 m (Blanchon et al., 2009) during the LIG to −130 to −120 m during the LGM (Gowan et al., 2021; Miller et al., 2020).

The major argument O’Dea et al. have is our reconstruction of sea levels around Panama during the LIG using the uncorrected values presented by Rohling and colleagues (2017), and that our divergence times do not match exactly with the LGM. However, in the context of our work, the absolute LIG sea-level height does not impact our results as our demographic models were only testing for the most recent discernible barrier to gene flow between each coast and subsequent potential secondary contacts. Any sea-level drop over current levels is expected to reduce gene flow across the Isthmus. The beginning of the LGM (19–26.5 ka) nearly overlaps with the recovered dates of population divergences between the allopatric adjacent populations (26.7–66.8 ka) and we suggested this was a driver in genetic isolation between the two coastal populations. Although these times do not completely overlap, the coincidence of the dramatic sea level drops during the LGM was discussed as a possible mechanism of gene flow restriction (Avila-Cervantes et al., 2021). Any earlier barriers to gene flow during older glacial maxima cannot be inferred, as these events would have been erased by subsequent pulses of genetic flow between the coasts.

However, the LGM low stand sea levels occur at least 2–40 ka after the estimated genetic divergence between the populations on the two coasts. This does not discount the possibility that during the drop in relative sea level (RSL) preceding the LGM, a threshold was passed where sea levels were low enough to cause the Isthmus to switch from a semi-permeable barrier to an impermeable barrier to gene flow in these animals. To further explore this, we present models of intermediate low RSLs that occurred throughout the time span between the LIG and LGM (Figure 1). These sea-level drops were nonlinear. Between approximately 20–70 ka were more than 60 m below current levels (Rohling et al., 2017). We hypothesize these lows were enough to significantly limit gene flow. This 50,000-year-long period preceding the LGM overlaps the divergence estimates we recovered and suggests this was the threshold in RSL drop that factored into their isolation (Figure 1B). Prior global glaciations and interglacial periods that occurred throughout the Ice Age may have had similar effects on crocodile gene flow across the Isthmus but are simply not recoverable due to the repeated periods of admixture between the two coasts. Only the last major restriction to gene flow event is recoverable with the present genomic data. There is no fossil record of the species to calibrate earlier divergences.

Figure 1.

(A) Relative sea level (RSL) of 6 m above present during the LIG, about 118–130 ka (Blanchon et al., 2009). (B) RSL of 60 m below present, about 50 ka preceding LGM (Rohling et al., 2017). Note that for most of this time RSL was even lower. (C) RSL of 130 m below present during the LGM, about 19–26.5 ka (Gowan et al., 2021; Miller et al., 2020). We used the GEBCO grids (2021) and the Panama province boundaries map from the STRI GIS data portal (https://stridata-si.opendata.arcgis.com/) to create the maps.

O’Dea et al. suggest paleoclimates may have played a more important role than RSL on C. acutus gene flow over the Isthmus. Below, we present data that challenge this hypothesis. Adaptations of the genus permit C. acutus to be the most widely distributed crocodile species in the Neotropics (Oaks, 2011). The species inhabit mangrove-lined coastal lagoons or estuaries, offshore cays and coral atolls, and rivers and reservoirs. These environments range from hypersaline to freshwater (Thorbjarnarson, 1989). Water salinities reported for C. acutus can vary from 0 ppt in coastal lagoons to >50 ppt in atolls, and offshore salinities range from 34 to 36 ppt (Platt et al., 2013). The species inhabit environments with contrasting temperatures, precipitation, and elevations over a vast range that spans from the northwestern Pacific coast of Mexico to Ecuador to the Caribbean islands. To explore the bioclimatic profile of C. acutus we used 400-point observations of the species (Thorbjarnarson et al., 2006) and extracted the values of the eight bioclimatic variables from WorldClim 2 (Fick & Hijmans, 2017) of each locality, at a resolution of 30 arc seconds. We used R V.3.6.1 (R Core Team, 2021) and the packages DISMO V.1.3 (Hijmans et al., 2020) and MAPTOOLS V.1.2 (Bivand & Lewin-Kohn, 2021). The results indicate a broad range of temperature and precipitation regimes for the species’ localities (Table 2). These extreme arid to tropical and cool to hot environments present a remarkably resilient species that subtle paleoclimatic changes are unlikely to have affected. At no point in the past 3.3 million years has the paleoclimate of Panama been estimated to have passed beyond the bioclimatic variables C. acutus encounters within its current distribution. Table 3 indicates the profile of 70 localities of C. acutus in Panama (Thorbjarnarson et al., 2006) for the bioclimatic variables from PaleoClim (Brown et al., 2018) with a resolution of 2.5 arc seconds during the last glacial maximum (21 ka), the LIG (130 ka), and the Pliocene (3.3 Ma). Temperature and precipitation values are well within the ranges of the species’ actual distribution today. The environmental variation and ranges among current populations of C. acutus are stunning and highlight the robustness of this species.

Table 2.

Bioclimatic profile of C. acutus based on 400 localities in its range of distribution (Thorbjarnarson et al., 2006) and the bioclimatic variables from WorldClim 2 (Fick & Hijmans, 2017).

Bioclimatic variable	Min	Mean	Max
Annual mean temperature (°C)	21.31	26.13	28.81
Max temperature of warmest month (°C)	29.1	32.55	35.9
Min temperature of coldest month (°C)	9.9	19.54	23.5
Temperature annual range (°C)	7.4	13.01	25.8
Annual precipitation (mm)	217	1,838	3,840
Precipitation of wettest month (mm)	70	307	734
Precipitation of driest month (mm)	0	29.1	132
Elevation (m)	0	91.38	1,265

Table 3.

Bioclimatic profile of C. acutus based on 70 localities in Panama (Thorbjarnarson et al., 2006) and the bioclimatic variables from PaleoClim (Brown et al., 2018) for the last glacial maximum (21 ka) (Karger et al., 2021), the last interglacial 130 ka (Otto-Bliesner et al., 2006), and the Pliocene (3.3 Ma) (Dolan et al., 2015).

	Last glacilal maximum (21 ka)			Last Intreglacial (130 ka)			Pliocene M2 (3.3 Ma)
Bioclimatic variable	Min	Mean	Max	Min	Mean	Max	Min	Mean	Max
Annual mean temperature (°C)	19.5	22.3	23.5	23.1	25.6	26.5	24.7	27.1	28.1
Max temperature of warmest month (°C)	22.4	26.6	30	26.9	29.5	30.5	NA	NA	NA
Min temperature of coldest month (°C)	16.2	18.9	20.4	20	22.4	24.1	NA	NA	NA
Temperature annual range (°C)	4	22.2	23.6	5.6	7.1	8.8	NA	NA	NA
Annual precipitation (mm)	1,358	2,121	3,214	1,380	2,381	3,340	1,908	2,766	3,753
Precipitation of wettest month (mm)	198	322.6	514	216	365	550	463	570.7	753
Precipitation of driest month (mm)	9	35.91	182	11	28.4	169	0	3.9	137

We agree that it is important to consider other drivers than changes in sea level to understand the biology and evolution of C. acutus in the Isthmus of Panama. However, the species biology and its capacity to inhabit a wide range of environments (Table 2) make it resilient to the increased salinity (discussed above) and climatic variation and rainfall changes related to the glacial–interglacial cycles (Table 3).

O’Dea et al. also suggest a series of other avenues to explore and interpret crocodile gene flow disruption across the Isthmus such as paleontology and archeology. Although there is a wealth of crocodyliform paleontological and archeological data we don’t see how the Miocene gavialoid or relatively recent archeological records O’Dea et al. suggest could help interpret Crocodylus acutus gene flow across the Isthmus during the Pliocene and Quaternary. The species of gavialoid suggested, Aktiogavialis caribesi, was recovered from the late Miocene Urumaco Formation (c. 10–9 Ma) (Salas-Gismondi et al., 2019), predating the formation of the Isthmus by up to seven million years and long before the immigration of Crocodylus from Africa, that occurred about 5 Ma (Avila-Cervantes & Larsson, 2018; Oaks, 2011). The archeological data postdates any of our inferred divergence dates by several thousand years and cannot offer clues to the factors acting on those divergences during the time span between the LIG and LGM.

In summary, we agree that we made an error in the reconstruction of the RSL high stand during the LIG. Sea-level high stands are likely to facilitate gene flow for these animals. The majority of terrestrial and marine systems studied in the region recover little to no influence of the Ice Age on population demographics but focus on the original formation of the Isthmus (Bacon et al., 2015; Cowman & Belwood, 2013; Lessios, 2008; Stange et al., 2018). However, more restricted coastal populations, such as mangroves, are revealing startling regional demographics that coincide with the LGM (Ceron-Souza et al., 2015).

We did find that the genetic divergence of Crocodylus acutus does not coincide with the formation of the Central American Isthmus, 3 Ma, but instead coincides with the drop in RSL leading up to and during the LGM, about 20 ka. Although paleoclimate most certainly played a role in crocodile local evolution, we demonstrate that the extreme tolerances of extant C. acutus make such claims for their Pacific-Caribbean divergence unlikely. Rather, the drop in RSL leading to the LGM seems like a plausible explanation that could have caused the Isthmus to become a barrier to gene flow for these amphibious taxa but further sampling and testing are needed to confirm this theory. The wide range of climate regimes for this species implies that past climates would have had little effect on its distribution; however, we expect that the relative elevation and width of the Isthmus should. We highlight that Ice Age Sea-level dynamics may have played larger roles in the evolution of Central American taxa and processes of genetic disruption than previously appreciated. Moreover, C. acutus may be a remarkable taxa to begin examining these effects due to its broad climate distribution and amphibious habitat.

Data availability

There is no data to be archived.

Author contributions

J.A.C. wrote the manuscript. H.C.E.L. wrote the manuscript. All authors gave final approval for publication.

Funding statement

Support for this project was provided by NSERC Discovery Grant to H.C.E.L.

Conflict of interest: The authors declare no conflict of interest.

Acknowledgments

This research was performed using the infrastructure of the Integrated Quantitative Biology Initiative, funded by the Quebec government, McGill University, and the Canadian Foundation of Innovation project 33122. We would like to thank the reviewers for their thoughtful comments toward improving our manuscript.