Two invasive Hieracium species’ potential distributions within the Greater Yellowstone Ecosystem were defined using invasion susceptibility models and habitat typing

Invasive Hieracium plant species are invading the Greater Yellowstone Ecosystem. The potential distribution of orange hawkweed (Hieracium aurantiacum) and meadow hawkweed (Hieracium caespitosum) were estimated using habitat susceptibility models to assist land managers’ management of these invasive plants. The objectives of this study were to: (1) develop models describing susceptibility of ecosystems to hawkweed invasion, (2) identify indicator species of orange hawkweed and meadow hawkweed, (3) determine habitat types where these invasive hawkweeds might occur, and (4) create habitat susceptibility maps for management planning and ground surveys. Models were developed using a Mahalanobis distance similarity technique from remotely sensed biotic and abiotic variables, as well as known location data for orange and meadow hawkweed. Ground validation was conducted to assess model weaknesses and subsequent model modification. Indicator plant species were identified as surrogates to determine the likelihood of hawkweed presence during ground survey. Transect data collected from areas susceptible to invasion also were used to determine habitat types where hawkweed might occur. The best model included eight variables: north–south aspect, east–west aspect, slope, NDVI, NDWI, blue spectral band, green spectral band, and precipitation. High susceptibility (65 + % likelihood of suitable habitat) consisted of 66,000 ha for meadow hawkweed and 35,000 ha for orange hawkweed, 5.0% and 2.7% of the study area, respectively. Meadow hawkweed and orange hawkweed had seven and three indicator plant species, respectively. Predicted hawkweed habitat susceptibility encompassed nine habitat types, ranging from xeric sagebrush steppe to wet forests and they overlapped except at the xeric habitat type. Habitat susceptibility models save costs and allow survey prioritization to those areas most susceptible to invasion.


