changeRangeR: An R package for reproducible biodiversity change metrics from species distribution estimates

Abstract

Conservation planning and decision-making rely on evaluations of biodiversity status and threats that are based upon species' distribution estimates. However, gaps exist regarding automated tools to delineate species' current ranges from distribution estimates and use those estimates to calculate both species- and community-level biodiversity metrics. Here, we introduce changeRangeR, an R package that facilitates workflows to reproducibly transform estimates of species' distributions into metrics relevant for conservation. For example, by combining predictions from species distribution models (SDMs) with other maps of environmental data (e.g., suitable forest cover), researchers can characterize the proportion of a species' range that is under protection, metrics used under the IUCN Criteria A and B guidelines (Area of Occupancy and Extent of Occurrence), and other more general metrics such as taxonomic and phylogenetic diversity and endemism. Further, changeRangeR facilitates temporal comparisons among biodiversity metrics to inform efforts toward complementarity and consideration of future scenarios in conservation decisions. changeRangeR also provides tools to determine the effects of modeling decisions through sensitivity tests. Transparent and repeatable workflows for calculating biodiversity change metrics from SDMs such as those provided by changeRangeR are essential to inform conservation decision-making efforts and represent key extensions for SDM methodology and associated metadata documentation.

INTRODUCTION

Reliable estimates of species' current ranges, as well as their recent and forecasted future trends, are essential to inform a wide range of biogeographical studies, and especially applications for conserving biodiversity (Whittaker et al., 2005) and related decision-making (Araújo et al., 2019). Species' geographic distributions are a well-defined Essential Biodiversity Variable (EBV; GEO, 2015; Pereira et al., 2013), which can be quantified through species distribution models (SDMs) and post-processed to better refine estimates of a species' current range (Merow et al., 2022). Species distribution estimates can potentially inform conservation status compilations, such as the Red List of Threatened Species, following assessment processes defined by the International Union for the Conservation of Nature (IUCN, 2022). However, the ability to transform a species' current estimated distribution into standardized conservation metrics in a transparent and repeatable way remains an outstanding issue in conservation biology (Anderson, 2022; IUCN, 2022; Kass, Meenan, et al., 2021).

In combination with expert information regarding range extents (Mainali et al., 2020), modeled species' range estimates (e.g., SDMs) provide several advantages for measuring biodiversity over expert-range maps, which are generally static and unable to be reproduced or transferred to other time periods or geographic extents. SDMs with high-quality data and following best practices (Araújo et al., 2019) generate distribution estimates that can provide fine-scale spatial information to support conservation planning (Di Febbraro et al., 2018). Further, SDMs can be transferred not only to unsampled areas but also to new time slices, allowing for calculations of trends in species' ranges (Veloz et al., 2012). SDM range predictions can be refined with other data (e.g., forest cover; Gavrutenko et al., 2021) or land conservation maps (e.g., protected or at-risk areas) to estimate recent fragmentation or habitat loss and determine estimates of the species' current range. Specifically, IUCN's Criteria A and B (2022) help advise practitioners of threats based on geographic ranges and their change through time (see Box 1). They employ two alternative metrics to measure a species' geographic range, extent of occurrence (EOO), and area of occupancy (AOO), which conservation practitioners can use to inform threat assessments, identify areas for priority assessment, and guide field surveys. The newly proposed area of habitat (AOH) adds to these by refining IUCN range maps using ancillary land-cover and remote-sensing data based on species' known habitat and elevational limits (Brooks et al., 2019).

1 BOXSummary of terminology pertaining to IUCN species extinction risk assessment criteria A and B as defined in Guidelines for Using the IUCN Red List Categories and Criteria (2022)

IUCN Criterion A pertains to population size reduction as measured through (a) direct observations, (b) taxon- appropriate indices of abundance, (c) declines in EOO, AOO (defined below), or habitat quality, (d) exploitation levels, or (e) the negative impacts of biotic interactors.

IUCN Criterion B pertains to geographic ranges, with the aim of identifying species with restricted distributions that are at least any two of (a) severely fragmented, (b) known from few locations, or (c) undergoing continuing decline or showing extreme fluctuations. Species' ranges are measured using either EOO or AOO.

Extent of occurrence (EOO) represents the geographic area currently occupied by a taxon. It is measured by the shortest continuous perimeter encompassing all of the “known, inferred, or projected sites of present occurrence of a taxon, excluding cases of vagrancy” (IUCN, 2022).

