Assessing Hybridization Risk Between ESA-Listed Native Bull Trout (Salvelinus confluentus) and Introduced Brook Trout (S. fontinalis) Using Habitat Modeling

The introduction of non-native species can negatively impact native species through reduced genetic fitness resulting from hybridization. The lack of spatiotemporal data on hybrid occurrences makes hybridization risk assessment difficult. Here, we developed a spatially-explicit Hybridization Risk Model (HRM) between native Oregon bull trout, an Endangered Species Act-listed Oregon species, and introduced brook trout by combining an intrinsic potential model (IPM) of brook trout spawning habitat and existing bull trout distribution and habitat use datasets in Oregon, United States. We created an expert-based brook trout IPM classification score (0–1) of streams based on the potential of geophysical attributes (i.e., temperature, discharge, gradient, and valley confinement) to sustain spawning habitats. The HRM included a risk matrix based on the presence/absence of both species as well as the type of habitat (spawning versus other) at 100-m stream segment resolution. We defined the hybridization risk as “extreme” when stream reaches contained bull trout spawning habitat and brook trout were present with IPM moderate or greater scores (IPM >0.5). Conversely, “low” risk reaches contained historic or non-spawning bull trout habitat, brook trout were absent, and IPM scores were low (IPM <0.25). Our HRM classified 34 km of streams with extreme risk of hybridization, 115 km with high risk, 178 km with moderate risk, and 6,023 km with low risk. Our HRM can identify a differential risk of hybridization at multiple spatial scales when either both species coexist in bull trout spawning habitat or are absent. The model can also identify stream reaches that would have higher risk of hybridization, but where brook trout are not currently present. Our modeling approach can be applied to other species, such as cutthroat trout and rainbow trout, Chinook and coho salmon, or similar species occurring elsewhere that potentially hybridize in freshwaters.