Introduction
Invasive plants, defined here as non-indigenous plants that change systems where they occur, are affecting areas of significant importance around the world as anthropogenic forces alter land use and transport seeds to new areas (Shiferaw et al. 2019;Szilassi et al. 2019;Pappgh et al. 2021). Invasive plants can reduce forage production, alter fire return intervals, reduce soil water content, and change native plant composition (Balch et al. 2013;Pfeiffer and Gorchov 2015;Shiferaw et al. 2019;Fenouillas et al. 2021;Papp et al. 2021). Land managers currently lack sufficient resources to reduce the ecological threat of invasive plants so it is imperative to prioritize areas with the greatest risk of invasion to receive focused management efforts.
Invasive plants cost billions of dollars each year from monitoring and control efforts, to higher food costs and reduced incomes on farms (Pyšek and Richardson 2010). Early detection of invasive plants, when infestations are small (< 1000 hectares), can mean the difference between efficient eradication and billions of dollars for ongoing control to stem biodiversity loss (Rejmánek and Pitcairn 2002;Shiferaw et al. 2019). Shiferaw et al. (2019) found the invasion of mesquite (Prosopis L.) in Ethiopia did not increase linearly when comparing Landsat data over time. They found that early on, the annual rate of spread of mesquite was over four-times slower compared to the spread in later years, highlighting the importance of early detection and early management strategies. Scientists are using remote sensing to better understand the distribution of invasive plants as they expand and move into more remote locations (e.g., Shiferaw et al. 2019;Masemola et al. 2020;Fenouillas et al. 2021;Papp et al. 2021;Pathak et al. 2021;Ghorbanian et al. 2022;Matas-Granados et al. 2022) because remote sensing reduces the cost associated with in situ sampling while providing greater spatial breadth and more frequent repeated measures than ground surveys could possibly achieve (Ahmed et al. 2020;Papp et al. 2021). Remote sensing has been essential in providing land managers with accurate data to manage invasive plants: Sentinel-2 multi-spectral data have been used to map Acacia dealbata (Link) and A. mearnsii (De Wild.) spread into grasslands and forests of South Africa (Masemola et al. 2020); Normalized Difference Vegetation Index (NDVI) was used to investigate long-term changes in land cover and identify vulnerability of threatened plants in Spain (Matas-Granados et al. 2022); and slope, aspect and elevation have been used to predict the likelihood of yellow star-thistle (Centaurea solstitialis L.) in the United States (Shafii et al. 2004).
Modeling has been used to predict habitat susceptibility to weed infestations using information of existing or previous infestations to prioritize survey areas (Shafii et al. 2003). Prioritizing efforts across expansive landscapes may be aided through use of Habitat Susceptibility Models (HSM) created within geographic information systems (GIS). GIS tools are used to combine spatial environmental data at known infestations and predict areas likely to be invaded based on similar ground conditions (Shafii et al. 2003(Shafii et al. , 2004Rew and Maxwell 2006). Areas that are similar to the known infestations of a target species are considered to be highly susceptible to invasion by that species, whereas areas that are different, are considered to have low susceptibility. HSM can efficiently direct ground survey efforts to locate invasive plant populations by identifying highly suitable habitat, saving land managers time and money.
The first step to developing an HSM is to identify abiotic and biotic factors (e.g., slope, aspect, elevation, vegetation indices) linked to the distribution of a target species (Lass et al. 2011;Adhikari et al. 2020). Remotely sensed digital elevation models can be used to interpret abiotic conditions such as aspect and slope, which have a direct impact on solar radiation and soil moisture levels (Holland and Steyn 1975). Some biotic conditions can be assessed using multispectral reflectance data (Lass et al. 2011;John et al. 2018;Ghorbanian et al. 2022). The reflectance bands encompass visible and near-infrared wavelengths; blue, green, red, and near-infrared (NIR) (central wavelengths: 490 nm, 560 nm, 665 nm, and 865 nm, respectively) often are used in ecological studies (Lass et al. 2011). The normalized difference vegetation index (NDVI, where NDVI = (NIR−Red) (NIR+Red) ) (Tucker 1979) is widely used because it differentiates between dense, living vegetation and sparse or senesced vegetation (e.g., Eitel et al. 2011;Lass et al. 2011;Ghorbanian et al. 2022). Recent studies have also incorporated the short-wave infrared band (SWIR, central wavelength = 1610 nm) to calculate additional indices such as the normalized difference water index (NDWI, where NDWI = (NIR−SWIR) (NIR+SWIR) ) which relates to plant water content (John et al. 2018). When models are unavailable, habitat types or indicator species may be used to aid on-the-ground surveys.
Habitat type is a term used to describe ecological characteristics of a location that act as a guide to distinguish between sites that can or cannot support certain plant species (Daubenmire 1984). This can be beneficial during field surveys for groundtruthing models or if maps cannot be loaded due to cellular service limitations or dead batteries in field tablets. Another option for determining habitat suitability when models are unavailable is to use an indicator species, a native or naturalized plant species with similar life-history characteristics of the target species, which may be considered an "indicator" of the potential presence of that target (Fleishman et al. 2005). An indicator species can inform managers on the susceptibility risk when the target species is absent or in low abundance.
In the United States as of 2017, 566,600 ha of National Parks were infested with invasive plants, with only 17,400 ha controlled (NPS 2019). Yellowstone National Park (YNP) is part of the Greater Yellowstone Ecosystem (GYE) which spans more than 8-million hectares across portions of Montana, Idaho, and Wyoming, USA. The GYE is a unique and intact ecosystem that is threatened by invasive species. It is characterized by a diverse flora with many vegetation types including: conifer forests, sagebrush steppes, and mountain meadows (Despain 1990). Invasive plants within the GYE challenge conservation goals, particularly meadow hawkweed (Hieracium caespitosum Dumort. = H. pratense Tausch) and orange hawkweed (H. aurantiacum L.). Meadow hawkweed and orange hawkweed were introduced into North America as ornamentals from Europe, as early as the nineteenth century (Wilson 2007). Hawkweed seeds are primarily wind-dispersed but can also be moved by animals. Once established, rhizomes and stolons allow plants to expand locally, at times excluding other plant species (Wilson and Callihan 1999).
Our study aimed to assist early detection of meadow hawkweed and orange hawkweed in the GYE by providing habitat susceptibility models for risk assessment and directing both plant survey and control activities for hawkweeds. The objectives were to: (1) develop susceptibility models for meadow hawkweed and orange hawkweed using known locations and environmental data within the GYE, (2) identify indicator species of meadow and orange hawkweed then assess the indicator power, (3) determine habitat types where meadow and orange hawkweed occur, and (4) create susceptibility maps that can be served on software platforms to allow invasive plant species mapping on mobile devices.

Spatial data acquisition
Ten environmental covariates were acquired for modeling: north-south aspect, east-west aspect, slope, hill shade, precipitation, blue and green spectral bands, NDVI, NDWI, and a moisture gradient index [(NDWI + 1) * (NDVI + 1)] (Fig. 2). Five spectral bands of Sentinel-2 data were downloaded from Earth Explorer USGS (European Space Agency 2018) (Fig. 2). All spectral data were analyzed at 20-m resolution. The NDVI and NDWI were calculated then combined to create a single variable to represent a moisture gradient of both. Each index was normalized by adding a value of one to remove negative values, then they were multiplied together.
Average annual precipitation data were acquired from Parameter-elevation Regressions on Independent Slopes Model (PRISM 2020) for 2012 through 2018 at 4-km spatial resolution. Data from each year were projected to 20-m resolution and then resampled using focal statistics with a 300 by 300-cell window to smooth the edges of 4-km data (TerrSet v. 18.31, Clark Labs, Worcester, MA; ArcMap v. 10.6.1, Esri, Redlands, CA). Average annual precipitation across the study area was 36 cm to 147 cm. A 30-m resolution digital elevation model (DEM) was acquired from Earth Explorer USGS (2018) and resampled to 20-m resolution using the nearest neighbor technique in ArcMap (Esri, Redlands, CA). Aspect, slope, and hillshade were derived from the DEM. In the northern hemisphere, south facing slopes receive more direct sunlight than north facing slopes and are typically drier. Because aspect is a circular variable, it was split into north-south (cosR) and east-west components (sinR) (Woodcock et al. 2008) where R = aspect in radians. Species presence data were obtained from the Greater Yellowstone Coordinating Committee (GYCC), Northern Rocky Mountain Exotic Plant Management Team, YNP Exotic Plant Management Team, and the US Forest Service. Locations were represented by polygons surrounding infestations. Ground surveys were conducted July 10-13, 2018 and July 22-26, 2019 to refine polygons and digitize new polygons using Collector (Esri, Redlands, CA). Location data were later converted from vectors to Labs, Worcester, MA) then approximately 60% of the known occurrence pixels were used as training data to develop the model (hereafter, training data) and 40% were used to validate the model (hereafter, validation data).

