Synergistic impacts of global warming and thermohaline circulation collapse on amphibians

Impacts on ecosystems and biodiversity are a prominent area of research in climate change. However, little is known about the effects of abrupt climate change and climate catastrophes on them. The probability of occurrence of such events is largely unknown but the associated risks could be large enough to influence global climate policy. Amphibians are indicators of ecosystems’ health and particularly sensitive to novel climate conditions. Using state-of-the-art climate model simulations, we present a global assessment of the effects of unabated global warming and a collapse of the Atlantic meridional overturning circulation (AMOC) on the distribution of 2509 amphibian species across six biogeographical realms and extinction risk categories. Global warming impacts are severe and strongly enhanced by additional and substantial AMOC weakening, showing tipping point behavior for many amphibian species. Further declines in climatically suitable areas are projected across multiple clades, and biogeographical regions. Species loss in regional assemblages is extensive across regions, with Neotropical, Nearctic and Palearctic regions being most affected. Results underline the need to expand existing knowledge about the consequences of climate catastrophes on human and natural systems to properly assess the risks of unabated warming and the benefits of active mitigation strategies.


Supplementary Information
Contents S1. Biogeographical realms S2. Description and analysis of climatic and bioclimatic variables S2.1 Description of control and hosing climate experiments S2.2 Differences between scenarios and time horizons S2.3 Climatic anomalies and novel climates in thermohaline circulation weakening scenarios S3. Ecological modelling approach and extended results S3.1 Selection of ecological niche modeling algorithms S3.2 Estimates of range contractions using ecological niche modelling algorithms S3.3 Range contractions across most diverse taxonomic groupings and extinction risk status S3.4 Potential effect of uncertainty in model algorithm, threshold criteria and grain size S3.5 Comparisons of projected range contractions between scenarios and time horizons S3.6 Geographical patterns of amphibian species losses S1. Biogeographical realms We selected endemic amphibian species from six biogeographical realms (Figure S1; see main text). These geographical delimitations are used by biogeographers and macroecologists to organize global diversity. Figure S1. Biogeographical realms (in color) used in this paper. Colors as follow: blue: Nearctic; yellow: Neotropical; green=Afrotropical; brown=Australian; cyan: Indomalayan; red: Palearctic.