INTRODUCTION
The introduction of non-native or invasive species poses several challenges to the conservation of native species. The main negative effects of species introductions on native species include predation (Brown and Moyle 1981;Ingeman 2017), competition for resources and habitat (Gurnell et al., 2004), and exposure to diseases carried by the introduced species (Chantrey et al., 2014;Young et al., 2017). There is also the risk of hybridization between native and non-native species, which has received less attention among researchers and natural resource managers (Dunham et al., 2004;Fausch et al., 2006;Gozlan et al., 2010;Al-Chokhachy et al., 2014;Kilshaw et al., 2016). Hybridization occurs naturally in wild populations as part of the evolutionary process of many taxa (Arnold 1997). Nonetheless, when hybridization occurs between native and nonnative species it can lead to reduced genetic fitness and/or genetic extinction of native species (Rhymer and Simberloff 1996). For example, the tiger salamander (Ambystoma tigrinum), a nonnative species introduced in central California in the 1950s for use as bait in recreational fishing (Johnson et al., 2011), now interbreeds with the rare native California tiger salamander (A. californiense), threatening its persistence (Riley et al., 2003). In the study by Riley et al. (2003), hybrid swarms were found in three of the six study sites and in only one of the remaining study sites did native California tiger salamanders comprise more than 8% of the population. In many western North American rivers, introgressive hybridization between salmonids (Bettles et al., 2005;Muhlfeld et al., 2009;Araujo et al., 2021) such as introduced rainbow trout (Oncorhynchus mykiss) and native cutthroat trout (O. clarkii lewisi and O. clarki clarki) has reduced the reproductive fitness of the latter (Bettles et al., 2005;Muhlfeld et al., 2009).
Salmonids are among the most introduced freshwater fish species worldwide (Dunham et al., 2004;Lowe et al., 2004;Arismendi et al., 2019). As transportation networks improved, greater access to more remote watersheds has made it easier for promoters of stocking programs to introduce salmonids in areas outside their native ranges (Krueger, 1991). The introduction of brook trout (Salvelinus fontinalis), a.k.a. brook char, native to eastern North America, has occurred in western North American rivers since at least the early 1900s (Brook trout Oregon Department of Fish and Wildlife, n.d). This introduction has increased the risk of hybridization with native bull trout (S. confluentus), a.k.a. bull char (Kanda et al., 1994(Kanda et al., , 2002Fausch et al., 2006;DeHaan et al., 2010;Howell 2018), a species that once was common throughout most major rivers in the Pacific Northwest of North America. By the 1990s, habitat loss and fragmentation resulted in local extirpations and significant range contractions of bull trout (Goetz 1994). As a result, bull trout have been listed as threatened under the U.S. Endangered Species Act since 1992 (U.S Fish and Wildlife Service, n.d), and now face additional threats from introduced brook trout, including the potential for introgressive hybridization. The combination of external fertilization and weak reproductive barriers are seen as significant contributors to the risk of hybridization (McKelvey et al., 2016).
Introgressive hybridization occurs when first-generation (F 1 ) progeny backcross with the parent species and transfer genetic material from one species to another. This has been documented between native bull trout and introduced brook trout in western Montana (Kanda et al., 2002). Kanda et al. (2002) found 30% of all specimens they collected were hybrids. In addition, 75% of the collected hybrids were first-generation (F 1 ) males whereas the rest were the result of backcrosses with the parent species. Further, most hybridization was observed to occur in one direction, involving bull trout females and brook trout males, representing wasted reproductive effort on the part of bull trout. While an understanding of hybridization risk is needed for its effective management, one of the many challenges faced by natural resource managers is the inability to quantify or categorize hybridization risk.
The ability to quantify hybridization risk at the stream reach level would be a useful tool when planning restoration projects or prioritizing conservation programs to remove non-native species. A simple method of quantifying risk in the context of emergency management is the Hazard/Risk Matrix (Pritchard 2010), where the probability of an incident occurring is combined with the potential impact if the event were to occur. Risk then becomes a function of Probability multiplied by Impact (NETC Emergency Management Institute 1998). This approach can be adapted to quantify the risk of hybridization between native and introduced species. If Impact is measured as the amount of spawning habitat used concurrently by introduced and native species, then the risk of hybridization can be defined as a function of the probability of simultaneous use of a given stream reach for spawning by the two species. For example, if brook trout and bull trout were sympatric in historical bull trout spawning habitat during the spawning season, there would be a greater risk of hybridization than if sympatric use only occurs in habitat used for migration or at some other time of the year. Unfortunately, there is a lack of information regarding seasonal variations on habitat use in relation to species life stages. This lack of related spatial and temporal distribution data has been identified as a gap that needs to be addressed in future work concerned with the structure and conservation of bull trout populations (Dunham and Rieman 1999). Filling this gap would require extensive sampling efforts over several years, an approach that would be labor-intensive and cost prohibitive. However, modeling methods of species distributions based on habitat intrinsic potential use (Burnett et al., 2007) may help address this gap and the quantification of hybridization risk between native and introduced species over large spatial extents in way that can be tested or validated.
This study was conceived with the following objectives: 1) develop a hybridization risk model based on the intrinsic potential of stream reaches to meet brook trout spawning habitat requirements, combined with existing distribution and habitat use data for bull trout, and 2) apply that model to Oregon bull trout habitat to quantify hybridization risk at different spatial scales. The Hybridization Risk Model (HRM) developed in this study combines the Intrinsic Potential Model (IPM) to identify potential brook trout spawning habitat, with empirical bull trout distribution and habitat use data in the state of Oregon, United States. The intrinsic potential modeling approach focuses on topographic attributes that influence stream habitat. This approach identifies stream reaches that have the greatest potential for spawning use based on life-stage, species-specific habitat needs. For this study, hybridization risk is then the product of the intrinsic potential of stream reach to be used for brook trout spawning, and the presence of bull trout and how bull trout use that same habitat.