Model building
Susceptibility models were developed using Mahalanobis distance (D 2 ) among landscape variables using Mahalanobis Typicality in the Habitat and Biodiversity Modeler: Habitat Suitability/Species Distribution Modeling tool (Clark Labs 2017). Mahalanobis Typicality is a semi-supervised classifier that assigns probabilities to pixels indicating their similarity to the training data pixels. The D 2 statistic is the standardized squared distance between habitat covariates for a given sample of training data and the mean vector of the covariates (Griffin et al. 2010). A D 2 value was calculated for each pixel in the study area based on the value of the covariates within that pixel, relative to the covariate averages in the training data: where ̂ is the vector of the mean values for habitat covariates and Σ is the variance-covariance matrix for habitat covariates at training data locations.
Seventeen different combinations of the environmental variables described above were tested in the Habitat and Biodiversity Modeler: Habitat Suitability/Species Distribution Modeling tool in TerrSet for predictions of meadow hawkweed and orange hawkweed susceptibility (Table 1). Each model included north-south aspect, east-west aspect, and slope, along with one to five additional variables. Models were created at 20-m resolution and assessed based on a 9% error of omission rate for the training data. The error rate was achieved by selecting a threshold for D 2 values that maintained 91% of training pixels in the "susceptible" category. A conservative threshold was chosen in hopes of finding invasions that may occur towards the ecological limits of the hawkweeds' ability to persist. D 2 values were set to create high, moderate, and low susceptibility categories (Table 2).
Steps for building a susceptibility model. Bold rectangles represent raw data and dashed rectangles represent environmental data layers created from the raw data. Ovals indicate soft or hard classification steps. Words in (italics) indicate the TerrSet tool used. Different combinations of environmental layers were tested with the species training data to find the best model. The final model was chosen to minimize susceptible area while maintaining a 9% error of omission for validation data. Refer to Table 1 for specific environmental variables used in each model iteration Table 1 Variables for meadow hawkweed (MeH) and orange hawkweed (OrH) susceptibility model selection All models include two aspect variables (north-south, east-west) and slope plus the variables indicated by x. Sun angle diff. represents the difference in hillshading in August subtracted from hillshading in May. Blue and green bands are Sentinel-2 spectral bands. Precipitation represents the 7-year average from 2012 to 2018. Validation error represents the portions of meadow hawkweed or orange hawkweed validation pixels that were classified as "not susceptible" based on the Mahalanobis cutoff (D 2 ) values (0.070-0.115 for MeH models; 0.140-0.220 for OrH models) that allowed for 8.5-9.0% training error (percent of training pixels categorized as "not susceptible") NDVI Sun angle diff NDWI (NDVI + 1) * (NDWI + 1) Reducing model error in lodgepole pine stands Positive model predictions in unsuitable habitat were compared to determine possible areas for model improvement. Although no meadow or orange hawkweed populations were documented or observed within dense lodgepole pine (Pinus contorta Douglas ex Loudon) stands, the hawkweed susceptibility models predicted varying levels of susceptibility in those areas. Therefore, we developed a method to identify and remove some of these lodgepole pine locations.
Fourteen polygons of dense lodgepole pine were digitized in summer 2019 and determined to be within fire perimeters that occurred in 1988. Fire polygons from 1982 to 1994 were extracted from the dataset and considered as potential dense lodgepole pine. Fire polygons from 1995 to 2016 (Gibson 2005(Gibson [years 1995(Gibson -2002; LANDFIRE 2019 [years 2003-2016]) were removed from the potential dense lodgepole pine perimeters to ensure recently burned locations were not considered as they could still be suitable to hawkweed invasion.
To predict locations of dense lodgepole pine, Habitat and Biodiversity Modeler: Habitat Suitability/Species Distribution Modeling (Clark Labs 2017) was used across the study area, following methods described above for hawkweed models. Sixty percent of dense lodgepole pine pixels were used to train five models using different combinations of blue and green spectral bands, NDVI, NDWI, and precipitation ( Table 3). The output was a map that assigned each pixel a probability from 0 to 1 based on similarity to the training data. The previously described fire region was then selected from the full study area model. Model accuracy was determined by the proportion of dense lodgepole pine validation pixels properly classified in "predicted dense lodgepole pine", and the proportion of open forest, meadow hawkweed occurrence pixels, and orange hawkweed occurrence pixels classified as "not dense lodgepole pine". The dense lodgepole pine model was based on Mahalanobis Typicality values.
Final model processing steps A hydrologic features map was digitized in Arc-GIS 10.6.1 (Esri, Redlands, CA) then imported into TerrSet as vector data and converted to a raster at 40-m resolution to remove water bodies. Because most paved roads in the study area were less than 15-m across, it did not seem appropriate to convert roads into 40-m pixels and lose susceptibility predictions adjacent to roads.

Plant community surveys
Plant community surveys were conducted July 10-13, 2018 and July 22-26, 2019. Due to a late spring in 2019, seasonal variation in sampling times was minimal. Surveys were less than 1.2 km from a road, though usually within 0.5 km. Each transect was 20-m long with five quadrats (i.e., plots), sized

Indicator species analysis
An indicator species analysis was conducted using a chi-squared test based on the historical occurrence of meadow hawkweed or orange hawkweed against each species found along each transect using the vegan package in R v. 3.5.1 (R Core Team 2018; Jones et al. 2018). Species found in fewer than three transects were removed from the statistical analysis, but were considered for habitat typing, leaving 66 species for statistical testing. A Benjamini and Hochberg (1995) p value adjustment was used to control the false discovery rate. A quantitative measure of the indicator power (IP) of an indicator species for the target invasive species was calculated from a presence-absence matrix on species with unadjusted chi-squared p values ≤ 0.05. A species is considered a strong indicator when it frequently co-occurs with the target invasive species while also infrequently occurring in the absence of the target. The IP equation: where O I is the frequency of the indicator species (I) occurrence, O T is the frequency of the target invasive species (T) occurrence, S is the frequency of shared occurrences of the indicator and the target species, and N is the total number of transects surveyed (Halme et al. 2009).

Habitat typing
A habitat vector layer for ArcMap and vegetation guide, both created by Despain (1990), were used to identify habitat types at transect locations. Due to the large spatial scale of Despain's (1990) habitat types, habitat types were further refined using Steele et al. (1983) forest descriptions and Mattson (1984) meadow or open riparian descriptions. Transect photographs and species lists were compared to descriptions to confirm habitat types.

Susceptibility models
Model 17 containing 8 variables: north-south aspect, east-west aspect, slope, NDVI, NDWI, blue spectral band, green spectral band, and precipitation had low error rates for both hawkweed species validation data while minimizing the total area  Figs. 3, 4). Generally, areas that received less than 40 cm (e.g., along the Yellowstone River north of YNP) or more than 95 cm (e.g. south-western portion of YNP and high elevations) of annual precipitation were not susceptible to meadow or orange hawkweed. Incorporating both blue and green spectral bands in model 17 greatly reduced the area predicted to be susceptible, compared to other models (Table 1). Model output was divided into four susceptibility categories based on training data occurrence rates.
The D 2 values were set to correspond with training data occurrence rates of approximately 9%, 20%, 30%, and 40% for the zero, low, moderate, and high categories, respectively (Table 2).
Approximately 614,000 ha (47%) of the study area were predicted to be susceptible to meadow hawkweed. Of the susceptible area, 66,000 ha were in the high susceptibility category; 192,000 ha in moderate susceptibility; and 356,000 ha in low susceptibility (Table 2). Approximately 314,000 ha (24%) of the study area were predicted to be susceptible to orange hawkweed. Of the area susceptible to  (Table 2).
Overlaying susceptibility models indicated 670,400 ha (51%) of the study area were not susceptible to either hawkweed (Table 4). Overlapping areas of low susceptibility encompassed 80,900 ha (6%), overlapping moderate susceptibility encompassed 36,500 ha (3%), and overlapping high susceptibility encompassed 20,800 ha (2%) of the study area (Table 4). Approximately 53% (324,000 ha) of area susceptible to meadow hawkweed was not susceptible to orange hawkweed and 7% (20,600 ha) of area susceptible to orange hawkweed was not susceptible to meadow hawkweed.
Approximately 413,000 ha (32%) of the study area were considered dense lodgepole pine based on wildfires between 1982 and 1994. The best model for dense lodgepole pine was DLP Model 4, which used blue and green spectral bands, NDVI, and precipitation (Table 3). The area selected for removal (considered to be dense lodgepole pine) had a D 2 value  (Table 3). DLP Model 4 was selected because it had D 2 ≤ 0.30, no error for meadow hawkweed or orange hawkweed location data, and only 8% error for open forest. Dense lodgepole pine exclusion was 9170 ha and 8640 ha for the meadow and orange hawkweed models, respectively.

Indicator species from plant community surveys
Forty-five transects were surveyed for plant cover in habitats predicted to be susceptible to meadow and orange hawkweeds based on model predictions ( Table 5). The meadow hawkweed model predicted 13, 10, 22, and zero transects were in the high, moderate, low, and no susceptibility categories, respectively. The orange hawkweed model predicted 11, 8, 17, and 9 transects were in the high, moderate, low, and no susceptibility categories, respectively.
A total of 209 plant species were recorded within the 45 transects, of which we identified 143 (Appendix A). Three genera were identified to genus only: strawberry (Fragaria L.), sedge (Carex L.), and horsetail (Equisetum L.). The chi-square contingency test revealed seven species were indicator species (p < 0.05) for meadow hawkweed and three species were indicator species for orange hawkweed ( Table 6). The unadjusted p values were used for selecting species for the indicator power analysis.
Two of the species significantly associated with orange hawkweed were also significantly associated with meadow hawkweed: Richard's geranium (Geranium richardsonii Fisch. & Trautv.) and fringed willowherb (Epilobium ciliatum Raf.) (Table 6). Species significantly associated only with meadow hawkweed were common selfheal (Prunella vulgaris L.), arrowleaf ragwort (Senecio triangularis Hook.), falsegold groundsel (Packera pseudaurea var. pseudaurea (Rydb.) W.A. Weber & Á. Löve), starry false lily of the valley (Maianthemum stellatum (L.) Link), horsetail (Equisetum spp.), sedge (Carex spp.) and ballhead ragwort (Senecio sphaerocephalus Greene). Alpine timothy (Phleum alpinum L.) was only significantly associated with orange hawkweed. Seven species significantly associated with meadow hawkweed had an IP of 0.998, indicating they frequently occurred with meadow hawkweed while also infrequently occurring in the absence meadow hawkweed ( Table 6). The strongest indicators of orange hawkweed were fringed willowherb and Richard's geranium (Table 6). Alpine timothy had moderate indicator power for orange hawkweed with an IP of 0.706.  Species that were important for habitat typing but either occurred in fewer than three transects or were not significant in the chi-square analysis, are listed in Table 7.

Habitat types
Eight habitat types were identified within predicted orange hawkweed-susceptible habitat: big sagebrush/Idaho fescue, sticky purple geranium phase; big sagebrush/bluebunch wheatgrass; Idaho fescue/bearded wheatgrass, sticky purple geranium phase; alpine timothy/water sedge; lodgepole pine/ Ross' sedge; subalpine fir/grouse whortleberry, grouse whortleberry phase; subalpine fir/bluejoint reedgrass; and subalpine fir/western meadowrue. All eight habitat types, plus big sagebrush/ Idaho fescue, were also identified within predicted meadow hawkweed-susceptible habitat (See Table 8 for habitat type descriptions; see Fig. 5 for habitat type examples). Habitat types ranged from xeric shrubland/grassland to wet forest. Historic infestations of orange and meadow hawkweed were found in fewer habitat types. Historic orange hawkweed infestations encompassed only four habitat types: Idaho Fescue/bearded wheatgrass, sticky purple geranium phase; subalpine fir/bluejoint reedgrass; alpine timothy/water sedge; and Lodgepole pine/ Ross' sedge. Historic infestations of meadow hawkweed were found in the same four habitat types as orange hawkweed infestations, plus two additional habitat types: big sagebrush/Idaho fescue and subalpine fir/grouse whortleberry.

Discussion
Remote sensing provides a cost-efficient method to monitor vegetation cover across a broad spatial extent and through time, something in situ sampling simply cannot achieve (Ahmed et al. 2020;Papp et al. 2021). Predictive modeling can help direct ground surveys to areas that are most likely to be invaded by a target species so new infestations can be treated before becoming large and unmanageable. Weed infestations can be a multi-billion dollar expense annually for large countries (Pyšek and Richardson 2010) so early detection and swift action are key.
In this study, we developed habitat susceptibility models for meadow hawkweed and orange hawkweed using the D 2 method with known locations and environmental data, determined indicator species, and classified plant communities at transect Habitat type descriptions from Despain (1990). Additional descriptions were included form Steele et al. (1983) and Mattson (1984). The elevation ranges described by Despain (1990) are included however, his GIS layer had some discrepancies. For example, the book listed big sagebrush/Idaho fescue occurring at 6800-9500 ft but a portion of the map labeled as big sagebrush/Idaho fescue was at 5400 ft locations into habitat types to assist land managers in their efforts to monitor and control these invasive plants in the GYE. Our models were based on presence-only data to avoid assumptions that a species' absence from an area indicated no susceptibility. A challenge of using the D 2 method to develop predictive habitat susceptibility models is the uncertainty if the invader is absent because the habitat is not suitable or simply because propagules have not yet arrived. It takes time for weed propagules to establish and occur in great enough abundance to be detected in a new location; the absence of the target species is not enough to confirm the location is not susceptible to invasion. Another limitation of the D 2 statistic is the lack of a significance test or established method to ascertain importance of the model variables (Griffin et al. 2010), in our case, which is why several models were compared. Once the best model was selected, the model prediction error rate was improved by excluding areas where neither hawkweed species were found within a habitat type. Dense lodgepole pine stands are not suitable habitat for meadow hawkweed or orange hawkweed, so refining models without dense lodgepole pine significantly reduced the amount of land where monitoring might otherwise have been considered. Even with these limitations, D 2 was the best method for presence-only data available for this study (Phillips et al. 2006). Our meadow hawkweed model indicated 20% of the study area was in the combined moderate and high categories and 27% was in the low category. Combined moderate and high categories for orange hawkweed susceptibility represent only 9% of the study area and the low category encompassed 15%. Prediction of habitat susceptibility for meadow hawkweed was nearly twice that of orange hawkweed which, is similar to the ratio of current mapped infestations of both species. The overlap of meadow and orange hawkweed susceptibility models suggests these species can co-occur across a wide range of habitats yet, portions of their ranges are unique. Survey crews should prioritize efforts along propagule corridors in areas predicted to be highly susceptible to both species and near known populations. Roads, streams, and trails allow propagule movement into new areas so, areas along these corridors should be scoured to maximize survey efficiency.
Outside of the region where the predictive model was developed, surrogate-based approaches may be used to estimate habitat susceptibility based on site conditions that can be evaluated on the ground, such as current habitat conditions or indicator species presence (Halme et al. 2009;Thuiller et al. 2012;Jones et al. 2018). Predicted meadow hawkweed susceptibility encompassed nine habitat types, and predicted orange hawkweed susceptibility encompassed eight of the same habitat types, excluding big sagebrush/ Idaho fescue (Table 8). However, historic hawkweed infestations were only recorded in six of these habitat types. The three newly identified potential habitat types: big sagebrush/Idaho fescue, sticky purple geranium phase; big sagebrush/bluebunch wheatgrass; and subalpine fir/western meadow-rue suggest more of GYE is susceptible to meadow and orange hawkweed than previously thought. These habitat types should be prioritized for future monitoring. Invasive hawkweeds are not typically found in dry shrub-steppe grasslands (Wilson 2007); however, meadow hawkweed was observed in a big sagebrush/ Idaho fescue habitat type. Orange hawkweed may be found in sites with higher moisture than described in this study as it is known to infest bogs and cranberry fields (Jorgensen and Nauman 1994). Wetter communities like willow and sedge dominated marshes were not surveyed in this study.
Indicator species are another surrogate-based approach we addressed in this study as training new field crews to identify key species may be more realistic than training them to identify broad habitat types. Alpine timothy was the only unique indicator of orange hawkweed. Thus, survey crews could mark locations of this plant to identify likely areas for invasion. Richard's geranium and fringed willowherb were strong indicators for both hawkweeds. If crews encounter those species they should methodically search the area.
Meadow hawkweed and orange hawkweed have the potential to impact a large portion of the GYE across a range of moisture conditions and habitat types. Early detection is vital because patch size of invaded areas grow at faster rates later in their establishment, reducing the likelihood of control (Rejmánek and Pitcairn 2002;Shiferaw et al. 2019). Prioritizing monitoring in habitats suitable for an invasive species can reduce the overall area that is surveyed and ultimately save land managers time and money on ground surveys and management efforts, while preserving the GYE in its most natural state. Using the susceptibility models to direct ground surveys within the predicted habitat types that have not yet had documented invasions could be instrumental in locating new infestations, avoiding further expansion, and refining models for these two invasive hawkweed species. Susceptibility models should be updated when new locations of invasive plants are discovered and when updated GIS data are available. When models are not available for an area, understanding the habitat types or indicator species associated with the species of interest can guide survey efforts. While our study focused on plant invasion susceptibility, the same process can be used for any species of concern.
Fig. 5 Meadow and orange hawkweed habitat types. a basin big sagebrush/Idaho fescue, b basin big sagebrush/Idaho fescue, sticky purple geranium phase, c Lodgepole pine/sedges, d subalpine fir/bluejoint grass, a is only susceptible to meadow hawkweed, b-d are susceptible to meadow hawkweed and orange hawkweed