Fine‐scale environmental DNA sampling reveals climate‐mediated interactions between native and invasive trout species

Abstract

It is widely recognized that biotic interactions may act as important mediators of species responses to climate change. However, collecting the abiotic and biotic covariates at the resolution and extent needed to reveal these interactions from species distribution models is often prohibitively expensive and labor‐intensive. Here we used crowd‐sourced environmental DNA sampling—the inference of species presence from genetic material in the environment—and high‐resolution habitat covariates across 630 sites over an area of nearly 10,000 km² to build an accurate species distribution model (AUC = 0.96; prediction accuracy = 0.90) for bull trout in cold‐water habitats that incorporates fine‐scale, context‐dependent interactions with invasive brook trout. We then used this model to project possible climate change and brook trout invasion scenarios for bull trout forward in time. Our environmental DNA sampling results were concordant with traditional electrofishing samples in the basin and revealed species patterns that were consistent with previous studies: Bull trout were positively associated with larger stream sizes and negatively associated with high brook trout abundances. However, our modeling also revealed an important nuance: At high abundance, brook trout appear to exclude bull trout from small streams, even those below the thermal optima for brook trout. Climate projections suggest a loss of suitable bull trout habitat as streams warm and summer flows decrease, which could make deleterious interactions with brook trout more common in the future. Where brook trout are invading bull trout habitats, streams that are both large and cold are most likely to provide native bull trout with long‐term refuges.

Introduction

Accurately predicting species distributions in response to climate change is key to conserving biodiversity this century (Guisan et al. ). Species distribution models (SDM) are often used to associate occurrence records with environmental conditions and then forecast future distributions based on predicted changes in those conditions (Guisan and Zimmermann ). Climate change will likely drive complex changes in the distribution of ecological communities and the interactions among their members (Hellmann et al. , Zarnetske et al. , Alexander et al. ). These factors are often overlooked or poorly addressed in SDMs (Davis et al. , Potter et al. ), though there has been increasing attention on the potential for covariates such as dispersal or species interactions to improve model accuracy (Guisan and Thuiller , Brooker et al. , Angert et al. ). Yet the greatest limitation to estimating those effects is often not conceptual or analytical, but logistical: It has been prohibitively expensive to collect spatially dense, multi‐species datasets across the broad scales needed to tease apart the ecological processes shaping species distributions (Guisan and Thuiller , Wisz et al. ).

New technologies are revolutionizing ecological data collection in ways that promise to overcome historical sampling limitations. In particular, environmental DNA (eDNA) sampling has emerged as a cost‐effective, reliable tool for surveying aquatic species in freshwater ecosystems. This sampling approach collects genetic material from the environment (e.g., water) to infer species presence (Jerde et al. ). Recent studies have found that eDNA sampling, when coupled with properly designed genetic markers, provides robust data on species occurrence with low error rates (e.g., Thomsen et al. , Pilliod et al. , Valentini et al. ). Because similar sampling and analysis protocols have been developed to detect eDNA from bacteria, fungi, plants, insects, mollusks, fishes, amphibians, mammals, birds, and reptiles (e.g., Schmidt et al. , Padgett‐Stewart et al. , Newton et al. , Valentini et al. ), individual environmental samples can be used to evaluate occupancy of many species, even those that were not targets of the original sampling effort (Dysthe et al. ). When those samples are collected at ecologically and demographically meaningful scales, patterns of species occurrence may reveal the influence of species interactions. A more nuanced understanding of these interactions may also be possible because, at least in some systems, DNA concentrations in environmental samples can be used as indices of local abundance (Takahara et al. , Wilcox et al. ). Finally, eDNA sampling is unlike most traditional sampling methods in that sample collection is rapid and requires little expertise if robust protocols are followed (e.g., Carim et al. ). This facilitates crowd‐sourced data collection and rapid inventories of large areas at relatively low cost (Biggs et al. , McKelvey et al. ). These attributes make eDNA sampling an efficient tool for developing datasets used in SDMs that enable assessments of multiple species and their interactions.

Salmonid fishes, which include trout, salmon, char, and whitefish, are socioeconomically important, cold‐water fishes distributed throughout the Northern Hemisphere. Many salmonid species have declined within their native ranges due to overharvest, habitat degradation and fragmentation, and interactions with nonnative species (e.g., Morita and Yamamoto , Fausch et al. ). The latter can be particularly troublesome because they often involve other salmonid species with extensive niche overlap that have been widely and repeatedly introduced and are difficult to remove once they become established. In western North America, the bull trout (Salvelinus confluentus) is a native species of conservation concern that declined dramatically during the 20th century and is now accorded federal protection in the USA and portions of Canada (COSEWIC , U.S. Fish and Wildlife Service ). Bull trout spawn exclusively in cold headwater habitats where their juveniles reside and grow for 2–3 yr, leading to populations that are spatially disjunct and fragmented among natal habitats separated by warmer main‐stem rivers that are inhospitable to juvenile fish (Rieman and Dunham , Isaak et al. ). Brook trout (S. fontinalis), a char native to eastern North America, have been widely introduced throughout the range of bull trout and, although adapted to somewhat higher stream temperatures, exhibit substantial niche overlap with bull trout (Isaak et al. ). This has led to widespread exploitation and interference competition and weakly introgressive hybridization (DeHaan et al. , Warnock and Rasmussen , b), and has contributed to the local displacement and extirpation of bull trout populations (Rieman et al. , McMahon et al. , Warnock and Rasmussen , b).