Target Species
Bull trout are a large, long-lived char native to western North America that are typically associated with clean, clear, cold water, and complex stream habitats that include large woody debris, riffles, and deep pools. In general, they prefer stable year-round stream flows, substrate with minimal fines, deep pools, and cover in the form of undercut banks and canopy (Rieman and Mclntyre 1993;Al-Chokhachy et al., 2010). Bull trout require several kilometers of interconnected river reaches with cold-water temperatures, preferably <12°C (Budy et al., 2019). Bull trout are considered one of the most thermally sensitive of salmonids (Dunham et al., 2003), restricting their distribution to colder streams with an upper-temperature threshold of 15°C (Rieman and Mclntyre 1993). Historically, this species has exhibited multiple life-history strategies. For example, after spending one to 3 years in natal or rearing habitats some individuals migrate downstream to larger rivers or lakes, not returning to spawn for another five to 7 years. Other individuals may remain in or near natal habitat for their entire lives . Bull trout reach maturity between four and 7 years of age, and spawn between late August through late November, depending on stream temperature conditions. Stream temperatures of <9°C are required for spawning, which typically limits spawning habitats to higher elevation headwaters (Ardren et al., 2011;Austin et al., 2018). Spawning occurs in reaches with some form of cover, primarily gravel substrate (Boag and Hvenegaard 1994) with minimal amounts of fine sediments, average depths of about 28 cm, and flows ranging from 1.0 to 39.0 cm s −1 (Guzevich and Thurow 2017).
Brook trout, a species of char native to eastern North America, have been introduced in many areas outside their native range since the late 1800s, including the Pacific Northwest (DeHaan et al., 2010). They mature at about 2 years, can tolerate a wider range of stream temperature conditions than bull trout, and have an uppertemperature threshold of approximately 20°C (Raleigh 1982). When both species are present, this results in a distribution pattern where bull trout occupy colder, higher elevation streams and brook trout occupy warmer, lower elevation streams, with an area of overlap between the two species . Brook trout typically spawn between September and November when stream temperatures are between 4.5 and 10°C, in substrate with little amount of fine sediments, and where there is upwelling groundwater (Raleigh 1982).

Study Area
The majority (55.34%) of all bull trout habitat types in Oregon are found in the northeastern part of the state within the John Day and Powder/Burnt River Basins, with additional bull trout habitats located in the Malheur/Owyhee, Klamath, Deschutes, To include all known Oregon bull trout habitat, our study area was selected at the eighth field hydrologic unit scale (HUC8), also known as a sub-basin, as defined by the U.S. Geological Survey (USGS Water Resources, n.d) (Supplementary Table S1 in Appendix). Sub-basins have an average area of 1,813 square kilometers. The sub-basins selected for the study included both current and historic habitat used by native bull trout populations, as well as habitat currently used by introduced brook trout. Habitat use is determined with the fhdUseType attribute in the Oregon Department of Fish and Wildlife (ODFW) 2021 Oregon Fish Habitat Distribution geodatabase. Data in the geodatabase is based on sampling, professional opinion, or modeling by staff biologist with ODFW and other natural resource agencies. Historical habitat use is defined as "habitat used, or potentially used, historically, but not currently." Current habitat use includes habitat used within the last five reproductive cycles for spawning, rearing, migration, foraging, overwintering, and resident multiple use.
Habitat distribution data for both species, as geographic information system (GIS) layers, data were obtained from the ODFW's Oregon Fish Habitat Distribution (FHD) dataset (Bowers 2019). The FHD dataset contained 818 records of bull trout distribution comprised of current life-stage specific habitat use information as well as areas of historical habitat use. There were 353 records of brook trout distribution that included known populations within the state, but did not include information about life-stage specific habitat use. Of the 47 sub-basins included in the study: 11 had only bull trout, 16 had only brook trout, and 20 had both species (Supplementary Table S2 in Appendix). The HRM was built using bull trout habitat within Oregon only, ODFW distribution data included bull trout habitat in three subbasins that straddle state lines between Oregon, Washington, and Idaho and contained streams in both Washington and Idaho. These were the Lower Grande Ronde, Walla Walla, and Hells Canyon sub-basins. Streams in these sub-basins that were outside Oregon were not included in the study.
The ODFW Fish Habitat Distribution Dataset was built on the National Hydrography Dataset (NHD) and included "habitat codes" that corresponded to bull trout habitat use by life-stage, such as spawning, rearing, migration, foraging, and overwintering. Because the distribution data were built on the NHD, we incorporated other NHD-based environmental data, such as flow and temperature, in the creation of GIS layers to develop stream profiles, which included brook trout spawning habitat intrinsic potential in relation to bull trout habitat and lifestage use.

