Skip to main content

The two chytrid pathogens of amphibians in Eurasia—climatic niches and future expansion

Abstract

Background

Climate affects the thermal adaptation and distribution of hosts, and drives the spread of Chytridiomycosis—a keratin-associated infectious disease of amphibians caused by the sister pathogens Batrachochytrium dendrobatidi (Bd) and B. salamandrivorans (Bsal). We focus on their climate-pathogen relationships in Eurasia, the only region where their geographical distributions overlap. Eurasia harbours invaded and native areas of both pathogens and the natural habitats where they co-exist, making it an ideal region to examine their environmental niche correlations. Our understanding of how climate change will affect their distribution is broadened by the differences in climate correlates and niche characteristics between Bd and Bsal in Asia and Europe. This knowledge has potential conservation implications, informing future spread of the disease in different regions.

Results

We quantified the environmental niche overlap between Bd and Bsal in Eurasia using niche analyses. Results revealed partial overlap in the niche with a unique 4% of non-overlapping values for Bsal, suggesting segregation along certain climate axes. Bd tolerates higher temperature fluctuations, while Bsal requires more stable, lower temperature and wetter conditions. Projections of their Realized Climatic Niches (RCNs) to future conditions show a larger expansion of suitable ranges (SRs) for Bd compared to Bsal in both Asia and Europe, with their centroids shifting in different directions. Notably, both pathogens' highly suitable areas in Asia are expected to shrink significantly, especially under the extreme climate scenarios. In Europe, they are expected to expand significantly.

Conclusions

Climate change will impact or increase disease risk to amphibian hosts, particularly in Europe. Given the shared niche space of the two pathogens across available climate gradients, as has already been witnessed in Eurasia with an increased range expansion and niche overlap due to climate change, we expect that regions where Bsal is currently absent but salamanders are present, and where Bd is already prevalent, may be conducive for the spread of Bsal.

Peer Review reports

Background

Global amphibian populations are continuing to decline, leading to species extinction, and the situation is predicted to worsen in the future [1,2,3]. These declines are the result of several factors, such as habitat change [4], pollutants and pesticides [2, 5], climate change [6, 7], and infectious diseases [8,9,10]. Climate has an effect on wildlife, both directly (e.g., thermal adaptation and species distribution), and indirectly; it is a key driver of disease spread, as it regulates temperature and humidity regimes [11, 12]. It is also indicated that the climate change may alter dynamics of species susceptibility by breaking the balance or long-term coexistence between hosts and pathogens, triggering disease shifts [13,14,15]. Recent studies have predicted that climate change will likely increase the prevalence of infectious pathogens, thus posing a risk of disease while affecting new hosts [16]. Therefore, it is essential to understand the ecological niche correlates of amphibian pathogens and their dispersal patterns under climate variations, to facilitate amphibian conservation.

Chytridiomycosis has caused the greatest amphibian population decline and species extinction in the world, posing an ongoing threat to biodiversity [8, 17]. The disease is caused by two sister chytrid pathogens, Batrachochytrium dendrobatidi (Bd) [18, 19] and B. salamandrivorans (Bsal) [20], which has been thought to originated in East Asia during the late Cretaceous or early Paleogene [17, 21]. Asian hosts are considered resistant to Chytridiomycosis due to their long interactions with the pathogens [22, 23], yet they have spread to other regions such as Europe, the Americas and Australia [21, 24], affecting over 500 species of anurans and urodelans [21, 25].

Bd has a broader worldwide distribution than Bsal, which has not been found in the Americas but is relatively restricted in Eurasia (Asia and Europe) [26, 27]. Notably, both pathogens co-occur in natural habitats of East and Southeast Asia and Western Europe [26, 28, 29]. In Asia, seventeen countries have reported the presence of Bd [30,31,32], with four also reporting Bsal [17, 28, 33]. In Europe, Bd has been detected in about 27 countries (amphibian disease web portal: https://amphibiandisease.org), and Bsal in four (http://bsaleurope.com/), with a more limited distribution.

In Eurasia, the pathogens have a considerable degree of spatial overlap, but there are also regions where they do not overlap [28, 29]. Complex biotic and abiotic factors and the differences in ecological requirements of the pathogens are the likely cause for these non-overlapping regions. Abiotic factors, in particular thermal parameters, likely play a role; laboratory studies suggest Bd prefers slightly higher temperatures (17–25℃) compared to much cooler conditions for Bsal (10–15℃) [20, 34]. Bd is a host-generalist that infects all orders of amphibians [30], while Bsal is a host-specialist, mainly infecting Urodela (salamanders and newts). However, the realized climatic niches (RCNs) of the two pathogens have not been comprehensively compared. Understanding the differences between their climate niches can help predict how they will respond to climate change, allowing mitigation actions to be taken to control their spread.

Climate warming is presumably shaping the potential suitable habitats of the two pathogens, Bd and Bsal. Studies for Bd based on previous versions of climate models suggest climate change will likely enhance expansion dynamics in temperate and cooler regions [16, 35], but Bsal's future range expansion is still poorly understood. Updated versions of the Coupled Model Inter-comparison Project Phase 6 (CMIP6) climate models putatively allow more accurate predictions of species distribution when compared to the previous version [36,37,38]. Given that the region of our interest includes both native and invaded areas, Eurasia is an ideal region to compare the shared environmental niche correlations between Bd and Bsal and evaluate how the two pathogens respond to climate change on a larger scale. In addition, the dissimilar climatic conditions in Asia and Europe, as well as the climatic niche requirements of Bd and Bsal, coupled with various climate scenarios, may lead to divergent dispersal patterns in the future. These variations could have critical implications for amphibian conservation, especially in areas with a high diversity of urodelans, where Bsal has not yet established itself.

We aimed to compare the ecological niches of Bd and Bsal in Eurasia and anticipate their range shifts in response to climate change. Patterns to formulate a hypothesis included: 1) differences of host features of both pathogens [8, 39], 2) differential environmental conditions in native and invaded areas in Eurasia [17, 21] and, 3) distributional overlap of pathogens in Asia and Europe [17, 28, 40]. Based on the observed patterns, we hypothesized that the two pathogens have subtle differences in their RCNs to occupy respective habitats, with Bd having a broader niche range. The hypothesis if true would predict that the two pathogens would respond differently to future climate changes and range shifts between them would differ under various climatic scenarios.

