Models combining multiple scales of inference capture hydrologic and climatic drivers of riparian tree distributions

Abstract

Predicting species geographic distributions is key to managing invasive species, conserving biodiversity, and understanding species' environmental requirements. Species distribution models (SDMs) commonly focus on climatic predictors, but other environmental factors can also be essential, particularly for species with specialized habitats defined by hydrologic, topographic, or edaphic conditions (e.g., riparian, wetland, alpine, coastal, serpentine). Here, we demonstrate a novel approach for capturing strong effects of both hydrologic and climatic predictors in SDMs for riparian plants, by merging analyses targeted at environmental drivers within riparian ecosystems and across the western USA (3.8 × 10⁶ km²). We developed presence-background SDMs from five algorithms for three invasive riparian trees (Tamarix ramossisima/chinensis [saltcedar], Elaeagnus angustifolia [Russian olive], and Ulmus pumila [Siberian elm]) and three native Populus spp. (cottonwoods). We used separate background datasets to develop models with different spatial scales of inference: (1) spatially filtered random points to represent available habitat across the study area and (2) target-group points from Salix (willow) occurrences to represent available riparian habitat. Random-background models captured hydrologic drivers of riparian tree distributions relative to the largely upland western USA, whereas Salix-background models captured climatic drivers within the context of riparian ecosystems. Combining predictions from the two backgrounds identified hydrologically suitable habitats within climatically suitable regions, resulting in fewer false “absences” than either background alone, improving predictions over previous SDMs, and providing more complete information to guide management decisions. Surprisingly, the predicted habitat for U. pumila, a newly recognized riparian invader, was as or more extensive than Populus deltoides/fremontii, T. ramossisima/chinensis, and E. angustifolia, the most common riparian tree complexes in the western USA. Watersheds constituting 20% of U. pumila predicted habitat contained no occurrence records, indicating high risk of future and unrecognized invasions. Combining models from random and ecosystem-specific target-group backgrounds may improve SDMs for species from many specialized habitats, providing a method to link predicted distributions to localized geographic features while capturing broad-scale climatic requirements.

INTRODUCTION

Globally, riparian ecosystems are often heavily invaded by non-native species (Hood & Naiman, 2000; Pysek et al., 2010; Stohlgren et al., 1998). Abundant resources, frequent disturbance from flooding and human activity, and long-distance dispersal by hydrochory and zoochory along riparian corridors all facilitate riparian invasion (Richardson et al., 2007). Riparian ecosystems occupy <5% of the landscape but provide important ecosystem services (Riis et al., 2020), supporting similar biodiversity to vastly more extensive uplands (Sabo et al., 2005). Riparian invasion can impact biodiversity, wildlife habitat, floodplain and channel geomorphology, water quality, and recreational value (Tickner et al., 2001).

The effective management of riparian invasive species requires accurate tools to predict where they are likely to occur and spread. Species distribution models (SDMs, also known as habitat suitability models, bioclimatic envelope models, and ecological niche models) are commonly used to predict geographic distributions based on relationships between known occurrences and environmental conditions (Elith & Leathwick, 2009). SDMs are useful tools for predicting invasive species distributions in their invaded range (Peterson, 2003; Thuiller et al., 2005), and can guide management decisions both for control of existing populations and for monitoring for early detection–rapid response to nascent populations (Crall et al., 2013).

Analytical methods for creating SDMs are well-developed, but most SDMs to date have focused on climatic predictors because climate is thought to be the dominant driver of large-scale species distributions, large-scale climate data are readily available, and many SDMs are generated to predict effects of climate change on species distributions (Austin & Van Niel, 2011; Pearson & Dawson, 2003). However, many factors other than climate also influence species' fundamental and realized environmental niches (Mod et al., 2016). In particular, climate-only SDMs are insufficient for predicting the local distributions of species adapted to hydrologically, topographically, or edaphically specialized habitats (e.g., riparian, wetland, coastal, alpine, serpentine), where other environmental factors are fundamental to habitat suitability.

For riparian plants, landscape-scale variation in surface and groundwater hydrology is key to species distributions (Aguiar et al., 2018; McShane et al., 2015). Riparian plants are largely restricted to small portions of the landscape with sufficient access to water to meet relatively high moisture requirements. Their distributions are strongly influenced by streamflow hydrology (timing, magnitude, and duration of high and low flows), which drives fluvial disturbance and moisture dynamics in riparian ecosystems (Poff et al., 1997). Altered streamflow hydrology due to water management (flow regulation, diversions, and groundwater pumping) alters natural disturbance and moisture regimes, and thus alters riparian community composition and facilitates invasion (Catford et al., 2011; Stromberg et al., 2007). Riparian plants vary considerably in hydrologic requirements, with interspecific differences in adaption to inundation, drought stress, and fluvial disturbance (Merritt et al., 2010). SDMs must capture these hydrologic requirements to predict where particular riparian species will occur at the scale of riparian ecosystems.

To date, few SDMs for riparian species have been developed, and most have included little hydrologic information. Approaches for including hydrologic information have varied, particularly between large-scale and small-scale riparian SDMs. Among large-scale SDMs intended to model distributions of widespread riparian taxa at regional or continental scales, most have included no hydrologic predictors, thus ignoring the fundamental roles of hydrology in riparian species distributions and artificially inflating the spatial extent of predicted suitable habitat across uplands. A few riparian SDMs, at both large and small scales, have included a single hydrologic predictor such as distance-to-surface water (e.g., Jarnevich et al., 2011) or length of shoreline (e.g., Collette & Pither, 2015). These simple hydrologic predictors are relatively easy to derive from available, large-scale hydrologic datasets and can capture hydrologic drivers of riparian species distributions to a limited extent, but they ignore the fundamental roles of topography and streamflow hydrology in determining riparian moisture dynamics and fluvial disturbance. Among small-scale riparian SDMs, some have defined the boundaries of hydrologically suitable habitat a priori, based on distance from the river (e.g., Collingham et al., 2000), aerial imagery of the floodplain (e.g., Nylen et al., 2019), or simulated inundation frequency (e.g., Murray et al., 2012). This approach requires making a priori assumptions about species' hydrologic requirements and reduces the utility of the resulting SDM for predicting and characterizing hydrologically suitable habitats. Further, it is often challenging to apply across large spatial scales because accurately and consistently delineating floodplain boundaries and/or deriving inundation frequencies is difficult and labor intensive. Some small-scale riparian SDMs have included multiple, mechanistic hydrologic predictors, but at such small scales that they did not include climatic drivers of species distributions (e.g., Fu & Guillaume, 2014; Nylen et al., 2019).

