1. Predictions of biodiversity trajectories under climate change are crucial in order to act effectively in maintaining the diversity of species. In many ecological applications, future predictions are made under various global warming scenarios as described by a range of different climate models. We propose a clustering methodology to synthesize and interpret the outputs of these various predictions. 2. We propose an interpretable and flexible two step methodology to measure the similarity between predicted species range maps and to cluster the future scenario predictions utilizing a spectral clustering technique. 3. We find that clustering based on predicted species range maps is mainly driven by the amount of warming rather than climate model or future scenario. We contrast this with clustering based only on predicted climate variables, which is driven primarily by climate models, i.e., scenarios of the same climate model are clustered together, even when the amount of warming input to the models is varied. 4. The differences between species-based and climate-based clusterings illustrate that it is crucial to incorporate ecological information to understand the relevant differences between climate models. Our findings can be used to better synthesize forecasts of biodiversity change under the wide spectrum of results that emerge when considering potential future scenarios. Natural and human-induced global climate change undeniably threatens the biodiversity of plants and animals and is expected to lead to large-scale range loss (Lovejoy, 2006;Mackey et al., 2008;Dantas-Torres, 2015;Radchuk et al., 2019). A common response of species to global climate change is dispersal or range map shifting (Pecl et al., 2017) or range map contraction. Understanding the extent of potential ecological responses to a rapidly changing climate is essential to enact effective conservation policies (Parry et al., 2007;Hannah et al., 2013Hannah et al., , 2020)). Forecasting species distributions or range maps under various scenarios of future climate provides a suite of potential outcomes (Burrows et al., 2014;Jones and Cheung, 2015;Molinos et al., 2016). Given a model for species range maps based on current climate, a suite of forecasts can be made using future climate projections such as those from the 6 (CMIP6). A comprehensive look at the forecasts under various scenarios for a large number of species can quickly become burdensome to summarize both within and among species. Clustering methods provide a mechanism for extracting lower dimensional information; clustering species range maps across scenarios elicits patterns of responses to climate change within a species and across species. Principal component analysis (PCA) is a widely used tool in ecological applications to interpret highdimensional data (Pearson, 1901;Jolliffe, 1986;Huettmann and Diamond, 2001;Janzekovic and Novak, 2012;Wu et al., 2019). However, PCA has significant drawbacks in the nonlinear, multivariate setting of the current application. PCA has a tight relationship with correlation which classically fails to capture nonlinear relationships. Secondly, it is not clear how to incorporate information from multiple species into a PCAbased approach, since different species have different sized range maps, and thus different species will naively vary in the principal component contributions due to size alone, which may not be desired. In contrast, we present a flexible alternative, spectral clustering, allowing for any similarity or distance between range maps. Spectral clustering is ideal for our dataset because of the high-dimensional range maps and the desire to amalgamate information from across species. Spectral clustering is a technique to cluster observations given only pairwise similarity or distance measures between the observations. It is well suited for this task because we can define a function that captures our desired notion of pairwise similarity measure between potential range maps, which will be discussed in Section 2.3.2. A popular alternative summarisation to clustering is a metric of change in species richness due to climate change known as climate velocity (Loarie et al., 2009) which can be used as a proxy for net effects on biodiversity change, particularly when detailed biodiversity information is not available (Burrows et al., 2014;Jones and Cheung, 2015;Molinos et al., 2016). However, these climate velocities rely entirely on climate information while ignoring ecological data. Such ecological information, such as species' exposure to climate conditions not found in the species' current niches that are described by species range maps models, are important factors to predict the species potential future range maps (Trisos et al., 2020). We propose to cluster the forecasted species range maps to interpret and investigate the similarities and differences among predictions across different warming scenarios, global climate models (GCMs), and species. We analyze a dataset of 1101 animals over 34 different future scenarios, for a total of over 37, 000 range maps. Ecological data informs the predicted range maps, and thus the spectral clustering results will make full use of the available ecological information. In order to assess the amount of information gained from the species range map model, we additionally apply spectral clustering to the scenarios based only on the predicted values of the climate variables. We hypothesize that clustering the scenarios based on these predicted range maps will produce significantly different results than clustering based only on the climatic variables of the scenarios. The differences will emphasize the value of monitoring species specific changes. Our methodology can be used to cluster based on individual species, on all species together, or on subsets of species such as those species most at risk. Clustering the predicted range maps for a single species across the suite of climate predictions is important to understand the climatic response for the single species, whereas clustering based on all species has the benefit of obtaining a broad understanding of the overarching patterns. Furthermore, clustering can be based only on a subset of the species such as those considered most at risk; summarizing predictions for species most at risk highlights climatic factors that drive high loss of biodiversity. We first describe prior modelling work to predict the future potential range maps given climate information (Sections 2.1 and 2.2). Then we discuss our methodology for clustering ecological and climate scenarios (Section 2.3); first based on predicted range maps, and secondly based on predicted climate variables.
Additionally, we provide an illustrating example. In Section 3, we discuss the scenario clusterings, with emphasis on the difference between the climate-based and ecological based clustering. We then demonstrate that we have detected meaningful clusters, with range maps within clusters similar to each other, but different than range maps in other clusters. Finally, we make some concluding remarks in Section 4. We used inhomogeneous Poisson point process models (PPPM) (Merow et al., 2017;Renner et al., 2015;Warton and Shepherd, 2010) to predict the future range maps of 1101 terrestrial mammal species in 2070 based on current day presence data and climatic covariates. For current day climatic conditions, five climate variables with minimum pairwise correlation were chosen from a set of 19 commonly available Bioclim variables in the world climate database, WorldCLIM2 (Fick and Hijmans, 2017). The mean values from 1970-2000 of the five chosen climate variables were used to fit a PPPM for the species occurrence. Occurrence data for each species were obtained from Miller (2020), and PPPMs were fit for species with at least 10 unique presence cells on a global 10km grid. The fitted model was used to predict future potential distributions based on the CMIP6 predictions of the five climate variables. Predicted relative occurrence rate was converted to binary range maps of presence/absence with a chosen threshold based on the 5% quantile of predicted relative occurrence rate values of training presences. This approach was used to make predictions for 1101 mammals with sufficient data using the 34 different sets of predicted climate variables. An example of predicted range maps for the Australian sandy inland mouse (Pseudomys hermannsburgensis) is shown in Figure 1b. Climate variables used for prediction are collected in the 6 (CMIP6). CMIP6 provides multiple climate predictions which vary by the underlying global climate model (GCM) and the representative concentration pathway (RCP) used. RCPs describe a predicted timeline of greenhouse gas concentrations (Eyring et al., 2016a). Four RCPs are included in CMIP6, which vary in the quantity of greenhouse gas emissions to capture the uncertainty of future emissions. We refer to the RCP trajectory of least emission as the "optimistic" trajectory, and refer to the most pessimistic trajectory as the "extreme" trajectory. See also Table 2 for more information about the different climate models. In this work we refer to a scenario as a (GCM, RCP) pair; a scenario represents uncertainty both in the evolution of climate, and in future greenhouse gas emissions. There is a variety of publications which studied the climate predictions considered in this paper; see Eyring et al. (2016b); Kageyama et al. (2018); Pascoe et al. (2020); Ma et al. (2020). First, we introduce some notation regarding the range maps (Section 2.3.1) followed by a way to quantify similarity between range maps (Section 2.3.2). Furthermore, we utilize spectral clustering to cluster scenarios based on the proposed similarity measure (Section 2.3.3). We conclude this section with a procedure to cluster scenarios based on different climate variables (Section 2.3.4), which can be compared to clustering based on the ecologically informed range maps which highlights the importance of utilizing ecological information. For notational simplicity we focus the presentation of the methodology for a single species. Furthermore, we suppose our range maps are on a r by c regular grid represented as (r, c) : r = 1, . . . , r; c = 1, . . . , c. We denote B s as the binary (presence/absence) predicted range map according to scenario s, to be more precise B s (r, c) = 1 if the cell at (r, c) represents a presence, and B s (r, c) = 0 otherwise. It is important to consider how these range maps differ from the present day; for this reason, we denote by P the present day map. Similarly as for B s (r, c), we define P (r, c) = 1 if the cell at (r, c) is presently occupied, and 0 otherwise. We define a candidate set A where presences and valid absences may occur, so that non-meaningful cells such as those in the ocean have no effect on quantifying similarity between range maps described in the next section. This binary candidate map A has size r by c, and could for example take the value 1 only when there is land, or alternatively the geographic area reachable from the species considering physical constraints on dispersal such as mountain ranges. We have chosen to take A as the union of the original range map and all scenarios, thus A(r, c) = 1 if P (r, c) = 1, or there is at least one scenario s such that B s (r, c) = 1. More generally, A is a mask for the r by c rectangular grid that is either occupied in the present or occupied in at least one future prediction. We choose to quantify the similarity (or differences) among the climate scenarios, s, with respect to changes in presence/absence from the present day, P . For each scenario s, we construct a corresponding weighted map, W s , that indicates the nature of the changes. The value of W s (r, c) is 0 if A(r, c) = 0, but when A(r, c) = 1, the value of W s (r, c) depends on the intersection of the value for the present day map and the value under scenario s at location (r, c). The four possible values of W s (r, c) are shown in The four possible cases represent unchanged (kept) presences, p keep , new presences, p new , new absences, a new , and unchanged absences, a keep , respectively, where "keep" and "new" are with respect to the present day range map, P . For example a grid cell (r, c) corresponding to p keep is a grid cell that is a presence in the current day range map, P (r, c) = 1, and remains a presence according to scenario s, B s (r, c) = 1.
The best choice of weightings, (p keep , p new , a new , a keep ), is dependent on the pairwise similarity metric used for clustering. We choose to compute the pairwise similarity between scenario range maps as the cosine similarity of their weighted range maps, CosineSimilarity(W s , W s ) or CS(W s , W s ) for short. We define the cosine similarity between two matrices W s , W s of dimension r x c as The product in the numerator in (1) suggests that when W s (r, c) and W s (r, c) are zero there is no contribution to the measure of similarity suggesting that the weighting scheme, (p keep , p new , a new , a keep ), should be constructed such that meaningful contributions to similarity have non-zero values. In our weighting scheme, presences are given positive weights, absences negative weights, and we choose to emphasize those cells whose ecological suitability for this species is vanishing. These "lost" cells corresponding to novel absences, {(r, c) : W s (r, c) = a new }, are particularly important; they represent cells where according to the range map model, the geographic location is no longer suitable given the change in climate. In addition, these cells should be weighted higher to emphasize their importance for conserving a species as they are known to have been occupied at some time point. By contrast the cells corresponding to p new represent regions that are predicted to become suitable, however such predictions do not imply that the locations will be occupied due to uncertainty in dispersal, which is not directly taken into account by the PPPM. Thus, we also make the choice to represent the fact we are more concerned about all cells in the region of present day occurrences. However, we emphasize the flexibility of our method as alternative choices can be made depending on the goals. A visualization of this weighting scheme is shown in Alternative weightings can be considered for various applications, for instance the choice of does not make use of the present day map, and does not differentiate cells that underwent a change from the present day map to the predicted scenarios. Alternative similarity metrics could have been considered. However, the cosine-similarity has certain advantages over other commonly used measures to quantify the similarities (or distances) between range maps. The Hellinger distance and Kullback-Leibler divergence (Hagen, 2002) rely on a probabilistic interpretation, and thus it is difficult to incorporate absences into these measures, however they might be useful when considering differences in range maps describing occurrence probability. The Kappa statistic (Hagen, 2002) can be used to measure similarities between categorical maps, however it is not as clear how to weight important cells, such as novel absences. In addition, the Kappa statistic does not directly incorporate the relative frequencies of absences and presences. For instance if two range maps both predict all presences except a single absence cell, these two range maps will have a negative Kappa statistic if their one absence cell differs in location. However, for our purposes we would like to consider such a pair of range maps to be very similar. The Wasserstein distance (Peyre et al., 2019) measures the distance a species moves, suggesting it is appropriate in our setting. However, the Wasserstein distance does not model the disappearance of regions, and attempts at computing the Wasserstein distance proved too slow or even non-convergent in all but the range maps with the smallest number of presence cells. The computation and resultant cosine similarity using these weighted range maps are readily interpretable; when two scenarios agree on the presence or absence of a cell, this cell has a positive contribution to the cosine similarity, whereas the cell has a negative contribution when the two scenarios disagree. For maps are identical and 1 if the two range maps completely disagree on presences. This interpretability allows for a simple method to combine information across species. Furthermore, cosine similarity can also be implemented very efficiently and quickly utilizing matrix multiplication. For each species m = 1, . . . , 1101, we compute the pairwise scenario similarity matrix by computing the cosine similarity (1) on each pair of scenario weighted range maps, (W m s , W m s ), where W m s is the weighted range map for species m under scenario s. That is for each mammal species, we construct the s by s matrix , where s = 34 is the number of scenarios. We cluster scenarios by spectral clustering on the cosine similarity measure (Von Luxburg, 2007). The properties of spectral clustering are understood from a graph theory perspective (Von Luxburg, 2007). The similarity matrix S m can be thought of as an undirected graph whose nodes are the scenarios and the edge weight between a pair of scenarios s and s is given by S m (s, s ). Spectral clustering has best performance on sparse graphs, thus for the dense similarity matrix S m , the first step is to sparsify it by taking the k-nearest neighbor graph, that is retaining an edge from s to s only if s is within the top k neighbors of s (i.e., it is within the top k scenarios of maximal similarity to s). When the graph is expressed as a matrix, sparsifying the graph corresponds to setting entries of the matrix to zero. Retaining the top k neighbors leads to a directed graph as this definition of nearest neighbor is not symmetric. Thus we retain the undirected edge from s to s if either s is within the top k neighbors of s, or vice-versa (Figures 2 A to B). The retained edges are weighted by the similarity of their endpoints. We denote by E this matrix of retained weights.
The main computational tool of spectral clustering is the graph Laplacian L = D  E, where D is a diagonal matrix of node degree, D(s, s) = ns s =1 E(s, s ) (i.e., the number of scenarios connected to scenario s). The graph Laplacian can be thought of as a matrix representation of a graph, with useful mathematical properties as will be discussed shortly. We used the random-walk normalized graph Laplacian, L rw = D 1 L, as suggested in Von Luxburg (2007), although this normalization choice is most significant when the node degrees vary significantly, which is not the case here. The normalization is used to prevent a single node from dominating the resulting spectral clustering when that node has many more connections, i.e., it is a "central node". However, such "central nodes" were not found in this application since each scenario is similar to only a few other scenarios (i.e. those with similar RCP and or GCM). If the graph corresponding to the normalized Laplacian is disconnected (i.e. there is at least one group of nodes that do not have any edges outside that group, Figure 2a.C), then there is at least one eigenvector of L rw which is sparse with non-zero entries only on nodes within that group, with an associated eigenvalue of 0. If instead the graph is "noisier" and fully connected as in our case, the first eigenvectors (corresponding to the smallest eigenvalues) are only approximately sparse; they have large magnitude on groups of nodes that are tightly connected, and small magnitude on nodes outside such groups. Thus these eigenvectors can be used to identify clusters, which is discussed next. The spectral embedding of the s by s matrix L rw is related to the eigendecomposition of L rw . The eigenvectors of L rw arranged by decreasing order of the corresponding eigenvalues form the columns of a s by s matrix U . The rows of U represent the coordinates of the corresponding observational unit (in our case observational units are scenarios). This embedding (transformation from a graph, L rw , to real valued coordinates, U ) is meaningful because of the properties of spectral clustering disucssed above. Illustrations of spectral clustering on both a 'clean' (perfectly separable clusters) and 'noisy' (one cross edge) examples are shown in Figures 2a and 2b. The images illustrate the steps of the spectral clustering described above. The clean example is 'cleanly' disconnected: the first four nodes form a cluster which share no edges to the last four nodes (i.e., the symmetrized kNN matrix, Figure 2a.B, is a block matrix). The spectral embedding process starting on the clean example can be summarized as follows: it starts with an input similarity matrix (Subfigure A)), where each row (and column) of the matrix represents an observation, and entry (i,j) of the matrix represents the similarity between observations i and j. The next step is to compute the nearest neighbor similarity matrix, which is shown in Subfigure B). Here we have chosen k=3 (the number of nearest neighbors is 3), so we set an entry in the nearest neighbor matrix to 1 if it is within the top 3 values in its row or column, else we set it to 0. The next step is to compute the graph Laplacian L = D E of the nearest neighbor graph, should in Subfigure C). Here L is the graph Laplacian, D is a diagonal matrix of the node degrees (in this case every node has degree 3) and E is the nearest neighbor graph. This is then normalized to random walk flavor of the graph Laplacian, L rw = D 1 L. The final step is to compute the spectral embedding from the graph Laplacian, which is done by taking the smallest two eigenvectors of the random walk graph Laplacian, shown in Figure D). The two columns serves as the two coordinates: the first column is the x-coordinate in subplot D), the second column the y-coordinate. We see that the spectral embedding in Figure D) faithfully represents the original graph, there are still two clusters of four nodes each, but now represented in a two dimensional plane rather than the original similarity matrix. The noisy example has a single cross edge entry of the matrix (i.e., there is one pair of scenarios from the original separate groups that are now considered neighbors, Figure 2b.B) that affects the final spectral embedding; the first and fifth point are closer together . Our case of 34 scenarios is more complex than this illustrated example of only 8 points, however the pattern that groups of scenarios which are all similar to each other will be near each other in the embedded space remains true. Once observations (in our case future climate scenarios) are transformed to real coordinates via spectral embedding as discussed above, a simple clustering algorithm can be applied such as k-means. However we chose the single-linkage clustering algorithm over k-means in this embedded space as it performs well according to the Davies-Bouldin criterion. The Davies-Bouldin criterion is a common clustering criterion quantifying the goodness of separation between clusters (Davies and Bouldin, 1979). The Davies-Bouldin criterion is interpreted as a ratio of the average intra-cluster distance and inter-cluster distance; a good clustering will have points that are close together within cluster (small intra-cluster distance), and points that are far apart in different clusters (large inter-cluster distance). Spectral clustering can be performed for a single species using only S m . Alternatively, one way to combine information across species is to capture group wide trends in scenarios. We cluster similarity matrices based on two aggregations. First, we average the similarity matrices across all m = 1101 species, as in S = (1/m) m m=1 S m . Second, we average over a subset of species whose fraction of area lost is among the highest 10%.
We average the similarity matrices across species in order to incorporate information from each species into a quantitative clustering approach, obtaining a broader understanding of the overarching (a) : the two groups are separable in that there are no shared members. These four plots show the steps of spectral embedding: Subfigure A) is the example input similarity matrix, where each row (or column) represents an observation. Subfigure B) shows the nearest neighbor step with the choice of k = 3, Subfigure C) shows the graph Laplacian (the random walk version of the Laplacian) of the nearest neighbor graph, and D) shows the final spectral embedding. We see that in this ideal example the spectral embedding perfectly separates the two groups of columns (the first four and second four observations are in separate clusters). The results on this clean example can be compared to a noisier example in Figure 2b. patterns between species range maps and climate scenarios. We contrast the clustering based on range maps with a clustering based only on predicted climate variables. This comparison will illustrate the importance of incorporating ecological information by contrasting clusters based only on climate to those using climate to first estimate species ranges before clustering. We discuss the process of clustering using only the five climate variables that were used to predict the range maps (see Section 2 for a description of the predictions). Each of these five globally distributed variables is predicted across the 34 scenarios. In order to directly compare climate-based clustering to our ecologically based clustering, we utilized the climate variables to perform a clustering in a similar fashion to the ecologically based clustering. However, for continuous data, the cosine similarity is not appropriate, as two maps that are shifted versions of each other would be considered very similar to each other. For instance, if one scenario predicted two degrees Celsius warmer everywhere than another scenario, the cosine similarity between these two maps would be very high, which is not desirable as these two maps represent significantly different predictions. Instead, we use the L 2 distance between maps, which will effectively use both the difference between the means of maps and differences in the spatial variation. To incorporate all five climate variables, each variable was normalized before applying the L 2 distance. We denote T f s (r, c) as the value of variable f according to scenario s at location (r, c). The variable scaled L 2 distance between a pair of scenarios s and s using all five variables is given by where  f is the standard deviation of variable f measured across all locations and scenarios, that is we calculate the mean of each variable across all locations and scenarios, We create a similarity matrix, S, from the distances by using the monotonically negative transformation S(s, s ) := 1/H(s, s ). Spectral clustering can be performed on this climate-informed matrix S which can be compared to the ecologically driven spectral clustering. In order to quantitatively compare cluster assignments, we used the adjusted Rand index (Hubert and Arabie, 1985). The adjusted Rand index measures how similar two cluster assignments are to each other. We make use of the adjusted Rand index by measuring the similarity of the clustering results to partition by either the representative concentration pathway (RCP) or global climate model (GCM), in order to quantitatively answer whether RCP or GCM differences is the driving factor behind the clustering results. Given two cluster assignments X = {X 1 , . . . , X n } and Y = {Y 1 , . . . , Y m }, where each X i and Y j are clusters (sets of scenarios). Define the cardinality of the set of overlaps between two clusters as The adjusted Rand index is given by where n 2 denotes the binomial coefficient and is calculated as n 2 = n(n1) . The adjusted Rand index measures how similar two cluster assignments are to each other compared to random assignments. Random assignments have an expected adjusted rand index of 0. 3 Results and Discussion One can get an overall sense of the changes predicted in the range maps from Figure 3, which shows how the diversity of mammals is spread spatially over the earth. We see from both Figures 3c and 3d that significant potential diversity losses are predicted around the equator in South America and Africa, whereas there is some potential diversity increases further north, consistent with previous findings (Chen et al., 2011). We note that these 'losses' do not necessarily imply a local extinction but rather that a species is exposed to climate beyond its current realized niche. Similarly species 'gains' do not necessarily imply that new species will occur in a given cell, but rather that the cell is newly suitable for a species climatically. A species with newly suitable locations may still not occur there due to dispersal limitation, biotic interactions, etc. By comparing the fraction of cells corresponding to a new and p new , we get a sense of potential range map shifts compared to present day range maps. This is shown in Figure 4, which shows potential range map shifts, as opposed to range map expansions or contractions; the number of new presences and new absences 2011) which finds that species are moving poleward and towards higher elevations, there is loss around the equator and some increase in diversity towards the northern pole. grow together, which would occur as range map shifts, instead of say absences growing as presences shrink, which would be indicative of an overall range size decrease. However, although novel absences and novel presences tend to occur together as shown in Figure 4, on average there are more novel absences (average 38% +/- 28% s.d. of current range map) than novel presences (average 25% +/- 26% s.d. ), which implies that the species ranges are both shifting and decreasing in size (average 13% net loss, +/- 31% s.d.)
due to the changing climate. The spectral clustering results using the species-specific similarity matrix, S m , is shown in Figure 5 for four species, illustrating the main types of patterns observed across all 1101 species. We found an interesting mixture of RCP and GCM dependence. RCP is a major driving factor of cluster composition; in most of the species-specific clusterings, the far left ("optimistic") cluster contains mainly scenarios with low RCP, and the far right ("extreme") cluster only scenarios with high RCP (Figures 5a-5c). However, the clustering for the slender treeshew (Tupaia gracilis; Figure 5d) predicted ranges is driven mainly by GCM. This clustering mainly by GCM was found in many species (30% species' cluster results had a higher Rand index with a GCM clustering than an RCP clustering), probably because these species niches had relatively weak dependence on mean annual temperature compared to other climate variabels. However, the most common trend is RCP dependence (70% of species). This variability among species is further evidence supporting the importance of the climate-ecology relationship; the most important difference between climate models (GCM or RCP) varies depending on the individual species. The spectral cluster results from using the similarity matrix averaged over all 1101 mammals considered, S, is shown in Figure 6a. We see in the bottom cluster of Figure 6a a mix of RCP and GCM. Similar to the individual species results, this suggests that RCP alone does not account for the variation, and the GCM is 5c), but a starkly different GCM-driven clustering for the treeshew . This variety was found among species, most species clusterings are strongly driven by RCP, but some are more reflective of GCM. These differences between species further demonstrates that ecological information is important to interpret the differences between scenarios. also important. We also performed clustering for only the species most at risk, defined by those species whose fraction of area lost is among the highest 10%. This loss in area can be used to approximate a loss in population abundance using the techniques in He (2012) and Che-Castaldo and Neel (2016). Performing spectral clustering by using the average of similarity matrices of this subset of species at risk is shown in This clustering puts a strong emphasis on RCP, that is, the species most at risk are more sensitive to RCP. One way to qualitatively measure the performance of the clustering is to look at the spatial overlap of For example, the overlap of all 34 scenario predictions for the Australian mouse is shown in Figure 8a. We see that the cluster associated with the lowest RCP scenarios ("optimistic" scenarios) accounts for most of the presences in the discrepant regions, whereas the cluster associated with the highest RCP ("extreme" scenarios) accounts for many of the absences . This demonstrates that we have discovered meaningful clusters; there is agreement within clusters but disagreements across clusters. We see that the climate-based clustering is meaningfully different than the ecological one, and the clustering is driven mainly by GCM instead of by RCP in the ecological based clustering. The climate-based clustering and ecological one have an adjusted Rand index of only 0.26, demonstrating their dissimilarity. This is evidence for the importance of considering the specific climate niche occupied by a species in relation to how those conditions are projected to change. Although these same five variables are used to predict the species' ranges, these predicted ranges paint a different picture of the scenario clustering because of how the species are influenced by the climate variables. In fact, we can get a sense of variable importance by clustering based on individual climate variables as opposed to their weighted average, shown in An interesting direction to explore further is to uncover subsets of species that respond differently than others. For instance it may be the case that rodents tend to fare worse under a specific climate model, compared to other mammals. A similar area of future research is to determine why some species like the slender treeshew cluster more by GCM instead of the more common pattern of RCP. Another extension is to consider how to combine information across species in a more holistic manner as opposed to the simple average. For instance, Dong et al. (2013) presents a methodology to cluster according to many graphs, which could be applied to the set of scenario graphs from each species. Financial support is gratefully acknowledged from a Xerox PARC 1455172, 1934985, 1940124, and 1940276, USAID, and for a Sustainable Future. CM acknowledges support from NSF Award HDR-1394790. Financial support is gratefully acknowledged from a Xerox PARC 1455172, 1934985, 1940124, and 1940276, USAID, and for a Sustainable Future. CM acknowledges support from NSF Award HDR-1394790. We provide here a table with information about the climate models considered in this work We provide here a table with information about the climate models considered in this work