Results

Realized niche overlaps between Bd and Bsal

The multivariate principal component analysis comparing the niches of Bd and Bsal in Eurasia (Fig. 1) revealed that the two pathogens exhibited environmental niche correlations. They occupied different niche spaces, and segregated based on factors such as mean diurnal range (Bio2), precipitation of the driest month (Bio14), precipitation seasonality (Bio15), precipitation of warmest quarter (Bio18) and elevation.

Fig. 1
figure 1

Distribution of Bd and Bsal and niche overlap analysis. Localities positive for Bd (green dots) and Bsal (red dots) on a background of altitudinal gradients are shown in a. The contributions of environmental variables to the first two PC axes are presented in b (standard definitions for each variable provided in Additional file 1: Table S1). A visualization of occupied and available niche space for the two pathogens is provided in c, based on the first two axes of a principal component analysis. Solid and dashed contour lines represent 100% and 50% of the available environment in Eurasia, respectively. The green, red and blue filled areas indicate the occupied niches of Bd, Bsal and overlap, respectively, with deeper colors representing higher overlap. Histograms d-e display the expected distribution of Schoener's D index values for the niche equivalence and niche similarity tests, with red diamonds indicating the actual observed niche overlap

However, the highest proportion of niche overlap for Bd and Bsal occurred along PC1 (Fig. 1c) with moderate values of D = 0.4365 and I = 0.6269. The niche equivalency and similarity between the pathogens were rejected due to subtle differences in the niches (Fig. 1d-e, Fig. 2). The niche stability value (environmental niche occupied by both Bsal and Bd) was high (0.9559) indicating that a large proportion of the niche space overlapping. Out of the total niche space, Bd occupied 0.2657 that did not overlap with Bsal, whereas Bsal occupied a distinct niche space of 0.0441.

Fig. 2
figure 2

Density distributions of Bd and Bsal occurrences along univariate axes. The hypothesis of niche equivalency is rejected when density distributions (Bd: red solid line and Bsal: green dotted line) exceed the estimated confidence limits (blue bands)

Univariate comparison of RCNs demonstrated varying degrees of niche overlap and differences between the two pathogens across all variables (Fig. 2). Notably, Bd had a broader diurnal range (Bio2) and mean temperature in driest quarter (Bio9) than Bsal. Distribution trends in temperature parameters suggested Bd could tolerate higher temperature fluctuations while Bsal had a relatively stable temperature regime. The comparison along the precipitation of the driest month variable highlighted Bsal's preference for slightly wetter conditions than Bd.

Performance of MaxEnt models and variable contributions

The MaxEnt models of the two pathogens recorded high AUC and TSS scores; Bd's AUC was 0.948 ± 0.002 and TSS was 0.795 ± 0.019, and Bsal's AUC and TSS were 0.984 ± 0.006 and 0.899 ± 0.035 respectively. This demonstrates the models' robustness and adequacy for predicting the climate suitability for Bd and Bsal.

Analysis of MaxEnt (Table 1; Additional file 2: Fig. S1) indicated variations in major bioclimatic variables between the two pathogens. Temperature-related parameters were the major contributors in the Bd model, with similar contributions from precipitation-related variables in the Bsal model.

Table 1 Contributions of bioclimatic variables to the MaxEnt models

Centroid distribution shifts of the two pathogens

The predicted future distribution changes reflected the differences in climatic niches of the two pathogens. The shifting of centroids in suitable ranges (SRs) was larger for Bd than Bsal, and larger centroid shifts occurred under SSP585 climate scenario rather than SSP245 (Fig. 3). These patterns held for both Asia and Europe. Additionally, the centroids of both pathogens in Europe shifted towards the east of their current ranges; the centroid of Bd in Asia moved towards northern latitudes (Fig. 3a), while that of Bsal was expected to move towards East Asia, suggesting they may behave differently in Asia in response to climate change (Fig. 3b).

Fig. 3
figure 3

Current climatic suitability and centroid shifts of Bd and Bsal. The current potential distribution (map) and shifts in centroid distribution (the boxes above maps) in Europe and Asia under future time periods (2041–2060, 2061–2080, 2081–2100) based on SSP245 and SSP585 climate scenarios for Bd a Bsal b

Changes in suitable habitats

Our analyses suggest that the overall geographic range size of Bd and Bsal would gradually increase in Eurasia during 2041–2060, 2061–2080 and 2081–2100 periods (Additional file 1: Table S2). In Europe, there was a marked expansion in suitable areas for both pathogens compared to their range size in Asia. In Asia, there was a substantial contraction in their range size in the southern regions across lower latitudes (Fig. 4).

Fig. 4
figure 4

Future range changes of Bd a and Bsal b in Europe and Asia. The changes in distributions are depicted for three future periods (2041–2060, 2061–2080, and 2081–2100) with two climate scenarios (SSP245 and SSP585)

In Europe, the current distributions of both pathogens are predicted to be concentrated in temperate and cold climate regions with weak dry seasons and warm summers (Fig. 3; Additional file 2: Fig. S2). The ensemble model identified more areas in Eastern Europe as potential suitable habitats for both pathogens compared to the single MaxEnt model (Additional file 2: Fig. S3). A distinct expansion in 2081–2100 period under SSP585 climate scenario was shown for both pathogens (Fig. 4, Additional file 2: Fig. S4). However, compared to Bsal, Bd is expected to have greater range shifts, and most regions would likely hold optimal climatic conditions for a long period (Fig. 4).