S2.1 Description of control and hosing climate experiments
The climate simulations used in this study were produced with the Institut Pierre Simon Laplace low-resolution coupled ocean-atmosphere model (IPSL-CM5-LR) 1 . The spatial resolution of the atmospheric component is 3.75ºx1.875º in longitude and latitude, respectively, and includes 39 vertical levels. The nominal resolution of the oceanic component is 2º with a higher latitudinal resolution of 0.5° in the equatorial ocean, and 31 vertical levels. The locations for the release of the freshwater are deep water formation regions in the North Atlantic (45ºN to 65ºN, 45ºW to 5ºE), which are classical regions of spread of the input of Greenland meltwater 2 . Since this freshwater input corresponds to Greenland ice sheet melting in terms of its long-term reservoir of freshwater, there is no reason to compensate for the water mass added, and no compensation is added in the model. Since this is a free surface model, this is not requested for conservation issues.
Five experiments are considered for the assessment presented here: a) the control run is the RCP8.5 emission scenario; b) four hosing experiments that are based on the RCP8.5 but include the addition of 0.11, 0.22, 0.34, and 0.68 Sv (1 Sv = 10 6 m 3 /s) of freshwater released in the North Atlantic from 2020 to 2070 3 . These experiments are labeled as A, B, C, and D, respectively. All simulations cover the period 2006-2100. Figure S2 illustrates the additional weakening of AMOC over the century that is produced by each of these hosing scenarios. The bioclimatic indices used in this paper are derived from the monthly temperature and precipitation data produced with the IPSL-CM5-LR.
Early studies about the effects of an AMOC collapse were based on climate models' simulations in which hosing experiments were imposed to pre-industrial climate conditions 4 (i.e., other forcing factors such greenhouse gases emissions were held constant). While this type of experiment can help to analyze effects such as climate feedbacks in isolation, a collapse of AMOC is not to be expected to occur under pre-industrial climate or without significant external forcing. In this paper, we use recent simulations in which the hosing experiments are imposed in the course of a high-emissions scenario (RCP8.5) 3 . A highwarming scenario such as the RCP8.5 provides a consistent baseline scenario for exploring the impacts of significant melting of Greenland that produces additional and substantial weakening or collapse of AMOC. Although the probability of a collapse of AMOC is not known, the literature strongly suggests that these probabilities are expected to increase with the level of warming 5,6 . The RCP8.5 represents a high-emission scenario and not a 'businessas-usual' scenario 7 . However, a recent study has shown that this scenario is in close agreement with historical cumulative CO 2 emissions, as well as with projections out to midcentury under current and stated policies, and is characterized by highly plausible levels of CO2 emissions in 2100 8 . We do not assign any probability of occurrence to the RCP8.5 scenario nor to the collapse of AMOC, instead we center on investigating the impacts on amphibian species this climate catastrophe could have.
The assessment of impacts in this paper focuses on climate and thus 30-year averages centered around 2030, 2050 and 2070 are used to represent climate conditions. The effects of climate variability including oscillations such as El Niño/Southern Oscillation (ENSO) and the North Atlantic Oscillation (NAO), or multidecadal climate variations are not addressed in the analysis and may require the use of large ensembles which are not available at the moment. Figure S2. Atlantic meridional overrunning circulation (AMOC) index at 26°N. The AMOC index is defined as the maximum of the Atlantic meridional stream function at 26°N below 500 m. The time evolutions of the indices are shown in percent, in reference to the mean state of the AMOC in the IPSL-CM5A-LR historical simulation over the period 1850-1950. The black line shows the control simulation based on the RCP8.5 projection, while the red, blue and light blue refer to the superimposed hosing experiments with the addition of 0.11, 0.22 and 0.68 Sv of freshwater from 2020 to 2070, respectively (1 Sv = 10⁶ m³/s). Thick and thin lines represent 30-year moving averages and annual frequency data, respectively.

