Every year, a substantial amount of academic research is published about, or with some connection to, the Pokémon franchise and/or fandom. While occasionally some of this research breaks through to the fandom and/or to mainstream news, much of this research sadly goes underreported, and unread by the wider community of Pokémon fans. We here at Bulbagarden think it's about time that changed. Over the next few months, Bulbagarden will be publishing a selection of recent research spanning an A to Z of academic disciplines, from agriculture and biological sciences, to zoology, and plenty of other topics in both the sciences and humanities in between. We don't have a particular schedule for how often these posts will go up, but we hope you'll look forward to them, and that you'll see something in these that catches your interest.
Our first featured article comes from a multinational team of researchers, whose study takes an extensive look at how climate change will impact the existence of suitable habitat for Kangaskhan. Dr. Dan L. Warren, the lead author of the study, noted the primary reason for choosing to use a Pokémon character in their research was to engage a broader audience with the issues conservation scientists regularly need to think about. We think that's a wonderful idea, and it's why we've chosen this article specifically to be our first featured article with this project.
Dr. Warren does offer one point of caution with respect to the culture of Pokémon. “It’s specifically based around over-exploitation with the tagline ‘gotta catch them all.’ The rarer they get, the more valuable they get. This is like some of the larger tuna species, which are in serious danger of going extinct. We might have focused on climate change in this study but, for many Pokémon (and many species around the world), overexploitation should also be a concern.” That's a sentiment that I'd like to think a lot of Pokémon fans can easily come to grips with right now, with the mind boggling levels of demand for Pokémon TCG cards that we're currently seeing worldwide.
The following article was originally published by the authors listed below as an open access article in Methods in Ecology and Evolution, a publication of the Britich Ecological Society (DOI: 10.1111/2041-210X.13591). It is republished here under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.
The effects of climate change on Australia’s only endemic Pokémon: Measuring bias in species distribution models
Authors
Dan L. Warren
Biodiversity and Biocomplexity Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan; Senckenberg Biodiversity and Climate Research Centre, Frankfurt, Germany
Alex Dornburg
Department of Bioinformatics and Genomics, University of North Carolina Charlotte, Charlotte, NC, USA
Katerina Zapfe
Department of Bioinformatics and Genomics, University of North Carolina Charlotte, Charlotte, NC, USA; Department of Biological Sciences, Clemson University, Clemson, SC, USA
Teresa L. Iglesias
Physics and Biology Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan; Animal Resource Section, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan
Abstract
1 | INTRODUCTION
Anthropogenic climate change has been demonstrated to negatively impact the distribution of suitable habitat for many species, leading to local extirpations, range shifts and extinctions (Bellard et al., 2012; Chen et al., 2011; Parmesan, 2006). Managers attempting to mitigate the effects of climate change often have to make decisions that require some information about the tolerance of species for given combinations of temperature, precipitation or other environmental variables in order to decide whether a given area is currently suitable, or will be so given expected environmental change (Hutchinson, 1978; Soberon & Nakamura, 2009). However, experimental approaches to estimating this ‘environmental niche’ are often not feasible for reasons of cost or due to other practical concerns. As such, a number of correlative approaches have been developed that attempt to derive estimates of species’ tolerances by correlating species occurrence data with the geographic distribution of environments. These species distribution models (SDMs, also known as environmental niche models or ENMs) are often one of the key lines of evidence used to make decisions to mitigate the effects of climate change. Unfortunately there are many methodological and conceptual issues involved in the construction and application of SDMs that remain poorly understood.
One major concern in SDM studies is that the methodological choices made in the SDM process may bias the predictions that the resulting models make (Barbet‐Massin et al., 2010; Beaumont et al., 2016; Datta et al., 2020). Using SDMs requires decisions to be made regarding the choice of algorithm, predictor variables, emissions scenario, climate model and study area among other parameters. In turn, these choices may affect the distribution of environments available for modelling as well as the distribution of environments that models will be projected onto, and as such each of these aspects of an SDM study's experimental design has the potential to introduce bias. Clearly, more approaches directed towards quantifying and mitigating biases in SDM studies are needed.
In order to better understand these issues, we use SDMs to examine the effects of climate change on Australia's only endemic Pokémon: kangaskhan (Garura kangaskhani, Figure 1). To estimate the extent to which our results might be driven by methodological biases, we compare our predictions of changing habitat suitability for G. kangaskhani to similar models built from virtual species (Warren et al., 2020), and to a recent study of 220 species of Australian mammals (Beaumont et al., 2016). Garura kangaskhani is the only Pokémon endemic to the geographic region used for these two studies, and is therefore the Pokémon most comparable to those studies when examining the effects of SDM study design. We also perform a post hoc analysis on models from Warren et al. (2020) to illuminate which aspects of SDM study design are responsible for inducing biases when predicting the effects of climate change. Finally, we develop a novel Monte Carlo method that can be used to estimate the bias in any given SDM study design, and compare our predictions for G. kangaskhani to those expected by bias alone.
1.1 | Study system
We provide a detailed discussion of the current (woefully inadequate) state of knowledge about G. kangaskhani in supplement S1, but for the sake of brevity we restrict ourselves here to those aspects of its biology most relevant to conservation. Garura kangaskhani is an omnivorous (Bulbapedia, 2020a), diurnal (PokéBase, 2020) vertebrate endemic to Australia. It is primarily associated with grasslands and savanna habitats, although frequent sightings by field researchers in highly disturbed metropolitan areas indicate that they are also fairly resilient to anthropogenic disturbance; one study of ‘Normal’ type Pokémon (which includes G. kangaskhani) found that parking lots and universities were the second and third most common sources of sightings, with sightings more frequent during partly cloudy weather (The Silph Road, 2016). Periodic spates of sightings across Europe, Asia and the Americas (Pokémon GO Wiki Community, n.d.) further suggests G. kangaskhani are capable of long range dispersal, but do not readily colonize new urban habitats outside of Australia.
1.2 | Threats
The vulnerability of animal and human populations is often a complex function of multiple interacting stressors (O’Brien et al., 2004). Garura kangaskhani were previously hunted to the brink of extinction (Bulbapedia, 2020b), and one of the primary threats to the persistence of G. kangaskhani populations continues to be the poaching of adults and eggs. While kangaskhan may sometimes be eaten (Gilbert, 2020), poaching is primarily motivated by the acquisition of specimens intended for blood sport; captured animals are used in contests in which the maternal defence of the young is exploited to motivate adult G. kangaskhani into fighting both conspecifics and heterospecifics (Molloy, 2013). Individuals deemed to be unfit for combat are frequently ground into ‘candy’ which is fed to combatants (D’Anastasio, 2016). This destructive practice is not limited to G. kangaskhani, nor is it geographically restricted to Australia. Current estimates indicate that tens of millions of ‘trainers’ and ‘breeders’ are participating in these matches worldwide (Wikipedia, 2020), with hundreds of millions of dollars being spent on supplies alone (Phillips, 2018). With the recent discovery of a rare ‘shiny’ morphotype of G. kangaskhani that can add a visual enhancement to a trainer's ‘battle league party’, we anticipate that collection pressure on this species will only increase in the future. We see little room for optimism in this scenario given that the culture around Pokémon harvesting so explicitly favours overexploitation, as seen in the common refrain of ‘gotta catch ‘em all’ (Prof. Y. Ohkido, pers. comm.).
Compounding the threats already posed by exploitation, anthropogenic climate change may constitute an additional challenge for G. kangaskhani. Little is known about the species’ climate niche beyond a marked preference for partly cloudy weather (PokemonGoHub, 2020; Pokémon GO Wiki Community, n.d.). It is therefore difficult to estimate a priori the response of G. kangaskhani to environmental change, or how this additional stressor might combine with the effects of poaching to impact the species’ long‐term survival. Moreover, their apparently broad habitat tolerances and the fact that most weather can be considered at least ‘partly’ cloudy, suggest that G. kangaskhani has few strong climate associations within Australia. This lack of strong associations may challenge our ability to make precise forecasts for this species, but may serve to make this species more useful for understanding the effects of bias; if a method shows a tendency to make strong predictions when there is no actual relationship between environmental predictors and the occurrence of the species, it suggests that those predictions may be driven by methodological biases.
In this study, we use a citizen science dataset of occurrences to estimate the effects of climate change on populations of G. kangaskhani. These data were initially recorded by hobbyists and professional trainers seeking out G. kangaskhani for exploitation, and as such the use of these data for conservation modelling does pose some ethical issues. However, given the uncertain conservation status of these noble creatures and the paucity of georeferenced localities of pokémon in natural history collection databases, few other options exist. We also acknowledge that the small sample size used here (37 points) may be cause for some concern, but note that this is not an unusually small sample size for SDM studies (Galante et al., 2018; Loe et al., 2012). Further, we compare our results to other SDM studies with N = 100 (Warren et al., 2020) and N ranging from 10 to 7,137 (Beaumont et al., 2016), and find broadly concordant patterns (see Section 4). We also note that seminal contributions have been made by investigators using data sources of similar quality (Hurlbert, 1990; Lozier et al., 2009; Michelot et al., 2016), and as such we feel that the use of these data are justified.
2 | MATERIALS AND METHODS
We obtained occurrence data for G. kangaskhani from online databases maintained by individuals seeking to collect them for use in blood sport or for personal gains that exploit the lack of international trade regulations for this species (PKMNGoTrading, 2016; KOPlayer, 2018; The Road and Community, 2016). Given how little is known about G. kangaskhani, many decisions in the modelling process were necessarily made based on common practice in the literature rather than a priori knowledge of the underlying biology (Araújo et al., 2019). In order to reduce the impact of spatial sampling bias we removed duplicate occurrences, as is common in SDM studies (Peterson et al., 2011). We also removed any occurrences falling in regions with no environmental data, resulting in 37 unique data points. Climate data were obtained from CliMond (Kriticos et al., 2012) at 10‐km resolution. We used the 19 Bioclim layers for both present and future projections. Future climate layers were obtained for four combinations of emissions scenario and climate model; A1B‐CS, A2‐CS, A1B‐MR and A2‐MR. Future climate projections were used for 2030, 2050, 2070, 2080, 2090 and 2100.
Using the ENMTools r package (Warren, Matzke, et al., 2021), we constructed SDMs using six different algorithms; generalized additive models (GAM), generalized linear models (GLM), maximum entropy (Phillips et al., 2006), random forests, Bioclim (Nix & Busby, 1986) and Domain (Carpenter et al., 1993). As this study is based on a citizen science dataset and limited research funds are available for systematic sampling of Pokemon distributions, true absence data are not available for G. kangaskhani. We therefore adopt a presence/background modelling approach that is common in SDM studies, in which environmental conditions at occurrence points are contrasted with the distribution of available environments within the species range (Peterson et al., 2011). This approach is common in SDM studies without true absence data (Araújo et al., 2019). Background points were extracted from a 250‐km circular buffer around each occurrence. Models were built using 10,000 background points. For each model, 30% of the presence data were randomly withheld for model validation. Code for model construction is given in Appendix S3 (Warren, Dornburg, et al., 2021).
In order to summarize the effects of climate change on the distribution of suitable habitat, we projected models onto each future climate scenario and then calculated the difference between the current and future suitability at each grid cell in the study area. We then average those projected changes to get the overall mean change in habitat suitability projected by a given combination of distribution model, emissions scenario and climate model. Change in habitat suitability was evaluated at two spatial scales; within the species’ current range, and at a continental scale including all of Australia and Tasmania.
To evaluate the extent to which predictions concerning the future of G. kangaskhani may have been driven by methodological bias, we compare our empirical predictions to predictions derived from two other sources. First, we examine the predicted effects of climate change on a set of models for 220 simulated species in Australia from Warren et al. (2020). Second, we compare our observations for G. kangaskhani to a distribution of predictions from a novel Monte Carlo approach that randomly chooses occurrence points from the study area.
2.1 | Post hoc analysis of simulations
In the simulations from Warren et al. (2020), virtual species’ ranges were derived from the availability of suitable habitat based on an underlying simulated niche, with a variety of niche breadths and range sizes. Within these ranges, species occurrences were sampled in proportion to the product of the suitability of habitat and a model of spatial sampling bias of varying strength (see Warren et al., 2020 for details). Models from these datasets were constructed using methods similar to those used to analyse G. kangaskhani, and were projected onto the same set of future scenarios. In order to understand the sources of bias and uncertainty in these predictions we also conducted a set of post hoc GLM analyses in which the average suitability change and associated errors from these models were regressed against year, with interaction terms for SDM algorithm, emissions scenario and climate model. Significant interactions from these regressions indicate that a given component of the study design modifies the direction or magnitude of the bias in predictions of future habitat suitability. See Appendix S2 (Warren, Dornburg, et al., 2021) for details and R code for these analyses.
2.2 | Monte Carlo estimates of methodological bias
To estimate biases in our predictions for kangaskhan, we developed a novel Monte Carlo approach that is derived from a test developed by Raes and ter Steege (2007). As originally designed, this approach is used to generate null distributions of predictive performance for species distribution models. The null distribution is constructed by selecting localities at random from within the study area until a dataset of the same sample size as the empirical dataset is obtained. These random localities are then used to generate SDMs using the same methods as those used for the empirical model, and the predictive performance of the model is measured. By iteratively repeating this process, the user can construct an approximate distribution of expected predictive performance for an SDM under a given set of modelling conditions under the null hypothesis that the occurrence of the species within the study area is uncorrelated with any of the predictor variables. Comparing the discrimination accuracy of the empirical model to this null distribution allows users to test the hypothesis that their model has predictive power beyond what can be explained by their study design alone.
Our Monte Carlo approach (Figure 2) modifies this method to examine the biases implicit in transferring species distribution models under a given set of analytical conditions. As in Raes and ter Steege (2007) and Bohl et al. (2019), we select N points at random from the study area (i.e. the area from which background data were taken), where N is the number of occurrence points used for modelling in the empirical dataset. Using these random points as ‘occurrences’ we build SDMs with all other modelling choices being identical to those for the empirical dataset. We then project our replicate models onto the different combinations of climate model and emissions scenario, and calculate the predicted changes in habitat suitability for those models in the same way as with the empirical data. This allows us to estimate a distribution of expected changes in habitat suitability given the same study area, SDM algorithm, emissions scenario and climate model as our empirical data, under the hypothesis that the actual occurrence of the species within the study area is unrelated to the environmental variables used for model construction.
2.3 | Interpretation
Our post hoc analyses of simulations from Warren et al. (2020) is intended to quantify the expected behaviour of our modelling approaches given a dataset in which there is an underlying (simulated) biological mechanism relating its probability of occurrence to a set of environmental predictors, while the novel Monte Carlo approach is intended to estimate the behaviour of a model built under identical conditions to our empirical model (including study region), but for which there is no biological mechanism underlying the distribution of occurrence points within that study area. The distributions of predictions produced by these two approaches are not necessarily intended to represent a null hypothesis that must be rejected in order to trust model results, but rather as estimates of the behaviour of our modelling approaches independent of our data. By comparing our empirical results to these distributions, we can determine to what extent our predictions are driven by our data, over and above the biases induced by our analytical framework.
3 | RESULTS
We find that the predicted change in suitability of habitat for G. kangaskhani is relatively consistent across spatial scales (continental vs. within species’ current range), but depends strongly on the method of analysis. Models built using Bioclim, Domain or GAM algorithms predict that habitat suitability will decline on average, but disagree substantially on the magnitude of that decline (Figure 3). GLM and Maxent models both make mixed predictions depending on the climate model used; both methods predict increased habitat suitability under the MIROC‐H climate model, but decreased suitability under the CSIRO model (Appendix S2, Warren, Dornburg, et al., 2021). Random forest models, by contrast, predict increasing habitat suitability regardless of emissions scenario or climate model.
The models from Warren et al., (2020) were constructed from occurrence data based on simulated organisms with randomly selected niche parameters. However, the predictions made from these models demonstrate many of the same behaviours as those seen in the empirical models for G. kangaskhani (Figure 4, see appendix S2 in Warren, Dornburg, et al., 2021 for details). In particular we notice a broad trend for Bioclim and Domain models to predict declining habitat suitability, and for random forests to predict increasing habitat suitability over time at the continental scale. In contrast to the results from G. kangaskhani and the Monte Carlo replicates, random forests models from simulated species predict lower‐than‐present suitability within species’ current ranges in 2030. However, they show a similar tendency towards increasing estimates of habitat suitability over subsequent years. GLM, GAM and Maxent models show much more mixed predictions across the set of simulated species, which in the case of GLM and Maxent mirror our predictions for G. kangaskhani.
As the Warren et al., (2020) data are based on simulated organisms for which the true suitability of habitat is known, we can project the true effects of climate change on those species and measure the extent to which our model predictions are inaccurate (Appendix S2, Warren, Dornburg, et al., 2021). At a continental scale, we find that the average true suitability of habitat changes little with time, but that the variance in suitability across simulated species increases substantially further in the future. Within simulated species’ current ranges, the modal true change in habitat suitability over time remains very close to zero, but the spread of the distribution of changes in true habitat suitability is increasingly biased downwards as predictions get further from the present.
Examining the distributions of errors in suitability predictions made by models shows that for the most part the general trends in behaviour of each modelling approach discussed above are in fact largely explained by biases in the kinds of errors each approach is making. GLM analyses of predictions of habitat suitability for these simulations and their accompanying errors demonstrate highly significant effects of every component of the modelling process (modelling algorithm, year, climate model, emissions scenario). This suggests that each one of these choices will affect the expected outcome of any modelling endeavour to a large (and typically un‐measured) degree.
The behaviour of models built from randomly chosen data points within the training area for G. kangaskhani showed patterns that were generally concordant with those from the simulated species, and in many cases the predictions for G. kangaskhani fall well within the distribution of changes from those Monte Carlo replicates (Figure 5). In most cases where the behaviour of the G. kangaskhani model did fall outside of the distribution estimated from random points, it was in the direction of making less extreme predictions for the effects of climate change on G. kangaskhani than expected.
4 | DISCUSSION
In this study, we used a citizen science dataset in conjunction with climate data to estimate the effects of climate change on Australia's only endemic Pokémon, and find that our conclusions about changes in habitat suitability are profoundly affected by methodological choices. By examining multiple simulated data sources, we highlight the important possibility that the results seen here might not be indicative of any direct effect of climate change on the suitability of habitat for G. kangaskhani, but instead represent the effects of bias that is imposed by the methods used to construct these models. The predicted effects of climate change on the suitability of habitat for G. kangaskhani generally reflect predictions we derive from models for simulated species as well as predictions from Monte Carlo tests based on randomly chosen data points within the range of G. kangaskhani.
It is interesting also to compare our results with Beaumont et al. (2016), a study that used models built with real data for 220 Australian mammal species to examine which modelling approaches were more prone to making predictions of range contractions or expansions (Figure 6). Those models were based on real data for species for which the climate niche is not known, and as a result it was not possible from that study to say which of those models better estimates the true impacts of climate change. However, we note that many of the model behaviours seen in Beaumont et al. (2016) echo those seen in our simulated data and Monte Carlo replicates; climate envelope methods show a strong tendency towards predictions of declining suitability and random forests show the strongest tendency to predict increasing suitability, while Maxent, GLM and GAM all produce predictions of change in habitat suitability that are relatively symmetric around zero (Figures 4–6-4–6). The similarity between the results for real species, simulated species and Monte Carlo replicates strongly suggests that the effects of the biases inherent in the design of SDM analyses are a ubiquitous feature of studies of this type.
In this study, our Monte Carlo tests focus on how a given set of modelling conditions may bias our estimates of the effects of climate change on average habitat suitability. However, we note that this framework can easily be extended. There are many different metrics that may be used to quantify the effects of climate change on species (e.g. change in habitable area based on thresholded suitability scores, geographic shifts in ranges or centroids, changes in climate niche breadth, etc.). Applying this approach to those metrics is simply a matter of applying the chosen metric to predictions from each Monte Carlo replicate to obtain the distribution of that metric given randomly sampled occurrence data. Similarly, this approach could easily be extended to obtain distributions for models being transferred in space, rather than in time; transferring models in time and space present very similar challenges for SDM studies, and it is likely that similar biases exist in any study of that type. Further, it is possible to calculate summary statistics on a per‐grid‐cell basis from the suitability maps for Monte Carlo replicates, producing a map of the geographic regions where model predictions might be driven by study design more than by the data used in modelling (Figure 7). Finally, we note that in the current study occurrence points for Monte Carlo replicates were sampled with uniform probability across the study with no spatial autocorrelation or sampling bias. While this approach has been widely used for constructing null distributions of model performance (Bohl et al., 2019; Raes & ter Steege, 2007; Warren, Dornburg, et al., 2021), other methods that attempt to maintain similar spatial autocorrelation to empirical data when sampling may yield different null distributions (Beale et al., 2008; Nunes & Pearson, 2017; Roxburgh & Chesson, 1998). We consider this a very promising area for future study, as it may help to better understand which aspects of data collection and curation are most instrumental in reducing bias in predictions.
Species distribution models are widely used with the goal of guiding management and policy. However, the fact that we produce consistent predictions from meaningless data indicates that there are some combinations of emissions scenario, climate model, study area and species distribution modelling method that are strongly biased towards a given prediction. Crucially, in some modelling conditions we find that all Monte Carlo replicates agreed on the direction of change in habitat suitability in the future, despite the fact that the points for these replicates were chosen at random from within Kangaskahan's range and therefore represented no underlying relationship between the environment and the probability of occurrence. At the very least this suggests that those study designs should be chosen only when there is a compelling reason to believe that they accurately represent the underlying biology. In the case of G. kangaskhani we have insufficient information to make such a decision; in addition to our scant knowledge of the species’ ecology we anticipate that there are likely to be serious issues with data quality due to spatial sampling bias and small sample size due to an unfortunate lack of scientific attention to wild Pokémon populations. In particular, we would like to point out that there may be some room for optimism given our level of ignorance of the distribution of G. kangaskhani; we do not currently know how far into the interior of Australia these animals might be found and so are unlikely to have fully characterized its environmental tolerances. Currently very little conservation funding is allocated for land acquisition and captive breeding programs for Pokémon, but increased investment in this area could be broadly beneficial; charismatic megavertebrates can often be used as flagship species to increase community involvement in conserving whole ecosystems (Cummins, 2020). In order to effectively manage this majestic species, we suggest that future funding efforts be directed towards sending the authors of this study into remote areas of Australia to collect more reliable data, preferably equipped with a large supply of the ‘incense’ (presumably some sort of pheromone) that Pokémon hunters use to attract animals for capture.
It is reasonable to view with heavy scepticism the results of any study in which the analytical design is only capable of producing one answer regardless of the data. This may not mean, however, that the predictions from such a study are necessarily wrong. To illustrate this point, it is worth considering an extreme case in which the entire study area becomes unsuitable for any organism currently living within it (e.g. in future scenarios all currently occupied habitat is constantly on fire). In this case, any reasonable model built from any set of random occurrence data in that region would show complete extirpation of the organism, which would also coincide with the results seen from any carefully vetted empirical dataset for a real species. The conclusions of those models would be correct, but they would be wholly unaffected by the data; they would be an inevitable result of the modelling conditions. Any dataset would likely support the same conclusions regardless of whether those data represent any underlying biological phenomenon.
In the real world, however, the effects of climate change are rarely so absolute, and are rarely known in advance. Finding that we are strongly biased towards making specific predictions irrespective of our data is therefore a cause for great concern; it implies that we may have essentially selected our conclusions by selecting our study design. In such a situation, the resulting outputs should at minimum be interpreted in the context of the biological plausibility of those methodological choices, the biases they induce, and their consequences for decisions made based on those models. In particular, if we find ourselves in a situation where different approaches show strong opposing biases, interpretation of those results must be made in the context of how biologically realistic the alternative modelling approaches are. Arguably one should go further than that, however, and simply not trust the outcome of a modelling approach that produces the same result regardless of the data provided to it. In these cases, additional lines of evidence are needed to adequately inform conservation decisions. Understanding when we are in such a situation is the key to making these decisions, and we believe that Monte Carlo tests such as the one presented here may be a valuable tool for measuring the biases inherent in the design of empirical studies.
ACKNOWLEDGEMENTS
Gotta thank ‘em all! The authors would like to thank the Westpunt Syndicate and acknowledge funding sources including an NSF GRFP grant to K.Z. (#0000‐0002‐5159‐641X). We are deeply grateful for many helpful comments from Erica Newman, Linda Beaumont, Megan Higgie, Jamie Kass, Nick Matzke, Peter Linder, Liam Langan, Jorge Lobo and an anonymous reviewer. We are also very grateful to Bob O'Hara and the rest of the editorial staff at MEE for indulging this foolishness. We would also like to thank all of the Pokémon trainers out there who made this possible. Finally the authors would like to acknowledge that this started out as a joke but actually got pretty dang interesting. The use of kangaskhan, Pokémon and related intellectual property here is intended as non‐commercial parody fair use, and does not in any way reflect an endorsement by Niantic, The Pokémon Company or Prof. Y. Ohkido. Open Access funding enabled and organized by Projekt DEAL. WOA Institution: SENCKENBERG GESELLSCHAFT FUR NATURFORSCHUNG; Blended DEAL: Projekt DEAL.
AUTHORS' CONTRIBUTIONS
Initial conception by D.L.W.; Study design by D.L.W., T.L.I., and A.D., analyses by D.W. The manuscript was written with equal contributions from all authors. Original artwork by K.Z.
REFERENCES
Appendix S1: Additional details about the life history of Garura kangaskhani and comments on how these may contribute to their response to climate change
Biology and life history
Despite being a charismatic megavertebrate (~220 cm adult height) that is frequently sighted near human population centers, G. kangaskhani was not formally described until 1997 (Bulbapedia, 2019). Information regarding its basic biology is scarce and much of what is available is based on anecdotal accounts from citizen scientists and unlicensed breeders, many of whom have no formal scientific training. However, what information is available suggests that the reproductive system and demography of G. kangaskhani may be unusual in a number of ways, and these factors may need to be taken into account when making conservation decisions for the species. Therefore, in addition to addressing the effects of climate change on G. kangaskhani we offer the following information as a summary of the current consensus on the biology of the species, with the hopes that it will inspire further investigation and guide conservation efforts.
Garura kangaskhani has never been sequenced for DNA, and therefore all available information about its genome is based on speculation (Butterfree, 2016). As a consequence of the lack of genetic information, the evolutionary history of G. kangaskhani remains contentious. This lineage is hypothesized to be either sister to a lineage leading to a clade of “fire type” Pokémon that includes Entei, Heatran, and Houndoom (Shelomi et al., 2012) or a member of a distinct lineage of ‘marsupial’ Pokémon that includes Cubone (Spaniard, 2016). While its breeding habits suggest that it is a marsupial, it possesses a series of dermal plates that are reminiscent of stegosaurs. This unusual feature may represent convergent evolution of a useful defense strategy or an ancient hybridization event between a stegosaur and a stem-group metatherian (Margulis & Sagan, 2008; Williamson, 2009). As G. kangaskhani has never been seen to evolve (PokéBase, 2020) it is further possible that this taxon represents a “living fossil”. However, disentangling these and other hypotheses of evolutionary history is stymied by this Pokémon’s enigmatic reproductive biology.
The reproductive biology of G. kangaskhani has received more attention than most other aspects of its biology, but what is known is perplexing and often contradictory. All sources agree that offspring are strongly altricial (Bulbapedia, 2020), living within the mother’s pouch for approximately three years. Mothers are said to form strong pair bonds with their young, and sleep standing upright in order to protect them (Figure 1). Interestingly, all observed individuals of G. kangaskhani have been reported as female. We suspect that this is largely due to the presence of young in the pouch, as detailed anatomical examinations have not been performed to our knowledge. We briefly explore two alternative hypotheses to explain this pattern: that male G. kangaskhani do exist but have not been correctly documented, and that G. kangaskhani are parthenogenic.
The lack of male G. kangaskhani observations may stem from males being rare, necessitating more extensive field environmental-genomic surveys to verify (Wilson Sayres, 2018). However, it is also possible that many male G. kangaskhani have simply been misidentified. For example, males may be identified as female based on the assumption of uniparental care. While in most marsupial species pouches are only present in females, there are known exceptions including the water opossum (or “yapok”; Chironectes minimus) and the recently extinct thylacine (Thylacinus cynocephalus), in which males have a pouch that is used to protect the genitals from snagging on vegetation (Nogueira et al., 2004; Paddle, 2002). It is therefore possible that male G. kangaskhani are using a highly modified intromittent organ as a nuptial display for females and that the rare ‘shiny’ forms are in fact males with display coloration. Alternatively, G. kangaskhani males may simply appear so different from females as to have been misidentified as another species. This is a mistake that has been made many times in other organisms (Moraes et al., 2017; Perez-Mellado & de la Riva, 1993). It is even possible that the male G. kangaskhani is in fact many times smaller than the female and either attached to, or embedded within, the female of the species (Pietsch, 2005) (Figure S1.1).
The other major possibility, and currently the favored hypothesis (Dyer, 2016), is that G. kangaskhani reproduce via parthenogenesis. While parthenogenesis is primarily associated with plant and invertebrate species, it has been observed in a number of vertebrate species (Booth et al., 2012; Fields et al., 2015; Fujita et al., 2007; Moritz & Heideman, 1993; Walker, 1986). Although the idea may seem somewhat farfetch’d, it is supported by the repeated observation that young G. kangaskhani are born with their own offspring in pouch. This observation would be hard to reconcile with the above hypotheses regarding cryptic males, but is quite compatible with a parthenogenic mating system, as a similar phenomenon is also seen in parthenogenic aphids (Stern, 2008).
The effects of multiple stressors
The results of this study indicate that climate change may have significant negative impacts on G. kangaskhani, as indicated by the decline in overall average habitat suitability. This is likely to compound the effects of existing threats to populations of G. kangaskhani, which are already under serious pressure from human exploitation. In particular, the presence of high density populations in areas of anthropogenic disturbance is cause for concern. While these areas of exceptional local abundance are often referred to as “nesting” habitat, there is no obvious reason to believe that they nest colonially under those conditions, and the current consensus is that the species may not nest communally at all (The Silph Road Community, 2016). It is not clear why these comparatively high densities are seen in human-modified environments, but we speculate that it is driven either by local food resources or by human activity uncovering deposits of kangaskhanite, a mineral they need for development. Given that these high density locations are also favored by hunters looking to capture G. kangaskhani for sport, these areas may in fact constitute an “ecological trap” (Hale & Swearer, 2016); extreme sink populations driven by a combination of attractive resources and high levels of unregulated exploitation. If this is the case, the effects of climate change on surrounding natural source populations may be compounded by emigration from those populations, leading to a potential for sudden and irreversible collapse. We see little room for optimism in this scenario given that the culture around Pokémon harvesting so explicitly favors exploitation, as seen in the common refrain of “gotta catch ‘em all” (Prof. Oak, pers. comm.).
Climate change places significant selective pressures on natural populations, and it is thought that some species may be able to adapt to changing climate regimes via natural selection (Catullo et al., 2019). The fact that G. kangaskhani are likely parthenogenic makes this significantly less likely; absence of sexual reproduction reduces effective population size and precludes recombination, leading to a curtailed ability to respond to selection (Rice & Friberg, 2009). However, we would also like to highlight again the similarity between what is known about G. kangaskhani reproduction and the reproductive system of the aphid (Stern, 2008), and in particular we note that aphids are capable of facultative sexual reproduction when the environment becomes unpredictable, increasing genetic variability of offspring in order to facilitate responses to selective pressures. If the similarities between G. kangaskhani and aphids extend to this response, it is possible that climate change will induce a wave of previously unseen sexual Kangaskhan morphs (Figure S1.2). The presence of large numbers of G. kangaskhani in human population centers makes this a somewhat alarming prospect, even more so if the sexual morphs are territorial.
Our first featured article comes from a multinational team of researchers, whose study takes an extensive look at how climate change will impact the existence of suitable habitat for Kangaskhan. Dr. Dan L. Warren, the lead author of the study, noted the primary reason for choosing to use a Pokémon character in their research was to engage a broader audience with the issues conservation scientists regularly need to think about. We think that's a wonderful idea, and it's why we've chosen this article specifically to be our first featured article with this project.
Dr. Warren does offer one point of caution with respect to the culture of Pokémon. “It’s specifically based around over-exploitation with the tagline ‘gotta catch them all.’ The rarer they get, the more valuable they get. This is like some of the larger tuna species, which are in serious danger of going extinct. We might have focused on climate change in this study but, for many Pokémon (and many species around the world), overexploitation should also be a concern.” That's a sentiment that I'd like to think a lot of Pokémon fans can easily come to grips with right now, with the mind boggling levels of demand for Pokémon TCG cards that we're currently seeing worldwide.
The following article was originally published by the authors listed below as an open access article in Methods in Ecology and Evolution, a publication of the Britich Ecological Society (DOI: 10.1111/2041-210X.13591). It is republished here under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.
The effects of climate change on Australia’s only endemic Pokémon: Measuring bias in species distribution models
FIGURE 1
Artist's rendition of adult Gangura kangaskhani with young in pouch
Artist's rendition of adult Gangura kangaskhani with young in pouch
Authors
Dan L. Warren
Biodiversity and Biocomplexity Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan; Senckenberg Biodiversity and Climate Research Centre, Frankfurt, Germany
Alex Dornburg
Department of Bioinformatics and Genomics, University of North Carolina Charlotte, Charlotte, NC, USA
Katerina Zapfe
Department of Bioinformatics and Genomics, University of North Carolina Charlotte, Charlotte, NC, USA; Department of Biological Sciences, Clemson University, Clemson, SC, USA
Teresa L. Iglesias
Physics and Biology Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan; Animal Resource Section, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan
Abstract
- Species distribution models (SDMs) are frequently used to predict the effects of climate change on species of conservation concern. Biases inherent in the process of constructing SDMs and transferring them to new climate scenarios may result in undesirable conservation outcomes. We explore these issues and demonstrate new methods to estimate biases induced by the design of SDM studies. We present these methods in the context of estimating the effects of climate change on Australia's only endemic Pokémon.
- Using a citizen science dataset, we build species distribution models for Garura kangaskhani to predict the effects of climate change on the suitability of habitat for the species. We demonstrate a novel Monte Carlo procedure for estimating the biases implicit in a given study design, and compare the results seen for Pokémon to those seen from our Monte Carlo tests as well as previous studies in the same region using both simulated and real data.
- Our models suggest that climate change will impact the suitability of habitat for G. kangaskhani, which may compound the effects of threats such as habitat loss and their use in blood sport. However, we also find that using SDMs to estimate the effects of climate change can be accompanied by biases so strong that the data themselves have minimal impact on modelling outcomes.
- We show that the direction and magnitude of bias in estimates of climate change impacts are affected by every aspect of the modelling process, and suggest that bias estimates should be included in future studies of this type. Given the widespread use of SDMs, systemic biases could have substantial financial and opportunity costs. By demonstrating these biases and presenting a novel statistical tool to estimate them, we hope to provide a more secure future for G. kangaskhani and the rest of the world's biodiversity.
1 | INTRODUCTION
Anthropogenic climate change has been demonstrated to negatively impact the distribution of suitable habitat for many species, leading to local extirpations, range shifts and extinctions (Bellard et al., 2012; Chen et al., 2011; Parmesan, 2006). Managers attempting to mitigate the effects of climate change often have to make decisions that require some information about the tolerance of species for given combinations of temperature, precipitation or other environmental variables in order to decide whether a given area is currently suitable, or will be so given expected environmental change (Hutchinson, 1978; Soberon & Nakamura, 2009). However, experimental approaches to estimating this ‘environmental niche’ are often not feasible for reasons of cost or due to other practical concerns. As such, a number of correlative approaches have been developed that attempt to derive estimates of species’ tolerances by correlating species occurrence data with the geographic distribution of environments. These species distribution models (SDMs, also known as environmental niche models or ENMs) are often one of the key lines of evidence used to make decisions to mitigate the effects of climate change. Unfortunately there are many methodological and conceptual issues involved in the construction and application of SDMs that remain poorly understood.
One major concern in SDM studies is that the methodological choices made in the SDM process may bias the predictions that the resulting models make (Barbet‐Massin et al., 2010; Beaumont et al., 2016; Datta et al., 2020). Using SDMs requires decisions to be made regarding the choice of algorithm, predictor variables, emissions scenario, climate model and study area among other parameters. In turn, these choices may affect the distribution of environments available for modelling as well as the distribution of environments that models will be projected onto, and as such each of these aspects of an SDM study's experimental design has the potential to introduce bias. Clearly, more approaches directed towards quantifying and mitigating biases in SDM studies are needed.
In order to better understand these issues, we use SDMs to examine the effects of climate change on Australia's only endemic Pokémon: kangaskhan (Garura kangaskhani, Figure 1). To estimate the extent to which our results might be driven by methodological biases, we compare our predictions of changing habitat suitability for G. kangaskhani to similar models built from virtual species (Warren et al., 2020), and to a recent study of 220 species of Australian mammals (Beaumont et al., 2016). Garura kangaskhani is the only Pokémon endemic to the geographic region used for these two studies, and is therefore the Pokémon most comparable to those studies when examining the effects of SDM study design. We also perform a post hoc analysis on models from Warren et al. (2020) to illuminate which aspects of SDM study design are responsible for inducing biases when predicting the effects of climate change. Finally, we develop a novel Monte Carlo method that can be used to estimate the bias in any given SDM study design, and compare our predictions for G. kangaskhani to those expected by bias alone.
1.1 | Study system
We provide a detailed discussion of the current (woefully inadequate) state of knowledge about G. kangaskhani in supplement S1, but for the sake of brevity we restrict ourselves here to those aspects of its biology most relevant to conservation. Garura kangaskhani is an omnivorous (Bulbapedia, 2020a), diurnal (PokéBase, 2020) vertebrate endemic to Australia. It is primarily associated with grasslands and savanna habitats, although frequent sightings by field researchers in highly disturbed metropolitan areas indicate that they are also fairly resilient to anthropogenic disturbance; one study of ‘Normal’ type Pokémon (which includes G. kangaskhani) found that parking lots and universities were the second and third most common sources of sightings, with sightings more frequent during partly cloudy weather (The Silph Road, 2016). Periodic spates of sightings across Europe, Asia and the Americas (Pokémon GO Wiki Community, n.d.) further suggests G. kangaskhani are capable of long range dispersal, but do not readily colonize new urban habitats outside of Australia.
1.2 | Threats
The vulnerability of animal and human populations is often a complex function of multiple interacting stressors (O’Brien et al., 2004). Garura kangaskhani were previously hunted to the brink of extinction (Bulbapedia, 2020b), and one of the primary threats to the persistence of G. kangaskhani populations continues to be the poaching of adults and eggs. While kangaskhan may sometimes be eaten (Gilbert, 2020), poaching is primarily motivated by the acquisition of specimens intended for blood sport; captured animals are used in contests in which the maternal defence of the young is exploited to motivate adult G. kangaskhani into fighting both conspecifics and heterospecifics (Molloy, 2013). Individuals deemed to be unfit for combat are frequently ground into ‘candy’ which is fed to combatants (D’Anastasio, 2016). This destructive practice is not limited to G. kangaskhani, nor is it geographically restricted to Australia. Current estimates indicate that tens of millions of ‘trainers’ and ‘breeders’ are participating in these matches worldwide (Wikipedia, 2020), with hundreds of millions of dollars being spent on supplies alone (Phillips, 2018). With the recent discovery of a rare ‘shiny’ morphotype of G. kangaskhani that can add a visual enhancement to a trainer's ‘battle league party’, we anticipate that collection pressure on this species will only increase in the future. We see little room for optimism in this scenario given that the culture around Pokémon harvesting so explicitly favours overexploitation, as seen in the common refrain of ‘gotta catch ‘em all’ (Prof. Y. Ohkido, pers. comm.).
Compounding the threats already posed by exploitation, anthropogenic climate change may constitute an additional challenge for G. kangaskhani. Little is known about the species’ climate niche beyond a marked preference for partly cloudy weather (PokemonGoHub, 2020; Pokémon GO Wiki Community, n.d.). It is therefore difficult to estimate a priori the response of G. kangaskhani to environmental change, or how this additional stressor might combine with the effects of poaching to impact the species’ long‐term survival. Moreover, their apparently broad habitat tolerances and the fact that most weather can be considered at least ‘partly’ cloudy, suggest that G. kangaskhani has few strong climate associations within Australia. This lack of strong associations may challenge our ability to make precise forecasts for this species, but may serve to make this species more useful for understanding the effects of bias; if a method shows a tendency to make strong predictions when there is no actual relationship between environmental predictors and the occurrence of the species, it suggests that those predictions may be driven by methodological biases.
In this study, we use a citizen science dataset of occurrences to estimate the effects of climate change on populations of G. kangaskhani. These data were initially recorded by hobbyists and professional trainers seeking out G. kangaskhani for exploitation, and as such the use of these data for conservation modelling does pose some ethical issues. However, given the uncertain conservation status of these noble creatures and the paucity of georeferenced localities of pokémon in natural history collection databases, few other options exist. We also acknowledge that the small sample size used here (37 points) may be cause for some concern, but note that this is not an unusually small sample size for SDM studies (Galante et al., 2018; Loe et al., 2012). Further, we compare our results to other SDM studies with N = 100 (Warren et al., 2020) and N ranging from 10 to 7,137 (Beaumont et al., 2016), and find broadly concordant patterns (see Section 4). We also note that seminal contributions have been made by investigators using data sources of similar quality (Hurlbert, 1990; Lozier et al., 2009; Michelot et al., 2016), and as such we feel that the use of these data are justified.
2 | MATERIALS AND METHODS
We obtained occurrence data for G. kangaskhani from online databases maintained by individuals seeking to collect them for use in blood sport or for personal gains that exploit the lack of international trade regulations for this species (PKMNGoTrading, 2016; KOPlayer, 2018; The Road and Community, 2016). Given how little is known about G. kangaskhani, many decisions in the modelling process were necessarily made based on common practice in the literature rather than a priori knowledge of the underlying biology (Araújo et al., 2019). In order to reduce the impact of spatial sampling bias we removed duplicate occurrences, as is common in SDM studies (Peterson et al., 2011). We also removed any occurrences falling in regions with no environmental data, resulting in 37 unique data points. Climate data were obtained from CliMond (Kriticos et al., 2012) at 10‐km resolution. We used the 19 Bioclim layers for both present and future projections. Future climate layers were obtained for four combinations of emissions scenario and climate model; A1B‐CS, A2‐CS, A1B‐MR and A2‐MR. Future climate projections were used for 2030, 2050, 2070, 2080, 2090 and 2100.
Using the ENMTools r package (Warren, Matzke, et al., 2021), we constructed SDMs using six different algorithms; generalized additive models (GAM), generalized linear models (GLM), maximum entropy (Phillips et al., 2006), random forests, Bioclim (Nix & Busby, 1986) and Domain (Carpenter et al., 1993). As this study is based on a citizen science dataset and limited research funds are available for systematic sampling of Pokemon distributions, true absence data are not available for G. kangaskhani. We therefore adopt a presence/background modelling approach that is common in SDM studies, in which environmental conditions at occurrence points are contrasted with the distribution of available environments within the species range (Peterson et al., 2011). This approach is common in SDM studies without true absence data (Araújo et al., 2019). Background points were extracted from a 250‐km circular buffer around each occurrence. Models were built using 10,000 background points. For each model, 30% of the presence data were randomly withheld for model validation. Code for model construction is given in Appendix S3 (Warren, Dornburg, et al., 2021).
In order to summarize the effects of climate change on the distribution of suitable habitat, we projected models onto each future climate scenario and then calculated the difference between the current and future suitability at each grid cell in the study area. We then average those projected changes to get the overall mean change in habitat suitability projected by a given combination of distribution model, emissions scenario and climate model. Change in habitat suitability was evaluated at two spatial scales; within the species’ current range, and at a continental scale including all of Australia and Tasmania.
To evaluate the extent to which predictions concerning the future of G. kangaskhani may have been driven by methodological bias, we compare our empirical predictions to predictions derived from two other sources. First, we examine the predicted effects of climate change on a set of models for 220 simulated species in Australia from Warren et al. (2020). Second, we compare our observations for G. kangaskhani to a distribution of predictions from a novel Monte Carlo approach that randomly chooses occurrence points from the study area.
2.1 | Post hoc analysis of simulations
In the simulations from Warren et al. (2020), virtual species’ ranges were derived from the availability of suitable habitat based on an underlying simulated niche, with a variety of niche breadths and range sizes. Within these ranges, species occurrences were sampled in proportion to the product of the suitability of habitat and a model of spatial sampling bias of varying strength (see Warren et al., 2020 for details). Models from these datasets were constructed using methods similar to those used to analyse G. kangaskhani, and were projected onto the same set of future scenarios. In order to understand the sources of bias and uncertainty in these predictions we also conducted a set of post hoc GLM analyses in which the average suitability change and associated errors from these models were regressed against year, with interaction terms for SDM algorithm, emissions scenario and climate model. Significant interactions from these regressions indicate that a given component of the study design modifies the direction or magnitude of the bias in predictions of future habitat suitability. See Appendix S2 (Warren, Dornburg, et al., 2021) for details and R code for these analyses.
2.2 | Monte Carlo estimates of methodological bias
To estimate biases in our predictions for kangaskhan, we developed a novel Monte Carlo approach that is derived from a test developed by Raes and ter Steege (2007). As originally designed, this approach is used to generate null distributions of predictive performance for species distribution models. The null distribution is constructed by selecting localities at random from within the study area until a dataset of the same sample size as the empirical dataset is obtained. These random localities are then used to generate SDMs using the same methods as those used for the empirical model, and the predictive performance of the model is measured. By iteratively repeating this process, the user can construct an approximate distribution of expected predictive performance for an SDM under a given set of modelling conditions under the null hypothesis that the occurrence of the species within the study area is uncorrelated with any of the predictor variables. Comparing the discrimination accuracy of the empirical model to this null distribution allows users to test the hypothesis that their model has predictive power beyond what can be explained by their study design alone.
Our Monte Carlo approach (Figure 2) modifies this method to examine the biases implicit in transferring species distribution models under a given set of analytical conditions. As in Raes and ter Steege (2007) and Bohl et al. (2019), we select N points at random from the study area (i.e. the area from which background data were taken), where N is the number of occurrence points used for modelling in the empirical dataset. Using these random points as ‘occurrences’ we build SDMs with all other modelling choices being identical to those for the empirical dataset. We then project our replicate models onto the different combinations of climate model and emissions scenario, and calculate the predicted changes in habitat suitability for those models in the same way as with the empirical data. This allows us to estimate a distribution of expected changes in habitat suitability given the same study area, SDM algorithm, emissions scenario and climate model as our empirical data, under the hypothesis that the actual occurrence of the species within the study area is unrelated to the environmental variables used for model construction.
FIGURE 2 (click to open thumbnail)
Graphical sketch of Monte Carlo method for estimating bias in model construction and transfer. For simplicity we demonstrate these effects using a simple bioclim model, but a similar argument holds for any algorithm. In panel (a), we choose our study area (dashed lines) based on the distribution of our occurrence points (Gangura kangaskhani silhouettes). This affects the distribution of environments available during model construction and evaluation (dashed area, panel b). Using the occurrence points, we build a model that parameterizes the distribution of the species’ occurrences in environment space (dark rectangle, panels b and c). Given an emissions scenario and climate model, the distribution of available environments may differ in a future time period (dashed area, panel c). Hence, almost any bioclim‐style model built using points from the available environment in panel (b) will estimate as suitable some set of environmental conditions that are no longer present in panel (c). Similarly, models will lack information to determine the suitability of the environmental conditions available in panel (c) that were not present in panel (b). Climate envelope methods such as bioclim do not tend to extrapolate much beyond the conditions they were trained in, so the change in available environments creates a bias towards predictions of declining habitat suitability and/or range contractions. Other modelling algorithms may have similar built‐in biases, but the direction and magnitude may differ with their individual tendencies to extrapolate. It is worth noting that the bias inherent in a given study design will be affected by the change in available environments (which is affected by choice of study area, emissions scenario and climate model) as well as the shape and extent of the region of environment space predicted to be suitable (which are affected by study area, sample size and modelling algorithm). As such, all of these choices may affect our baseline expectation for what models will predict in a given study design independent of the actual locations of the species occurrences. To measure this bias, we repeatedly permute the locations of our data points within the study area and construct models using these random data (panel d), with every other aspect of the modelling process identical to the study design used for our empirical data. This eliminates any biological connection between the occurrence of our species and the environmental gradients, allowing us to estimate a distribution (panel e) of expected changes under the hypothesis that our occurrence data contain no information, and our predictions are consequently driven entirely by the bias implicit in the other choices made in our study design (sample size, study area, emissions scenario, climate model, etc.)
Graphical sketch of Monte Carlo method for estimating bias in model construction and transfer. For simplicity we demonstrate these effects using a simple bioclim model, but a similar argument holds for any algorithm. In panel (a), we choose our study area (dashed lines) based on the distribution of our occurrence points (Gangura kangaskhani silhouettes). This affects the distribution of environments available during model construction and evaluation (dashed area, panel b). Using the occurrence points, we build a model that parameterizes the distribution of the species’ occurrences in environment space (dark rectangle, panels b and c). Given an emissions scenario and climate model, the distribution of available environments may differ in a future time period (dashed area, panel c). Hence, almost any bioclim‐style model built using points from the available environment in panel (b) will estimate as suitable some set of environmental conditions that are no longer present in panel (c). Similarly, models will lack information to determine the suitability of the environmental conditions available in panel (c) that were not present in panel (b). Climate envelope methods such as bioclim do not tend to extrapolate much beyond the conditions they were trained in, so the change in available environments creates a bias towards predictions of declining habitat suitability and/or range contractions. Other modelling algorithms may have similar built‐in biases, but the direction and magnitude may differ with their individual tendencies to extrapolate. It is worth noting that the bias inherent in a given study design will be affected by the change in available environments (which is affected by choice of study area, emissions scenario and climate model) as well as the shape and extent of the region of environment space predicted to be suitable (which are affected by study area, sample size and modelling algorithm). As such, all of these choices may affect our baseline expectation for what models will predict in a given study design independent of the actual locations of the species occurrences. To measure this bias, we repeatedly permute the locations of our data points within the study area and construct models using these random data (panel d), with every other aspect of the modelling process identical to the study design used for our empirical data. This eliminates any biological connection between the occurrence of our species and the environmental gradients, allowing us to estimate a distribution (panel e) of expected changes under the hypothesis that our occurrence data contain no information, and our predictions are consequently driven entirely by the bias implicit in the other choices made in our study design (sample size, study area, emissions scenario, climate model, etc.)
2.3 | Interpretation
Our post hoc analyses of simulations from Warren et al. (2020) is intended to quantify the expected behaviour of our modelling approaches given a dataset in which there is an underlying (simulated) biological mechanism relating its probability of occurrence to a set of environmental predictors, while the novel Monte Carlo approach is intended to estimate the behaviour of a model built under identical conditions to our empirical model (including study region), but for which there is no biological mechanism underlying the distribution of occurrence points within that study area. The distributions of predictions produced by these two approaches are not necessarily intended to represent a null hypothesis that must be rejected in order to trust model results, but rather as estimates of the behaviour of our modelling approaches independent of our data. By comparing our empirical results to these distributions, we can determine to what extent our predictions are driven by our data, over and above the biases induced by our analytical framework.
3 | RESULTS
We find that the predicted change in suitability of habitat for G. kangaskhani is relatively consistent across spatial scales (continental vs. within species’ current range), but depends strongly on the method of analysis. Models built using Bioclim, Domain or GAM algorithms predict that habitat suitability will decline on average, but disagree substantially on the magnitude of that decline (Figure 3). GLM and Maxent models both make mixed predictions depending on the climate model used; both methods predict increased habitat suitability under the MIROC‐H climate model, but decreased suitability under the CSIRO model (Appendix S2, Warren, Dornburg, et al., 2021). Random forest models, by contrast, predict increasing habitat suitability regardless of emissions scenario or climate model.
FIGURE 3 (click to open thumbnail)
Projected effects of climate change on the suitability of habitat over time for Gangura kangaskhani. Results are shown for models projected to the continental scale (top) and within the species current range as defined by the 250‐km buffer around occurrence points (bottom). Each panel represents projections onto all combinations of emissions scenario and climate model. For similar plots separated by scenario and climate model, see Appendix S2 (Warren, Dornburg, et al., 2021)
Projected effects of climate change on the suitability of habitat over time for Gangura kangaskhani. Results are shown for models projected to the continental scale (top) and within the species current range as defined by the 250‐km buffer around occurrence points (bottom). Each panel represents projections onto all combinations of emissions scenario and climate model. For similar plots separated by scenario and climate model, see Appendix S2 (Warren, Dornburg, et al., 2021)
The models from Warren et al., (2020) were constructed from occurrence data based on simulated organisms with randomly selected niche parameters. However, the predictions made from these models demonstrate many of the same behaviours as those seen in the empirical models for G. kangaskhani (Figure 4, see appendix S2 in Warren, Dornburg, et al., 2021 for details). In particular we notice a broad trend for Bioclim and Domain models to predict declining habitat suitability, and for random forests to predict increasing habitat suitability over time at the continental scale. In contrast to the results from G. kangaskhani and the Monte Carlo replicates, random forests models from simulated species predict lower‐than‐present suitability within species’ current ranges in 2030. However, they show a similar tendency towards increasing estimates of habitat suitability over subsequent years. GLM, GAM and Maxent models show much more mixed predictions across the set of simulated species, which in the case of GLM and Maxent mirror our predictions for G. kangaskhani.
FIGURE 4 (click to open thumbnail)
Predicted change in suitability based on simulated species from Warren et al. (2020)
Predicted change in suitability based on simulated species from Warren et al. (2020)
As the Warren et al., (2020) data are based on simulated organisms for which the true suitability of habitat is known, we can project the true effects of climate change on those species and measure the extent to which our model predictions are inaccurate (Appendix S2, Warren, Dornburg, et al., 2021). At a continental scale, we find that the average true suitability of habitat changes little with time, but that the variance in suitability across simulated species increases substantially further in the future. Within simulated species’ current ranges, the modal true change in habitat suitability over time remains very close to zero, but the spread of the distribution of changes in true habitat suitability is increasingly biased downwards as predictions get further from the present.
Examining the distributions of errors in suitability predictions made by models shows that for the most part the general trends in behaviour of each modelling approach discussed above are in fact largely explained by biases in the kinds of errors each approach is making. GLM analyses of predictions of habitat suitability for these simulations and their accompanying errors demonstrate highly significant effects of every component of the modelling process (modelling algorithm, year, climate model, emissions scenario). This suggests that each one of these choices will affect the expected outcome of any modelling endeavour to a large (and typically un‐measured) degree.
The behaviour of models built from randomly chosen data points within the training area for G. kangaskhani showed patterns that were generally concordant with those from the simulated species, and in many cases the predictions for G. kangaskhani fall well within the distribution of changes from those Monte Carlo replicates (Figure 5). In most cases where the behaviour of the G. kangaskhani model did fall outside of the distribution estimated from random points, it was in the direction of making less extreme predictions for the effects of climate change on G. kangaskhani than expected.
FIGURE 5 (click to open thumbnail)
Distribution of predicted changes in suitability values from Monte Carlo replicates based on the study area for Gangura kangaskhani. Black circles represent predictions of change in habitat suitability from empirical models
Distribution of predicted changes in suitability values from Monte Carlo replicates based on the study area for Gangura kangaskhani. Black circles represent predictions of change in habitat suitability from empirical models
4 | DISCUSSION
In this study, we used a citizen science dataset in conjunction with climate data to estimate the effects of climate change on Australia's only endemic Pokémon, and find that our conclusions about changes in habitat suitability are profoundly affected by methodological choices. By examining multiple simulated data sources, we highlight the important possibility that the results seen here might not be indicative of any direct effect of climate change on the suitability of habitat for G. kangaskhani, but instead represent the effects of bias that is imposed by the methods used to construct these models. The predicted effects of climate change on the suitability of habitat for G. kangaskhani generally reflect predictions we derive from models for simulated species as well as predictions from Monte Carlo tests based on randomly chosen data points within the range of G. kangaskhani.
It is interesting also to compare our results with Beaumont et al. (2016), a study that used models built with real data for 220 Australian mammal species to examine which modelling approaches were more prone to making predictions of range contractions or expansions (Figure 6). Those models were based on real data for species for which the climate niche is not known, and as a result it was not possible from that study to say which of those models better estimates the true impacts of climate change. However, we note that many of the model behaviours seen in Beaumont et al. (2016) echo those seen in our simulated data and Monte Carlo replicates; climate envelope methods show a strong tendency towards predictions of declining suitability and random forests show the strongest tendency to predict increasing suitability, while Maxent, GLM and GAM all produce predictions of change in habitat suitability that are relatively symmetric around zero (Figures 4–6-4–6). The similarity between the results for real species, simulated species and Monte Carlo replicates strongly suggests that the effects of the biases inherent in the design of SDM analyses are a ubiquitous feature of studies of this type.
FIGURE 6 (click to open thumbnail)
Projected proportional range change for 220 Australian mammal species under two different climate change scenarios, modified from Beaumont et al. (2016) Figure 3. In that study models were built using fourteen different algorithms, but here the figure is modified to include only the algorithms used for Gangura kangaskhani to facilitate comparison with simulation and Monte Carlo results. Proportional range change is calculated as the future area of suitable habitat divided by the present area of suitable habitat based on binary predictions made from thresholded models. As such it is not directly comparable to the change in average suitability presented here for the other analyses, but the two are expected to be highly correlated (i.e. declining suitability scores will typically lead to a predicted reduction in the area of suitable habitat)
Projected proportional range change for 220 Australian mammal species under two different climate change scenarios, modified from Beaumont et al. (2016) Figure 3. In that study models were built using fourteen different algorithms, but here the figure is modified to include only the algorithms used for Gangura kangaskhani to facilitate comparison with simulation and Monte Carlo results. Proportional range change is calculated as the future area of suitable habitat divided by the present area of suitable habitat based on binary predictions made from thresholded models. As such it is not directly comparable to the change in average suitability presented here for the other analyses, but the two are expected to be highly correlated (i.e. declining suitability scores will typically lead to a predicted reduction in the area of suitable habitat)
In this study, our Monte Carlo tests focus on how a given set of modelling conditions may bias our estimates of the effects of climate change on average habitat suitability. However, we note that this framework can easily be extended. There are many different metrics that may be used to quantify the effects of climate change on species (e.g. change in habitable area based on thresholded suitability scores, geographic shifts in ranges or centroids, changes in climate niche breadth, etc.). Applying this approach to those metrics is simply a matter of applying the chosen metric to predictions from each Monte Carlo replicate to obtain the distribution of that metric given randomly sampled occurrence data. Similarly, this approach could easily be extended to obtain distributions for models being transferred in space, rather than in time; transferring models in time and space present very similar challenges for SDM studies, and it is likely that similar biases exist in any study of that type. Further, it is possible to calculate summary statistics on a per‐grid‐cell basis from the suitability maps for Monte Carlo replicates, producing a map of the geographic regions where model predictions might be driven by study design more than by the data used in modelling (Figure 7). Finally, we note that in the current study occurrence points for Monte Carlo replicates were sampled with uniform probability across the study with no spatial autocorrelation or sampling bias. While this approach has been widely used for constructing null distributions of model performance (Bohl et al., 2019; Raes & ter Steege, 2007; Warren, Dornburg, et al., 2021), other methods that attempt to maintain similar spatial autocorrelation to empirical data when sampling may yield different null distributions (Beale et al., 2008; Nunes & Pearson, 2017; Roxburgh & Chesson, 1998). We consider this a very promising area for future study, as it may help to better understand which aspects of data collection and curation are most instrumental in reducing bias in predictions.
FIGURE 7 (click to open thumbnail)
Using the 100 Monte Carlo replicates for the Maxent models for Gangura kangaskhani, we recorded the mean change in suitability (panel a) and the standard deviations of those predictions (panel b) when predicted to all four combinations of emissions scenario and climate model for the year 2100. As these models represent the behaviour of the study design given uninformative occurrence data, the desired behaviour is for the mean predicted change (panel a) to be near zero (indicating no directional bias), and the standard deviation around that mean (panel b) to be high (indicating that model predictions are sensitive to changes in the data). Although there are some regions that approximately fit this set of desired outcomes, we note that there are areas along the northern and southern coasts where predicted changes are largely insensitive to the data (SD is low, panel b), and these areas include predictions of both increased (northern coast and Tasmania) and decreased (southwestern coastal areas) habitat suitability from the present day (panel a). These indicate regions where Maxent models are biased towards predicting a specific direction of change in habitat suitability, and where the data have little leverage with which to counteract these biases
Using the 100 Monte Carlo replicates for the Maxent models for Gangura kangaskhani, we recorded the mean change in suitability (panel a) and the standard deviations of those predictions (panel b) when predicted to all four combinations of emissions scenario and climate model for the year 2100. As these models represent the behaviour of the study design given uninformative occurrence data, the desired behaviour is for the mean predicted change (panel a) to be near zero (indicating no directional bias), and the standard deviation around that mean (panel b) to be high (indicating that model predictions are sensitive to changes in the data). Although there are some regions that approximately fit this set of desired outcomes, we note that there are areas along the northern and southern coasts where predicted changes are largely insensitive to the data (SD is low, panel b), and these areas include predictions of both increased (northern coast and Tasmania) and decreased (southwestern coastal areas) habitat suitability from the present day (panel a). These indicate regions where Maxent models are biased towards predicting a specific direction of change in habitat suitability, and where the data have little leverage with which to counteract these biases
Species distribution models are widely used with the goal of guiding management and policy. However, the fact that we produce consistent predictions from meaningless data indicates that there are some combinations of emissions scenario, climate model, study area and species distribution modelling method that are strongly biased towards a given prediction. Crucially, in some modelling conditions we find that all Monte Carlo replicates agreed on the direction of change in habitat suitability in the future, despite the fact that the points for these replicates were chosen at random from within Kangaskahan's range and therefore represented no underlying relationship between the environment and the probability of occurrence. At the very least this suggests that those study designs should be chosen only when there is a compelling reason to believe that they accurately represent the underlying biology. In the case of G. kangaskhani we have insufficient information to make such a decision; in addition to our scant knowledge of the species’ ecology we anticipate that there are likely to be serious issues with data quality due to spatial sampling bias and small sample size due to an unfortunate lack of scientific attention to wild Pokémon populations. In particular, we would like to point out that there may be some room for optimism given our level of ignorance of the distribution of G. kangaskhani; we do not currently know how far into the interior of Australia these animals might be found and so are unlikely to have fully characterized its environmental tolerances. Currently very little conservation funding is allocated for land acquisition and captive breeding programs for Pokémon, but increased investment in this area could be broadly beneficial; charismatic megavertebrates can often be used as flagship species to increase community involvement in conserving whole ecosystems (Cummins, 2020). In order to effectively manage this majestic species, we suggest that future funding efforts be directed towards sending the authors of this study into remote areas of Australia to collect more reliable data, preferably equipped with a large supply of the ‘incense’ (presumably some sort of pheromone) that Pokémon hunters use to attract animals for capture.
It is reasonable to view with heavy scepticism the results of any study in which the analytical design is only capable of producing one answer regardless of the data. This may not mean, however, that the predictions from such a study are necessarily wrong. To illustrate this point, it is worth considering an extreme case in which the entire study area becomes unsuitable for any organism currently living within it (e.g. in future scenarios all currently occupied habitat is constantly on fire). In this case, any reasonable model built from any set of random occurrence data in that region would show complete extirpation of the organism, which would also coincide with the results seen from any carefully vetted empirical dataset for a real species. The conclusions of those models would be correct, but they would be wholly unaffected by the data; they would be an inevitable result of the modelling conditions. Any dataset would likely support the same conclusions regardless of whether those data represent any underlying biological phenomenon.
In the real world, however, the effects of climate change are rarely so absolute, and are rarely known in advance. Finding that we are strongly biased towards making specific predictions irrespective of our data is therefore a cause for great concern; it implies that we may have essentially selected our conclusions by selecting our study design. In such a situation, the resulting outputs should at minimum be interpreted in the context of the biological plausibility of those methodological choices, the biases they induce, and their consequences for decisions made based on those models. In particular, if we find ourselves in a situation where different approaches show strong opposing biases, interpretation of those results must be made in the context of how biologically realistic the alternative modelling approaches are. Arguably one should go further than that, however, and simply not trust the outcome of a modelling approach that produces the same result regardless of the data provided to it. In these cases, additional lines of evidence are needed to adequately inform conservation decisions. Understanding when we are in such a situation is the key to making these decisions, and we believe that Monte Carlo tests such as the one presented here may be a valuable tool for measuring the biases inherent in the design of empirical studies.
ACKNOWLEDGEMENTS
Gotta thank ‘em all! The authors would like to thank the Westpunt Syndicate and acknowledge funding sources including an NSF GRFP grant to K.Z. (#0000‐0002‐5159‐641X). We are deeply grateful for many helpful comments from Erica Newman, Linda Beaumont, Megan Higgie, Jamie Kass, Nick Matzke, Peter Linder, Liam Langan, Jorge Lobo and an anonymous reviewer. We are also very grateful to Bob O'Hara and the rest of the editorial staff at MEE for indulging this foolishness. We would also like to thank all of the Pokémon trainers out there who made this possible. Finally the authors would like to acknowledge that this started out as a joke but actually got pretty dang interesting. The use of kangaskhan, Pokémon and related intellectual property here is intended as non‐commercial parody fair use, and does not in any way reflect an endorsement by Niantic, The Pokémon Company or Prof. Y. Ohkido. Open Access funding enabled and organized by Projekt DEAL. WOA Institution: SENCKENBERG GESELLSCHAFT FUR NATURFORSCHUNG; Blended DEAL: Projekt DEAL.
AUTHORS' CONTRIBUTIONS
Initial conception by D.L.W.; Study design by D.L.W., T.L.I., and A.D., analyses by D.W. The manuscript was written with equal contributions from all authors. Original artwork by K.Z.
REFERENCES
- Araújo, M. B., Anderson, R. P., Barbosa, A. M., Beale, C. M., Dormann, C. F., Early, R., Garcia, R. A., Guisan, A., Maiorano, L., Naimi, B., O’Hara, R. B., Zimmermann, N. E., & Rahbek, C. (2019). Standards for distribution models in biodiversity assessments. Science Advances, 5(1), eaat4858.
- Barbet‐Massin, M., Thuiller, W., & Jiguet, F. (2010). How much do we overestimate future local extinction rates when restricting the range of occurrence data in climate suitability models? Ecography, 33(5), 878– 886.
Beale, C. M., Lennon, J. J., & Gimona, A. (2008). Opening the climate envelope reveals no macroscale associations with climate in European birds. Proceedings of the National Academy of Sciences of the United States of America, 105(39), 14908– 14912. - Beaumont, L. J., Graham, E., Duursma, D. E., Wilson, P. D., Cabrelli, A., Baumgartner, J. B., Hallgren, W., Esperón‐Rodríguez, M., Nipperess, D. A., Warren, D. L., Laffan, S. W., & VanDerWal, J. (2016). Which species distribution models are more (or less) likely to project broad‐scale, climate‐induced shifts in species ranges? Ecological Modelling, 342(December), 135– 146.
- Bellard, C., Bertelsmeier, C., Leadley, P., Thuiller, W., & Courchamp, F. (2012). Impacts of climate change on the future of biodiversity. Ecology Letters, 15(4), 365– 377.
- Bohl, C. L., Kass, J. M., & Anderson, R. P. (2019). A new null model approach to quantify performance and significance for ecological niche models of species distributions. Journal of Biogeography, 46(6), 1101– 1111.
- Bulbapedia. (2020a). Pokémon Food. Bulbapedia. Retrieved from Pokémon food - Bulbapedia, the community-driven Pokémon encyclopedia
- Bulbapedia. (2020b). Kangaskhan (Pokémon). Bulbapedia. Retrieved from Kangaskhan (Pokémon) - Bulbapedia, the community-driven Pokémon encyclopedia
- Carpenter, G., Gillison, A. N., & Winter, J. (1993). DOMAIN: A flexible modelling procedure for mapping potential distributions of plants and animals. Biodiversity & Conservation, 2(6), 667– 680.
- Chen, I. C., Hill, J. K., Ohlemüller, R., Roy, D. B., & Thomas, C. D. (2011). Rapid range shifts of species associated with high levels of climate warming. Science, 333(6045), 1024– 1026.
- Cummins, E. (2020). What Pokémon can teach us about conservation and climate change. The Verge. Retrieved from https://www.theverge.com/2020/1/30/21115034/pokemon‐conservation‐climate‐change‐detective‐pikachu
- D’Anastasio, C. (2016). A dark theory on what happens to transferred Pokemon. Kotaku. Retrieved from https://www.kotaku.com.au/2016/07/a‐dark‐theory‐on‐what‐happens‐to‐transferred‐pokmon/
- Datta, A., Schweiger, O., & Kühn, I. (2020). Origin of climatic data can determine the transferability of species distribution models. NeoBiota, 59(July), 61– 76.
- Galante, P. J., Alade, B., Muscarella, R., Jansa, S. A., Goodman, S. M., & Anderson, R. P. (2018). The challenge of modeling niches and distributions for data‐poor species: A comprehensive approach to model complexity. Ecography, 41(5), 726– 736.
- Gilbert, B. D. (2020). Pokémon Edibility | Unraveled. Polygon. Retrieved from Polygon YouTube channel
- Hurlbert, S. H. (1990). Spatial distribution of the montane unicorn. Oikos, 58(3), 257– 271.
- Hutchinson, G. E. (1978). An introduction to population biology (p. 260). Yale University Press.
- KOPlayer. (2018). Pokemon Catch. KO Player. Retrieved from http://www.koplayer.com/act/PokemonCatch/
- Kriticos, D. J., Webber, B. L., Leriche, A., Ota, N., Macadam, I., Bathols, J., & Scott, J. K. (2012). CliMond: Global high‐resolution historical and future scenario climate surfaces for bioclimatic modelling. Methods in Ecology and Evolution/British Ecological Society, 3(1), 53– 64.
- Loe, L. E., Bonenfant, C., Meisingset, E. L., & Mysterud, A. (2012). Effects of spatial scale and sample size in GPS‐based species distribution models: Are the best models trivial for red deer management? European Journal of Wildlife Research. https://doi.org/10.1007/s10344‐011‐0563‐5
- Lozier, J. D., Aniello, P., & Hickerson, M. J. (2009). Predicting the distribution of sasquatch in western North America: Anything goes with ecological niche modelling. Journal of Biogeography, 36(9), 1623– 1627.
- Michelot, T., Langrock, R., & Patterson, T. A. (2016). moveHMM: an R package for the statistical modelling of animal movement data using hidden Markov models. Methods in Ecology and Evolutio, 7(11): 1308– 1315.
- Molloy, D. (2013). Pokemon is actually a pretty brutal blood sport. The Irish Examiner. Retrieved from https://www.irishexaminer.com/breakingnews/discover/pokemon‐is‐actually‐a‐pretty‐brutal‐blood‐sport‐608939.html
- Nix, H. A., & Busby, J. (1986). BIOCLIM, a bioclimatic analysis and prediction system. Annual Report CSIRO, CSIRO Division of Water and Land Resources.
- Nunes, L. A., & Pearson, R. G. (2017). A null biogeographical test for assessing ecological niche evolution. Journal of Biogeography, 44(6), 1331– 1343.
- O’Brien, K., Leichenko, R., Kelkar, U., Venema, H., Aandahl, G., Tompkins, H., Javed, A., Bhadwal, S., Barg, S., Nygaard, L., & West, J. (2004). Mapping vulnerability to multiple stressors: Climate change and globalization in India. Global Environmental Change. Redirecting
- Parmesan, C. (2006). Ecological and evolutionary responses to recent climate change. Retrieved from https://doi.org/10.1146/annurev.ecolsys.37.091305.110100
- Peterson, A. T., Soberón, J., Pearson, R. G., Anderson, R. P., Martínez‐Meyer, E., Nakamura, M., & Araújo, M. B. (2011). Ecological niches and geographic distributions (MPB‐49). Princeton University Press.
- Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190(3), 231–259.
- Phillips, T. (2018). Pokémon go active player count highest since 2016 summer launch. Eurogamer. Retrieved from https://www.eurogamer.net/articles/2018‐06‐27‐pokemon‐go‐player‐count‐at‐highest‐since‐2016‐summer‐launch
- PKMNGoTrading. (2016). Where to Find and Catch Kangaskhan. PKMN Go Trading. Retrieved from Where to Find and Catch Kangaskhan - Pokemon Go Wiki
- PokéBase. (2020). Kangaskhan. Pokémon Database. Retrieved from Kangaskhan Pokédex: stats, moves, evolution & locations
- Pokémon GO Wiki Community. (n.d.). Pokémon GO Safari Zone. Pokémon GO Wiki. Retrieved from Pokémon GO Safari Zone
- PokemonGoHub. (2020). Pokemon Go Hub: Kangaskhan. Pokemon Go Hub. Retrieved from https://db.pokemongohub.net/pokemon/115
- Raes, N., & ter Steege, H. (2007). A null‐model for significance testing of presence‐only species distribution models. Ecography, 30(5), 727– 736.
- Roxburgh, S. H., & Chesson, P. (1998). A new method for detecting species associations with spatially autocorrelated data. Ecology, 79(6), 2180– 2192.
- Soberon, J., & Nakamura, M. (2009). Niches and distributional areas: Concepts, methods, and assumptions. Proceedings of the National Academy of Sciences of the United States of America. https://doi.org/10.1073/pnas.0901637106
- The Silph Road. (2016). How Spawns Work in Pokemon GO ‐ Research from The Silph Road. Youtube. June 30, 2016. Retrieved from The Silph Road YouTube Channel
- The Silph Road Community. (2016). The Silph Road Pokémon Go Nest Atlas. The Silph Road. January 2016. Retrieved from https://thesilphroad.com/atlas/#3.64/‐24.68/137.16
- Warren, D. L., Dornburg, A., Zapfe, K., & Iglesias, T. L. (2021). Data and code for analysis of effects of climate change on kangaskhan and summary of simulations from Warren et al 2020. Dryad Digital Repository, Dryad Data -- Data and code for analysis of effects of climate change on kangaskhan and summary of simulations from Warren et al. 2020
- Warren, D. L., Matzke, N. J., Cardillo, M., & Baumgartner, J. B., Beaumont, L. J., Turelli, M., Glor, R. E., Huron, N. A., Simões, M., Iglesias, T. L., Piquet, J. C., & Dinnage, R. (2021). ENMTools 1.0: an R package for comparative ecological biogeography. Ecography (early view). https://doi.org/10.1111/ecog.05485
- Warren, D. L., Matzke, N. J., & Iglesias, T. L. (2020). Evaluating presence‐only species distribution models with discrimination accuracy is uninformative for many applications. Journal of Biogeography, 47(1), 167– 180.
- Wikipedia. (2020). Pokémon Go. Wikipedia. Retrieved from Pokémon Go - Wikipedia
Appendix S1: Additional details about the life history of Garura kangaskhani and comments on how these may contribute to their response to climate change
Biology and life history
Despite being a charismatic megavertebrate (~220 cm adult height) that is frequently sighted near human population centers, G. kangaskhani was not formally described until 1997 (Bulbapedia, 2019). Information regarding its basic biology is scarce and much of what is available is based on anecdotal accounts from citizen scientists and unlicensed breeders, many of whom have no formal scientific training. However, what information is available suggests that the reproductive system and demography of G. kangaskhani may be unusual in a number of ways, and these factors may need to be taken into account when making conservation decisions for the species. Therefore, in addition to addressing the effects of climate change on G. kangaskhani we offer the following information as a summary of the current consensus on the biology of the species, with the hopes that it will inspire further investigation and guide conservation efforts.
Garura kangaskhani has never been sequenced for DNA, and therefore all available information about its genome is based on speculation (Butterfree, 2016). As a consequence of the lack of genetic information, the evolutionary history of G. kangaskhani remains contentious. This lineage is hypothesized to be either sister to a lineage leading to a clade of “fire type” Pokémon that includes Entei, Heatran, and Houndoom (Shelomi et al., 2012) or a member of a distinct lineage of ‘marsupial’ Pokémon that includes Cubone (Spaniard, 2016). While its breeding habits suggest that it is a marsupial, it possesses a series of dermal plates that are reminiscent of stegosaurs. This unusual feature may represent convergent evolution of a useful defense strategy or an ancient hybridization event between a stegosaur and a stem-group metatherian (Margulis & Sagan, 2008; Williamson, 2009). As G. kangaskhani has never been seen to evolve (PokéBase, 2020) it is further possible that this taxon represents a “living fossil”. However, disentangling these and other hypotheses of evolutionary history is stymied by this Pokémon’s enigmatic reproductive biology.
The reproductive biology of G. kangaskhani has received more attention than most other aspects of its biology, but what is known is perplexing and often contradictory. All sources agree that offspring are strongly altricial (Bulbapedia, 2020), living within the mother’s pouch for approximately three years. Mothers are said to form strong pair bonds with their young, and sleep standing upright in order to protect them (Figure 1). Interestingly, all observed individuals of G. kangaskhani have been reported as female. We suspect that this is largely due to the presence of young in the pouch, as detailed anatomical examinations have not been performed to our knowledge. We briefly explore two alternative hypotheses to explain this pattern: that male G. kangaskhani do exist but have not been correctly documented, and that G. kangaskhani are parthenogenic.
The lack of male G. kangaskhani observations may stem from males being rare, necessitating more extensive field environmental-genomic surveys to verify (Wilson Sayres, 2018). However, it is also possible that many male G. kangaskhani have simply been misidentified. For example, males may be identified as female based on the assumption of uniparental care. While in most marsupial species pouches are only present in females, there are known exceptions including the water opossum (or “yapok”; Chironectes minimus) and the recently extinct thylacine (Thylacinus cynocephalus), in which males have a pouch that is used to protect the genitals from snagging on vegetation (Nogueira et al., 2004; Paddle, 2002). It is therefore possible that male G. kangaskhani are using a highly modified intromittent organ as a nuptial display for females and that the rare ‘shiny’ forms are in fact males with display coloration. Alternatively, G. kangaskhani males may simply appear so different from females as to have been misidentified as another species. This is a mistake that has been made many times in other organisms (Moraes et al., 2017; Perez-Mellado & de la Riva, 1993). It is even possible that the male G. kangaskhani is in fact many times smaller than the female and either attached to, or embedded within, the female of the species (Pietsch, 2005) (Figure S1.1).
Figure S1.1.
G. kangaskhani male embedding under dorsal dermal plates. Females may carry a fused dispersal-morph male (see Figure S1.2 (C) and (D)) under each dermal plate at differing fusion stages, and accept another once the previous male has become fully absorbed. Over time, scar tissue builds up at the fusion site, raising the plates.
G. kangaskhani male embedding under dorsal dermal plates. Females may carry a fused dispersal-morph male (see Figure S1.2 (C) and (D)) under each dermal plate at differing fusion stages, and accept another once the previous male has become fully absorbed. Over time, scar tissue builds up at the fusion site, raising the plates.
The other major possibility, and currently the favored hypothesis (Dyer, 2016), is that G. kangaskhani reproduce via parthenogenesis. While parthenogenesis is primarily associated with plant and invertebrate species, it has been observed in a number of vertebrate species (Booth et al., 2012; Fields et al., 2015; Fujita et al., 2007; Moritz & Heideman, 1993; Walker, 1986). Although the idea may seem somewhat farfetch’d, it is supported by the repeated observation that young G. kangaskhani are born with their own offspring in pouch. This observation would be hard to reconcile with the above hypotheses regarding cryptic males, but is quite compatible with a parthenogenic mating system, as a similar phenomenon is also seen in parthenogenic aphids (Stern, 2008).
The effects of multiple stressors
The results of this study indicate that climate change may have significant negative impacts on G. kangaskhani, as indicated by the decline in overall average habitat suitability. This is likely to compound the effects of existing threats to populations of G. kangaskhani, which are already under serious pressure from human exploitation. In particular, the presence of high density populations in areas of anthropogenic disturbance is cause for concern. While these areas of exceptional local abundance are often referred to as “nesting” habitat, there is no obvious reason to believe that they nest colonially under those conditions, and the current consensus is that the species may not nest communally at all (The Silph Road Community, 2016). It is not clear why these comparatively high densities are seen in human-modified environments, but we speculate that it is driven either by local food resources or by human activity uncovering deposits of kangaskhanite, a mineral they need for development. Given that these high density locations are also favored by hunters looking to capture G. kangaskhani for sport, these areas may in fact constitute an “ecological trap” (Hale & Swearer, 2016); extreme sink populations driven by a combination of attractive resources and high levels of unregulated exploitation. If this is the case, the effects of climate change on surrounding natural source populations may be compounded by emigration from those populations, leading to a potential for sudden and irreversible collapse. We see little room for optimism in this scenario given that the culture around Pokémon harvesting so explicitly favors exploitation, as seen in the common refrain of “gotta catch ‘em all” (Prof. Oak, pers. comm.).
Climate change places significant selective pressures on natural populations, and it is thought that some species may be able to adapt to changing climate regimes via natural selection (Catullo et al., 2019). The fact that G. kangaskhani are likely parthenogenic makes this significantly less likely; absence of sexual reproduction reduces effective population size and precludes recombination, leading to a curtailed ability to respond to selection (Rice & Friberg, 2009). However, we would also like to highlight again the similarity between what is known about G. kangaskhani reproduction and the reproductive system of the aphid (Stern, 2008), and in particular we note that aphids are capable of facultative sexual reproduction when the environment becomes unpredictable, increasing genetic variability of offspring in order to facilitate responses to selective pressures. If the similarities between G. kangaskhani and aphids extend to this response, it is possible that climate change will induce a wave of previously unseen sexual Kangaskhan morphs (Figure S1.2). The presence of large numbers of G. kangaskhani in human population centers makes this a somewhat alarming prospect, even more so if the sexual morphs are territorial.
Figure S1.2.
G. kangaskhani adult females (A) may produce sexual morphs, including territorial males (B), or considerably smaller males capable of dispersing over greater distances by flight (C) or hopping similar to Dipodomys spp. (D).
G. kangaskhani adult females (A) may produce sexual morphs, including territorial males (B), or considerably smaller males capable of dispersing over greater distances by flight (C) or hopping similar to Dipodomys spp. (D).
- Booth, W., Smith, C. F., Eskridge, P. H., Hoss, S. K., Mendelson, J. R., 3rd, & Schuett, G. W. (2012). Facultative parthenogenesis discovered in wild vertebrates. Biology Letters, 8(6), 983–985.
- Bulbapedia. (2019, June 9). EP034. Bulbapedia. EP034 - Bulbapedia, the community-driven Pokémon encyclopedia
- Bulbapedia. (2020, March 8). Kangaskhan (Pokémon). Bulbapedia. https://bulbapedia.bulbagarden.net/wiki/Kangaskhan_(Pokêmon)
- Butterfree. (2016). Pokémon Genetics. The Cave of Dragonflies. Pokémon Genetics
- Catullo, R. A., Llewelyn, J., Phillips, B. L., & Moritz, C. C. (2019). The Potential for Rapid Evolution under Anthropogenic Climate Change. Current Biology: CB, 29(19), R996–R1007.
- Dyer, J. (2016, April 1). Kangaskhan – Parental Bond and Parthenogenesis. The Biology of Pokémon. https://pokemonbiology.com/2016/04/01/kangaskhan-parental-bond-and-parthenogenesis/
- Fields, A. T., Feldheim, K. A., Poulakis, G. R., & Chapman, D. D. (2015). Facultative parthenogenesis in a critically endangered wild vertebrate. Current Biology: CB, 25(11), R446–R447.
- Fujita, M. K., Boore, J. L., & Moritz, C. (2007). Multiple origins and rapid evolution of duplicated mitochondrial genes in parthenogenetic geckos (Heteronotia binoei; Squamata, Gekkonidae). Molecular Biology and Evolution, 24(12), 2775–2786.
- Hale, R., & Swearer, S. E. (2016). Ecological traps: current evidence and future directions. Proceedings. Biological Sciences / The Royal Society, 283(1824). https://doi.org/10.1098/rspb.2015.2647
- Margulis, L., & Sagan, D. (2008). Acquiring genomes: A theory of the origin of species. Basic books.
- Moraes, S. D. S., Cardoso, L. W., Silva-brandÃo, K. L., & Duarte, M. (2017). Extreme sexual dimorphism and polymorphism in two species of the tiger moth genus Dysschema (Lepidoptera: Erebidae): association between males and females, sexual mimicry and melanism revealed by integrative taxonomy. Systematics and Biodiversity, 15(3), 259–273.
- Moritz, C., & Heideman, A. (1993). The Origin and Evolution of Parthenogenesis in Heteronotia Binoei (Gekkonidae): Reciprocal Origins and Diverse Mitochondrial DNA in Western Populations. In Systematic Biology (Vol. 42, Issue 3, pp. 293–306). Origin and Evolution of Parthenogenesis in Heteronotia Binoei (Gekkonidae): Reciprocal Origins and Diverse Mitochondrial DNA in Western Populations
- Nogueira, J. C., Castro, A. C. S., Câmara, E. V. C., & Câmara, B. G. O. (2004). Morphology of the Male Genital system of Chironectes minimus and Comparison to other didelphid marsupials. Journal of Mammalogy, 85(5), 834–841.
- Paddle, R. (2002). The Last Tasmanian Tiger: The History and Extinction of the Thylacine. Cambridge University Press.
- Perez-Mellado, V., & de la Riva, I. (1993). Sexual Size Dimorphism and Ecology: The Case of a Tropical Lizard, Tropidurus melanopleurus (Sauria: Tropiduridae). Copeia, 1993(4), 969–976.
- Pietsch, T. W. (2005). Dimorphism, parasitism, and sex revisited: modes of reproduction among deep-sea ceratioid anglerfishes (Teleostei: Lophiiformes). Ichthyological Research, 52(3), 207–236.
- PokéBase. (2020, March 12). Kangaskhan. Pokémon Database. Kangaskhan Pokédex: stats, moves, evolution & locations
- Rice, W. R., & Friberg, U. (2009). A Graphical Approach to Lineage Selection Between Clonals and Sexuals. In I. Schön, K. Martens, & P. Dijk (Eds.), Lost Sex: The Evolutionary Biology of Parthenogenesis (pp. 75–97). Springer Netherlands.
- Shelomi, M., Richards, A., Li, I., & Okido, Y. (2012). A phylogeny and evolutionary history of the Pokémon. Annals of Improbable Research, 18(4), 15.
- Spaniard, I. (2016, December 15). Pokemon Tree of life for generation 7. Reddit. r/gaming - Pokemon Tree of life for generation 7
- Stern, D. L. (2008). Aphids. Current Biology: CB, 18(12), R504–R505.
- Walker, J. M. (1986). The Taxonomy of Parthenogenetic Species of Hybrid Origin: Cloned Hybrid Populations of Cnemidophorus (Sauria: Teiidae). In Systematic Zoology (Vol. 35, Issue 3, p. 427). The Taxonomy of Parthenogenetic Species of Hybrid Origin: Cloned Hybrid Populations of Cnemidophorus (Sauria: Teiidae) on JSTOR
- Williamson, D. I. (2009). Caterpillars evolved from onychophorans by hybridogenesis. Proceedings of the National Academy of Sciences of the United States of America, 106(47), 19901–19905.
- Wilson Sayres, M. A. (2018). Genetic Diversity on the Sex Chromosomes. Genome Biology and Evolution, 10(4), 1064–1078.
- Sources
Unique macro-vertebrate at risk from blood sport and climate change
A two-meter-tall Pokémon, endemic to Australia, was the focus of a new study, which has revealed that some species distribution models are inherently biased.www.oist.jp