In Asia, current distributions for both pathogens were across eastern mountainous coastal areas and high altitudes of temperate climate regions (Fig. 3; Additional file 2: Fig. S2). The ensemble models predicted more suitable habitats in Southeast Asia (Additional file 2: Fig. S3). Future predictions suggest that their ranges would expand into the Hengduan Mountains, North China Plain, and coastal mountainous regions in China. Bsal's future distribution was fragmented and concentrated in higher altitudes, with Southern and Central China becoming more suitable than Bd (Fig. 4).

A major proportion of range expansion for both pathogens (Fig. 4) was reflected more in areas currently classified as low and/or moderately suitable. We therefore examined changes in ‘highly suitable areas’ in Asia and Europe (Fig. 5). Our results showed that the highly suitable areas for Bd in Asia would shrink (Additional file 2: Fig. S5) despite shifting towards East Asia (Fig. 5a). Bsal's extension of highly suitable habitats would become less suitable (Fig. 5b) despite the overall expansion shown in Fig. 4. In Europe, Bd was projected to have low and moderate suitability regions in Eastern Europe, while Western and Central Europe were projected to maintain high suitability (Fig. 5a). Bsal was expected to have a wider extension in low suitability habitats in the future, while highly and moderately suitable areas were projected to decrease (Fig. 5b, Additional file 2: Fig. S5).

Fig. 5
figure 5

Future climatic suitability of Bd a and Bsal b. Different colors represent classes of future habitat suitability; low suitability, moderate suitability and high suitability. Suitability categories are looked under two climate scenarios (SSP245 and SSP585) across future three periods (2041–2060, 2061–2080, 2081–2100)

Future niche variations between Bd and Bsal

The niche breadth (B value) of both pathogens is likely to expand under future climate scenarios (Table 2). Specifically, Bd's niche breadth is expected to be larger than that of Bsal. The niche overlap between the two pathogens was projected to increase during the three future periods considered. Range overlap was also predicted to gradually increase under SSP245 climate scenario.

Table 2 Estimated niche breadth (B), niche overlap (D and I) and range overlap. The niche values of Bd and Bsal were calculated based on the maps of predicted habitat suitability under SSP245 and SSP585 climate scenarios. Range overlap was estimated based on the threshold value (0.6) of high suitability

Discussion

Eurasia is an important region for distribution of Bd and Bsal [17, 21], as they co-exist there under suitable climatic conditions [28, 29]. We studied their RCNs and used future climate models to assess their habitat suitability and niche variations. We also considered how range shifts might differ between Asia and Europe given the different climatic factors in each region.

Niche overlap and climate preferences between Bd and Bsal

Our results suggest that Bd and Bsal share a similar environmental niche space. Their overlapping climatic space is associated with mean temperatures of wettest and driest quarters (Fig. 1b and c). Bsal showed highest occurrence densities than Bd in all variables except for mean temperature of wettest quarter (Fig. 2), indicating that its optimum conditions for Bsal may be suitable for Bd as well. Notably, Bsal is able to occupy under highly humid conditions as reflected by higher occurrence densities under precipitation of driest month and precipitation of warmest quarter (Fig. 2). Therefore, Bsal's environmental associations and high virulence [20, 41] suggest it may spread to Bd-present habitats with salamanders in the future. Overall, wet conditions appear to shape Bsal's occurrence, as shown by other studies [39, 42].

Current suitable habitats and future range shifts under scenarios of climate change

MaxEnt-based models agree with most prior studies, but differ slightly in Bd and Bsal's distributions. Bd is predicted to spread across Europe, while past models suggested less suitability in Northern-Central-Eastern Europe [43,44,45,46], the increased number of occurrence records added in the current analysis may have resulted in this difference. Bsal's model predicts a higher distribution in Western-Central Europe, while region-scale studies predicted a more restricted area in southern Germany [40, 47, 48]. For Asia, both Bd and Bsal have potential distributions concentrated in Southeast-East Asia (Fig. 3) but previous studies suggested Central Asia, South Korea [35] and Kamchatka [43,44,45] as different SRs for Bd. Under future climate scenarios, range contractions of Bd are predicted in southern China and Southeast Asia due to warmer temperatures and excessive rainfall [49]. Previous estimates suggested Bd would decrease and be restricted in Europe, but our study found a broader distribution in Europe. These discrepancies in future predictions could be attributed to the use of varying spatial ranges used in different models [30, 50]. Variations in species distribution models are expected due to several factors, such as differences in spatial resolutions, predictor variables used in the analyses, and the number of occurrence records utilized. Furthermore, some studies utilize more predictors, such as amphibian species richness, which can also contribute to differences in predictions [35, 45].

Bd is likely to track climate change better than Bsal, due to its generalist nature [51], which could explain the differences in range expansion between the two pathogens (Figs. 4 and 5). Many of Bsal's habitats will become less suitable due to climate change, but habitats for Bsal will remain highly suitable in Germany and Belgium. The centroid distribution shifts of Bd and Bsal in Asia (Fig. 3) suggest that Bd will move northward towards optimal temperature conditions and Bsal towards wetter and cooler conditions in East Asia, as will their hosts.

Implications in amphibian conservation

Climate change could impact the range expansion and niche breadth of Bd and Bsal, potentially threatening vulnerable amphibians across Eurasia (Fig. 4 and Table 2). In Asia, amphibians may shift towards northern regions [52, 53], and Bd may also expand its range. As global warming continues, amphibian hosts adapting to cooler climates may be at increased infection risk [54] as temperature can influence both the kinetics of host–pathogen interactions and act as a host stressor [11]. In case where certain areas with suitable climates remain unaffected by chytrid pathogens, such as Switzerland in Europe and Himalaya highlands in Asia, should be monitored on a long-term basis. If they continue to remain free of pathogen infiltration over time [55], it would call for further studies to understand why such patterns are being observed.