Intrinsic Potential Model for Brook Trout Spawning Habitat
The IPM approach uses geomorphic and hydrologic characteristics to classify the "intrinsic potential" of a given stream reach to define suitable habitat requirements based on species-specific life-stage habitat needs. Physical characteristics such as stream gradient, discharge, and valley constraint have been used to model these interactions and derive intrinsic potential for habitat use in salmonids at the reach level (Agrawal et al., 2005;Burnett et al., 2007;Busch et al., 2013). When considering hybridization risk between the two thermally sensitive species in this study, temperature conditions over space are added to the model to account for species-specific thermal limitations to habitat use. The intrinsic potential of a stream reach is then modeled using landscape variables that influence channel characteristics. Burnett et al. (2007) developed a model for coho salmon using the geometric mean of channel gradient, valley width index, and mean annual discharge. It has been used to assess stream habitats relative to land use and ownership and has been useful in identifying salmonid habitat in need of restoration. To the best of our knowledge, it has not been used for assessing interactions between multiple species, nor to assess hybridization risk.
The IPM for brook trout spawning habitat was developed in a two-step process. First, a two-part survey of species experts was used to construct brook trout spawning habitat suitability indices based on temperature, discharge, gradient, and valley width index. Species-experts were identified through a review of published brook trout research and from the International Charr Symposium participant list (Muir 2020). In the first survey, respondents were asked to identify threshold values for four specific habitat variables in relation to brook trout spawning requirements: minimum/maximum stream temperature, maximum gradient, minimum/maximum discharge, and minimum/maximum valley width index (See survey details in Supplementary Appendix C). Results were then used to plot suitability indices based on the mean response values for each variable. The suitability index scores use a zero value to represent unsuitable habitat and a value of one to represent optimum values for each variable. In the second survey, sent to the same general group of experts, respondents were asked to select the appropriate value for four inflection points associated with temperature, gradient, discharge, and valley width index. The inflection points represented the absolute minimum, minimum optimal, maximum optimal, and absolute maximum values of each variable in relation to brook trout spawning use. Final suitability curves were then developed based on the response means for each of the inflection points for each variable.
The second step in the IPM process was applying the suitability indices for the classification of environmental data using GIS. Spatial data can be represented in GIS so that specific features, such as species distributions and habitat types, can be analyzed in both space and time (Fisher and Rahel 2004). The attributes of different feature classes representing the same physical location can be combined to model the interaction of different characteristics associated with that location. Interactions between riverscapes, defined as the "... spatially heterogeneous scene of the river environment, ... unfolding through time." (Fausch et al., 2002), and organisms have both spatial and temporal components that are complex and have the potential to strongly affect freshwater salmonid population dynamics (Neville et al., 2006).
A spatially explicit geohydrologic dataset was constructed in ESRI's ArcGIS Pro (2.7) program (ESRI 2020). All bull trout habitats from the ODFW dataset were used to define the extent of  (Isaak et al., 2017). The historical stream temperature models from 1993-2015 were used for IPM model input.
For each variable in the IPM model (temperature, discharge, gradient, valley width index), values below the Absolute Minimum (AbsMin) or above the Absolute Maximum (AbsMax), had an index score of zero. Values that were greater than or equal to the Minimum Optimal (MinOpt) and less than the Maximum Optimal (MaxOpt) were given a score of one. For simplicity, we assumed linear relationships between the AbsMin and MinOpt breakpoints, and between MaxOpt and AbsMax. We then used the geometric mean of the variables as the IPM final score for that reach. Final IPM scores of less than 0.25 were assigned as "Low intrinsic potential" for brook trout spawning. Scores between 0.25-0.50 mean were assigned as "Moderate intrinsic potential". Scores between 0.50-0.75 represented "High intrinsic potential", and scores between 0.75-1.0 indicated "Very High intrinsic potential".

Hybridization Risk Model
The hybridization risk model (HRM) was developed by combining 1) the results of a novel Intrinsic Potential Model (IPM) of brook trout spawning habitat requirements and empirical distribution data, with 2) Oregon bull trout empirical distribution and specific life-stage habitat use data (Figure 2). Hybridization risk was quantified as the overlap between the IPM score of reaches suitable for brook trout spawning and bull trout habitat use by life stage (Table 1). This approach allowed the evaluation of both current risk and the potential risk of hybridization based on the existing and potential presence of brook trout. The highest hybridization risk was identified when 1) a stream reach had brook trout present, 2) the IPM score was moderate or greater (see section below for details), and 3) the stream reach overlapped with bull trout spawning habitat. Conversely, a stream reach with the lowest hybridization risk occurred when 1) brook trout were absent, 2) the IPM score was classified as low, and 3) the stream reach only overlapped historic or non-spawning bull trout habitat ( Table 1).