The likelihood of bull trout persistence is further jeopardized by climate change, which is altering thermal and hydrological regimes of aquatic ecosystems across western North America (Luce and Holden , Leppi et al. , Isaak et al. , Young et al. ). In response, bull trout populations are contracting into smaller and more isolated headwater habitats (Eby et al. , Al‐Chokhachy et al. ) where the probability of extirpation may be higher (Rieman et al. , Peterson et al. , Isaak et al. ). The presence of brook trout has been shown to decrease bull trout occupancy and is projected to do so under climate change scenarios (Rieman et al. , Isaak et al. ), but a precise understanding of the ecological mechanisms responsible for this decline—and their spatial predictability—has been missing.

Environmental DNA sampling is particularly useful for studies of bull trout because the species naturally occurs at low population densities and detecting populations via conventional methods is both labor‐intensive and expensive (Peterson et al. ). Previous studies have found that eDNA sampling is more sensitive than traditional sampling for detecting rare salmonids in streams, is a reliable indicator of brook trout abundance, and accurately reflects bull trout distributions previously determined by traditional sampling (McKelvey et al. , Wilcox et al. ; Appendix S1). Here, we demonstrate the utility of repurposing a systematically distributed and ecologically guided set of eDNA samples originally collected to evaluate the presence of bull trout in natal habitats across two river basins (Isaak et al. , Young et al. ) to develop an accurate SDM for bull trout that incorporates the effects of sympatry with different densities of brook trout. First, we re‐analyzed eDNA samples to estimate brook trout relative abundance. We then linked the biological observations to geospatial stream habitat and climate covariates and used mixed‐effects SDMs that predicted the probability of bull trout occurrence. Finally, we used this model in conjunction with high‐resolution climate change scenarios to predict distributional shifts in bull trout due to climate change and brook trout invasion.

Methods

Bull trout and brook trout detection using eDNA sampling