Area of occupancy (AOO) is a more conservative subset of EOO that identifies areas of suitable habitat currently occupied by a taxon. It is measured as the summed area of 2 km×2 km pixels that capture all of the areas where the species is known or very likely to occur.

In addition to those from single species, range estimates from multiple species can be combined to calculate assemblage-level macroecological and conservation metrics. For example, overlaying species' ranges can be used to estimate maps of potential species richness and endemism (Calabrese et al., 2014; Zurell, Zimmermann, et al., 2020). The benefits and disadvantages of using binary range estimates (whether expert-drawn or post-processed SDM-derived) for macroecological applications have been well documented (Graham & Hijmans, 2006). Several additional metrics can be assessed when combined with phylogenetic data: phylogenetic diversity can elucidate the spatial distribution of evolutionary history (Kembel et al., 2010), while phylogenetic endemism determines the relative spatial restriction of evolutionary history in an area (Rosauer et al., 2009).

Several R packages already exist for calculating biodiversity metrics such as IUCN's range metrics EOO and AOO (e.g., ConR; Dauby et al., 2017, RED; Cardoso, 2017, speciesgeocodeR; Töpel et al., 2016; rCAT; Moat & Bachman, 2020; See Table 1), but their workflows operate only on occurrence datasets, precluding estimates of these metrics based on SDMs. Additionally, existing tools lack functionalities to harness remotely sensed data (e.g., satellite imagery) that provide high-resolution information that can be included to improve further insights—including trends over time—for Red Listing and other biodiversity measurements (Duncan et al., 2021). We built the package changeRangeR to address the need for reproducible quantification of biodiversity metrics from both occurrence datasets and species' range estimates, to enable related analyses at both the species and community levels, and to include easy integration of remotely sensed information within a single workflow. Additionally, changeRangeR helps researchers report on methodological uncertainty and temporal trends in range estimates. It can also incorporate relevant multispecies metrics.

TABLE 1 Comparison of changeRangeR's features with those of other software packages for calculating biodiversity metrics such as IUCN's range metrics EOO and AOO (adapted from Chang, 2022)

	Software tools
changeRangeR's features	speciesgeocodeR	rCAT	ConR	Red
Range proportion
Calculates the proportion of species range binary map overlap with categorical features (e.g., shapefiles), or any user-defined georeferenced categorical spatial vector layer.
Estimates current species- and community-level metrics by integrating SDMs and selective masking	[IMAGE OMITTED. SEE PDF.]
Range correlation
Assesses the proportion of overlap between a species' range map (binary values) and a continuous environmental layer broken into quantiles, or user-specified values (i.e., values between 0 and 1)	[IMAGE OMITTED. SEE PDF.]
Returns the average Pearson correlation coefficient between the species' range and a continuous layer.
Multispecies estimates
Calculations relating to trends in species' and communities' ranges over time			[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]
Compares ordinal raster values (e.g., species richness) with regions of interest and returns proportions of the raster that fall within and outside of these regions	[IMAGE OMITTED. SEE PDF.]		[IMAGE OMITTED. SEE PDF.]
Calculates phylogenetic endemism by dividing each branch length of a phylogenetic tree by the total range area of the species that descend from that branch, then summing
Extent of occurrence (EOO) calculation
MCP around the occurrence localities	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]
MCP around modeled continuous range estimates thresholded using a range of user-defined values
Adheres to model metadata and reporting standards through logging of modeling decisions and outputs		[IMAGE OMITTED. SEE PDF.]
Area of occupancy (AOO)
Calculates area of occupied cells (2 km width/height) based on species occurrence data		[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]
Calculates area of occupied cells (2 km width/height) based on binary grid of species' range				[IMAGE OMITTED. SEE PDF.]
Other Biodiversity Evaluations
Calculates metrics relevant to conservation evaluations by assessing species' geographic ranges' proximity and overlap with geographic features such as land cover and sensitivity to climatic changes				[IMAGE OMITTED. SEE PDF.]
Provides statistics relevant for monitoring biodiversity at the community level	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]