S2.2 Differences between scenarios and time horizons
Bioclimatic variables are more biologically meaningful than precipitation and temperature variables alone and thus are commonly used in ecological modelling. We constructed five bioclimatic variables used to model amphibian species' distributions and that have been considered important for amphibian's biology in previous studies 9 : Annual mean temperature (AMT; bio1); warmest month temperature (WMT; bio5); coldest month temperature (CMT; bio6); total annual precipitation (AP; bio12); and precipitation seasonality (PS; bio15). These bioclimatic variables have a strong influence on species' distributions and functional traits [9][10][11][12] and have been used to model species' distributional areas both under current 9 and future climate change conditions 12,13 . Part I of Figures S3 to S7 shows maps of these bioclimatic variables under the control scenario and all the hosing experiments described in S2.1 for the 2070 period, while Part II shows the differences between the control (RCP8.5) and each one of the hosing experiments. Figure S8 shows radial graphs of the median change for each region, bioclimatic variable and horizon (2030,2050,2070), with respect to the reference climatology . Departures from symmetry in each axis denote differences between the RCP8.5 and the RCP8.5 plus additional and substantial weakening of AMOC (experiment D) (e. Some salient features of this figure illustrate the effects of additional and substantial weakening over the climate of the selected regions. For instance, under the RCP8.5 scenario median changes in AMT, CMT and WMT tend to be similar in magnitude and sign for all regions. However, with the additional weakening of the AMOC important contrasts between changes in CMT and WMT appear in all regions and large differences in precipitation with respect to the RCP8.5 control simulation tend to occur. Moreover, the variation of deltas across scenarios illustrates important nonstationarities across geography in the departure of current climate conditions in the hosing experiments ( Figure S9-S14). For instance, CMT and WMT exhibit large departures in all freshwater discharge scenarios with respect to the control simulation. These figures illustrate the degree of climatic departure from baseline conditions for each of the five bioclimatic variables used in the ecological niche models ( Figure S9-S14). The variation of these departures is context-dependent and some regions would be more affected by temperature whereas other by changes in rainfall regimes. As denoted by the boxplots in these figures, all bioclimatic variables show very large variability within regions. The delta values from these figures ( Figure S9-S14) represent the differences between a given future scenario (e.g., control and hosing experiments) and baseline conditions . The non-stationarity across time and space shows the complexity of these abrupt climate changes and of their potential impacts on ecosystems. For instance, the largest departures were observed in precipitation variables (AP and PS) across all six realms ( Figure S9-S14). The departure of extreme monthly temperatures (CMT and WMT) from baseline conditions is contextdependent and in some regions is larger (e.g., Palearctic) than in others (e.g., Neotropical). As it is discussed in the following subsection, large heterogeneities in climatic departures from current conditions are also present within regions which also contribute to the large spatial variability in the projected impacts on biodiversity.

S2.3 Climatic anomalies and novel climates in thermohaline circulation weakening scenarios
We also analyzed the spatial emergence of novel climates (i.e., non-analog climates) 14 and climatic departures based on the set of five bioclimatic variables described above. These analyses were performed for all experiments: the control reference scenario (RCP 8.5) and each of the four hosing experiments (0.11 Sv; 0.22 Sv; 0.34 Sv; 0.68 Sv). Climatic departures were calculated as the Euclidean distances between future climatic conditions and current climatic conditions ( Figure S15) 15 . Some regions experience larger departures in the short term (2030; e.g., Australasian, Nearctic and Palearctic) than others ( Figure S15). However, Euclidean distances reveal that more extreme climatic conditions could emerge in several regions across time under an AMOC collapse scenario (see outliers in Figure S15). Figure

Selection of ecological niche modeling algorithms
A set of 15 species were selected randomly for each biogeographical realm and species distribution models were fitted using five algorithms: MaxEnt, MARS, CART, ANN, GLM and GBM to explore which algorithms have the best predictive performance (e.g., omission rate, AUC, TSS, and Kappa; Table S1; see also 16,17 for full description of these validation metrics). A full description of each model algorithm can be found in Peterson et al. 18 and Guisan et al. 19 . We selected only three algorithms (MaxEnt; CART; and BRT) that have a good model performance based on the lowest omission rate, high AUC, and TSS values and a low omission rate. Species distribution models for 2509 species were estimated with these three algorithms. All models show a high predictive performance (high TSS, specificity and sensitivity values; Figure S17). We excluded from subsequent analyses those species with a poor performance denoted by low values in validation metrics (e.g., TSS <= 0.4, sensitivity and specificity <= 0.5).

S3.2 Estimates of range contractions using ecological niche modelling algorithms
We performed ecological niche modeling for 2509 endemic species to six biogeographical realms across the world (Table S1). The biogeographical realm with the highest number of modelled species was Neotropical (1211) followed by Afrotropical (466 species). The Palearctic realm had the lowest number of modelled species (108 species). We found an extensive variation in the response of amphibian species to future climate change scenarios. This suggests that amphibian species can respond either positively (i.e, expand its ranges) or negatively (i.e., contract its ranges) to future climate conditions. More species will experience range contractions than range expansions under all climate change scenarios. This result is relatively similar across the three model algorithms (Table  S2). Furthermore, the percentage of range contractions was relatively similar across biogeographical realms (Table S3). In Nearctic and Palearctic realms, we found that more species contract their ranges across all the hosing scenarios in comparison with a control scenario (IPLS-RCP 8.5). This suggests that endemic amphibians in extra-tropical realms seem to be more sensitive to potential contractions in their distributional areas with climate change scenarios simulating a thermohaline weakening. Table S3. Percentage of species suffering range contractions across biogeographical realms in a high-warming scenario (RCP 8.5) and four hosing experiments (0.11, 0.22, 0.34, 0.68 Sv; 1 Sv = 10 6 m 3 /s) of freshwater release in the North Atlantic from 2020 to 2100), modelled using a maximum entropy algorithm.