Estimating Risk of Hybridization Between Brook Trout and Bull Trout
The results of the IPM and the species distribution/habitat data were then used to calculate the hybridization potential at two spatial extents: 1) for all bull trout habitat in the study area and 2) only for bull trout spawning habitat. Using ArcGIS, a polyline layer was created that mapped Oregon bull trout habitat in 100 m reaches. NHDPlus HR flow lines were used to calculate the total length of the stream network of each HUC8 sub-basin corresponding to each spatial extent. The total length of bull trout habitat, brook trout habitat, sympatric habitat, and bull trout spawning habitat for each sub-basin was then calculated using the stream network in combination with ODFW species/habitat use distribution data. Bull trout spawning habitat was defined as reaches within the ODFW dataset that had habitat use codes of "Primarily spawning with some rearing" or "Resident species only, multiple uses".
Hybridization risk was quantified by multiplying the IPM score for each stream reach by a risk category multiplier based on the presence or absence of brook trout in that reach and bull trout habitat use ( Table 1). Bull trout spawning habitats belonged to the highest risk category (multiplier = 1.0) if brook trout were present, and to the second highest (multiplier = 0.75) if brook trout were absent. All other bull trout habitat had risk category multipliers of 0.5 when brook trout were present, and 0.25 if brook trout were absent. We then created hybridization risk maps at the state, sub-basin (HUC8), and selected watershed (HUC10) spatial scales in the study area (Supplementary Appendix D). Oregon has limited information available for use in model validation using actual hybrid occurrences. Therefore, in the discussion section, we contrast model output with two published studies and data collected over several years as part of a brook trout removal effort in Oregon.

Expert Opinion Surveys
The first survey (Supplementary Appendix C) was sent to 220 researchers of which 49 responded, representing a 22.3% response rate. The response means for the lower and upper stream temperature limiting brook trout spawning were 4.7 and 13.8°C, respectively. The response mean for maximum gradient was 6.2%. The response means for monthly minimum and maximum discharge during the spawning season were 1.4 m 3 /s and 17.3 m 3 /s, respectively. The means for minimum and maximum valley width index (VWI) were 2.6 and 13.1 respectively.
The responses from the first survey were used to develop habitat suitability curves that were tested during the second survey. There were several researchers from the first survey that did not participate in further surveys, along with the inclusion of additional researchers recommended by respondents to the first survey, which resulted in the second survey (Supplementary Appendix C) being sent to 210 researchers, of which 50 responded, representing a 23.8% response rate.
Habitat suitability indices (Figure 3) revealed the geographic distribution of preferred spawning habitats for brook trout based on interactions among stream temperature, stream gradient, mean monthly discharge, and valley confinement (measured as valley width index). Stream reaches with mean monthly temperatures between 3.16 and 13.10°C had potential for spawning use, with reaches having temperature conditions between 5.5 and 10.3°C classified as of greatest potential. Stream gradient ranged from 0.88% (Rosgen G stream type) to 7.55% (Rosgen A stream type), with optimal spawning habitat between 1.89% (Rosgen C) and 5.38% (Rosgen A). Mean monthly discharge rates, measured in cubic meters per second, ranged from 0.77 m 3 /s to 12.07 m 3 /s, with optimal discharge rates between 1.86 and 7.99 m 3 /s. Valley confinement, measured as valley width index, ranged from 1.37 to 12.69, with 2.72-9.45 being the optimal values ( Table 2).

Bull Trout and Brook Trout Habitat Distributions
Bull trout habitat occurred in 31 of the 47 sub-basins in the study area (Supplementary Table S2 in Appendix). Brook trout habitat existed in 36 of the 47 sub-basins, and habitat with both species   (Figure 4).