To be most useful, SDMs for widespread riparian invasive species must integrate fundamental effects of climate and hydrology at large spatial scales (representing the species' full potential range) and at fine resolution (from 10s to 100s of meters, at the scale of gradients from riparian to upland habitats) (Collingham et al., 2000). Some SDMs have employed a “hierarchical” approach to address the problem of modeling fine-resolution environmental drivers across large spatial scales by first modeling climatically suitable habitats at a broad spatial scale, and then modeling effects of finer-resolution predictors (e.g., land use) at a smaller scale within areas identified as climatically suitable (Bradley, 2010; Pearson et al., 2004). However, modeling climate and other environmental drivers separately in a hierarchical framework prevents explicit analysis of interactions between climate and other predictors. Such interactions are particularly likely to be important when both climate and other predictors influence the same environmental condition, such as moisture or temperature. For example, proximity to surface water in drier climates, or south-facing aspects in cooler climates (Austin & Van Niel, 2011), may support occurrence in areas that otherwise would be climatically unsuitable.

In the present study, we demonstrate an alternative approach to developing riparian SDMs that can capture strong effects of both climate and hydrology at fine resolution across large spatial scales. This approach involves combining predictions from separate models developed using different background points (i.e., pseudo-absences), which provide different spatial scales of inference for identifying environmental drivers of habitat suitability.

Many SDMs are developed based on statistical comparisons between occurrence locations and background points, which are intended to represent the full range of available habitats when reliable absence data are unavailable (Phillips et al., 2009). Background points are often selected at random, but occurrence sampling is typically non-random and correlated with environmental conditions, so the resulting SDMs may represent spatial variation in sampling intensity rather than habitat suitability. One common solution is to spatially filter (i.e., thin) occurrence records to reduce clustering due to survey bias. Alternatively, background points are sometimes selected from occurrence records of other taxa (i.e., “target-group” background points), in an attempt to match, and thereby account for, survey bias in occurrence records of the study taxa (Phillips et al., 2009).

Because random and target-group background points differ in placement and spatial extent, the spatial scale of inference differs for the resulting SDMs (Jarnevich et al., 2015). With random background points, the scale of inference spans the full geographic extent of the sampled area, assuming survey bias is fully removed by spatial filtering. The resulting SDMs reflect differences in environmental conditions between the study taxon habitat and the broader landscape. In contrast, with target-group background points, the scale of inference is restricted to environmental conditions where the study taxa and/or target group can occur. The resulting SDMs reflect differences in suitable habitat between the study taxa and target group (Vollering et al., 2019). For example, for riparian SDMs, target-group background points could be selected from occurrence records of other riparian taxa, making the scale of inference specific to riparian ecosystems and reflecting environmental conditions that distinguish study taxon habitat from other riparian habitats. For riparian species, the climatic and hydrologic conditions that distinguish the study taxon habitat from the broader landscape may differ considerably from conditions that distinguish the study taxon habitat from other riparian habitats. Thus, random and target-group background datasets may result in very different SDMs, emphasizing different environmental drivers that influence habitat suitability within their different scales of inference. Combining their predictions may result in SDMs that can capture the fundamental roles of both climate and other, specialized habitat conditions in defining species distributions.

Across western USA riparian ecosystems, three of the most common woody taxa are non-native: Tamarix ramosissima/chinensis (T. ramosissima Ledeb., T. chinensis Lour., hybrids; saltcedar), Elaeagnus angustifolia L. (Russian olive), and Ulmus pumila L. (Siberian elm) (Friedman et al., 2005; Perry et al., 2018). T. ramosissima/chinensis and E. angustifolia invasions have received considerable scientific and management attention (Shafroth et al., 2010; Sher & Quigley, 2013). Although U. pumila invasion has received less attention, a recent survey revealed that U. pumila was as or more common than T. ramosissima/chinensis and E. angustifolia across 282,000 km² of the western USA (Perry et al., 2018). Accurate SDMs for these taxa are needed to prioritize monitoring and management efforts, especially for U. pumila, as no SDMs for it have yet been developed and little is known about its distribution and potential spread. SDMs that include hydrologic information would be particularly useful because these invasive species differ from one another, and from native riparian trees, in their hydrologic requirements. For example, all three invasive species are more drought tolerant than native riparian Populus (L., cottonwoods) and Salix (L., willows) (Lite & Stromberg, 2005; Perry et al., 2013; Reynolds & Cooper, 2010), and E. angustifolia and U. pumila are larger seeded, which make them less reliant on fluvial disturbance to create bare, moist surfaces for germination and seedling establishment than T. ramosissima/chinensis, Populus, and Salix (Nagler et al., 2011). As a result, E. angustifolia and U. pumila can occur in drier riparian areas and uplands (Olsen & Knopf, 1986; Park et al., 2016). Several large-scale SDMs have been developed for T. ramosissima/chinensis and E. angustifolia, but none of them included hydrologic predictors beyond distance-to-water or shoreline length (Bradley et al., 2009; Collette & Pither, 2015; Jarnevich et al., 2011; Jarnevich & Reynolds, 2011; Kerns et al., 2009; Liu et al., 2014; Morisette et al., 2006). Most of these SDMs used random background datasets; one used a target-group background dataset, with marginal model performance (Cord et al., 2010).

In this study, we developed SDMs across the western USA for U. pumila, T. ramosissima/chinensis, E. angustifolia, and several native Populus that dominate western riparian forests (Friedman et al., 2005). Our primary objective was to predict the potential distributions of these widespread invasive and native riparian trees. To accomplish this, our second objective was to evaluate the utility of combining models from random and riparian-specific, target-group background datasets to capture climatic and hydrologic drivers in riparian SDMs. Our third objective was to use the resulting predictions to help guide management and monitoring decisions, by characterizing environmental drivers of species occurrence, comparing species' potential ecoregional distributions, and identifying areas at risk of undetected and future spread by the invasive taxa.

METHODS

We developed SDMs at 90-m resolution across the conterminous USA west of the 100th Meridian (hereafter, western USA) for T. ramosissima/chinensis, E. angustifolia, U. pumila, and the following three Populus taxonomic groups: (1) “Populus deltoides/fremontii” (section Aigeiros; P. deltoides spp. monilifera [Aiton] Eckenwalder [plains cottonwood], P. deltoides spp. wislizeni [S.Watson] Eckenwalder [Rio Grande cottonwood], P. fremontii S.Watson [Fremont cottonwood]), (2) “P. balsamifera” (section Tacamahaca; P. balsamifera spp. trichocarpa [Torr. & A.Gray ex Hook.] Brayshaw [black cottonwood], P. balsamifera L. ssp. balsamifera [balsam poplar]), and (3) “P. angustifolia” (section Tacamahaca; James, narrow-leaf cottonwood). All of these Populus can hybridize. Following McShane et al. (2015), we grouped Populus for which even non-hybrids are difficult to distinguish reliably where their ranges overlap.

Occurrence records

We obtained occurrence records from regional, national, and international databases and large-scale ecological datasets (Appendix S1: Table S1). We searched for all taxonomic synonyms listed in the United States Department of Agriculture (USDA) plants database (https://plants.sc.egov.usda.gov/). To reduce location errors, we removed all records collected prior to 1995 or for which the geographic coordinates were missing, had less than three decimal places, or had zero degree-seconds. We thinned to one record per taxon per 90-m pixel (Figure 1). For model training, we used occurrence records obtained from all sources in January 2020. For taxa with above 10,000 records, we randomly selected 10,000 for model training. For model evaluation, we used the remaining occurrence records (for taxa with >10,000 records), plus new records that were added to regional, national, and international databases between January 2020 and May 2021, as independent test data.

FIGURE 1. Occurrence record locations for (a) Ulmus pumila, (b) Tamarix ramosissima/chinensis, (c) Elaeagnus angustifolia, (d) Populus deltoides/fremontii, (e) Populus balsamifera, and (f) Populus angustifolia used in species distribution model training and assessment for the western USA. Spatially filtered training data were used for model training with random background points (minimum distance of 50 km between points). Unfiltered training data were used for model training with Salix background points. Independent test data were used for model assessment and not for model training. The total number of occurrence records in Salix (and random)-background model training was 1184 (192) U. pumila, 10,000 (504) T. ramosissima/chinensis, 10,000 (341) E. angustifolia, 3079 (381) P. deltoides/fremontii, 1021 (198) P. balsamifera, and 1279 (183) P. angustifolia.

Environmental predictors

We examined 48 environmental predictors describing climate, topography, hydrology, and land use (Appendix S1: Table S2). The climate predictors included bioclimatic variables, land surface temperature/emissivity, and seasonal mean and minimum temperatures, precipitation, potential water deficits, and evapotranspiration. The topographic predictors included elevation, slope, and normalized topographic position index (TPI) to indicate low landscape position (De Reu et al., 2013). Hydrologic predictors included distance-to-water, stream order, annual streamflow, upstream drainage area, upstream ditch length, and upstream reservoir storage volume. Land-use predictors included the global human modification index and percent cover by impervious surfaces, developed land uses, and agricultural land uses. To reduce collinearity, we selected subsets of predictors for analysis such that no included predictors were strongly correlated (Pearson |r|, Spearman |ρ|, Kendall |τ| ≤ 0.7; Dormann et al., 2013; Appendix S1: Table S2).

Species distribution modeling

We used the United States Geological Survey (USGS) Software for Assisted Habitat Modeling (SAHM, v2.1.2; Morisette et al., 2013), a package in VisTrails (Freire et al., 2006). The modeling workflow is depicted in Appendix S2: Figure S1. We developed model ensembles for each taxon using all five statistical algorithms available in SAHM, which include a range of algorithm types commonly used in SDMs: boosted regression trees (Elith et al., 2008), generalized linear models (McCullagh & Nelder, 1989), multivariate adaptive regression splines (Elith & Leathwick, 2007), maximum entropy modeling (Phillips et al., 2006), and random forest analyses (Friedman, 2001).

We developed separate models using the following two background datasets: (1) spatially filtered random background points and (2) target-group occurrence records of Salix species (Appendix S2: Figure S2). For random-background models, we randomly selected 10,000 points across the western USA, and spatially filtered occurrences and background points to a minimum distance of 50 km. For the Salix-background models, we randomly selected 10,000 points from 21,962 Salix occurrence records (Appendix S1: Table S1). Seventy Salix species occur in the western USA (Appendix S1: Table S3), including some of the most common woody riparian taxa (Friedman et al., 2005). All are riparian and/or wetland species, and range in wetland indicator status from facultative to obligate, making their collective distribution broadly representative of the range of riparian and wetland conditions in the study region. As woody riparian species, they are also relatively likely to have been sampled with similar survey bias to the study taxa.

We performed 10-fold cross-validation for each model, by splitting occurrence records into 10 equal-sized random subsets (see Figure 1 for sample sizes) and training the model 10 times, each time withholding one subset for testing. We assessed model performance using the area under the receiver operating characteristic curve (AUC) and the true skill statistic (TSS). We used SAHM, v2.1.2, default parameter values for all models (Appendix S1: Table S4), revisiting these when overfitting was indicated by differences between training and cross-validation AUC > 0.05 or by implausibly complex response curves. Visual examination of model residuals indicated no evidence of remaining spatial autocorrelation.

We computed multivariate environmental similarity surfaces (MESS) to quantify the degree of extrapolation necessary to project species distributions across the western USA based on the training data. These indicated that the full parameter space across the western USA was well-represented during model training; areas requiring extrapolation to conditions outside the range of occurrence and background points were few, small, and scattered.

To compare the importance of different environmental predictors in models with different background datasets, we calculated predictor permutation importance for each model as the change in AUC when values for that predictor were randomly assigned to the occurrence and background designations. We then averaged predictor importance values across statistical algorithms for each background dataset.

To assess the consistency of habitat suitability across statistical algorithms for each background dataset, we created discretized versions (binary maps) of the continuous output maps for each model, with locations assigned to a suitable or unsuitable category. We discretized the maps using the following three suitability thresholds (in order of increasing conservatism): (1) minimum relative suitability for 99% of occurrence records (hereafter, X1st), (2) minimum relative suitability for 90% of occurrence records (hereafter, X10th), and (3) maximized sum of sensitivity and specificity (MaxSS) (Appendix S1: Table S5). The X1st and X10th thresholds served as presence-only counterparts to MaxSS, which assumes background points are absences. Discretizing the maps made it possible to compare predicted species distributions across models from the different algorithms; by contrast, the continuous output from presence-background models indicates relative habitat suitability (i.e., compared with the other values within each model), leading to different meanings of the same values from different models. To generate ensemble maps for each taxon with each background dataset, we summed the discretized maps across the five statistical algorithms (maximum = 5; hereafter, random-background and Salix-background maps). Further, we summed random- and Salix-background maps to generate ensemble maps across the two backgrounds (maximum = 10; hereafter, random + Salix maps).

To compare model performance among random-background, Salix-background, and random-Salix ensemble maps, we examined the accuracy of predicted habitat suitability within each ensemble map for occurrence records not used in model training (independent test data).

Calculations for ecological and management implications

To characterize areas at risk of invasion by different taxa, we quantified the ensemble suitable habitat for each taxon within US Environmental Protection Agency ecoregions, by summing habitat suitability across all pixels in each Level I and Level III ecoregion and then dividing by the maximum possible value (total number of pixels in the ecoregion × total number of models). Maps of the US EPA ecoregions are available at https://www.epa.gov/eco-research/ecoregions and in Omernik and Griffith (2014).

To identify areas at higher risk of undetected or future invasion, we identified major watersheds (USGS hydrologic unit code 6 [HUC6]) where each taxon was under-represented in occurrence records relative to predicted habitat suitability. For each HUC6 watershed within the study area, we calculated the difference between actual occurrence record density and the density of occurrence records that would be expected if they were distributed among watersheds in proportion to habitat suitability in MaxSS random + Salix maps, using the following equation:[Image Omitted. See PDF]

Here, n_HUC (the number of occurrence records within the watershed) divided by km²_HUC (area of the watershed) calculated the actual occurrence record density, and ∑ens_suitability_HUC (sum of ensemble habitat suitability across pixels within the watershed) divided by ∑ens_suitability_total (sum of ensemble habitat suitability across the study area) calculated the proportion of total suitable habitat that occurred within the watershed. Multiplying this proportion by n_total (total number of occurrence records across the study area) apportioned occurrence records to the watershed in proportion to habitat suitability. Dividing by km²_HUC transformed this expected number of occurrences into density units. To assess current knowledge of species occurrences as fully as possible, we included all occurrence records used for model training and independent test data, as well as records from 1980 to 1994 from Appendix S1: Table S1 sources. All records were spatially filtered to a minimum distance of 50 km within each watershed to reduce the effects of survey bias.

RESULTS Effects of background dataset on model performance and predictions Model performance

All models performed well for all taxa, with high discrimination capacity in cross-validations (Appendix S1: Table S6). The average cross-validation AUC was >0.85 for 88% of models (minimum = 0.79). Average cross-validation TSS was >0.55 for 73% of models and >0.4 for all but two models (minimum = 0.36).

Random + Salix maps were consistently superior to random- and Salix-background maps for reducing false “absences” in independent test data (Figure 2; Appendix S2: Figures S3 and S4). In MaxSS maps, there were 1%–9% fewer occurrences where all models predicted unsuitable habitat in random + Salix maps than in random- or Salix-background maps, and 2%–18% fewer occurrences where >30% of models predicted unsuitable habitat, indicating that random- and Salix-background maps alone each missed different components of suitable habitat. Accordingly, more occurrences (by 1%–16%) were successfully predicted as suitable habitat by ≥50% of models in random + Salix maps (excepting for U. pumila random-background maps), suggesting that agreement between random- and Salix-background maps strengthened predictions for areas with suitable habitat. However, a complete agreement among models was less common in random + Salix maps for T. ramosissima/chinensis and P. deltoides/fremontii, with 6%–16% fewer occurrences predicted as suitable habitat by 90%–100% of models in random + Salix maps than in random-background maps.

FIGURE 2. Ensemble map validation: habitat suitability at occurrence record locations for independent test data, for (a–c) Ulmus pumila, (d–f) Tamarix ramosissima/chinensis, (g–i) Elaeagnus angustifolia, (j–l) Populus deltoides/fremontii, (m–o) Populus balsamifera, and (p–r) Populus angustifolia. Bars show the number of occurrence records that received each possible value for ensemble habitat suitability based on maximized sum of sensitivity and specificity suitability thresholds: 0–5 models for (a, d, g, j, m, p) random-background maps and (b, e, h, k, n, q) Salix-background maps and 0–10 models for (c, f, i, l, o, r) random + Salix maps.

Predictor importance

Predictor importance differed strongly between models from different background datasets. In random-background models, TPI₁₆ and/or distance-to-water were among the most important predictors (Table 1; Appendix S1: Table S7), reflecting the riparian nature of the study taxa relative to the predominately upland surrounding landscape. This result was strongest for P. deltoides/fremontii and P. angustifolia, and weakest for U. pumila. Developed land cover (dev_90-m) and human modification index were important for U. pumila and to a lesser extent E. angustifolia and T. ramosissima/chinensis (Appendix S2: Figures S5–S7). In addition, low annual precipitation was important for U. pumila, T. ramosissima/chinensis, and E. angustifolia (Appendix S2: Figures S5–S7).

TABLE 1 Estimated importance of each predictor variable (in percentage) in species distribution models for the six study taxa across the western USA.

Predictor variablea by type	Random background						Salix background
	U	T	E	Pd/f	Pb	Pa	U	T	E	Pd/f	Pb	Pa
Climate
maxTwarm_mo	5	3	3	8	3	8	4	57	38	60	18	27
minTwinter	4	1	9	1	1	7	40	3	21	4	11	20
meanTwet_qtr	5	1	1	3	21	4	22	2	6	9	24	5
meanPannual	13	20	15	9	4	2	9	19	8	2	4	3
meanPwarm_qtr	3	2	3	1	1	6	5	3	9	3	6	17
Tdiurnal_range	2	4	2	2	5	1	3	2	3	3	19	13
Isothermality	1	1	0	…	1	2	3	6	4	…	4	7
Pseasonality	…	1	1	…	…	…	…	4	5	…	…	…
ETa_summer	4	3	5	3	2	3	1	0	0	2	1	1
Topography
TPI16	10	28	15	29	18	27	2	0	0	1	2	1
TPI0.4	1	4	1	3	3	2	0	0	0	1	1	0
slope	1	1	1	1	6	3	1	0	0	1	1	3
Hydrology
dist_water	5	11	12	15	6	27	2	0	0	0	0	1
annual_flow	2	2	2	5	16	2	1	1	1	1	6	1
reservoirsupstrm	1	2	1	5	2	0	0	0	1	1	1	0
ditchesupstrm	1	1	1	1	1	1	1	0	1	1	0	0
reservoirsHUC12	0	0	0	0	0	0	0	0	0	0	1	0
Land use
HMc	12	1	5	1	2	1	6	1	3	10	1	1
dev90-m	21	11	15	8	2	1	1	1	0	0	0	0
agr90-m	9	2	4	3	5	2	0	0	0	0	0	0
agrHUC12	0	1	1	2	0	0	0	0	0	1	0	0
impervious	0	0	0	0	0	0	0	0	0	0	0	0

Note: Values shown are ensemble means of relativized permutation importance from five models developed with different statistical algorithms. Importances ≥10 are shown in boldface. Ellipsis (…) indicates a variable was excluded from analysis because of strong correlations (|r| > 0.7) with other variables. Abbreviations: E, Elaeagnus angustifolia; Pa, Populus angustifolia; Pb, Populus balsamifera; Pd/f, Populus deltoides/fremontii; T, Tamarix ramosissima/chinensis; U, Ulmus pumila. a For detailed information on predictor variables, see Appendix S1: Table S1.

In Salix-background models, topographic, hydrologic, and land-use variables were less important (Table 1; Appendix S1: Table S8). Instead, the maximum temperature of the warmest month (maxT_{warm_mo}) was the most frequently important predictor, with higher relative suitability for T. ramosissima/chinensis, E. angustifolia, and P. deltoides/fremontii with hot summers, and P. balsamifera and P. angustifolia with moderate summers (Figure 3). In contrast, for U. pumila, minimum winter temperature (minT_winter) was most important, with higher relative suitability with moderate winters. MinT_winter was also important for E. angustifolia and P. angustifolia, with higher relative suitability with cold winters. Low annual precipitation remained important for T. ramosissima/chinensis (Appendix S2: Figure S6).

FIGURE 3. Predictor response curves for (a) maximum temperature of the warmest month (maxTwarm_mo) and (b) minimum winter (December–February) temperature (minTwinter), from species distribution models trained with random (brown) and Salix (blue) background points for Ulmus pumila, Tamarix ramosissima/chinensis, Elaeagnus angustifolia, Populus deltoides/fremontii, Populus balsamifera, and Populus angustifolia across the western USA. Models were developed using each of five statistical algorithms, namely, boosted regression trees (BRT), generalized linear models (GLM), multivariate adaptive regression splines (MARS), maximum entropy modeling (MAXENT), and random forest analyses (RF). Values below species names are relativized predictor importance: mean (minimum, maximum). Red hash marks along the x-axes indicate predictor values at occurrence record locations.

Predicted distributions

In random-background maps, predicted distributions mainly followed the distribution of surface water across the western USA, in accord with the importance of TPI₁₆ and distance-to-water (Figures 4 and 5; Appendix S3: Figures S1–S6; Appendix S4: Figures S1–S6). U. pumila, T. ramosissima/chinensis, E. angustifolia, and P. deltoides/fremontii predicted distributions also followed roads and other human development. In contrast, in Salix-background maps, suitable habitat was sometimes concentrated along large rivers but the importance of hydrology to these riparian taxa was not generally apparent. Instead, predicted distributions mainly reflected broad-scale climatic variation. Combining the predicted distributions from the two backgrounds resulted in random + Salix maps that reflected both broad-scale climatic variation and local-scale hydrologic variation, with suitable habitats concentrated in low topographic positions near surface water (and/or human development) within larger, climatically suitable regions (Figures 4 and 5; Appendix S4: Figures S1–S6).

FIGURE 4. Ensemble maps of suitable habitat for the invasive study taxa (Ulmus pumila, Tamarix ramosissima/chinensis, and Elaeagnus angustifolia) across the western USA. Ensembles depicted are from (a–c) five statistical algorithms trained with spatially filtered random background points, (d–f) five statistical algorithms trained with Salix background points, and (g–i) all 10 models (five statistical algorithms with each background dataset). To create these ensemble maps, we transformed continuous relative habitat suitability (0%–100%) from each model (Appendix S3: Figures S1–S3) into binary predictions of suitable versus unsuitable habitat based on the maximized sum of sensitivity and specificity (MaxSS), and counted models with suitable habitat within each 90-m pixel for a maximum of 5 in (a–f) and 10 in (g–i).

FIGURE 5. Ensemble maps of suitable habitat for the native study taxa (Populus deltoides/fremontii, Populus balsamifera, and Populus angustifolia) across the western USA. (a–c) Ensembles depicted are from (a–c) five statistical algorithms trained with spatially filtered random background points, (d–f) five statistical algorithms trained with Salix background points, and (g–i) all 10 models (five statistical algorithms with each background dataset). To create these ensemble maps, we transformed continuous relative habitat suitability (0%–100%) from each model (Appendix S3: Figures S4–S6) into binary predictions of suitable versus unsuitable habitat based on the maximized sum of sensitivity and specificity (MaxSS), and counted models with suitable habitat within each 90-m pixel for a maximum of 5 in (a–f) and 10 in (g–i).

Ecological and management implications of predicted distributions Differences among taxa in extent of suitable habitat

Predicted habitats for P. deltoides/fremontii, U. pumila, and T. ramosissima/chinensis were extensive, each constituting 25%–50% of maximum possible habitat suitability across the western USA in X10th and MaxSS maps (Figure 6; Appendix S4: Figures S9–S10). P. angustifolia habitat was slightly less extensive (23%–28%), and E. angustifolia and P. balsamifera habitats were considerably less extensive (16%–22%).

FIGURE 6. Percent of total area across the western USA and within US Environmental Protection Agency Level I ecoregions that was identified as suitable habitat for each study taxon in random + Salix maps, based on each of three thresholds for identifying suitable habitat: (1) minimum relative suitability for 99% of occurrence records (X1st), (2) minimum relative suitability for 90% of occurrence records (X10th), and (3) the maximized sum of sensitivity and specificity (MaxSS). Values are the proportion of 90-m pixels containing predicted suitable habitat, weighted by the number of models (maximum = 10) that predicted suitable habitat in each pixel. E, Elaeagnus angustifolia; Pa, Populus angustifolia; Pb, Populus balsamifera; Pd/f, Populus deltoides/fremontii; T, Tamarix ramosissima/chinensis; U, Ulmus pumila.

In North American Deserts, which constitute 37% of the western United States, T. ramosissima/chinensis and P. deltoides/fremontii had more extensive habitats than other study taxa (Figure 6; Appendix S4: Figures S9 and S10). In X10th and MaxSS maps, T. ramosissima/chinensis habitat was also more extensive than P. deltoides/fremontii. However, within different Level III sub-ecoregions of the North American Desert (Appendix S1: Tables S9–S11), other taxa had similar or more extensive habitat than T. ramosissima/chinensis, including P. deltoides/fremontii in Warm Deserts and southeastern Cold Deserts, E. angustifolia in southern and central Cold Deserts, P. balsamifera in northern Cold Deserts, P. angustifolia in eastern Cold Deserts, and U. pumila in southeastern and northern Cold Deserts. T. ramosissima/chinensis habitat was least extensive in northern Cold Deserts.

In contrast, in the Great Plains, which constitute 31% of the western USA, U. pumila and P. deltoides/fremontii had more extensive habitats than other study taxa across the ecoregion and within most Level III sub-ecoregions. U. pumila and P. deltoides/fremontii habitats were less extensive in the northern Great Plains, where they were only moderately more extensive than P. angustifolia, T. ramosissima/chinensis, and E. angustifolia. U. pumila habitat was also less extensive in the southernmost Great Plains, where only P. deltoides/fremontii habitat was common.

These ecoregional patterns were apparent in all ensemble maps. However, differences among taxa in the extent of suitable habitat were smaller in random- than Salix-background maps, because the importance of topographic, hydrologic, and land-use predictors in random-background models resulted in relatively similar predicted distributions across taxa, at low topographic positions near surface water and/or human development.

Risk of undetected or future invasion

We identified major watersheds at higher risk of undetected or future invasion as those with low occurrence record densities relative to what would be expected if occurrence records were distributed in proportion to habitat suitability (Figure 7; Appendix S4: Figure S11). U. pumila was most under-represented in occurrence records relative to habitat suitability in portions of the Missouri, Arkansas, and Red River Basins in the northwestern and central Great Plains. Strikingly, 25 watersheds within these basins contained 20% of U. pumila predicted suitable habitat but no occurrence records.

FIGURE 7. Under- and over-representation in occurrence records relative to predicted habitat suitability for the invasive study taxa (Ulmus pumila, Tamarix ramosissima/chinensis, and Elaeagnus angustifolia) in United States Geological Survey hydrologic unit code 6 watersheds (HUC6) across the western USA. (a–c) Density of occurrence records within each watershed (in number per square kilometer). (d–f) Proportion of total habitat suitability across the study region contained within each watershed, adjusted for the watershed area (in percentage per square kilometer; from random + Salix maps generated using the maximized sum of sensitivity and specificity). (g–i) Difference between the actual occurrence record density and the occurrence record density that would be expected if occurrence records were distributed in proportion to habitat suitability within each watershed (in number per square kilometer). Negative (orange) values in (g–i) indicate watersheds with few occurrence records relative to predicted habitat suitability (i.e., dark shades in (a–c) but light shades in (d–f)), suggesting potential undetected or future invasion.

Like U. pumila, E. angustifolia was under-represented in occurrence records relative to habitat suitability in portions of the Missouri, Arkansas, and Red River Basins, although often not in the same portions as U. pumila. In contrast, T. ramosissima/chinensis was more often under-represented in portions of the lower Colorado River and Rio Grande Basins, in Warm Desert and southern Cold Desert. Few watersheds with substantial E. angustifolia or T. ramosissima/chinensis under-representation contained no occurrence records (1%–2% of predicted suitable habitat). In several large Cold Desert watersheds, E. angustifolia and/or T. ramosissima/chinensis were under-represented relative to abundant suitable habitat despite having many occurrence records.

DISCUSSION Modeling riparian species distributions using multiple background datasets

The SDMs developed here illustrate an effective approach for capturing the fundamental roles of both hydrologic and climatic predictors in riparian SDMs, by combining predictions from separate models based on random and riparian-specific, target-group background datasets. Because the random- and Salix-background datasets differed in placement and spatial extent, the resulting SDMs emphasized distinct environmental drivers and each alone overlooked key limits to riparian plant distributions. The random background points provided a broad spatial extent for comparison across the largely upland and undeveloped western USA, resulting in models that emphasized topographic, hydrologic, and land-use drivers of riparian plant distributions. However, the random-background models did a poor job of discerning climatic drivers, perhaps because large differences in topographic, hydrologic, and land-use predictors between occurrence and random background points obscured the role of climate. In contrast, the Salix background points restricted the spatial extent for comparison to riparian woody plant habitats. The resulting models emphasized climatic drivers, identifying climatic conditions that distinguished the study taxa habitat from Salix habitat within the context of riparian ecosystems. However, the Salix-background models did a poor job of discerning hydrologic drivers, both here and in a previous T. ramosissima/chinensis SDM (Cord et al., 2010), perhaps because hydrologic predictors differed relatively little between Salix and the study taxa. The novel approach of combining predictions from random and target-group backgrounds corrected their separate weaknesses, generating maps that identified hydrologically suitable habitats within climatically suitable regions (Figures 4 and 5).

The resulting random + Salix maps led to fewer false “absences” in independent test data than model ensembles from either background alone (Figure 2). Further, they were an improvement over most previous SDMs for E. angustifolia and T. ramosissima/chinensis, and effectively predicted well-known Populus distributions. For E. angustifolia, two previous large-scale SDMs generated maps that closely matched our Salix-background maps but did not identify hydrologically suitable habitats, because hydrologic predictors were not examined (Liu et al., 2014) or the hydrologic metric (shoreline length) was insufficient to characterize access to water (Collette & Pither, 2015). A third previous SDM included distance-to-water but failed to predict the southern limit to E. angustifolia's distribution (Guilbault et al., 2012; Jarnevich & Reynolds, 2011), perhaps because the importance of distance-to-water relative to random background points partially obscured effects of climate.

For T. ramosissima/chinensis, several previous large-scale SDMs did not attempt to model hydrologic drivers, had weaker model performance than our SDMs, and predicted substantially broader suitable habitat than random + Salix maps (Bradley et al., 2009; Cord et al., 2010; Ikeda et al., 2014; Liu et al., 2014). Two others included distance-to-water and discerned strong effects of both temperature and distance-to-water despite using only random background points (Jarnevich et al., 2011; Kerns et al., 2009). In Jarnevich et al. (2011), distance-to-water was twice as important as any climatic predictor and predicted habitat extended farther into the Great Plains and Mediterranean California than in random + Salix maps. In Kerns et al. (2009), the SDM spanned only a portion of the northwestern USA and used different suitability thresholds, so the predicted distributions are difficult to compare.

For Populus, two previous large-scale SDMs did not attempt to capture hydrologic drivers and predicted broader suitable habitat for P. fremontii (Ikeda et al., 2014) and narrower suitable habitat for P. angustifolia (Bothwell et al., 2021) than random + Salix maps. Random + Salix maps largely matched Populus spp. distribution maps based on historical occurrence records (Little, 1971), with P. deltoides/fremontii in the eastern and southwestern portions of the study area, P. balsamifera in the northwest, and P. angustifolia at higher elevations. However, random + Salix maps suggested broader suitable habitats, probably both because they incorporated more thorough, current occurrence information (Figure 1) and because they included areas with potential habitats similar to known occurrences. Compared to Little's maps, P. deltoides/fremontii predicted habitat extended farther south and southwest, P. balsamifera habitat extended farther southeast, and P. angustifolia habitat extended farther from the Northwestern Forested Mountains into the Great Plains and Cold Deserts.

Combining models from random and target-group background datasets could also improve SDMs for species from other hydrologically, topographically, or edaphically specialized habitats, such as wetland, coastal, alpine, and serpentine species. For specialized habitats, accurate SDMs must capture the environmental conditions that restrict species to specific geographic features, while also capturing other, broad-scale climatic and environmental requirements. Comparisons with the full study area based on random background points are likely to capture the relevant geographic features, but the overwhelming importance of those features to species distributions may partially obscure other environmental drivers. Target-group background points based on taxa from similar specialized habitats can reveal environmental conditions that distinguish study taxa habitats from other, similar taxa, but cannot capture geographic features that drive both study taxa and target-group distributions (Vollering et al., 2019). Combining information from both can provide more complete characterizations of species distributions for specialized habitats.

Environmental drivers of riparian tree distributions

SDMs can be useful for generating hypotheses about mechanistic drivers of species distributions (Guisan & Thuiller, 2005). Understanding invasive and native species' ecological niches in turn can help to inform invasive species management. Bringing together the SDMs from different background datasets provided a more complete understanding of environmental drivers of riparian woody species distributions than either would have alone. Random-background SDMs revealed differences among native and invasive taxa in the importance and roles of topographic, hydrologic, and land-use drivers, while Salix-background SDMs revealed differences among taxa in climatic drivers and potential responses to climate change.

Hydrology and topography

The combined importance of low topographic position (TPI₁₆) and distance-to-water in most random-background models across study taxa suggests that these two variables captured different components of hydrologically suitable habitat, with TPI₁₆ identifying low areas with potentially shallow water tables, including those relatively far from surface water, and distance-to-water identifying areas close to surface water, including those not much lower than the surrounding landscape. TPI has not been used commonly in SDMs, but its consistent importance here suggests it may often improve riparian SDMs. Predictions based on both TPI and distance-to-water may provide more complete information for riparian invasive species risk assessments, especially for areas that are farther from major surface waters.

The lower importance of TPI₁₆ and distance-to-water for U. pumila supports the hypothesis that U. pumila is more drought-tolerant than the other taxa. In its native range, U. pumila occupies both riparian forests and sparse, upland woodlands in semi-arid savannas, sandlands, and montane steppe (Dulamsuren et al., 2009; Park et al., 2016; Su et al., 2014). In its introduced range, water stress reduces seedling growth of U. pumila less than common Populus and Salix species (Perry et al., 2013), and adult U. pumila are more likely to dominate riparian communities along intermittent than perennial streams (Perry et al., 2018). Shallow U. pumila response curves for TPI₁₆ and distance-to-water in our models predicted a broader hydrologic range of habitats at risk of U. pumila invasion, with suitable habitats less tightly linked to river corridors (Appendix S4: Figures S7 and S8). TPI₁₆ and/or distance-to-water were also less important for E. angustifolia and P. balsamifera, in accord with E. angustifolia greater drought-tolerance and lower reliance on fluvial disturbance for establishment than Populus spp. and T. ramosissima/chinensis (Reynolds & Cooper, 2010) and with P. balsamifera upland occurrences in wetter portions of its range (e.g., Pacific Northwest).

Contrary to our expectations, annual streamflow and water management metrics were unimportant predictors. At the broad spatial scale of random background points, this suggests that proximity to water was more informative than streamflow hydrology for identifying suitable riparian habitat relative to uplands. At the narrower scale of comparison with Salix background points, this suggests that annual streamflow and water management differed relatively little between study taxa and Salix occurrences. This is surprising because the study taxa are known to differ widely in drought-tolerance and response to water management, with low streamflow and altered flood timing and magnitude benefitting T. ramosissima/chinensis and E. angustifolia over Populus and Salix species (McShane et al., 2015; Merritt & Poff, 2010; Nagler et al., 2011; Stromberg et al., 2007). Large differences in climate between study taxa and Salix occurrences may have obscured comparatively smaller effects of hydrologic predictors. Also, discerning hydrological differences in suitable habitat between the study taxa and Salix species may require a finer resolution than the 90-m pixel size of our SDMs, as the study taxa and Salix often coexist at scales of 10–100 s of meters within riparian forest mosaics varying in composition, successional stage, and hydrologic adaptation (Richter & Richter, 2000; Stromberg et al., 2007). In addition, annual streamflow, ditch length, and reservoir storage were likely less accurate than TPI and distance-to-water metrics and may have been poor indicators of the flow regime attributes that influence riparian tree occurrence (e.g., stream power, flood timing, magnitude; McShane et al., 2015).

Land use

The importance of land use in random- but not Salix-background models illustrates the challenges of discerning land-use drivers of species distributions when land use is associated with both study taxa and target-group background points. Like most study taxa, Salix background points occurred more frequently in developed areas than did random background points (mean dev_90-m = 10% ± 22% vs. 2% ± 11%). As a result, land use was unimportant to study taxa occurrence relative to Salix background points. Associations between land use and both Salix and study taxa could be causal or due to survey bias, because developed areas have likely been surveyed more frequently. However, survey bias is unlikely to be the full explanation for relationships between dev_90-m and U. pumila and E. angustifolia, because multiple studies with absence data also have reported positive relationships between human development and U. pumila and E. angustifolia (McShane et al., 2015; Perry et al., 2018; Reynolds et al., 2022; Ringold et al., 2008). Streambank and channel stabilization, soil and vegetation disturbance, and propagule pressure from upland plantings may facilitate U. pumila and E. angustifolia invasions in urban landscapes (Perry et al., 2018; Reynolds et al., 2022).

Climate

The importance of low annual precipitation for the invasive taxa (and not native taxa) in both random- and Salix-background models supports the hypothesis that drought-tolerance allows them to invade riparian habitats in drier climates (Glenn & Nagler, 2005; Perry et al., 2013; Reynolds & Cooper, 2010). Although moisture in riparian ecosystems is largely derived from streamflow fed by upstream precipitation, local precipitation may support riparian plants during low streamflows and on higher geomorphic surfaces. Lower local precipitation may reduce biotic resistance from mesic native plants and thus favor drought-tolerant invasive taxa. In particular, the consistent importance of low precipitation for T. ramosissima/chinensis is supported by previous large-scale analyses (Jarnevich et al., 2011; McShane et al., 2015; Perry et al., 2018) and suggests that adaptation to low precipitation contributes to T. ramosissima/chinensis invasion.

Maximum and/or minimum temperatures (maxT_{warm_mo}, minT_winter) were important for all taxa in Salix-background models, but with distinct roles for different taxa, even among closely related Populus. The importance of low minT_winter for E. angustifolia supports the hypothesis that unmet phenological chilling requirements restrict E. angustifolia's southern distribution (Collette & Pither, 2015; Friedman et al., 2005; Guilbault et al., 2012). The importance of high maxT_{warm_mo}, and low importance of minT_winter, for T. ramosissima/chinensis matches previous T. ramosissima/chinensis SDMs (Jarnevich et al., 2011; Kerns et al., 2009) and suggests that low warm-season temperatures, not just winter frosts, restrict T. ramosissima/chinensis's northern distribution (Friedman et al., 2005; Lesica & Miles, 2001; Sexton et al., 2002). Response curves for Populus species support the hypotheses that warm temperatures drive P. deltoides/fremontii's distribution (Ikeda et al., 2014) and moderate warm-season temperatures support P. angustifolia (Bothwell et al., 2021). They also suggest that unmet chilling requirements restrict P. angustifolia's distribution and moderate warm-season temperatures support P. balsamifera's distribution. Together, these differences resulted in distinct predicted climatic niches for the three Populus (Figure 5; Appendix S4: Figures S7 and S8). For U. pumila, the dominant, parabolic relationship with minT_winter suggests that U. pumila's distribution may be restricted by both insufficient cold tolerance in cold regions and unmet chilling requirements in warm regions. Cold intolerance in U. pumila would be surprising given its northern distribution relative to T. ramosissima/chinensis and E. angustifolia in Eurasia (Perry et al., 2018), but could result from U. pumila originally introduced to North America having been collected near the southern edge of its native distribution (Leopold, 1980). In accord with previous SDMs, these patterns suggest that predicted warming in the western USA may expand suitable habitat for T. ramosissima/chinensis (Bradley et al., 2009; Ikeda et al., 2014; Kerns et al., 2009) and P. deltoides/fremontii (Ikeda et al., 2014), while shifting suitable habitat northward for U. pumila, E. angustifolia, P. balsamifera, and P. angustifolia (Bothwell et al., 2021).

Management implications for invasive riparian trees

Habitat suitability predictions from our SDMs can help guide monitoring and management decisions for woody riparian invasion in different ecoregions and watersheds of the western USA. In particular, our models suggest that potentially widespread U. pumila invasion deserves greater scientific and management attention. U. pumila suitable habitat was similarly extensive to P. deltoides/fremontii and T. ramosissima/chinensis, which are the most common native and invasive riparian tree complexes across the western USA, and was more extensive than E. angustifolia, which is the second most common invasive woody riparian species across the study area (Friedman et al., 2005). Further, U. pumila could spread within Great Plains riparian ecosystems that have been of less concern for T. ramosissima/chinensis and E. angustifolia invasion, including southern Great Plains sub-ecoregions where T. ramosissima/chinensis and E. angustifolia had little predicted suitable habitat, central Great Plains sub-ecoregions where they had moderate suitable habitat but considerably less than U. pumila, and northeastern Great Plains sub-ecoregions where T. ramosissima/chinensis had little predicted habitat and E. angustifolia had considerably less than U. pumila.

Invasive species monitoring is particularly important for suitable habitat where the species has not yet been detected, because increased monitoring in these areas can enable early detection and rapid response to nascent invasions, identify overlooked invasions, and improve datasets for understanding invasive species distributions. We identified watersheds where the study taxa have been less sampled or have not yet invaded, by comparing occurrence record densities with what would be expected if occurrence records were distributed in proportion to habitat suitability. For U. pumila, this approach highlighted watersheds in the northwestern and central Great Plains that contained 20% of predicted suitable habitat yet contained no occurrence records. These areas could be surveyed to determine whether they contain well-established but unsampled populations, nascent populations, or suitable habitat at risk of invasion, as each circumstance would demand a different management approach.

For T. ramosissima/chinensis and E. angustifolia, which have been more intensively sampled than U. pumila, most watersheds with substantial suitable habitat contained at least some occurrence records. Watersheds with fewer-than-expected occurrences in portions of the Great Plains and North American Deserts may highlight populations that have not been fully documented or have greater capacity for expansion. However, especially for T. ramosissima/chinensis, which is thought to have reached its full potential range in the western USA (Friedman et al., 2005), watersheds with fewer-than-expected occurrences also might indicate areas that have been well-sampled but are still much less sampled than the most intensively monitored watersheds, such as near major human population centers.

It is also possible that some predicted suitable habitats with few occurrence records are not actually suitable. SDMs are useful tools for generating hypotheses about potential distributions based on available data, but they cannot predict distributions or causality with certainty. Our modeling choices allowed for at least acceptable levels of model credibility and utility for all criteria in the SDM assessment rubric developed by Sofaer et al. (2019) (Appendix S1: Table S12). In particular, absence data or targeted surveys of under-sampled predicted suitable habitat may improve the accuracy of future models.

ACKNOWLEDGMENTS

We thank Peder Engelstad for downloading and processing occurrence records, Tim Assal for coding TPI calculations, Pairsa Belamaric for processing environmental predictor data, and Bradley Butterfield for reviewing a previous version of the manuscript. This research was funded by the United States Geological Survey Biological Threats and Invasive Species Research Program. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government. L. G. Perry worked in cooperation with the U.S. Geological Survey Fort Collins Science Center.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Data (Perry et al., 2022) are available from the USGS ScienceBase-Catalog: https://doi.org/10.5066/P9LIB2TF.