Specifically, changeRangeR allows users to report current species- and community-level metrics estimates by integrating species' range estimates (rasters) and selective masking. Additionally, changeRangeR provides a workflow for users to calculate several aspects of trends over time in metrics related to species' geographic distributions, for example, the change in the area of a species' distribution or the change in percentage of forest within a species' range over time. Users can also make comparisons with calculations made from occurrence localities only. The package also allows users to assess the sensitivity of range estimates to differences in methodological choices (e.g., the threshold used to make a continuous SDM prediction into a binary range estimate). changeRangeR is designed to be an independent R package and provides a workflow to perform comprehensive spatial and temporal range map refinements and evaluations. changeRangeR can receive as input binary maps that estimate species' realized current distributions, such as continuous suitability predictions, AOH maps, expert-drawn maps, or post-processed SDM outputs from the R package maskRangeR (Merow et al., 2022). This workflow can facilitate the measurement of biodiversity conservation metrics made using range maps and enables both single and multiple species analyses. The following sections outline how changeRangeR fills gaps for assessing range characterizations that build upon range maps for quantifying species- and community-level biodiversity metrics associated with the geographic range.

PACKAGE DESCRIPTION

changeRangeR meets the needs of researchers and conservation practitioners in four main ways. First, changeRangeR allows users to calculate metrics relevant to conservation evaluations by assessing the proximity and overlap of species' ranges with geographic features such as land cover or sensitivity to climatic change. Second, using any form of occurrence estimate, users can help assess species' risk levels using IUCN's criterion B guidelines regarding geographic ranges. Third, changeRangeR provides calculations relating to temporal trends in species ranges and communities for proactive planning that can also be useful information for IUCN's criterion A. Finally, changeRangeR helps users adhere to model metadata and reporting standards through recording of modeling decisions and outputs.

Range protections and threats

As species' ranges shift in response to environmental changes, the importance of geographic features such as parks or deforested areas will increase, resulting in strong landscape implications for species' ranges (Krauss et al., 2014). The ratioOverlap function calculates the geographic overlap between a species' range and geographic features of interest.

Range proportion

Using ratioOverlap, users can calculate the proportion of overlap between a species' range with categorical features (e.g., shapefiles) such as land-cover classes, habitat types, ecoregions, political borders (country, state/province, etc.), protected areas (e.g., national parks), or any user-defined and georeferenced categorical spatial vector layer. For example, users can calculate the proportion of a species' total range that falls within different types of protected areas, such as national parks, national wildlife refuges, and local parks. The species' range map should be any binary range estimate (e.g., expert-refined model, AOO, binary SDM, etc.).

Range correlation

Users may also use ratioOverlap to assess the proportion of overlap between a species' range map (binary values) and a continuous environmental layer separated into quantiles. By default, the function will determine the proportion of overlap that is within the first (<25%), second (25%–50%), third (51%–75%), and fourth (>75%) quartiles, respectively. The option for user-specified values (i.e., values between 0 and 1) is also available. In addition, the function will return the average Pearson correlation coefficient between the species' range and the continuous layer. Examples of environmental layers for calculating overlap with species' range maps are the Humboldt's Human Footprint Indices (Correa Ayram et al., 2018) and hunting pressure level datasets (Deith & Brodie, 2020).

Inputs to assess species' threatened status

changeRangeR also provides metrics needed to assess species' threatened status using IUCN's criterion B guidelines regarding geographic range. Specifically, the package can calculate estimates of EOO, AOO, and sensitivity to threshold decisions. Recommendations of the relevant IUCN Species Survival Commission Red List Authority should be consulted when contributing EOO or AOO measurements (IUCN, 2022).

EOO–IUCN's extent of occurrence

IUCN's EOO calculation is defined as the area contained within the shortest imaginary, continuous boundary drawn to encompass all the regions currently occupied by a taxon, excluding vagrant localities (e.g., those that do not represent the species' core habitat area; Section 4.9, IUCN, 2022). This measure is intended to include discontinuities or disjunctions within the overall distribution of a taxon (e.g., large areas of unsuitable conditions, in contrast, see AOO, below, IUCN, 2022). EOO is typically measured by drawing a minimum convex polygon (MCP, also called convex hull) around current localities and may exclude large areas that are actually occupied by the species but remain unsampled (Joppa et al., 2016; Figure 1). Alternatively, a convex hull around a thresholded SDM prediction that predicts suitability for unsampled areas can be a better approach for data-poor species (Kass, Meenan, et al., 2021; Syfert et al., 2014). Such an estimate should be processed to exclude areas beyond likely dispersal barriers and those the species is unlikely to occupy due to unfavorable biotic conditions, and it should still include all species' presence localities in keeping with the definition of EOO (IUCN, 2022).