Conclusions

Our analysis across Eurasia revealed shared environmental conditions between the two chytrid pathogens and moderate niche overlap, with Bd having a broader RCN than Bsal. Our study identifies areas of potential high-priority for monitoring the spread of the chytrid fungi.

Methods

Study area

We delimited a large area covering Eurasia (Xmin = 14.28S and Xmax = 180.00; Ymin = 13.00W and Ymin = 90.00E, WGS84) to capture all Bd and Bsal records in the region. This includes Asia (42,968,070.93 km2) and Europe (10,631,747.53 km2). We downloaded pathogen range shapefiles from http://tapiquen-sig.jimdofree.com for use in subsequent analyses.

Occurrence records

To model the distribution of Bd and Bsal, we obtained occurrence records of their amphibian hosts in the wild (excluding those in captivity). Bd coordinates were from the Bd-Maps Legacy Project (https://amphibiandisease.org/projects/, accessed Jun 15 2021) and Bsal records from relevant literature [17, 24, 28, 33, 40, 56, 57]. We applied the nearest neighbor joining index to reduce the effect of spatial autocorrelation due to sampling bias [58]. The final positive records consisted of 446 Bd and 57 Bsal observations (Fig. 1; Additional file 1: Table S3).

Environmental variables

We initially downloaded 19 bioclimatic variables and the elevation variable from WorldClim (http://www.worldclim.org/) at a spatial resolution of 2.5 arc minutes (~ 5 × 5 km2) for the period 1970–2000 [59]. To account for multicollinearity [60], we calculated correlations between variables in ENMTools [61] and visualized them in a Pearson correlation matrix (Additional file 2: Fig. S6). We then selected six non-correlated bioclimatic variables with Pearson’s r >|0.7|, plus elevation, for further analysis [62] (Additional file 1: Table S1). These variables were used to characterize climate niches, analyze niche overlap and predict current distribution.

To reduce uncertainty in predictions of future distribution patterns, we downloaded six bioclimatic layers representative of three future time periods (2041–2060, 2061–2080, and 2081–2020). We used the latest CMIP6 models with shared socioeconomic pathways (SSPs) to enhance confidence of the outcomes [36]. Different Global Climate Models (GCMs) can affect the predicted distributions of species differently [63, 64], so we used multiple GCMs (BCC-CSM2-MR, CNRM-CW6-1, CNRM-ESM2-1, CanESM5, IPSL-CM6A-LR, MIROC-ES2L, MIROC and MRI-ESM2-0) for our models. We employed two representative concentration pathways (SSP245 and SSP585, resulting in 2100 radiative forcing levels of 4.5 W/m2 and 8.5 W/m2, respectively) [37] to represent moderate and pessimistic climate change scenarios [38].

Niche overlap analyses

We implemented a quantitative approach to evaluate environmental niche correlates between Bd and Bsal based on ecological factors associated with pathogen presence. We performed the realized niche analysis of the two pathogens using “ecospat” package [65], where the first two axes of ordination framework were gridded into 100 × 100 cells depicting the climatic space. This permitted us to test whether differences of observed niches result from partial niche filling or niche expansion, as environmental availability is accounted for during ordination. To quantitatively compare the niches of the two pathogens, we used the principal component analysis to extract environmental variables (Additional file 1: Table S1) of each pathogen and used a smooth kernel density function to determine the ‘smoothed’ density of occurrences in each cell in the niche space [66]. We used Schoener's D metric and Hellinger's-based I index to quantify niche overlap and 1000 repetitions for generating a null distribution to analyze niche equivalence and similarity [66,67,68]. The scores for D and I index range from 0 (no overlap) to 1 (complete overlap). Niche stability, expansion, and unfilling between Bd and Bsal were calculated based on each cell in the ordination using “ecospat” package [65]. We performed a bootstrap hypothesis test of equality, based on a single environmental variable occupied by each pathogen using “sm” package [69], to compare their realized niches in a univariate space.

Distribution modeling

We used MaxEnt-based species distribution models (SDMs) with a presence-only method, ensemble techniques, and an analysis of centroid shift and climatic suitability change for each pathogen in both Europe and Asia.

We predicted and quantitatively calculated the extent and changes of habitat suitability for the two pathogens, Bd and Bsal, under current climate and two future climatic scenarios (SSP245 and SSP585). The maximum entropy algorithm (MaxEnt) was selected to create habitat suitability models based on their occurrence records in Eurasia, as it is particularly suited for presence-background algorithms by comparing the available environmental conditions in the background with the conditions favorable to species occurrences [70,71,72].

To reduce model over-fitting and model complexity in trained current and future models [60], we performed relevant analysis using the “kuenm” package [73] and generated all possible combinations of regularization multipliers (RM: ranging from 0.0 to 4.0 in 0.1 steps) with all combinations by including the linear (L), quadratic (Q), product (P), threshold (T), and hinge (H) features. Among all created models, the model with the lowest AICc score was selected as the optimal model. Therefore, the model with the combination of QPTH feature and regularization multiplier of 1.2 was used to predict Bd distribution and the model with the combination of T feature and regularization multiplier of 1.0 was used to predict Bsal distribution. We conducted a MaxEnt analysis by using the maximum training sensitivity plus the specificity threshold (MTSS), with a training set consisting of 25% occurrence coordinates, ten replicates, and cloglog outputs [74, 75]. The Area Under the receiver operating characteristic Curve (AUC) and true skill statistic (TSS) were used to evaluate model performance [70, 76]. The model performance was evaluated as ‘good’ when the value of AUC was above 0.9 [77] and the value range of TSS between 0.6 and 0.8 was regarded as useful while it was considered good or excellent when value was over 0.8 [76].

We converted the above results from MaxEnt models into binary habitat suitability (suitable habitats and unsuitable habitats) to reclassify it into four categories of climatic suitability based on the values of MTSS [78]: unsuitability (< MTSS), low suitability (MTSS-0.4), moderate suitability (0.4–0.6), and high suitability (> 0.6) based on continuous suitability from the outputs under the background of current climate and future climatic scenarios. To examine the distribution trends for both pathogen in Eurasia, the size changes in SRs were calculated as follows [79]:

$$\mathrm{Range}\;\mathrm{size}\;\mathrm{change}=\frac{\mathrm{future}\;\mathrm{suitable}\;\mathrm{area}-\mathrm{current}\;\mathrm{suitable}\;\mathrm{area}}{\mathrm{current}\;\mathrm{suitable}\;\mathrm{area}}\times100\%$$

To compare the differences and changes in suitable habitats of Bd and Bsal between Asia and Europe, we divided the maps of habitat suitability and binary suitable habitats of Eurasia into the two regions (Europe and Asia). Binary suitable habitats in each region were used to evaluate changes in centroid positions and distributions, and these processes were completed using SDMtoolbox [80]. The area changes in reclassified habitat suitability for each region were calculated using spatial analysis tools in ArcMap v10.4 (ESRI).

To better understand the changes in the ecological niches of Bd and Bsal pathogens over time, we measured the extent of overlap between their niches and ranges. For this, we estimated their niche breadth by utilizing the inverse concentration metric of mean 'B' in ENMTools, which was based on the suitability space (cells) generated from Maxent models for both pathogens [61, 81]. Levin’s metric range was the same as D metric and I index, which ranged from 0 to 1, where 0 represents minimum niche breadth and 1 represents maximum niche breadth [82].

While MaxEnt is generally a solid approach, some authors criticize the single-algorithm strategy [83, 84]. Therefore in order to verify the robustness of our predictions we also performed an ensemble of multiple techniques in “Biomod2” package [85].We used eight algorithms, Generalized Linear Models (GLM), Generalized Boosting Models (GBM), Classification Tree Analysis (CTA), Surface Range Envelop (SRE), Flexible Discriminant Analysis (FDA), Multivariate Adaptive Regression Splines (MARS), Random Forest (RF) and Maximum Entropy (MAXENT.phillips), to evaluate their current distributions. In the modelling process, pathogen records were divided into the training data (75%) and test data (25%) and repeated 10 times; three different sets of pseudo-absences were created by randomly sampling from climatic space available within Eurasia. Out of the established 240 different models for each pathogen, 158 models of Bd had TSS scores > 0.6 (KAPPAaverage = 0.66, TSSaverage = 0.80, ROCaverage = 0.91), 181 models of Bsal had TSS scores > 0.6 (KAPPAaverage = 0.83, TSSaverage = 0.85, ROCaverage = 0.93). The final ensemble model was constructed to simulate and predict the suitable climatic areas at present based on the metric threshold (TSS > 0.6) [84].

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article and its additional files. All codes used in this study are publicly available in specific R packages as indicated in our paper.

Abbreviations

Bd :

Batrachochytrium dendrobatidis

Bsal :

Batrachochytrium salamandrivoransin

RCNs:

Realized Climatic Niches

SRs:

Suitable ranges

CMIP6:

Coupled Model Intercomparison Project Phase 6

Bio2:

Mean diurnal range

Bio8:

Mean temperature in wettest quarter

Bio9:

Mean temperature in driest quarter

Bio14:

Precipitation of driest month

Bio15:

Precipitation seasonality

Bio18:

Precipitation of warmest quarter

References

  1. International Union for Conservation of Nature (IUCN). The IUCN red list of threatened species. Version 2021–3. Available at: http://www.iucnredlist.org. 2021.

  2. Green DM, Lannoo MJ, Lesbarrères D, Muths E. Amphibian population declines: 30 years of progress in confronting a complex problem. Herpetologica. 2020;76(2):97–100.

    Article  Google Scholar 

  3. Ceballos G, Ehrlich PR, Raven PH. Vertebrates on the brink as indicators of biological annihilation and the sixth mass extinction. Proc Natl Acad Sci U S A. 2020;117(24):13596–602.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Cordier JM, Aguilar R, Lescano JN, Leynaud GC, Bonino A, Miloch D, et al. A global assessment of amphibian and reptile responses to land-use changes. Biol Conserv. 2021;253: 108863.

    Article  Google Scholar 

  5. Saaristo M, Brodin T, Balshine S, Bertram MG, Brooks BW, Ehlman SM, et al. Direct and indirect effects of chemical contaminants on the behaviour, ecology and evolution of wildlife. Proc Biol Sci. 1885;2018(285):20181297.

    Google Scholar 

  6. Li Y, Cohen JM, Rohr JR. Review and synthesis of the effects of climate change on amphibians. Integr Zool. 2013;8(2):145–61.

    Article  PubMed  Google Scholar 

  7. Araujo MB, Thuiller W, Pearson RG. Climate warming and the decline of amphibians and reptiles in Europe. J Biogeogr. 2006;33(10):1712–28.

    Article  Google Scholar 

  8. Scheele BC, Pasmans F, Skerratt LF, Berger L, Martel A, Beukema W, et al. Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity. Science. 2019;363(6434):1459.

    Article  CAS  PubMed  Google Scholar 

  9. Fisher MC, Henk DA, Briggs CJ, Brownstein JS, Madoff LC, McCraw SL, et al. Emerging fungal threats to animal, plant and ecosystem health. Nature. 2012;484(7393):186–94.

    Article  CAS  PubMed  Google Scholar 

  10. Daszak P, Cunningham AA, Hyatt AD. Emerging infectious diseases of wildlife-threats to biodiversity and human health. Science. 2000;287(5452):443–9.

    Article  CAS  PubMed  Google Scholar 

  11. Altizer S, Ostfeld RS, Johnson PTJ, Kutz S, Harvell CD. Climate change and infectious diseases: from evidence to a predictive framework. Science. 2013;341(6145):514–9.

    Article  CAS  PubMed  Google Scholar 

  12. Rohr JR, Raffel TR. Linking global climate and temperature variability to widespread amphibian declines putatively caused by disease. Proc Natl Acad Sci USA. 2010;107(18):8269–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Moura-Campos D, Greenspan SE, DiRenzo GV, Neely WJ, Toledo LF, Becker CG. Fungal disease cluster in tropical terrestrial frogs predicted by low rainfall. Biol Conserv. 2021;261: 109246.

    Article  Google Scholar 

  14. Cohen JM, Civitello DJ, Venesky MD, McMahon TA, Rohr JR. An interaction between climate change and infectious disease drove widespread amphibian declines. Glob Chang Biol. 2018;25(3):927–37.

    Article  PubMed  Google Scholar 

  15. Spitzen-van der Sluijs A, Canessa S, Martel A, Pasmans F. Fragile coexistence of a global chytrid pathogen with amphibian populations is mediated by environment and demography. Proc R Soc B. 1864;2017(284):20171444.

    Google Scholar 

  16. Cohen JM, Sauer EL, Santiago O, Spencer S, Rohr JR. Divergent impacts of warming weather on wildlife disease risk across climates. Science. 2020;370(6519):eabb1702.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Martel A, Blooi M, Adriaensen C, Van Rooij P, Beukema W, Fisher MC, et al. Recent introduction of a chytrid fungus endangers Western Palearctic salamanders. Science. 2014;346(6209):630.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Berger L, Speare R, Daszak P, Green DE, Cunningham AA, Goggin CL, et al. Chytridiomycosis causes amphibian mortality associated with population declines in the rain forests of Australia and Central America. Proc Natl Acad Sci USA. 1998;95(15):9031–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Longcore JE, Pessier AP, Nichols DK. Batrachochytrium dendrobatidis gen. et sp. nov., a chytrid pathogenic to amphibians. Mycologia. 1999;91(2):219–27.

    Article  Google Scholar 

  20. Martel A, Spitzen-van der Sluijs A, Blooi M, Bert W, Ducatelle R, Fisher MC, et al. Batrachochytrium salamandrivorans sp. nov. causes lethal chytridiomycosis in amphibians. Proc Natl Acad Sci USA. 2013;110(38):15325–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. O’Hanlon SJ, Rieux A, Farrer RA, Rosa GM, Waldman B, Bataille A, et al. Recent Asian origin of chytrid fungi causing global amphibian declines. Science. 2018;360(6389):621.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Swei A, Rowley JJ, Rodder D, Diesmos ML, Diesmos AC, Briggs CJ, et al. Is chytridiomycosis an emerging infectious disease in Asia? PLoS ONE. 2011;6(8): e23179.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Fu M, Waldman B. Ancestral chytrid pathogen remains hypervirulent following its long coevolution with amphibian hosts. Proc R Soc B. 1904;2019(286):20190833.

    Google Scholar 

  24. Spitzen-van der Sluijs A, Martel A, Asselberghs J, Bales EK, Beukema W, Bletz MC, et al. Expanding distribution of lethal amphibian fungus Batrachochytrium salamandrivoransin Europe. Emerg Infect Dis. 2016;22(7):1286–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Schmeller DS, Utzel R, Pasmans F, Martel A. Batrachochytrium salamandrivorans kills alpine newts (Ichthyosaura alpestris) in southernmost Germany. Salamandra. 2020;56(3):230–2.

    Google Scholar 

  26. Fisher MC, Garner TWJ. Chytrid fungi and global amphibian declines. Nat Rev Microbiol. 2020;18(6):332–43.

    Article  CAS  PubMed  Google Scholar 

  27. Castro Monzon F, Rödel M-O, Ruland F, Parra-Olea G, Jeschke JM. Batrachochytrium salamandrivorans’ amphibian host species and invasion range. EcoHealth. 2023;19(4):475–86.

    Article  PubMed Central  Google Scholar 

  28. Laking AE, Ngo HN, Pasmans F, Martel A, Nguyen TT. Batrachochytrium salamandrivorans is the predominant chytrid fungus in Vietnamese salamanders. Sci Rep. 2017;7(1):44443.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Lötters S, Wagner N, Kerres A, Vences M, Steinfartz S, Sabino-Pinto J, et al. First report of host co-infection of parasitic amphibian chytrid fungi. Salamandra. 2018;54(4):287–90.

    Google Scholar 

  30. Olson DH, Ronnenberg KL, Glidden CK, Christiansen KR, Blaustein AR. Global patterns of the fungal pathogen Batrachochytrium dendrobatidis support conservation Urgency. Front Vet Sci. 2021;8: 685877.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Rahman MM, Badhon MK, Salauddin M, Rabbe MF, Islam MS. Chytrid infection in Asia: how much do we know and what else do we need to know? Herpetol J. 2020;30(2):99–111.

    Article  Google Scholar 

  32. Rahman MM, Jahan H, Rabbe MF, Chakraborty M, Salauddin M. First detection of Batrachochytrium dendrobatidis in wild frogs from Bangladesh. EcoHealth. 2021;18(1):31–43.

    Article  PubMed  Google Scholar 

  33. Yuan ZY, Martel A, Wu J, Van Praet S, Canessa S, Pasmans F. Widespread occurrence of an emerging fungal pathogen in heavily traded Chinese urodelan species. Conserv Lett. 2018;11(4): e12436.

    Article  Google Scholar 

  34. Piotrowski JS, Annis SL, Longcore JE. Physiology of Batrachochytrium dendrobatidis, a chytrid pathogen of amphibians. Mycologia. 2004;96(1):9–15.

    Article  PubMed  Google Scholar 

  35. Xie GY, Olson DH, Blaustein AR. Projecting the global distribution of the emerging amphibian fungal pathogen, Batrachochytrium dendrobatidis, based on IPCC climate futures. PLoS ONE. 2016;11(8): e0160746.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Eyring V, Bony S, Meehl GA, Senior CA, Stevens B, Stouffer RJ, et al. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. Geosci Model Dev. 2016;9(5):1937–58.

    Article  Google Scholar 

  37. O’Neill BC, Tebaldi C, van Vuuren DP, Eyring V, Friedlingstein P, Hurtt G, et al. The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6. Geosci Model Dev. 2016;9(9):3461–82.

    Article  Google Scholar 

  38. Hausfather Z, Peters GP. Emissions-the ‘business as usual’ story is misleading. Nature. 2020;577(7792):618–20.

    Article  CAS  PubMed  Google Scholar 

  39. Stegen G, Pasmans F, Schmidt BR, Rouffaer LO, Van Praet S, Schaub M, et al. Drivers of salamander extirpation mediated by Batrachochytrium salamandrivorans. Nature. 2017;544(7650):353–6.

    Article  CAS  PubMed  Google Scholar 

  40. Lötters S, Wagner N, Albaladejo G, Boning P, Dalbeck L, Dussel H, et al. The amphibian pathogen Batrachochytrium salamandrivorans in the hotspot of its European invasive range: past-present-future. Salamandra. 2020;56(3):173–88.

    Google Scholar 

  41. Farrer RA, Martel A, Verbrugghe E, Abouelleil A, Ducatelle R, Longcore JE, et al. Genomic innovations linked to infection strategies across emerging pathogenic chytrid fungi. Nat Commun. 2017;8(1):14742.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Feldmeier S, Schefczyk L, Wagner N, Heinemann G, Veith M, Lötters S. Exploring the distribution of the spreading lethal salamander chytrid fungus in its invasive range in Europe – a macroecological approach. PLoS ONE. 2016;11(10): e0165682.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Ron SR. Predicting the distribution of the amphibian pathogen Batrachochytrium dendrobatidis in the New World. Biotropica. 2005;37(2):209–21.

    Article  Google Scholar 

  44. Rödder D, Kielgast J, Lötters S. Future potential distribution of the emerging amphibian chytrid fungus under anthropogenic climate change. Dis Aquat Org. 2010;92(3):201–7.

    Article  Google Scholar 

  45. Liu X, Rohr JR, Li YM. Climate, vegetation, introduced hosts and trade shape a global wildlife pandemic. Proc R Soc B. 2013;280(1753):20122506.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Olson DH, Aanensen DM, Ronnenberg KL, Powell CI, Walker SF, Bielby J, et al. Mapping the global emergence of Batrachochytrium dendrobatidis, the amphibian chytrid fungus. PLoS ONE. 2013;8(2): e56802.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Katz TS, Zellmer AJ. Comparison of model selection technique performance in predicting the spread of newly invasive species: a case study with Batrachochytrium salamandrivorans. Biol Invasions. 2018;20(8):2107–19.

    Article  Google Scholar 

  48. Beukema W, Martel A, Nguyen TT, Goka K, Schmeller DS, Yuan ZY, et al. Environmental context and differences between native and invasive observed niches of Batrachochytrium salamandrivorans affect invasion risk assessments in the Western Palaearctic. Divers Distrib. 2018;24(12):1788–801.

    Article  Google Scholar 

  49. Mahlstein I, Knutti R, Solomon S, Portmann RW. Early onset of significant local warming in low latitude countries. Environ Res Lett. 2011;6(3): 034009.

    Article  Google Scholar 

  50. Soberón J, Nakamura M. Niches and distributional areas: concepts, methods, and assumptions. Proc Natl Acad Sci USA. 2009;106(supplement_2):19644–50.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Rohr JR, Dobson AP, Johnson PTJ, Kilpatrick AM, Paull SH, Raffel TR, et al. Frontiers in climate change–disease research. Trends Ecol Evol. 2011;26(6):270–7.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Duan RY, Kong XQ, Huang MY, Varela S, Ji X. The potential effects of climate change on amphibian distribution, range fragmentation and turnover in China. PeerJ. 2016;4: e2185.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Mi C, Huettmann F, Li X, Jiang Z, Du W, Sun B. Effects of climate and human activity on the current distribution of amphibians in China. Conserv Biol. 2022:e13964.

  54. Cohen JM, Civitello DJ, Brace AJ, Feichtinger EM, Ortega CN, Richardson JC, et al. Spatial scale modulates the strength of ecological processes driving disease distributions. Proc Natl Acad Sci USA. 2016;113(24):E3359–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Ellepola G, Herath J, Dan S, Pie MR, Murray KA, Pethiyagoda R, et al. Climatic niche evolution of infectious diseases driving amphibian declines. bioRxiv. 2022.

  56. Dalbeck L, Düssel-Siebert H, Kerres A, Kirst K, Koch A, Lötters S, et al. Die salamanderpest und ihr Erreger Batrachochytrium salamandrivorans (Bsal): Aktueller Stand in Deutschland. Z Feldherpetol. 2018;25:1–22.

    Google Scholar 

  57. Martel A, Vila-Escale M, Fernández-Giberteau D, Martinez-Silvestre A, Canessa S, Van Praet S, et al. Integral chain management of wildlife diseases. Conserv Lett. 2020;13(2): e12707.

    Article  Google Scholar 

  58. Dormann CF, McPherson JM, Araújo MB, Bivand R, Bolliger J, Carl G, et al. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography. 2007;30(5):609–28.

    Article  Google Scholar 

  59. Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37(12):4302–15.

    Article  Google Scholar 

  60. Warren DL, Seifert SN. Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol Appl. 2011;21(2):335–42.

    Article  PubMed  Google Scholar 

  61. Warren DL, Glor RE, Turelli M. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography. 2010;33:607–11.

    Google Scholar 

  62. Elith J, Phillips SJ, Hastie T, Dudík M, Chee YE, Yates CJ. A statistical explanation of MaxEnt for ecologists. Divers Distrib. 2011;17(1):43–57.

    Article  Google Scholar 

  63. Heikkinen RK, Luoto M, Araújo MB, Virkkala R, Thuiller W, Sykes MT. Methods and uncertainties in bioclimatic envelope modelling under climate change. Prog Phys Geogr. 2006;30(6):751–77.

    Article  Google Scholar 

  64. Goberville E, Beaugrand G, Hautekèete NC, Piquot Y, Luczak C. Uncertainties in the projection of species distributions related to general circulation models. Ecol Evol. 2015;5(5):1100–16.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Di Cola V, Broennimann O, Petitpierre B, Breiner FT, D’Amen M, Randin C, et al. Ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography. 2017;40(6):774–87.

    Article  Google Scholar 

  66. Broennimann O, Fitzpatrick MC, Pearman PB, Petitpierre B, Pellissier L, Yoccoz NG, et al. Measuring ecological niche overlap from occurrence and spatial environmental data. Global Ecol Biogeogr. 2012;21(4):481–97.

    Article  Google Scholar 

  67. Schoener TW. Sizes of feeding territories among birds. Ecology. 1968;49(1):123–41.

    Article  Google Scholar 

  68. Warren DL, Glor RE, Turelli M. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution. 2008;62(11):2868–83.

    Article  PubMed  Google Scholar 

  69. Bowman AW, Azzalini A. R package sm: nonparametric smoothing methods (version 2.2–4). 2010.

  70. Elith J, H. Graham C, P. Anderson R, Dudík M, Ferrier S, Guisan A, et al. Novel methods improve prediction of species' distributions from occurrence data. Ecography. 2006;29(2):129–51.

  71. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Modell. 2006;190(3–4):231–59.

    Article  Google Scholar 

  72. Phillips SJ, Anderson RP, Dudík M, Schapire RE, Blair ME. Opening the black box: an open-source release of Maxent. Ecography. 2017;40(7):887–93.

    Article  Google Scholar 

  73. Cobos ME, Peterson AT, Barve N, Osorio-Olvera L. Kuenm: an R package for detailed development of ecological niche models using Maxent. PeerJ. 2019;7: e6281.

    Article  PubMed  PubMed Central  Google Scholar 

  74. Guisan A, Zimmermann NE. Predictive habitat distribution models in ecology. Ecol Modell. 2000;135(2–3):147–86.

    Article  Google Scholar 

  75. Liu C, White M, Newell G, Griffioen P. Species distribution modelling for conservation planning in Victoria. Australia Ecol Modell. 2013;249:68–74.

    Article  Google Scholar 

  76. Allouche O, Tsoar A, Kadmon R. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). J Appl Ecol. 2006;43(6):1223–32.

    Article  Google Scholar 

  77. Swets JA. Measuring the accuracy of diagnostic systems. Science. 1988;240(4857):1285–93.

    Article  CAS  PubMed  Google Scholar 

  78. Aidoo OF, Souza PGC, Silva RS, Santana PA, Picanço MC, Kyerematen R, et al. Climate-induced range shifts of invasive species (Diaphorina citri Kuwayama). Pest Manag Sci. 2022;78(6):2534–49.

    Article  CAS  PubMed  Google Scholar 

  79. Hu ZM, Zhang QS, Zhang J, Kass JM, Mammola S, Fresia P, et al. Intraspecific genetic variation matters when predicting seagrass distribution under climate change. Mol Ecol. 2021;30(15):3840–55.

    Article  PubMed  Google Scholar 

  80. Brown JL, Anderson B. SDMtoolbox: a python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. Methods Ecol Evol. 2014;5(7):694–700.

    Article  Google Scholar 

  81. Levins R. Evolution in changing environments: some theoretical explorations. United States: Princeton University Press; 1968.

    Book  Google Scholar 

  82. Mandle L, Warren DL, Hoffmann MH, Peterson AT, Schmitt J, von Wettberg EJ. Conclusions about niche expansion in introduced Impatiens walleriana populations depend on method of analysis. PLoS ONE. 2010;5(12): e15297.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Thuiller W. BIOMOD - optimizing predictions of species distributions and projecting potential future shifts under global change. Glob Chang Biol. 2003;9(10):1353–62.

    Article  Google Scholar 

  84. James TY, Toledo LF, Rödder D, Silva Leite D, Belasen AM, Betancourt-Román CM, et al. Disentangling host, pathogen, and environmental determinants of a recently emerged wildlife disease: lessons from the first 15 years of amphibian chytridiomycosis research. Ecol Evol. 2015;5(18):4079–97.

    Article  PubMed  PubMed Central  Google Scholar 

  85. Thuiller W GD, Gueguen M, Engler R, Breiner F, Lafourcade B. biomod2: ensemble platform for species distribution modeling. R package version 4.1–1. 2022.

Download references

Acknowledgements

We thank the Editor Dr. Duarte Goncalves and two anonymous reviewers for providing feedback to improve the clarity of the paper. We also thank Tai Gao and Christos Mammides for advice on data analyses.

Funding

Guangxi University Laboratory Startup Funding (MM), Chinese Scholarship Council (CSC) for graduate studies (JH, GE) for financial assistance. The funding bodies played no role in the design of the study and collection, analysis, interpretation of data, and in writing the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

DS, GE, JH and MM conceptualized the study. DS and GE designed methodology. All authors contributed in writing the original draft, reviewing and editing. MM provided resources and supervised the study.

Corresponding author

Correspondence to Madhava Meegaskumbura.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

Madhava Meegaskumbura is an editorial board member for BMC Ecology and Evolution.

All other authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, D., Ellepola, G., Herath, J. et al. The two chytrid pathogens of amphibians in Eurasia—climatic niches and future expansion. BMC Ecol Evo 23, 26 (2023). https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-023-02132-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-023-02132-y

Keywords