2.1 Model calibration

The basic information to carry out the ecological niche modeling are geographic records and climatic variables. We obtained the geographic records of A. phalloides between 1950-2019 from the Global Biodiversity Information Facility (GBIF) (^{Wolfe et al., 2010}; ^{Mushroom Observer, 2019}). With these data, firstly we visualized the change in the distribution of this toxic mushroom across time and geographic space: 1950-1975, 1975-2000, and 2000-2019. Afterwards, we used records from 1970-2019 to carry out the modeling and identify climatic suitability in North America under a baseline scenario and two different time-series: 2050 (2041-2060), and 2070 (2061-2080). We used records from 1970-2019 to calibrate the model given that in this time there was more certainty about the taxonomic identification of A. phalloides (^{Pringle and Vellinga, 2006}). We obtained a final number of records of 435 in Europe and 339 in North America. To avoid sampling bias and to achieve a better modeling result, we tested different minimum records distances and used 30 km in Europe to avoid aggregation. In North America the minimum distance between records was 10 km.

Climatic variables were obtained from WorldClim (https://www.worldclim.org/), where 19 bioclimatic variables have been created. These variables have more biological sense than only using annual tendencies; they also represent limiting values and variability (^{Fick and Hijmans, 2017}). From these 19 climatic variables we chose eight using biological and statistical criteria (Table I). The number of variables were reduced through a correlation matrix (eliminating those highly correlated variables; R > 0.8) and choosing the most important variables in A. phalloides biology. The variables had a 5-min resolution (~10 km), which has proven to give sound results for biological systems at continental scales (^{Ureta et al., 2022}).

We used the CanESM2 general circulation model, whose last generation (CMPI5) has a good performance for Mexico and the Northern Hemisphere in general (^{Altamirano et al., 2021}). Given that there is no independence between generations of general circulation models, we can assume that CanESM5 (CMIP6) still has a good performance for Mexico and the Northern hemisphere. General circulation models are mathematical models of the atmospheric and oceanic circulation that allow the projection of climatic conditions at specific times. These models have different alternatives of scenarios or socioeconomic hypothesis that range from optimistic to pessimistic (^{O’Neill et al., 2020}). In this study we used the pessimistic scenario SSP5-8.5 in the time-series 2050 (2041-2060) and 2070 (2061-2080) to have an insight of the consequences through a precautionary principle approach (^{UN, 1992}).

2.2 Climatic potential distribution of Amanita phalloides

To choose the final modeling algorithm, we firstly carried out a pilot analysis using seven different algorithms (MaxEnt, MARS, CTA, SVM, RF, GBM, and ANN) (^{Thuiller et al 2009}) Table II and using different predictive power metrics (ROC tests, Kappa, and TSS) to evaluate their performance (^{Allouche et al., 2006}; ^{Peterson et al., 2008}). These three metrics are based on omission and commission errors. For the algorithm evaluation we used the multimodel platform BIOMOD in R (^{Thuiller et al., 2009}). The algorithm with best performance was MaxEnt (^{Phillips et al., 2006}) (Table II).

To improve the MaxEnt performance, different parameter combinations were evaluated through the ENMeval package proposed by ^{Muscarella et al. (2014)}. The ENMeval tool helps identifying the best feature combinations to optimize the modeling. We used the 10 best combinations. The criteria to select the combinations were based on the ROC test, omission rate, and AICc (^{Akaike, 1974}). We compared the validation metrics between the default tuning and the suggested by default and the models obtained through tuning the algorithm with best performance. To obtain the variables that contributed more to the model we used a Jacknife test. We used 70% of the European records to calibrate the model and 30% to validate it. We carried out 10 replicates and used the cloglog suitability output, since it has been suggested as the most adequate to estimate presence probability (^{Phillips et al., 2017}).

Once we got the best algorithm tuning we carried out three different experiments under the baseline scenario (1970-2000): (1) using 100% of the European records to calibrate the model and the North American records for validation; (2) using 70% of the South European records to calibrate the model, 30% of the South European records for validation and projecting suitability in North America, and (3) using the North American records to calibrate (70%) and validate the model (30%), projected also to North America. For the experiment (1), we used 100% of the European records to calibrate the model, given that we had an independent database for its validation in North America. To convert all of our suitability maps into binary maps (presence/absence), we used the minimum suitability value in which there was at least one North American presence. With these three experiments we evaluated the potential risk for the species finding climatic suitability in Mexico.

We carried out two different experiments under climate change scenarios: (1) using all European records, and (2) using only North American records. These two experiments allowed us to have a better understanding of the potential risk of A. phalloides reaching Mexico under climate change scenarios with the populations that already exist in our continent vs. the potential risk if all European populations reach the continent.

To identify the areas under greater risk of intoxication in Mexico, we overlaid two different suitability maps with a map of the indigenous languages of Mexico (^{INEGI, 2010}): (1) climatic suitability (representing lower risk), and (2) climatic suitability + biological suitability map (representing the highest risk). A. phalloides establishes ectomycorrhizal relationships; consequently, its distribution not only depends on the climatic niche, but it strongly depends on the distribution of plant hosts. Even when biotic variables were not directly incorporated into the modeling, the climatic suitability map was overlaid with the potential distribution map of ecosystems in which this mushroom develops: boreal forest, coniferous forest, temperate mountain forest, temperate forest, wet subtropical forest, and subtropical mountain forest. In this way and under the assumption that these plant communities will be static across time, it is possible to identify geographic areas where climatic and biotic suitability occur under a baseline scenario and climate change scenarios. The map of indigenous languages is a way to infer where the mycophylic communities are. The consumption and commercialization of WEM is not exclusive to indigenous communities, but in Mexico all indigenous communities are mycophylic and use this resource (^{Garibay-Orijel and Ruan-Soto, 2014}). Also, indigenous populations are at greater risk in terms of food security (^{González-Martell et al., 2022}); consequently, the overlay of indigenous populations with the suitability of A. phalloides under a baseline and a climate change scenario helps identifying areas of greater intoxication risk. Even though there are high emigration rates of indigenous people in Mexico, their ancestral territories are fixed and indigenous population always remains there.

3.1 Historic, current, and future distribution of Amanita phalloides

The first time-period visualized was from 1950 to 1975, given that before this time records have greater taxonomic uncertainty. In this period, the great majority of GBIF records were localized in the north of the European continent and some observations occurred in United States and Canada (Fig. 1). In the next time-period (1976-2000), GBIF records start to appear in the north of Italy, France and even Spain. Finally, in the last time-period (2000-2018) European and North American records became more abundant in warm areas such as the south of Italy, Spain and Portugal (Fig. 1). In the United States this mushroom distribution reaches the Mexican border in the east coast.

Fig. 1 Observed distribution of Amanita phalloides across time. 

The first climatic suitability experiment (baseline scenario, 1970-2000, using all European records), shows that A. phalloides can expand its distribution down to Mexico (Fig. 2). In this scenario Mexico has climatic suitability in 34% of its area, while this percentage in the United States and Canada accounts to 56 and 15%, respectively (Fig. 3). When we overlaid the map of indigenous communities with the climatic suitability map we obtained that about 70% of the indigenous communities’ areas have climatic conditions suitable for the development of this toxic mushroom. This area comprises the original territories of at least 17 Mexican indigenous groups that consume traditionally wild mushrooms: Tarahumaras in Chihuahua; Tepehuanos, Coras, and Huicholes in Nayarit and Jalisco; Huastecos, Tepehuas, Nahuas, and Totonacos in Veracruz; Chinantecos, Mazahuas, Mazatecos, Mixes, and Zapotecos in northern Oaxaca; and Choles, Tzeltales, Tzotziles, and Tojolabales in Chiapas (Fig. 3). The percentage of indigenous communities’ areas that overlay with climatic and biological suitability (potential distribution of known hosts) under the baseline scenario was of 5% (Table III). The main affected communities were Tarahumaras, Tepehuanes and Nahuas, who can be considered at highest risk.

Fig. 2 Climatic suitability of Amanita phalloides under the baseline scenario. (a) Using European records to calibrate the model, (b) using south-European records to calibrate the model, and (c) using North American records to calibrate the model. 

Fig. 3 Potential distribution map of Amanita phalloides using European records. (a) Overlap between A. phalloides potential distribution and the distribution of its vegetation hosts in North America. (b) Overlap between A. phalloides potential distribution, host distribution and indigenous territories in Mexico. Grey: potential distribution of A. phalloides; green: potential distribution of its hosts; pink: indigenous territories. 

Climatic suitability was also found in very warm areas of Mexico such as the Yucatan peninsula, where the natural hosts of A. phalloides do not occur. This result is the main reason why we decided to carry out the second experiment using the South European records to project into North America. The second experiment helped us to confirm that the projected climatic suitability in warm areas of Mexico and the United States originates from populations distributed in the south of Europe (Fig. 2).

To have a better insight of a more realistic risk of A. phalloides expanding its distribution to Mexico we carried out a third experiment that used only North American records to model climatic suitability (Fig. 2). In this experiment, the coincidence of climatic suitability with pine and oak forests is very clear, but there are areas with this vegetation type in which the species has not yet been registered. Very scarce climatic suitability could be found in Mexico with this experiment, significantly reducing the risk to local communities except for Cochimíes and Kumiais, which inhabit a small area in the north of Baja California and hardly consume WEM.

3.2 Amanita phalloides under climate change scenarios

When using all European records to carry out the modeling, for the 2050 scenario (SSP5-8.5) the suitability decreases only -4% (52% of climatic suitability) in the United States, while in Canada the suitability is reduced -11.3% (3.7% of climatic suitability) (Fig. 4). In contrast, climatic suitability in Mexico increases up to +46% (80% of climatic suitability). In other words, the climatic suitability of A. phalloides is expected to reach further south (warmer areas). Under scenario 2070, A. phalloides will further increase its climatic suitability in the Unites States and Mexico by +2 (54% of climatic suitability) and +4% (84% of climatic suitability), respectively. In Canada the climatic suitability remains the same in comparison to 2050. When overlaying the distribution of indigenous communities in Mexico with the future climatic suitability of A. phalloides, we found that in comparison to the baseline scenario there is a reduction in co-occurrence in scenario 2050 and an increase of co-occurrence in 2070 (64 and 70%, respectively). These climatic suitability increases under climate change scenarios (2050 and 2070) is not reflected in the percentages of area with low (climatic suitability) and high risk (climatic and biological suitability) (Table III). The percentage of low-risk area is lower in 2050 than under the baseline scenario and higher under the 2070 scenario. A lower percentage of high-risk area is expected to co-occur with indigenous communities (~3 and 2% for 2050 and 2070, respectively) than under the baseline scenario. Still, a high number of communities are expected to be affected: Tarahumaras, Tepehuanes, Otomíes, Nahuas, and Pames.

Fig. 4 Climatic suitability of Amanita phalloides under the baseline scenario and climate change scenarios (2050 and 2070) using European records. Overlap between A. phalloides potential distribution and the distribution of its vegetation hosts in North America and indigenous territories in Mexico. Grey: potential distribution of A. phalloides; green: Potential distribution of its hosts; pink: indigenous territories. 

When using North American records to carry out the modeling, we found that under climate change scenarios there is almost no climatic suitability in Mexico (Fig. 5), while under the baseline scenario climatic suitability remains only in a small portion of the Baja Californian Peninsula, where indigenous communities are scarce but still exist (Cochimíes and Kumiais).

Fig. 5 Amanita phalloides climatic suitability maps under climate change scenarios using just North American records. The suitability in these maps shows that current Amanita phalloides’ populations in the United States are expected to have low climatic suitability in Mexico under the climate change scenarios used. 

4.1 Historic and observed distribution of Amanita phalloides

Amanita phalloides records, which were revised and visualized from 1950 to 2019, show that with time this mushroom has been recorded in a wider area through Europe and North America. In Europe, A. phalloides is strongly associated with oaks and species from the Fagaceae family. Between 1950 and 1975, records were concentrated in the northern region of the European continent, where mushrooms have been intensively studied (^{Wollan et al., 2008}). In more recent years (1976-2000), A. phalloides was registered in GBIF for the first time on the east and south regions of Europe. This lack of A. phalloides registers in GBIF prior to 1980 from southern Europe might be a database bias since the species has been distributed in Spain, Portugal, and Italy for decades (F. Martínez, A. Rinaldi, S. Goncalvez com. pers.). Experts also agree that the species has become more abundant in southern Europe, suggesting that thermophilic genotypes adapted to the Mediterranean climate.

In the United States, during the period from 2001 to 2019, there was a great increment of A. phalloides distribution across the western coast, and associations with new endemic hosts were registered (^{Wolfe and Pringle, 2012}). The fact that this mushroom can establish ectomycorrhizal associations with new hosts unveils that its distribution is not necessarily restricted to its known hosts. The association with new hosts in introduction areas has already been reported in other species from the same genus such as Amanita muscaria (^{Vellinga et al., 2009}). In areas where oaks are frequent, they remain as the main host of the species (^{Wolfe and Pringle, 2012}). In the east coast, A. phalloides is restricted to parks and pine-forest plantations. In this geographic area, the species is more associated to pines than to oaks. Even when changes of hosts have been registered in the United States, there is no genetic isolation from European populations, confirming that genetic divergence is recent (^{Pringle et al., 2009}).

4.2 Potential distribution of Amanita phalloides

The modeling suggests important differences between the observed distribution of A. phalloides in the North American range and the potential distribution under baseline and climate change scenarios. Our results reveal that, when all European records are considered, the species finds climatic suitability in Canada, the United States, and Mexico. Our results also show that there are geographic areas with climatic suitability that do not co-occur with vegetation types with which the species has been historically associated.

After obtaining the results of the second experiment (modeling with European southern records) it is evident that these southern populations mainly determine climatic suitability in Mexico. In contrast, the third experiment (using North American records) suggests that, until now, populations present in the United States do not find climatic suitability in most of Mexico under current and climate change scenarios. One explanation is that populations in North America are locally adapted and, consequently, they have their own climatic niche that cannot be found in Mexico. The ecological niche, and consequently the climatic niche, is a species’ ecological characteristic, but there are studies that suggest different ecological niches for different taxonomic levels beneath the species (^{Nakazawa et al. 2010}; ^{Ureta et al. 2012}).

The other explanation would be that the territorial expansion of the species in North America is due to a slow dispersal over short distances. Therefore, populations in the USA have not yet faced the climates of Mexico, although they are already very close to the border. However, it is a fact that the species niche is quite broad if all geographic records are considered and, until now, there is no genetic isolation between North American and European populations (^{Pringle et al., 2009}). So, even when the short-term risk might be low, there is still a risk of A. phalloides intoxication for Mexican indigenous communities because of North American approaching populations. If the precautionary principle is considered, we should warn local authorities and communities in sites where climatic suitability (using European records), vegetation hosts, and indigenous communities overlay. A priority area for this warning campaign should be Chihuahua in northern Mexico. In this region, the original Tarahumara group subsists in precarious conditions and during the rainy season they supplement their diet with WEM (^{Moreno-Fuentes et al., 2004}).

As an ectomycorrhizal species, it could be thought that A. phalloides distribution should be restricted to the distribution of its historical hosts, but it has already been demonstrated that this species is capable of adapting to new hosts in a few years (^{Vellinga et al., 2009}). Our results show that the species has an enormous climatic potential if we consider its entire distribution (ecological niche), but also that there are populations with different climatic thresholds (i.e., European southern populations or populations in the United States). If populations from southern Europe reach the American continent, about 70% of the indigenous communities could be at risk of intoxication.

4.3 Potential distribution of Amanita phalloides under climate change scenarios

Amanita phalloides suitability under climate change scenarios is expected to decrease in the east and west coasts of the United States and Canada, at the same time that the climatic suitability increases in the northwest, southeast and center of Mexico. This increase in Mexico’s climatic suitability represents a higher risk in terms of food security for local communities given that A. phalloides is morphological similar to other WEM. Up to 23 species of Amanita are consumed traditionally in Mexico (^{Garibay-Orijel and Ruan-Soto, 2014}); of these, the most similar to A. phalloides are the Amanita caesarea complex (A. arkansana, A. bassi, A. caesarea s.l., A. calyptratoides, A. calyptroderma, A. hayalyuy, A. hemibapha, A. jacksonii, A. laurae, A. tecomate, and A. yema), consumed by the 17 indigenous groups; and the Amanita rubescens complex (A. novinupta, A. amerirubescens, and A. rubescens s.l.), consumed by Nahuas, Totonacos, Tarahumaras, Tzeltales, Tzotziles, and Tojolabales. Amanita phalloides can also be confused with other edible mushrooms as Armillaria mellea s.l., consumed by Huicholes, Mazatecos, Tzeltales, and Nahuas (Martínez Mart et al., 1983; ^{Villarreal and Pérez-Moreno, 1989}; ^{Moreno-Fuentes et al., 2004}; ^{Garibay-Orijel et al., 2006}; ^{Robles-Porras et al., 2006}; ^{Lampman, 2007}; ^{Sheppard et al., 2008}; ^{Haro-Luna et al., 2019}; ^{López-García et al., 2020}; ^{Ramírez-Terrazo et al., 2021}).

Even when mushroom intoxications are not one of the main causes of deaths in rural communities in Mexico, they represent a health problem that affects ancestral cultural practices and local economies (^{Ruan-Soto et al., 2007}). Mushroom intoxications can happen by a lack of experience in recognizing edible from toxic mushrooms (^{Guzmán, 1987}). Although people that traditionally eat them have a strong mycological knowledge, it is true that species are changing their distribution due to climate change and it is possible that new toxic species appear in areas where people is not used to them.

The expected shift in A. phalloides’ distribution is towards southern areas in contrast to the general tendency of temperate species moving north to cooler areas. In most studies, the pattern is that temperate species will shift their distribution to temperate areas in northern latitudes and higher altitudes (^{Araújo and Rahbeck, 2006}; ^{Urban, 2018}; ^{Mamet et al., 2019}; ^{Osland et al., 2021}). Our results are based on observed records in which the species has recently been recorded abundantly in southern Europe. However, if only the North American populations are used in the modeling, the risk for local communities is lower. This result makes us highlight that when exotic populations are successful, they adapt locally and change their ecological requirements (^{Fitzpatrick and Keller, 2015}). This gives a strong argument to carry out ecological niche modeling with taxonomic entities below the species level (^{Ureta et al., 2012}). It was very helpful to identify in this study populations that represent a higher risk in Mexico under current and climate change scenarios.

Finally, it is important to recognize that model transferring to other climatic scenarios and geographic areas comes with uncertainty (see Table SI and Figs. S1-S4 in the supplementary material).The variables in the baseline scenario exhibit a range of values that can be exceeded in other scenarios or periods. This situation diminishes the model performance. Also, new variables’ combinations can be present, also diminishing the predictive power of the model. Consequently, it is important to identify geographic areas where uncertainty is greater and be cautious about the results (^{Mesgaran et al., 2014}).