FIGURE 1. Bassaricyon neblina is a small, arboreal, carnivoran mammal found in the Andean cloud forests of Ecuador and Colombia. The map shows the extent of occurrence (EOO) resulting from a minimum convex polygon surrounding all occurrence records (red polygon) and the EOO resulting from a minimum convex polygon surrounding the binary SDM prediction (black polygon) following Gerstner et al. (2018). The upper-bound estimate of the area of occupancy (AOO; light brown areas) is defined as the area of this binary SDM prediction, resampled to the 2 km × 2 km resolution suitable for IUCN Red List assessments. The inset map shows an example lower-bound AOO represented as 2 km × 2 km red pixels calculated from occurrence localities.

changeRangeR allows users to calculate EOO via two options. The first option uses the mcp function to define the EOO by drawing an MCP around the occurrence localities, while the second option, mcpSDM, emulates the methodology described by Syfert et al. (2014): it draws an MCP around modeled continuous range estimates that are thresholded using a range of user-defined values, then finds which of these MCPs is most similar to the original MCP delineated around the occurrence localities, resulting in an EOO estimate that takes into account predictions for unsampled areas. An important caveat is that thresholds greater than the minimum suitability across all occurrence localities will result in some occurrences being excluded from the EOO, which does not match the definition from IUCN (IUCN, 2022).

AOO–IUCN's area of occupancy upper bound

Users can calculate the total area of 2 km×2 km grid cells to determine bounds of IUCN's area of occupancy (AOO). AOO is intended to represent only areas that are likely to be presently occupied, inferred through habitat information or predicted via modeling, and thus is meant to exclude unsuitable or unoccupied habitat patches (IUCN, 2022). To ensure consistency and comparability of results in IUCN assessments, for most applications AOO should be calculated at a reference scale with grids of 2 km×2 km (a cell area of 4 km²). In changeRangeR, users can calculate AOO one of two ways: using occurrence localities as a lower-bound estimate, or alternatively using a binary range estimate derived from a distribution model or an area of habitat map as an upper-bound estimate—the latter should ideally be processed to exclude areas that do not currently offer habitat for the species (Figure 1). Lower-bound estimates of a species' range show the minimum area that the species most likely inhabits, whereas upper-bound estimates are more likely to include more of the current range but may overpredict it. For accurate areal calculations appropriate for IUCN measurements, users need to use grids that are defined in a projected coordinate system, ideally one that uses an equal area projection, which best preserves area relationships and has units in meters instead of degrees (IUCN, 2021).

The function AOOarea calculates AOO based on either grid cells containing occurrence records (lower bound) or a binary range estimate (upper bound). The lowest estimate of AOO is defined as the total area of 2 km × 2 km grid cells that contain occurrence records. Lower-bound estimates like this are strongly affected by sampling biases that can result in the omission of large, unsampled areas that may be part of the species' current range (Anderson, 2022). An alternative is to estimate a higher bound of AOO based on a gridded binary range estimate from an SDM. Upper bounds for AOO can include unsampled grid cells predicted as suitable by the SDM and are therefore more robust to sampling bias, though they also have more uncertainty as they include model predictions. If the function AOOarea is supplied with a different spatial resolution, it will resample the grid to a 2 km × 2 km resolution before calculating the area (the input raster must be projected in meters). Users are advised to assess the sensitivity of several thresholds used to create the binary map by calculating areas based on different thresholds (also see IUCN, 2022, which describes an approach to calculate AOO first at the highest resolution possible and then upscale successively to determine the scaling correction). Examples include different suitability thresholds for an SDM prediction or different levels of current forest cover (Figure 2). As estimates of AOO, by definition, must include all grid cells that contain occurrence records that still represent extant populations (IUCN, 2022), some higher thresholds that exclude such grid cells would be inappropriate for calculating AOO (Kass, Meenan, et al., 2021).

FIGURE 2. Area of occupancy (AOO) can be calculated from an SDM prediction. The sensitivity of thresholds used to create a binary modeled range estimate from a continuous one can be assessed. AOO calculated from these SDMs is affected by the threshold used (left panel; darker colors representing a more constrained distribution as the forest threshold increases) and the area becomes reduced as the threshold increases (center panel). The calculated area is also affected by the year of the forest cover data used to apply the forest threshold to the SDM, allowing a comparison of area change over time (right panel).

Range trends through time

Given the appropriate temporal data, changeRangeR can leverage EOO and AOO predictions through time to calculate IUCN Red List criteria pertaining to population and range size reduction. Predictions of range changes over time have varying levels of associated uncertainty, yet they represent a data-driven methodology for forecasting a species' conservation status. Climatic and land-cover information from different time periods can be combined with range information to determine trends in several biodiversity metrics including continuing decline, extreme fluctuations, and proximity to important landscape features (e.g., protected areas). Calculating these trends allows researchers to predict how species may respond to current and future environmental changes.