S3.3 Range contractions across most diverse taxonomic groupings and extinction risk status
Range contractions for the most diverse taxonomic groupings and extinction risk status were also estimated for all climate experiments described in section S2.1. Estimated range contractions are large for all taxonomic groupings and risk status even for the 2030 horizon. For the control simulation, the median reduction in distribution range for all taxonomic groupings and risk status is close to 50%, increasing to about 75% by the end of the century (Figures S18, S19). "Critically endangered" and "endangered" show lower range contractions than other categories, particularly when full dispersal is allowed. The range contraction of species is much larger compared to the control scenario across taxonomic groupings ( Figure S18), extinction risk status ( Figure S19), but not between freshwater discharge levels (Figures 1, S18, S19).

S3.4 Potential effect of uncertainty in model algorithm, threshold criteria and grain size
Our results are robust to different kinds of uncertainty in projected geographical distributions across a high-warming scenario and four hosing experiments. First, we found that the three algorithms show that, on average, species are projected to experience extensive losses of suitable areas across hosing experiments ( Figure S20). We found similar results for the percentage of range contraction, independently of the threshold criteria used ( Figure S21). Finally, we did not find differences using a coarser-grain scale (1°) in the patterns of range contraction ( Figure S22).   Our results suggest that impacts of a thermohaline weakening are severe across the three time periods examined (i.e., 2030, 2050, 2070), levels of freshwater discharge (i.e., 0.11 to 0.68 Sv; 1 Sv = 10 6 m 3 /s), biogeographical realms, and two dispersal scenarios (Figure 1). Comparisons between projected range losses from the Greenland ice sheet melting simulations and a control simulation reveal substantial increases in projected range losses are expected from these catastrophic scenarios (Figure 1). In addition, we note that endemic species from some biogeographical realms were more susceptible to range contractions. These regions include Palearctic, Nearctic and Australasian, although with a slight variation between model algorithms. The severity of the hosing experiments on a high-warming scenario was evaluated across the conservation status of amphibian species using the extinction risk ranking from IUCN 20 . Some studies have provided evidence that life history and geographical range traits that predict extinction risk are useful as indicators to anticipate vulnerability to climate change 21 . Many of these categories for amphibian species were defined based on geographical range size, population density or relative abundance and potential and direct human stressors. We evaluated the percentage of range contractions across IUCN categories from critically endangered (CR) to least concern (LC) (Figure 1). Threatened species (i.e., species classified as critically endangered, endangered, and vulnerable) will suffer more substantial contractions across hosing experiment scenarios than non-threatened species (i.e., least concern; Figure 1). These results suggest that species with small geographical ranges and less abundant in the wild will be more affected by a rapid ice sheet melting affecting global climate patterns.