Samples were collected from 2015 to 2016 across two, 4th‐code hydrologic units (Flint‐Rock and Upper Clark Fork; 9569 km² in extent) in west‐central Montana, USA, as part of the Rangewide Bull Trout eDNA Project (Young et al. ; https://www.fs.fed.us/rm/boise/AWAE/projects/BullTrout_eDNA.html) led by the National Genomics Center for Wildlife and Fish Conservation (Fig. ). For our analysis, we used 630 samples collected at 1‐km intervals throughout 54 potential juvenile bull trout cold‐water habitat patches identified using a rangewide juvenile bull trout occurrence model (Isaak et al. ). Potential habitat patches consisted of contiguous stream segments in which individual reaches had mean August water temperatures ≤11°C, summer flows ≥0.0057 m³/s, and stream reach gradients ≤15% (Isaak et al. ). In a few habitat patches, short (<1 km) segments of stream exceeded the 11°C threshold, but these were retained as part of the larger surrounding natal area when delineating patches (Isaak et al. ). Field sampling maps and spatial coordinates for sampling sites were downloaded from the Rangewide Bull Trout eDNA Project website, and sampling was conducted by personnel from the University of Montana, U.S. Forest Service, Trout Unlimited, the Clark Fork Coalition, and Montana Fish, Wildlife, and Parks. For this analysis, we removed from consideration any eDNA samples that were collected above a known fish passage barrier. Samples at the 630 field locations were collected using a standardized protocol (Carim et al. ). Briefly, this protocol involves filtering 5 L of stream water through a 1.5‐micron pore size microfiber glass filter, which is stored in silica desiccant until arrival at the laboratory where samples are held at −20°C prior to analysis. All sampling equipment was sterilized with a 50% household bleach solution in a dedicated laboratory space, then provided to field crews in single‐use sample packets to minimize contamination risk.

Genetic analysis

In the laboratory, we extracted DNA from one half of each filter using the DNeasy Blood and Tissue Kit and QIAshredder columns (QIAGEN) with modifications from the manufacturer's protocol as described in Carim et al. (). Each extraction was assayed for mitochondrial DNA from bull trout and brook trout using species‐specific quantitative PCR (qPCR) markers described in Wilcox et al. (); BUT1 and BRK2, respectively. Fifteen‐microliter reactions were analyzed in triplicate for each sample and composed of 7.5 μL 2 × TaqMan Environmental Mastermix 2.0 (Life Technologies, Carlsbad, California, USA), 0.75 μL of 20 × assay mix (primers and FAM‐labeled, MGB‐NFQ hydrolysis probe; Integrated DNA Technologies and Life Technologies), 4 μL template DNA, and laboratory‐grade sterile water. Thermocycling conditions for both assays were 95°C 10 min, [95°C 15 s, 60°C 60 s] × 45 cycles. Each sample was first analyzed with the bull trout qPCR assay, which was multiplexed with an internal amplification control (IPC; TaqMan Exogenous Internal Positive Control from Life Technologies) to test for the presence of PCR inhibitors. Any samples with evidence of PCR inhibition (≥1 C_t shift in the IPC amplification curve) were treated using an inhibitor removal kit (Zymo Research; McKee et al. ) and re‐analyzed (n = 44). Only, samples with complete PCR inhibition removal were assayed for brook trout DNA and used for analysis. Each PCR plate included at least one set of triplicate positive control and no‐template control wells. Plates for brook trout also included a five‐point, triplicate, 1:5 dilution sequence (standard curve) from 10 to 6250 DNA copies/reaction constructed from a synthetic gene containing the brook trout assay amplicon sequence (Integrated DNA Technologies; see Wilcox et al. , for details on standard curve preparation). We used this standard curve to estimate the concentration of brook trout eDNA in each sample. We considered linear amplification in any of the three replicate reactions for an assay to be a positive detection. For quantification, we averaged across all three replicates (assigning wells without amplification a quantity of zero; Ellison et al. ). All extractions and PCR setup were done in separate, dedicated laboratory spaces. We regarded species as present at all sites with a positive detection, acknowledging that on occasion, eDNA may be collected below the downstream most individual in a population (Jane et al. , Pont et al. ).

Habitat covariates

Based on previous research, we hypothesized that the distribution and abundance of bull trout and brook trout would be influenced by the large gradients in stream temperature, discharge, and channel slope (gradient) that occur across mountainous landscapes. To describe those covariates, our sample sites were linked to the 1:100,000‐scale National Hydrography Dataset (http://www.horizon-systems.com/NHDPlus/index.php), which provides a nationally consistent reach vector representation of stream networks (McKay et al. ). Channel gradient for each reach in the study area was estimated as the drop in channel elevation per unit length (m/m × 100) using a custom Python script in ArcGIS (GRAD). Summer flow estimates (FLOW) were derived for each reach for a recent historical period (1993–2011) and two future climate scenarios that were downloaded from the Western U.S. Flow Metrics website (http://www.fs.fed.us/rm/boise/AWAE/projects/modeled_stream_flow_metrics.shtml; Wenger et al. ). The two future scenarios were based on the A1B emissions scenario for the 2040s (2030–2059) and 2080s (2070–2099; IPCC 2007) and global circulation model (GCM) projections of air temperature and precipitation changes from a ten‐climate model ensemble (Hamlet et al. ). Those projections were used to force the Variable Infiltration Capacity (VIC) hydrologic model and obtain future summer flow deltas (Wenger et al. ). Summer temperature estimates (TEMP) were based on mean August water temperatures for the same recent historical period and two A1B future periods, and were downloaded from the NorWeST website (http://www.fs.fed.us/rm/boise/AWAE/projects/NorWeST.html). These future scenarios were based on the same ten‐climate model ensemble used for the flow scenarios. Within the study area, stream temperature scenarios were interpolated from 340 observations at 152 sites from 1993 to 2011 (mean number of years data per site = 2.24, range = 1–16) and were part of an accurate underlying temperature model (r² = 0.91; RMSPE = 0.98°C; MAE = 0.62°C) fit to a much larger regional dataset (Isaak et al. ). Although the A1B emissions scenario was originally run for the third phase of the Coupled Model Intercomparison Project (CMIP3) and has been superseded by Representative Concentration Pathways (RCP) in CMIP5, there is a strong similarity between A1B and RCP 6.0 (Wright et al. ).

Bull trout models and scenarios

Samples within each habitat patch were spatially autocorrelated, so we modeled detection of bull trout eDNA using binomial generalized linear mixed models (GLMMs), treating habitat patch as a random effect group. This modeling framework allowed us to consider independent regression intercepts and coefficient slopes for each habitat patch. We used a binomial model because bull trout eDNA analyses were conducted without use of a standard curve, and so inference about bull trout eDNA concentration was not possible. Main effects in the model included GRAD, FLOW, TEMP, and brook trout eDNA concentration (BRK; estimated mtDNA copies/reaction). We also considered interactions between BRK, FLOW, and TEMP because the latter two covariates are thought to be important abiotic controls on competitive interactions between brook trout and bull trout (Rieman et al. , Wenger et al. , Warnock and Rasmussen ). Because local brook trout abundance is correlated with brook trout eDNA concentration in mountain streams in this region (r = 0.727 in Wilcox et al. , Appendix S1: Fig. S1), we used estimated brook trout eDNA concentration as an index of brook trout abundance in our analyses. Brook trout eDNA concentration and mean summer stream discharge were right‐skewed and so were log‐transformed to improve model fit (1.0 was added to brook trout eDNA concentration prior to log‐transformation to accommodate values of 0.0; Table ). A quadratic term (BRK²) was also included as a candidate covariate because initial data exploration suggested that the response of bull trout occupancy to brook trout abundance was non‐linear, with a small positive relationship at low brook trout abundances and a negative effect at high brook trout abundances. Stream temperature relationships with species occurrence are also usually non‐linear when examined across a broad range of thermal conditions (Isaak et al. ), but that pattern did not appear in the truncated range of conditions sampled here (i.e., sites <11°C) and so a linear parameter was used. Some covariates were significantly correlated (Table ), but correlations were lower than the 0.7 level considered problematic with regard to collinearity (Dormann et al. ).

Summaries of candidate covariates used to build models predicting probability of bull trout eDNA presence at a site

Covariate	Min	Max	Mean	Median
TEMP	5.56	12.72	9.06	9.10
FLOW	0.21	98.70	8.27	3.79
GRAD	0.00	0.13	0.05	0.05
BRK	0	841	34	0

Note

Mean August stream temperature (°C; TEMP), mean summer stream discharge (cfs; FLOW), stream gradient (GRAD), and brook trout eDNA concentration (mtDNA copies/PCR; BRK).

Correlation matrix of candidate covariates used to build models predicting probability of bull trout eDNA presence at a site

Covariate	TEMP	FLOW	GRAD	BRK
FLOW	0.235
GRAD	−0.543	−0.441
BRK	0.184	0.011	−0.150
BULL	0.169	0.457	−0.310	−0.145

Models were fit using the package lme4 (Bates et al. ) in the statistical environment R (R Core Development Team ) using the glmer function which fits maximum likelihood models using a Laplace approximation. Following the recommendations of Zuur et al. (), we conducted model selection in two stages. First, we selected the random effect structure using a full fixed‐effects model containing all of the terms under consideration, and then, we selected the fixed‐effects terms. Our full model contained all candidate main‐effect and interaction terms (TEMP, FLOW, BRK, GRAD, BRK × TEMP, BRK × FLOW, and BRK²). We used the Akaike Information Criterion (AIC) to compare random‐effect structures with a random intercept only versus a random intercept and all seven possible combinations of random slopes for FLOW, TEMP, and BRK (Bozdogan , Burnham and Anderson ). We did not consider a random‐slope effect for GRAD as it did not appear to be an important term in any models in initial data exploration. We note that Zuur et al. () recommend not comparing information criteria such as AIC to select random effects from maximum likelihood model fits because maximum likelihood estimates of variance components are biased. However, this bias is small when sample sizes are large (Zuur et al. ), as in this study. For calculating AIC, we counted one degree of freedom for each random effect variance term and correlation term (i.e., one random effect = one degree of freedom, two random effects = three degrees of freedom). The lowest AIC model contained random slopes for BRK and FLOW (Table ).

Random effect structure selection for a generalized linear mixed model to predict probability of bull trout eDNA presence at a site with random effects by habitat patch

Random effect structure	AIC	ΔAIC
Intercept + FLOW + BRK	467.0
Intercept + BRK	469.4	2.4
Intercept + FLOW + TEMP + BRK	471.7	4.7
Intercept + TEMP + BRK	472.9	5.9
Intercept + FLOW	476.6	9.6
Intercept	477.7	10.7
Intercept + TEMP	478.8	11.8
Intercept + FLOW + TEMP	478.8	11.8

Notes

We considered a random intercept and random slopes for each covariate. Each model was fit with the same full fixed‐effect structure (all considered main effects and interaction terms). Models are ranked by AIC.

Using this random effect structure, we then used AIC to compare fixed‐effect structures for both models. All candidate models included the main effects TEMP, FLOW, and BRK, plus all possible combinations of the interactions BRK × TEMP and BRK × FLOW and main effects for GRAD and BRK² (n = 16 models total). For each model, we estimated the area under the curve (AUC) of the receiver operating characteristic (ROC) using the pROC package (Xavier et al. ) and model accuracy with a 0.5 threshold (custom R script), then ranked all models by AIC (Table ).

Fixed‐effect structure selection for a generalized linear mixed model to predict probability of bull trout eDNA presence at a site with random effects by habitat patch

Model	AIC	ΔAIC	AUC	Acc.
FLOW + TEMP + BRK + BRK2 + BRK × FLOW	464.0		0.964	0.897
FLOW + TEMP + BRK + BRK2	465.2	1.2	0.968	0.898
FLOW + TEMP + GRAD + BRK + BRK2 + BRK × FLOW	465.2	1.2	0.965	0.898
FLOW + TEMP + BRK + BRK2 + BRK × TEMP + BRK × FLOW	465.6	1.6	0.964	0.900
FLOW + TEMP + BRK + BRK2 + BRK × TEMP	466.5	2.5	0.967	0.900
FLOW + TEMP + GRAD + BRK + BRK2 + BRK × TEMP + BRK × FLOW	467.0	3.0	0.965	0.900
FLOW + TEMP + GRAD + BRK + BRK2	467.1	3.1	0.968	0.898
FLOW + TEMP + BRK + BRK × FLOW	467.8	3.8	0.963	0.903
FLOW + TEMP + BRK	468.3	4.3	0.966	0.897
FLOW + TEMP + BRK + BRK × TEMP + BRK × FLOW	468.5	4.4	0.963	0.900
FLOW + TEMP + GRAD + BRK + BRK2 + BRK × TEMP	468.5	4.5	0.966	0.894
FLOW + TEMP + GRAD + BRK + BRK × FLOW	468.7	4.7	0.963	0.898
FLOW + TEMP + BRK + BRK × TEMP	468.7	4.7	0.966	0.894
FLOW + TEMP + GRAD + BRK + BRK × TEMP + BRK × FLOW	469.7	5.7	0.964	0.892
FLOW + TEMP + GRAD + BRK	470.2	6.1	0.966	0.894
FLOW + TEMP + GRAD + BRK + BRK × TEMP	470.4	6.4	0.964	0.894

Notes

All models used the same random effect structure (random intercept and slopes for FLOW and BRK). Models are ranked by AIC. Also shown is the estimate model AUC and accuracy (threshold = 0.5).

We used the package influence.ME (Nieuwenhuis and Grotenhuis ) to test for highly influential habitat patches that could be driving parameter estimates in our model. We sequentially held out each habitat patch from the data and then refit the model to derive parameter estimates with and without each patch and to estimate patch influence (Cook's distance; Cook ). The lme4 package used for model fitting approximates fixed‐effect P‐values using a Wald's Z‐test. For the final model, we also used a parametric bootstrap to estimate a 95% confidence interval around each coefficient estimate (500 replicates using the boot.mer function in lme4; Table ).

Fixed‐effect estimates for the top model predicting probability of bull trout eDNA presence at a site

Parameter	Estimate	95% CI
Intercept	−8.896	−12.130 to −6.420
FLOW	2.370	1.850–3.192
TEMP	0.425	0.126–0.690
BRK	1.241	0.444–2.231
BRK2	−0.187	−0.379 to −0.061
BRK × FLOW	−0.292	−0.510 to −0.075

Note

The 95% CI around each estimate was determined from 500 parametric bootstraps.

The best bull trout model was used to forecast the probability of bull trout occupancy at each of the 630 sites using covariate values from the 2040 and 2080 climate scenarios. Potential future BRK abundances under these climate scenarios were derived from a model described in the section below. For the 2040 and 2080 scenarios, we extracted the mean and 0.025 and 0.975 percentiles (95% confidence interval) of the bull trout occupancy predictions at each site based on the BRK scenarios. Sites in the future scenarios with TEMP >12°C were assigned values of zero bull trout occurrence probability because earlier empirical assessments based on regional fish surveys found these habitats to be thermally unsuitable for juvenile bull trout (Isaak et al. , , b). To assess the impact of brook trout on future bull trout occupancy probabilities, we compared the predicted probability of bull trout under scenarios that allowed brook trout abundance to increase in response to changing environmental conditions with predicted probability of bull trout presence assuming no change in brook trout abundance.

Brook trout abundance model

To forecast brook trout abundance for use with climate projections, we used the R package lme4 to build a linear mixed model that predicted log‐transformed BRK + 1 as a function of TEMP, FLOW, and GRAD because these covariates are also known to affect the distribution of this species (Wenger et al. , Isaak et al. ). For future projections, the model was applied only to sites within cold‐water habitat patches where brook trout were detected at least once (n = 476 sites within 32 habitat patches; other sites held to ≤5 mtDNA copies/reaction) because locations in other patches were often inaccessible to brook trout. We used 500 replicates of a parametric bootstrap to estimate a 95% CI around each fixed‐effect parameter (the boot.mer function in lme4). We then used covariate values associated with recent, 2040, and 2080 stream climate scenarios to predict BRK + 1 concentrations at each site. We assessed the suitability of this brook trout model by assessing the performance (classification accuracy and AUC) of our bull trout model when using predicted brook trout eDNA concentrations instead of observed brook trout eDNA concentrations. Further, to address high uncertainty in our brook trout model, under the 2040 and 2080 scenarios for each site we took 500 random draws from a Gaussian distribution based on the fixed‐effect prediction interval (estimated using the function predictInterval in the MuMin package).

Results

Genetic analysis

All positive control samples amplified DNA of the targeted species. One plate for bull trout was re‐run because of low‐level amplification in a no‐template control, whereas there was no amplification in the remaining negative control samples. Brook trout standard curve efficiency was 85–106% (mean = 94%; r² > 0.98). We detected bull trout at 262 of the 630 sites, brook trout at 293 sites, both species at 121 sites, and neither species at 196 sites. Brook trout eDNA concentrations were low at most sites (76% of positive sites had <100 copies/reaction), but the distribution of eDNA concentrations was skewed (range = 0–841 copies/reaction) toward a small number of sites with high values.

Habitat relationships

Bull trout were detected across a range of stream sizes, but disproportionately so in larger streams (e.g., median FLOW for all sites was 3.79 cfs, but 9.32 cfs for sites where bull trout were detected; Fig. ). Sympatry with brook trout was associated with bull trout being found in even larger streams (median FLOW = 15.98 cfs; Fig. ). Brook trout were also found across the range of stream sizes, but most often in smaller channels; 50.0% of brook trout detections were at sites <5.00 cfs. At sites where brook trout were allopatric, this pattern was even more pronounced (median FLOW = 2.59 cfs; Fig. ). Neither taxon was detected at sites <6.69°C, but one or both were present across the range of warmer stream temperatures. Overall, sites lacking both species tended to be smaller, colder, and steeper than sites with either or both of them (median FLOW, 1.88 vs. 6.30 cfs; median TEMP, 8.61 vs. 9.24°C; median GRAD, 6.2 vs. 4.2%). Finally, brook trout eDNA concentrations were positively correlated with smaller, warmer, and less steep streams (Fig. ).

Trout models

Our top fixed‐effects model for bull trout was FLOW + TEMP + BRK + BRK² + BRK × FLOW (random effect structure: intercept + FLOW + BRK). This model had excellent fit (AUC = 0.964, accuracy = 0.897 at a 0.5 threshold), and all of the coefficient estimates were statistically significant (Walds Z‐test P < 0.05 and parametric bootstrap 95% CIs did not include zero; Table ). The BRK² and BRK × FLOW interaction terms were also included in two of three of the other models within 2 ΔAIC points of the top model (Table ). The influence of individual habitat patches on model parameter estimates was minor and largely a function of the number of sites per patch (Fig. ; Appendix S2: Figs. S2, S3).

Bull trout occupancy probability was positively related to FLOW and TEMP (Fig. ). Occupancy probability was positively correlated with BRK but negatively correlated with BRK² (Table ). The result was that bull trout occupancy probability increased at low brook trout eDNA concentrations (<10 copies/PCR) but declined at higher concentrations (Fig. ; Appendix S2: Fig. S1). However, the negative interaction between brook trout eDNA concentration and stream discharge indicated that bull trout occupancy probability only increased at sites with low discharge when brook trout abundance was low (Fig. ).

The model predicting brook trout eDNA concentration had reasonable fit (RMSE estimated using merTools = 1.464, conditional r² estimated using MuMin = 0.521). Our model predicted that brook trout eDNA concentrations increased with mean August stream temperature (coefficient mean = 0.486, 95% CI = 0.310–0.673) and declined with mean summer discharge (coefficient mean = −0.003, 95% CI = −0.017 to 0.011) and reach gradient (coefficient mean = −12.036, 95% CI = −19.833 to 4.496). Using predicted brook trout eDNA concentrations instead of observed values to predict current bull trout distributions substantially reduced predictive power, but it still retained good accuracy (accuracy at a 0.5 threshold = 0.759, AUC = 0.801).

Bull trout scenario responses

Relative to recent conditions, the 630 sites are predicted to increase in mean August stream temperature (mean change = 1.18 and 2.04°C for 2040 and 2080, respectively), decline in mean summer flow (mean change = −2.76 to −3.73 cfs [27.0–35.8%] for 2040 and 2080, respectively), and increase in brook trout eDNA concentration (mean change = 7.46 and 16.39 mtDNA copies/PCR for 2040 and 2080, respectively). Consequently, the bull trout model predicted future declines in occurrence probabilities that averaged 0.101 by the 2040s and 0.197 by the 2080s (Fig. ). The probability of bull trout occupancy was predicted to increase by small amounts at about a quarter of the sites by 2040 or 2080, but only among sites with a relatively low initial probability of occupancy (<0.5). Conversely, large declines in occupancy probability (≥0.25) were predicted for 6.7–27.9% of sites (2040–2080, respectively; Fig. ). Most of the large predicted declines in occupancy probability were due to future stream temperatures exceeding 12°C and becoming thermally unsuitable for juvenile bull trout. These same sites, however, also tended to be larger and to have above‐average probabilities of bull trout occupancy under historic conditions (Figs. , ). Unsurprisingly, removal of brook trout was predicted to increase the likelihood of bull trout occupancy at most sites in all future scenarios, but these improvements were highly variable. The greatest improvements were in those streams which currently have the highest probabilities (>0.75) of bull trout presence.

Discussion

Systematically distributed eDNA samples coupled with high‐resolution geospatial covariates throughout stream networks revealed fine‐scale biotic interactions that shape potential bull trout responses to future climate change. Earlier studies have indicated that bull trout were likely to undergo broad‐scale declines as a consequence of warming stream temperatures (Rieman et al. ), but lacked patch‐scale specificity. Further refinements resolved how biotic and abiotic characteristics influenced patterns of occupancy in individual cold‐water habitats (i.e., patches) and identified which habitats were likely to remain suitable in the future (Isaak et al. ). The present study, albeit of more limited geographic scope, adds greater spatial precision and ecological realism by revealing the reach‐specific characteristics within cold‐water habitats that dictate the outcome of interactions between brook trout. Overall, when brook trout are abundant, they appear capable of excluding bull trout from small streams across a range of temperatures, but rarely occur at high densities in larger streams. Flows in these streams, however, are expected to decline in the future, rendering them more suitable for brook trout and less suitable for bull trout. This shift, and the warming temperatures expected this century (Isaak et al. ), is expected to cause bull trout distributions to contract (Eby et al. , Al‐Chokhachy et al. ) and may eventually lead to extirpations of local populations, particularly where the only remaining habitats are small streams.

The relationships revealed in this study, though more ecologically nuanced and spatially explicit, are consistent with previous research on these species, providing further support for the robustness of ecological inference that can be made sight‐unseen about aquatic taxa with eDNA sampling. That brook trout are associated with smaller, warmer, and lower‐gradient streams, and that bull trout are associated with larger, colder streams, but may be competitively excluded by brook trout from some habitats, are all expectations based on previous work (e.g., Rieman et al. , Wenger et al. , Benjamin et al. ), that were rendered more spatially explicit by comprehensive eDNA sampling. One non‐intuitive finding, however, was that bull trout occupancy was positively related to low numbers of brook trout in small streams. We suspect that this was indicative of streams that were suitable for and accessible to both species, but in which bull trout were largely unaffected by the relatively small numbers of brook trout. For example, the relation between bull trout presence and brook trout abundance was only positive when concentrations of brook trout eDNA were low (Fig. , likely equivalent to <3 fish/100 m in small streams; Wilcox et al. ).

Landscape consequences of brook trout invasion and climate change

Even given temperatures below their thermal optimum, brook trout appear to be better at exploiting small streams than are bull trout. This may be particularly problematic for bull trout in small or isolated streams that adopt resident life histories in which adults mature at small sizes and are non‐migratory (McPhail and Baxter ) and therefore exhibit greater niche overlap with, and weaker biotic resistance to, brook trout (Warnock and Rasmussen ). Resident bull trout may also suffer greater demographic losses, in the form of wasted reproductive effort, as a consequence of hybridization with brook trout. In contrast, bull trout populations that include migratory individuals may be more resistant to brook trout invasions. These individuals mature at much larger adult sizes and are vastly more fecund, and the size‐selective mate choice typical of salmonids would tend to favor pairing of migratory adult bull trout during spawning. Their large size, moreover, may preclude them from occupying very small stream channels where interactions with brook trout might be the strongest, and might contribute to the positive influence of flow on bull trout habitat occupancy in our model.

The susceptibility of bull trout to brook trout invasion in small streams has important population‐scale consequences. Displacement of bull trout from small streams by brook trout effectively reduces cold‐water habitat size, which is a key driver of bull trout occupancy probability (Dunham and Rieman , Isaak et al. ). This explanation ties our reach‐scale observations with those at the patch‐scale by Isaak et al. (), who found that brook trout invasion of individual cold‐water habitats had the same effect on probability of occupancy by juvenile bull trout as did halving the size of the habitat patch. This is of particular concern because even though small streams represent lower quality habitat for bull trout, most thermally suitable habitat resides in small streams. For example, among our study sites, about three‐fourths have a predicted mean summer flow <10 cfs and a third are <2 cfs. In some cases, this might result in brook trout displacing bull trout from headwater reaches into larger downstream habitats that are at risk of exceeding the thermal tolerance of juvenile bull trout. An immediate implication of our analysis is that the few large streams that continue to be thermally suitable for bull trout in coming decades could play a critical role in shaping their distribution on the landscape. Thus, conservation investments within headwater streams that might have the greatest benefit for future bull trout populations may include efforts to cool temperatures of larger streams, enhance flows of colder ones, and to suppress or eradicate brook trout populations.

Sources of uncertainty

Models predicting species responses to climate change are vulnerable to a host of assumptions. In our case, this is particularly acute with respect to predictions about future climate, brook trout invasions, and interactions between the two. There is uncertainty about future climatic conditions because of uncertainty about future emissions and because of variation among global circulation models, even under identical emission scenarios (Stainforth et al. ), which is compounded by uncertainty induced by model scale downsizing (Ahmadalipour et al. ). For bull trout, this uncertainty is most important in larger streams that are near the thermal threshold of becoming unsuitable habitat.

As with previous efforts to model the distribution of brook trout in their invaded range, our simple brook trout model left a large proportion of the variance in brook trout eDNA concentration unexplained. Besides unmeasured, fine‐scale covariates within sampled areas, it may be that proximity to very‐low‐gradient, warm reaches that serve as springboards for brook trout invasions is an important factor explaining their distribution and abundance, but this is something that our sampling did not address (Adams et al. , Benjamin et al. , Wenger et al. ). When predicting future scenarios, we restricted our analyses and projections to cold‐water habitats where brook trout were detected at least once (i.e., brook trout were known to be present). This assumption that brook trout have achieved equilibrium within currently invaded habitat is consistent with the length of time that brook trout have been in many of these systems (50–150 yr for many locations in the western USA; Dunham et al. ) and with observations of brook trout distributions over time in similar invaded systems (Adams et al. , Howell ). We note that our future projections are based on the assumption that human translocation will not lead to further brook trout invasions in this area, which may not be realistic. Conversely, although brook trout abundances were predicted to increase across invaded habitats, our modeling of brook trout does not account for other impacts of climate change on brook trout populations. For example, Wenger et al. () suggested that climate change may result in an increased frequency of high winter flows, which are detrimental to survival of the eggs and fry of juvenile brook trout (Wenger et al. ). Thus, it is possible that climate change could also lead to declines in brook trout abundance in some habitats.

Finally, we recognize that eDNA sampling is limited in the sense that demographic and population genetic characteristics of the target species are either unknown or must be inferred. As a result, we had to make some assumptions in interpreting our data. First, we assumed that the presence of mtDNA indicated the presence of the corresponding parental taxa. Brook trout and bull trout can hybridize, and our eDNA markers are unable to distinguish between hybrids and parental taxa with the same mitochondrial haplotype. However, high concordance with traditional sampling (Appendix S1) indicates that hybrids had little or no effect on inference, which is consistent with extremely low fitness of brook‐bull trout hybrids (Kanda et al. ). Second, we assumed for model interpretation that bull trout eDNA presence indicated the presence of juvenile bull trout. To meet this assumption, we intentionally restricted sampling to thermally suitable potential natal habitats for bull trout. In a previous study of 81 sites in streams with thermal environments similar to those sampled here, 74 (91%) were occupied by bull trout <150 mm (i.e., putative juvenile fish; Rieman et al. ). Thus, when bull trout are detected within these habitats, there is a high probability that juveniles are present. However, occasional detection of individual adult bull trout where there are no juvenile bull trout could cause us to exaggerate the suitability of warmer habitats for juvenile animals.

Species distribution modeling via crowd‐sourced sampling

This study was made possible by drawing on data and samples from the Rangewide Bull Trout eDNA project—a collaborative effort that has leveraged crowd‐sourced field sampling by dozens of collaborators across more than 5000 sites. Our ability to engage with multiple groups in this effort was enabled by outreach efforts directed at groups of biologists interested in bull trout conservation as well as by ensuring the results of their sampling would be made available in a user‐friendly, open‐access database that ties species detections to relevant geospatial data (i.e., the environmental covariates used in this study). At the same time, we maintained data quality by controlling where and how samples were collected; that is, all samples were collected based on an ecologically relevant sampling design and used a standard, robust field protocol (Carim et al. ). Thus, this study represents an example of how existing, broad‐scale ecological data, when accessible, can be leveraged to answer new questions (Hampton et al. ). This is particularly relevant to eDNA sampling because samples collected with one taxon in mind are readily archived for re‐analysis to make inferences about the presence of other aquatic taxa (Dysthe et al. ). Here, we did so to determine the distribution of brook trout in the study area at a fraction of the cost and time that would have been required for a new sampling effort. We anticipate that archived eDNA samples will only become more valuable for understanding ecological communities as additional samples accumulate, new species‐specific markers are validated, and improved multi‐taxa analytical approaches are developed (Deiner et al. ).

Acknowledgments

The Rangewide Bull Trout eDNA Project relies on field collection by dozens of collaborators including personnel from the University of Montana, U.S. Forest Service, Trout Unlimited, Clark Fork Coalition, and Montana Fish, Wildlife, and Parks, along with logistical support and genetic analysis at the National Genomics Center for Wildlife and Fish Conservation by staff including Kellie Carim, Tommy Franklin, Caleb Dysthe, and Samuel Greaves. The project was in part funded by Region 1 of the U.S. Forest Service and by the U.S. Fish and Wildlife Service. TMW received support from the NSF (GRFP Grant No. DGE‐1313190) and the University of Montana W.A. Franke College of Forestry and Conservation (Wesley M Dixon Graduate Fellowship and the George & Mildred Cirica Graduate Student Support Fund). This manuscript also benefited from feedback on study design, analysis, and writing by Winsor Lowe.