Users can also calculate the overlap between SDMs and projected environmental variables through time. The function futureOverlap assesses the overlap of future binary range projections from SDMs with future environmental vector layers for the same time periods (e.g., annual forest cover). For example, a user can project an SDM to several future decades (2021–2040, 2041–2060 from the CMIP6 climate projections) and compare the overlap with projected deforestation available for those decades (for example, through 2050 using Clark Lab's Land Change Model; Clark Labs, 2009).

The function envChange determines how the area of a binary range estimate for the present can change due to the dynamics of an environmental variable that changes through time. This is especially useful for short-term estimates where climate change projections are either unavailable or not expected to shift the range within the time period. For example, monitoring how a current species' range is affected by changes to forest cover over time could help to identify trends and zones for implementing regional conservation actions (Figure 2). With this function, users supply a time series of relevant environmental rasters, an appropriate suitability threshold, and a binary range estimate in a matching coordinate reference system projection. Users can then decide how the range estimate is thresholded, with options for upper (e.g., maximum forest cover), lower (e.g., minimum human footprint), neither (e.g., environmental rasters are already binary), or both (e.g., two values that define a temperature range). Users can also determine the relationship between the range estimate and environmental variables with the Pearson correlation coefficient.

Community statistics

Community-level metrics can aid conservation planning and decision-making by providing information on species' representativeness and complementarity across conservation areas (e.g. Astudillo-Scalia & Albuquerque, 2020). These indicators can reflect different dimensions of biological diversity such as species richness, spatial and temporal turnover, endemism, and phylogenetic diversity (e.g. see Guisan & Rahbek, 2011; Varzinczak et al., 2020). The complementarity function specifically compares ordinal raster values (e.g., species richness) with regions of interest and returns proportions of the raster that fall within and outside of these regions.

Areas with many endemic species represent unique combinations of phylogeny and geography (Varzinczak et al., 2020) and are important targets for conservation efforts. Such areas may be left out of conservation priorities because of a lack of data or methods available to calculate information about them (Rosauer et al., 2009). We provide the function speciesEndemism for calculating an estimate of species endemism using a rasterStack of binary range estimates. This function calculates endemism as the number of species in a given area divided by the total area in which those species are found (this is also referred to as “range rarity”; Roberts et al., 2002). Phylogenetic endemism, computed based on species' phylogenies and spatial range data (Rosauer et al., 2009), is a measure of areas with phylogenetically distinct and geographically range-restricted species that can help identify potentially important conservation areas. The calc_PE function calculates phylogenetic endemism by dividing each branch length of a phylogenetic tree by the total range area of the species that descend from that branch, then summing.

Modeling metadata and reporting standards

Accurate recording and reporting of analytical decision-making are of critical importance, especially for science-backed policy. The R package rangeModelMetadata (rmms; Merow et al., 2019) defines an rmm R object representing metadata that can be populated with modeling decisions and results related to SDMs. This facilitates reproducibility of outputs and enables an easily accessible assessment of model quality. To further increase transparency, we made a custom function to document analytical steps in changeRangeR and the resulting outputs of post-processed SDMs. The buildRMM function catalogs inputs and results of other changeRangeR functions, adding this information to either a pre-existing rmm object, or creating a new one if not supplied.

PACKAGE DEVELOPMENT AND TESTING

changeRangeR was developed in partnership between the American Museum of Natural History's Center for Biodiversity and Conservation and user-consultation workshops by The Colombia Biodiversity Observation Network (BON), hosted by the Alexander von Humboldt Biological Resources Research Institute (<http://www.humboldt.org.co/en>). In particular, we consulted closely with developers and users of the Colombia BON application BioModelos (<http://biomodelos.humboldt.org.co/en>), a system for vetting and hosting species' distribution estimates (Velásquez-Tibatá et al., 2019). This software benefited greatly from the involvement of these conservation practitioner end users at the design stage and through iterative beta-testing of software components during the development process. Our engagement with end users helped us create a product that was deemed effective for the needs of a global audience of stakeholders and users.

CONCLUSIONS

Transparent and repeatable workflows for SDMs and associated biodiversity change metrics are essential to inform evidence-based conservation decision-making such as Red List assessments and conservation offsets (Moilanen et al., 2020). They even have potential to inform new proposed workflows for the IUCN Green Status of Species, which measures the success of a species conservation plan and the potential for recovery (Grace et al., 2021). Because of this, uncertainty related to SDMs should be better communicated to decision makers in a meaningful and informative way (Costa et al., 2018).

changeRangeR allows users to measure several biodiversity metrics for species and communities and report on temporal range trends and maps. It facilitates the inclusion of biodiversity metrics such as species richness alongside other considerations such as complementarity in conservation decisions (Astudillo-Scalia & Albuquerque, 2020; Varzinczak et al., 2020). Further, this package gives users the tools to determine the effects of their analytical decisions through sensitivity tests. SDM post-processing toolboxes such as changeRangeR are an important step for measuring and monitoring biodiversity and represent key extensions for SDM methodology and associated documentation and metadata protocols such as ODMAP (Fitzpatrick et al., 2021; Zurell, Franklin, et al., 2020).

changeRangeR can work in tandem with other SDM software by sharing inputs and outputs. The package is designed to receive inputs from another recent package, maskRangeR (Merow et al., 2022), which provides current species range estimates using remotely sensed data and expert knowledge. Both packages are specifically designed to be included as new modules in Wallace, an R-based interactive platform and workflow for reproducible modeling of species niches and distributions (Kass et al., 2018; Kass, et al. In Press); however, they can also be effectively used independently from Wallace to provide users with flexible tools to measure biodiversity indicators. Wallace's workflow includes other R packages designed to help users follow best practices for SDMs, including ENMeval (Kass, Muscarella, et al., 2021) for tuning model complexity and rmms for model metadata generation (Merow et al., 2019).

The translation of research outputs into effective policy-making can be bolstered by strengthening the collaborative ties between taxonomic experts, modeling experts, and policy makers to progress the transition of science into policy (Urbina-Cardona et al., 2019). These metrics must also be considered alongside other policy recommendations such as sociopolitical acceptability and economic feasibility to result in effective and sustainable conservation impacts (Brown, 2017).

AUTHOR CONTRIBUTIONS

Peter J. Galante and Mary E. Blair led the conceptual design and writing of the manuscript. Peter J. Galante led code development. Samuel Chang Triguero, Andrea Paz, Matthew Aiello-Lammens, Beth E. Gerstner, Jamie M. Kass, and Cory Merow contributed to code development, conceptual design, vignette development, and writing the manuscript. Bethany A. Johnson, Elkin A. Noguera-Urbano, and Gonzalo E. Pinilla-Buitrago contributed to conceptual design, code testing, vignette development, and writing the manuscript.

ACKNOWLEDGMENTS

We thank Erica Johnson, Jorge Velásquez-Tibatá, Dan Rosauer, and Ned Horning for helping to craft the ideas presented in this work and Robert P. Anderson for providing feedback on previous versions of this manuscript. This research was made possible through funding to all authors from the National Aeronautics and Space Administration (Grant No. 80NSSC18K0406 to MEB). AP, BAJ, GEPB, CM, JMK, SC, and MAL acknowledge support from the US National Science Foundation (NSF; DBI-1661510 to Robert P. Anderson). JMK was supported by a Japan Society for the Promotion of Science (JSPS) Postdoctoral Fellowship for Foreign Researchers. BEG was supported by NASA FINESST Award 80NSSC19K133. Finally, we are grateful for feedback provided by all of the participants of our end-user consultation and testing workshops, especially to Cristina Lopez-Gallego, Hector Manuel Arango, Nicolas Urbina, and Leonardo Lemus Mejia, and all of the institutions that our participants represent: Asociación Colombiana de Herpetología, National Audubon Society, The Nature Conservancy—Colombia, Wildlife Conservation Society-Colombia, Comisión Nacional para el Conocimiento y Uso de la Biodiversidad-México, Instituto de Ecología- México, A.C. INECOL, Parques Nacionales Naturales de Colombia, Sociedad Colombiana de Mastozoología, Universidad Nacional-Colombia, Pontificia Universidad Javeriana, Universidad de Antioquia, Universidad del Valle, Universidad EAFIT, the Sistema de Información sobre la Biodiversidad de Colombia, Texas A&M University, McGill University, Temple University, among others. We also acknowledge workshop support from the Alexander von Humboldt Biological Resources Research Institute and the Center for Biodiversity and Conservation at the American Museum of Natural History.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are openly available on GitHub at https://github.com/cmerow/changeRangeR. Updated software vignettes will be maintained on GitHub at https://cmerow.github.io/changeRangeR/. Stable updates to the changeRangeR package will be maintained on CRAN (https://CRAN.R-project.org/package=changeRangeR).