Model Results
The corresponding IPM scores were calculated at the 100-m reach scale and classified as Low (<0.25), Moderate (0.25-0.50), High (0.50-0.75), and Very High (>0.75) intrinsic potential for brook trout spawning use ( Figure 5). Bull trout habitat was found in three broader areas within Oregon including the northern Cascades ( Figure 5A), the northeastern part of the state in the Blue Mountains ( Figure 5B), and the southwestern area near Sprague River ( Figure 5C). Major areas of overlap between bull trout habitat and brook trout spawning habitat were in the northeastern region-279 km ( Figure 5B), followed by the northern Cascades-258 km ( Figure 5A), and then the Sprague River area-177 km ( Figure 5C).
Model results for brook trout spawning habitat in bull trout spawning habitat (Figure 6) showed the greatest concentration of Moderate (0.25-0.50) to Very High (>0.75) values in the northeastern part of the state-257 km ( Figure 6B), followed by the northern Cascades region-22 km ( Figure 6A), and then the Sprague River area-7 km ( Figure 6C). There were 466 km of bull trout spawning habitat with Low (<0.25) brook trout spawning IPM in the northern Cascades region, 1,265 km in the northeastern part of the state, and 57 km in the Sprague River area.
Moderate or higher brook trout IPM scores were found in 286 km (13%) of bull trout spawning habitat (Figure 7). At the HUC8 sub-basin level, the overlap of Moderate to Very High brook trout spawning IPM with bull trout spawning habitat ranged from  no overlap of the 142 km of bull trout spawning habitat in the Imnaha sub-basin, to 34% of the 119 km of bull trout spawning habitat in the Upper John Day sub-basin (Figure 7). Other areas with an overlap of 20% or more of bull trout spawning habitat by Moderate to Very High brook trout spawning IPM were the Upper Malheur-45 km (20%), Powder-19 km (24%), Middle Fork John Day 21 km (26%), and North Fork John Day 71 km (28%). Stream reaches with high IPM scores that overlapped bull trout spawning habitat had Extreme hybridization risk (Table 1). There were 344 reaches covering 33.9 km categorized as having Extreme hybridization risk. Most of these reaches were in the North Fork John Day sub-basin (51%), followed by the Upper Malheur sub-basin (19.8%), and Middle Fork John Day and Powder sub-basins (12.3% each). Conversely, bull trout spawning habitat where brook trout were not currently present were given a risk value of 0.75 and considered to have High hybridization risk based on the potential for hybridization occurrence if brook trout became established in the future. There were 1,190 reaches in the study area, covering 115.3 km, categorized as having High hybridization risk. Most of these stream reaches were in the Upper Malheur sub-basin (22.2%), followed by the Upper John Day sub-basin (21.1%), and the North Fork John Day sub-basin (16.1%). Other bull trout habitat, such as overwintering or foraging habitat, were given a moderate risk value of 0.5 if brook trout were present, and a low-risk value of 0.25 if brook trout were absent. There were 1,819 reaches covering 177.5 km that were categorized as Moderate risk, and 60,977 reaches covering 6,023.2 km that were categorized as Low risk. Most of these stream reaches were concentrated in the northeastern part of the state in the John Day basin (21.8%), followed by the Lower Snake basin (20.1%), and the Deschutes basin (15.6%) on the eastern slope of the Cascades Range. An example of a risk map at multiple scales is illustrated in Figure 8 and Supplementary Appendix D.   (Bowers 2019). The value of the HRM approach lies in its ability to identify stream reaches with a differential risk of hybridization depending on whether brook trout potentially co-occur in bull trout spawning habitat or are currently absent. The HRM model provides risk values from the 100-m reach scale up to the HUC8 sub-basin, giving managers a high-level overview for use in discussions with policymakers, and a detailed view of reach-level risk for planning future fieldwork or specific management actions. For example, the HRM can be used at multiple spatial scales to identify individual reaches that might be candidates for brook trout removal, to model potential hybridization risk if bull trout are reintroduced, or to incorporate scenarios of brook trout introductions or expansion. Moreover, the HRM has the potential to incorporate different climate model scenarios to assess hybridization risk under different hydroclimatic conditions. Furthermore, to build the HRM, we created an IPM for brook trout that provides a relatively simple procedure to identify potentially suitable brook trout spawning habitats using publicly available data. This procedure is based on the IPM developed by Burnett et al. (2003), which has been applied to several salmonid species in the Pacific Northwest (Agrawal et al., 2005;Cooney and Holzer 2007;Busch et al., 2013;Bidlack et al., 2014;Flitcroft et al., 2014). Our IPM for brook trout spawning habitat incorporates habitat suitability indices for three variables used by Burnett et al. (2003): discharge, gradient, and valley confinement, and a fourth variable, temperature, to account for the thermal sensitivity of both species. Previous habitat suitability indices for brook trout have been used on their native range including the Southern Blue Ridge (Schmitt et al., 1993) and the Mid-Atlantic region of the eastern United states (Smith and Sklarew 2012). However, these indices do not account for specific ontogenic stages (i.e., spawning habitat). Other habitat suitability indices for brook trout (Raleigh 1982) do not account for effects that landscape characteristics may have on stream habitats (Smith and Sklarew 2012). Here, we used expert opinion surveys to define habitat suitability index curves specifically for brook trout spawning habitat. To our knowledge, this is the first attempt to create a comprehensive IPM for brook trout outside its native range.
We included temperature as an additional variable in the IPM given the thermal sensitivity of both bull trout and brook trout.
Temperature has been identified as one of the primary factors limiting the distribution of these species. For spawning, bull trout require stream temperatures below 9°C (Ardren et al., 2011;Austin et al., 2018), whereas brook trout has a higher upper thermal threshold (up to 10°C). Therefore, projected warming trends could result in the loss of 18-92% of suitable bull trout habitat in some areas (Rieman et al., 2007). As water temperatures increase, available suitable bull trout spawning habitat will likely decrease, while brook trout spawning habitat may increase. This may lead to greater areas of habitat overlap and thus, increased risk of hybridization between the two species.
The HRM approach does have some limitations, however. First, it does not model the actual rate of hybridization, rather it focuses on the potential risk of hybridization. Second, other habitat or landscape characteristics (e.g., groundwater input) not included in this IPM could play an important role in the selection of redd sites (e.g., Snyder et al., 2015). These characteristics may be able to be incorporated in the future to further refine the IPM for brook trout spawning. Third, uncertainties on survey responses to construct brook trout habitat suitability curves are not incorporated, instead we used the mean values of the responses. Survey results may be improved using a more focused group of participants. Fourth, monthly temperature data for specific streams during spawning season are not always available, requiring the HRM to use estimated temperature values from NorWeST . This is an illustrative example of the HRM applied to multiple spatial scales (spatial resolution includes 100 m river segments).
Frontiers in Environmental Science | www.frontiersin.org March 2022 | Volume 10 | Article 834860 stream temperature models (Isaak et al., 2017). In the future, collection of field data in specific streams may address some of these limitations. The HRM model is still unvalidated. However, sub-basins classified as of Extreme or High hybridization risk show an overall agreement with documented presence of hybrids in the Upper Malheur headwaters (DeHaan et al., 2010;Haslick 2021) and the Powder River (Howell 2018). Future efforts to validate this modelling approach could include intensive sampling in areas where species are already known to be hybridizing. For example, the use of environmental DNA (eDNA) can provide a non-invasive method to establish presence of both parent species across stream reaches during the spawning season. In addition, field data collected at the reach scale, such as the presence of parent species and/or hybrid offspring along with physical stream characteristics, may be useful to validate and improve this modelling approach in the future.
Managing threatened and/or endangered species is challenging under the best circumstances. Data collection for specific streams is time-consuming and cost prohibitive. Changing climatic conditions are presenting additional challenges as well as an increased urgency in identifying streams and populations that can benefit most from management efforts. Emergency managers realized that assessing and communicating risk can be subjective and developed methods to quantify risk in a defensible manner (Waugh and Tierney 2007). A similar problem exists when assessing hybridization risk in freshwater species. By applying some of the techniques used in emergency management we can begin to quantify hybridization risk in a way that is easier to communicate to policymakers, funding agencies, and the public. The HMR approach used here provides a tool for natural resource managers to use when evaluating stream reaches with species at risk for future stream restoration and protection efforts.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
MM, IA, and GG contributed to conception and design of the study. MM organized the database. MM and JAO developed the intrinsic potential model. MM wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.