S3.5 Comparisons of projected range contractions between scenarios and time horizons
We found significant differences between all four hosing experiments and the control one. The differences were broadly similar either under a full dispersal or no dispersal scenario (Table S4). Full dispersal refers to a scenario where species can occupy suitable climates across a biogeographical realm. In contrast, no dispersal refers to the lack of ability to colonize suitable climates outside its distributional area (Table S4). Species richness maps were generated based on the stacking of presence-absence maps for each species in each biogeographical realm under current climatic conditions and for each climate change scenario. We estimated the percentage of species loss and species gains by each pixel using these species richness maps.
The percentage of species loss were calculated based on the total number of current species richness and how many species were lost in each pixel under each of the climate change scenarios ( Figure S23-S29). The geographical patterns of species richness were calculated separately for each ecological niche modeling algorithm (Figures S23-S27) and then were averaged across algorithms to explicitly incorporate this uncertainty source (Figures 2, S28-S29). Similarly, the percentage of species gains (Figures S34-S35) were calculated based on the total number of current species richness and how many species were gained in each pixel under each of the climate change scenarios.
Under a full dispersal scenario, the percentage of projected species losses varies extensively across geography but there are strong similarities across the hosing experiments and temporal horizons, except for 2070 and the most extreme hosing experiment (0.68 sv; Figures S24-S27) in which species loss was much more severe. All realms would be affected by several local extinctions. These results suggest that the potential impacts of climate change plus a thermohaline shutdown would be non-stationary and widespread across regions. We note that the projected local extinctions would be more severe in the Palearctic, northwestern Nearctic, north of India, east of Africa, and southeastern Australia (Figure 2 and Figure S23-S29). The amphibian assemblages inhabiting these regions would be more impacted by severe range contractions driven by a catastrophic climate change event. To further illustrate the effects of additional weakening of AMOC, figures S30 to S32 show the difference in the geographical patterns of percentage of species losses between RCP8.5 and each one of the hosing experiments. These figures reveal the highly non-linear response of the amphibian assemblages to these freshwater discharge scenarios and suggest the possibility of a tipping point behavior.
The estimated percentage of species loss is the result of the interplay between changes in climate and the species response to such changes. As discussed in section S2 and in previous paragraphs, changes in climate are spatially very heterogeneous, particularly under the hosing experiments. Similarly, the responses to changes in climate are idiosyncratic and the same climatic variation can produce substantially different impacts on the 2,509 amphibian species considered in this paper. Figure S30 illustrates the spatial heterogeneity of changes in climate and of impacts on amphibian species, as well as their association at the grid scale for experiment D (0.68 sv) and the control scenario (RCP 8.5). As is shown in this figure, large changes in climate are typically associated with high percentage of species loss. However, due to the diversity of species and their sensitivity, in some cases large (small) losses can occur for relatively small (large) changes in climate. The diversity in changes in climate at the grid cell level and the idiosyncratic responses of the amphibian species modelled makes it difficult to extract general conclusions about what combinations of the bioclimatic variables drive the overall impacts in each realm.
Under a no dispersal scenario, the percentage of species losses are consequently more drastic ( Figure S34-S36). The species losses under a high-emission scenario for the first time horizon are less severe than for the hosing scenarios ( Fig. S34-a vs. Fig. S34-c). The differences are not very large across time horizons for experiment D ( Figure S34-b vs. Fig. S34-d). In addition, the species losses were highly similar across hosing scenarios and time periods (Fig.  S36). However, we consider that these no dispersal scenarios likely are unrealistic because they reduce to zero the ability of species to track their climate requirements through geography and the projected species losses might be extremely inflated. Figure S37 shows the percentage of range expansions in each one of the scenarios. Range expansions tend to be reduced and homogeneous across all climate change scenarios.
Finally, the geographical patterns of percentage of species gains (figures S38-S39) show that species gains in regional assemblages are more frequent and pervasive in a high-emissions scenario (RCP8.5) than in hosing experiments ( Figure S34). Although the Neotropical region tends to exhibit some gains of species in comparison to other regions, the highest gains were toward temperate regions (e.g., the United States and Canada). Although these maps show the opposite patterns of species losses, these gains implicate that any species can colonize a given cell from any part of each realm and this can be unrealistic for many species due to limited dispersal. A further study will be necessary to evaluate properly the turnover in species richness across time horizons and scenarios using recently developed metrics of temporal turnover 22 .          Red values indicates that species losses were higher in a high-emission scenario (RCP8.5) than a given hosing experiment. By contrast, blue values indicate that species losses were higher in a given hosing experiment than the control high-emission scenario (RCP8.5). Figure S33. Association between percentage of species loss and climatic departures in 2070. Panels a) and b) show the association between percentage of species loss and climatic departures for the hosing experiment D and RCP8.5 control reference simulation, respectively. Color gradient from: 1) grey to blue denotes increasing levels of climatic departure; 2) grey to yellow indicates increasing percentage of species loss; 3) grey to red denotes both increasing levels of climate departure and percentage of species loss.