Riparian Bird Occupancy in a Mountain Watershed in the Colorado Mineral Belt Appears Resilient to Climate-Change-Driven Increases in Metals and Rare Earth Elements in Water and Aquatic Macroinvertebrates

Abstract

Acid rock drainage (ARD) impacts species composition in mountain streams. The potential impact for riparian birds experiencing elevated metal uptake by consumption of benthic invertebrates is concerning but not well studied. We investigated the influence of metal and rare earth element (REE) content in benthic invertebrates on the presence of breeding birds in an ARD-impacted watershed in Colorado, USA, where tree swallows in nest boxes had previously been found to have elevated metal concentrations at some sites. The concentrations of particular REEs in invertebrates were higher than those for cadmium or lead. Avian point counts indicated that most bird species detected were present at most sites, and that tree swallows were rarely found. Occupancy models showed that the availability of shrub or forest habitat was a good predictor for a few habitat-specialists, but metal and REE concentrations in water and invertebrates were not good predictors of avian presence. For other species, neither habitat type nor water quality were important predictors. Overall, this study indicates that the climate-change-driven increases in metals and REEs may not influence the presence of riparian birds in ARD-impacted streams.

1. Introduction

Acidic drainage can occur naturally wherever geology puts sulfidic minerals in contact with air and water and is referred to as acid rock drainage (ARD). Major water quality concerns arise in regions where mining activities greatly accelerate this weathering process, which is referred to as acid mine drainage (AMD) [1]. In the late 1800s, during the Colorado Gold Rush, extensive hard-rock mining took place throughout the Colorado Mineral Belt in the Rocky Mountains (USA). Many mining claims were abandoned just as quickly as they were excavated and piles of waste rock or mine tailings are common throughout the Colorado Rocky Mountains to this day. The US Government Accountability Office has estimated that there are hundreds of thousands of abandoned mines scattered throughout the Western United States [2]. These abandoned mines continue to contribute acidic, metal-enriched water to mountain streams, impairing the aquatic ecosystems and acting as a source of trace metals to drinking-water supplies [3]. Recent studies have found that rare earth elements (REEs) are also elevated in ARD-impacted streams in the Colorado Mineral Belt and that summertime concentrations of both metals and REEs are increasing due to climate change [4].

A pervasive impact of ARD streams is the reduction in species richness, occurring across microbial, algal, and invertebrate communities [5,6]. Genomics methods have shown that microbial taxa in ARD streams are dominated by specialized acidophile organisms [7,8]. Fish and amphibians are often absent due to metal or acid toxicity [9,10]. The drivers for ARD impacts are not only the acidity itself and metal toxicity but also the precipitation of iron and aluminum oxides [11]. In streams with high amounts of metal oxide precipitates covering the streambed, benthic organisms such as primary producers, periphyton, and invertebrates can be smothered completely [12,13]. Ecosystem processes, such as litter decomposition, are also suppressed [14,15]. However, in streams with a low pH and high dissolved metal concentrations, but without precipitates, a few tolerant species of benthic invertebrates may be as abundant as in pristine streams. For example, McKnight and Feder [12] found no difference in the benthic invertebrate abundance in an ARD stream compared to a nearby pristine stream, and the community was dominated by the stonefly Zapada haysii which is common in ARD streams.

Benthic macroinvertebrates are often used as bioindicators of stream health, and studies of acidic drainages have shown that invertebrates accumulate metals in their bodies [16,17]. There are few studies of REE content in benthic invertebrates and no studies of REEs’ accumulation in benthic invertebrates in ARD streams. In a study of temperate lakes in Canada, Amyot et al. [18] found that REE concentrations ranged from 0.1 to 1 ng/g and that REEs were subject to trophic dilution, with lower concentrations at higher trophic levels. The accumulation of metals and REEs in benthic invertebrates is of concern from an ecosystem perspective because of the potential for transport up the food web through bioaccumulation to fish and other predators [19,20].

There has been little investigation into whether elevated metal or REE concentrations in streams affect nearby terrestrial organisms. One potential pathway for these elements to enter terrestrial food chains is the emergence of aquatic insects that spend their larval stages in streams and can be consumed by birds as adult flying insects. Studies have shown that birds acquire elevated metals from prey in a variety of environments, including from aquatic sources [21,22,23]. Environmental acidification and the resulting increased mobility of metals, especially cadmium, mercury, and lead, has been a long-term concern for birds [24,25].

The few studies of ARD effects on terrestrial organisms have addressed concerns related to ingesting elevated levels of toxic metals such as lead (Pb) and cadmium (Cd), which may lead to a lower survival in breeding birds. There is some evidence for reduced density, delayed breeding, and fewer offspring produced by Eurasian dippers (Cinclus cinclus) and Louisiana waterthrushes (Seiurus motacilla) breeding along acidified streams in Wales and Pennsylvania, respectively [26,27]. The most acidic streams in the dipper and waterthrush studies had pH values of 4.5, which is higher than many ARD streams. In a study of tree swallows (Tachycineta bicolor) occupying nest boxes adjacent to ARD streams in the Snake River watershed in Colorado (USA), Custer et al. [28] found that at sites where invertebrates had elevated metal concentrations, nestlings also had elevated, but not toxic, Pb concentrations in liver tissue. It has not yet been determined whether ARD affects the behavior, survival, reproduction, or population dynamics of free-living avian species living adjacent to ARD streams. In a recent study of 35 first- and second-order stream ecosystems impacted by ARD in the Colorado Mineral Belt, Kraus et al. [29] found that adult aquatic invertebrates had lower contents of Cd and Pb than the larval stages living in the streams, reducing the impact of metal transfer to insectivorous riparian spiders. These lower metal contents could potentially also limit the impact of metal transfer to insectivorous riparian birds.

The present study was conducted in the same watershed as the study by Custer et al. [28] and assessed the potential impacts to all riparian birds in a large montane watershed with both pristine and metal-impacted sub-watersheds. The two objectives were to determine (i) both the trace metal and REE contents of benthic invertebrates in the watershed to obtain more complete information on ARD impacts, and (ii) whether birds breeding in riparian forests and wetlands in the watershed are less likely to occupy sites adjacent to the ARD-impacted streams compared to the pristine streams. Given that climate change is expected to continue to exacerbate the impacts of ARD [30], the results of this study help to establish the degree of impact of ARD on both aquatic and terrestrial wildlife.

2. Methods

2.1. Study Area

This study builds on the extensive and ongoing record of water-quality monitoring and hydrologic studies in the Snake River watershed (e.g., [4,12,31,32,33]). For this study, the lowest point of the watershed was defined as the confluence of the Snake River and Peru Creek, forming a watershed with an area of 80 km² (Figure 1). Within each sub-watershed, valleys are dominated by spruce-fir forest, with riparian willow shrublands near the streams. At higher elevations, the dominant land-cover types are alpine turf and exposed bedrock or scree (based on personal observation and the USGS GAP/LANDFIRE National Terrestrial Ecosystems dataset). Runoff in the watershed is dominated by spring snowmelt, as recorded by a United States Geological Survey (USGS) stream gage (09047500) seven kilometers downstream of the study area (Figure S1). Flow in the Upper Snake River is well correlated with flow at the USGS gage [3].

Snake River and Peru Creek are identified on Colorado’s 2022 303(d) list of impaired water bodies due to high concentrations of zinc, cadmium, copper, and lead [34]. The Upper Snake River contains naturally-occurring ARD due to underlying acid-producing schists and gneisses [35]. The pH has ranged between 3.4 and 5.3 since 1980 [3,36]. Monitoring records since 1980 for the Upper Snake indicate that warmer summer temperatures have driven increasing concentrations of sulfate, zinc, manganese, other trace metals, and REEs [3,4]. A “tributary of interest” contributes significantly to the trace metal and REE load in the Upper Snake [4,33]. The elevated concentrations of REEs have made REEs an emerging concern for the watershed [4,37].

In contrast to the Upper Snake, the water quality in Deer Creek is good, despite the presence of abandoned mines [38]. The pH values since 1980 have ranged between 6.2 and 7.6.

At the confluence of the Upper Snake with Deer Creek, the pH increases rapidly and there is subsequent precipitation of aluminum and iron oxides. The pH below the confluence has ranged between 4.8 and 6.8 since 1980 [39,40,41]. Peru Creek, which joins the Lower Snake 4.3 km downstream of the confluence of the Upper Snake and Deer Creek, is heavily impacted by AMD from abandoned mines. The main AMD source in Peru Creek is the Pennsylvania Mine, which was recently remediated to improve water quality.

The Rocky Mountains of Colorado include a wide array of avifauna. In riparian habitats above 3000 m (10,000 feet), common species are the American robin (Turdus migratorius), Wilson’s warbler (Cardellina pusilla), Lincoln’s sparrow (Melospiza lincolnii), and the white-crowned sparrow (Zonotrichia leucophrys) [42]. In spruce-fir forest above 3000 m, common species again include the American robin, as well as the mountain chickadee (Poecile gambeli), golden-crowned kinglet (Regulus satrapa), ruby-crowned kinglet (Regulus calendula), yellow-rumped warbler (Setophaga coronata), dark-eyed junco (Junco hyemalis), and Cassin’s finch (Haemorhous cassinii) [42].

Twenty study sites were established with five sites along each of four main streams in the Snake River watershed: Deer Creek, the Upper Snake, the Lower Snake, and Peru Creek (Figure 1). Sites were numbered from upstream to downstream (i.e., D1 was the highest site on Deer Creek, LS5 was the lowest site along the Lower Snake). Six of the sites corresponded to sites used by Custer et al. [28] and Yang [36]: D1 = Upper Deer (Custer et al. site), D3 = Middle Deer, D4 = Lower Deer, US1 = Upper Snake, US4 = Middle Snake, and LS4 = Montezuma. Using satellite imagery, additional sites were selected to ensure at least 600 m between sites. Given that the home ranges of the birds commonly found in the area are much less than 600 m, this buffer meant that the double-counting of birds was unlikely. A GPS unit (Garmin GPSmap 64) was used for navigation in the field. Final sites were determined based on access and proximity to the stream. The elevations of the sites ranged from 3047 to 3465 m (10,000 to 11,370 feet) above sea level. Photos of select sites are provided in Figure S1. For each site, bird counts were conducted between 23 and 140 m away from the stream. The bird counts could not be conducted directly adjacent to the stream because the sound of running water negatively impacted the observer’s ability to detect birds.

2.2. Sample Collection, Observations, and Analytical Methods

In order to address the first objective of this study, water and benthic macroinvertebrate samples were collected from both acidic and relatively pristine sites and analyzed for concentrations of trace metals and rare earth elements (REEs). Water-quality data in the Snake River, Deer Creek, and Peru Creek were collected at 14 sites (D1, D2, D3, D4, D5, US1, US3, US4, LS1, LS3, LS5, P1, P2, P4; see Figure 1) between 17 July and 26 July 2017. The mean daily discharge at USGS gage 09047500 was between 96 and 124 cfs on the four sampling days, with higher flows during the first week (Figure S1). Ten sites were visited twice and four sites were only visited once due to access constraints. Specific conductance, water temperature, and pH were measured with YSI meters. Water samples (150 mL volume) were collected in acid-washed HDPE bottles, filtered on the day of collection through 0.45 μm glass-fiber filters, and acidified to below pH 2 with Fisher Scientific Trace Metal Grade nitric acid.

Benthic macroinvertebrates were collected with a Surber sampler at all sites except P2, where the stream was too deep and fast for the sampler [43]. The contents of the net were placed in a white tray and macroinvertebrates were picked out, placed in a zip-top bag, and frozen immediately. If the number of invertebrates was low, the Surber sampler was deployed a second time to ensure sufficient mass for analysis, thereby allowing for comparison with previous studies in this watershed that used similar methods. Invertebrate samples were stored in a −10 °C freezer before being freeze dried. The dried invertebrates were separated from any detritus with forceps and the dry mass of each sample was recorded. On two occasions (the second visit to P4 and the only visit to LS5), no benthic macroinvertebrates were detected. The smallest mass of benthic macroinvertebrates was 2.2 milligrams at LS1; the largest was 84.6 milligrams at D4. Due to the low mass of some samples, invertebrates from the two sampling dates for D1, D4, LS1, P1, US1, US3, and US4 were combined, for a total of 13 samples. These samples were then digested in a microwave digestor with Fisher Scientific (Waltham, MA, USA) Trace Metal Grade nitric acid.

Both macroinvertebrate and water samples were analyzed for concentrations of 57 cations, trace metals, and REEs using inductively coupled plasma mass spectrometry (ICP-MS). The concentrations of eight commonly measured trace metals (aluminum [Al], cadmium [Cd], copper [Cu], iron [Fe], manganese [Mn], nickel [Ni], lead [Pb], and zinc [Zn]) and 16 rare earth elements (REEs) are reported here. One water sample and one macroinvertebrate sample were run as duplicates to confirm the analytical error of the ICP-MS, which was also confirmed by comparison to standards. The maximum percent error for the water sample was 4.17% (europium), with an average error of 1.60%. The maximum percent error for the macroinvertebrate sample was 11.11% (lutetium), with an average error of 2.01%. The error values were consistently low; thus, metal concentrations are reported without adjustments for error.

In order to address the second objective of this study, the entire suite of bird species in the watershed was observed to gain a broad picture of the community-wide effects of ARD in the terrestrial system. Avian point counts were conducted during three weeks of the avian breeding season in the summer of 2017. Point counts took place between 9 June and 30 June, 2017, based on nesting dates reported in Yang [36]. Water sampling did not take place concurrently with the avian point counts; during the avian breeding season, the water levels were too high to use the Surber sampler. Each site was visited on three separate days for a total of three counts per site. Surveys took place between ten minutes after sunrise and 10 a.m. During each eight-minute count, all individual birds of all species were recorded that were seen or heard. The avian abundance data were converted to presence–absence detection histories for all species that were detected at more than 4 of the 20 sites (>20% naïve occupancy).

The survey-specific covariates recorded for each count included the Julian date, start time, start temperature (from a Garmin tempe wireless temperature sensor), and categorical indices for wind (0 = smoke rises vertically; 1 = smoke drifts; 2 = leaves rustle, wind felt on face; 3 = leaves and small twigs in constant motion; 4 = raises dust, small branches are moved; 5 = small trees in leaf begin to sway), sky (0 = few clouds; 1 = partly cloudy or variable sky; 2 = cloudy), and noise (1 = quiet; 2 = wind rustles leaves; 3 = distant roar of water; 4 = nearby roar of water). In order to reduce the number of modeled parameters, wind was converted to a binary variable with values of 0 and 1 comprising one category and values of 2 or greater comprising the other. Noise was similarly divided into two categories: values of 1 or 2 comprised the first category and values of 3 or 4 comprised the second.

The site-specific covariates were metrics of either habitat or water quality. Habitat variables were derived from remote-sensing data, specifically the USGS GAP/LANDFIRE National Terrestrial Ecosystems dataset (https://gapanalysis.usgs.gov/gaplandcover/data/, accessed on 1 March 2018). Landcover for the state of Colorado was imported into ArcMap along with the study site locations. The first habitat covariate was the landcover type, reclassified as either forest, shrub, or grass. The percent of forest and the percent of wetland shrub were determined around each point at 100-, 200-, and 300-m scales. Finally, elevation (as recorded on a Garmin GPSmap) was included as a habitat variable, because species richness usually decreases with altitude [44]. The Snake River watershed may be at the altitudinal limit of some species. Water-quality covariates were based on the metal and REE concentrations in water or benthic macroinvertebrates.

2.3. Data Analysis

The effects of dissolved metals on biota can be additive [45], and an additive measure of toxicity is more conservative in that it likely overestimates organisms’ exposure to harmful conditions. For this study, the “cumulative criterion unit” (CCU) was used and was calculated as the ratio of the measured metal concentration to the U.S. Environmental Protection Agency (EPA) criterion value for that metal, summed for all metals [46]:

$$\mathrm{CCU}={\displaystyle \sum }\frac{m_i}{c_i}$$

where mi is the total dissolved metal concentration and ci is the criterion value for the ith metal. CCU values are unitless. Values above 1.0 are likely to cause harm to aquatic organisms [46]. The Colorado Water Quality Control Commission (Regulation No. 33) chronic criterion values were used as the criterion values for Al, Cd, Cu, Mn, Ni, Pb, and Zn, which are hardness dependent and varied with the calcium and magnesium concentrations. The Al standard is for total recoverable metal, but the measurements were made of the dissolved phase and the Al criterion is likely too high. Colorado does not have an Fe standard; thus, the national EPA Criteria Continuous Concentration (CCC), was used—this standard is not hardness dependent and is equal to 1000 μg/L [47]. Chronic values were used instead of acute standards because the chronic values represent a higher threshold, making the CCU a more conservative estimate of toxicity. Only the trace metals were included in the CCU because criterion values for REEs do not exist. Metals that were below the ICP-MS detection limit were not included in the CCU calculation.

In order to evaluate potential spatial autocorrelation, Moran I values were calculated for the dissolved metal and REE concentrations at the stream sites and in the benthic macroinvertebrates, as well as for bird species. The coordinates of the sample sites were used to create a matrix of Euclidean distances for each pair of sites. The inverse of this matrix was used with the Moran I function to determine the Moran I value for each analyte.

A univariate regression was used to determine whether the dissolved metal and REE concentrations were a significant predictor of the invertebrate metal concentration. Separate regressions were conducted for the 9 sites in metal-impacted sub-watersheds and the 5 Deer Creek sites with neutral pH, because these two groups of sites were dominated by invertebrate species in different functional feeding groups. At sites where benthic macroinvertebrate samples from different weeks were combined, the concentrations were averaged across weeks. Data were log-transformed prior to regression. Only regressions with significant p-values (alpha = 0.05) are reported. Regressions, and all other statistical analyses, were conducted in R version 3.4.2 [48].

Occupancy models were used to relate habitat and trace metal and REE exposure variables to the presence of various avian species. In order to avoid the multicollinearity of covariates in occupancy models, aggregate metrics were chosen as site-specific covariates. Principal component analyses (PCAs) were conducted on centered and scaled metal concentrations in water and in invertebrates to reduce the dimensionality of the multivariate dataset while retaining most of the dataset’s original variation. The first principal components (PC1) of dissolved metals and REEs in water samples and all metals and REEs in benthic macroinvertebrates were used as site-specific covariates. The CCU (the additive measure of trace-metal concentrations in water) and the concentrations of Pb in invertebrates were also used as predictor variables. Pb was chosen because Pb was the only metal shown to have some cellular-level effects on nesting tree swallows in the Snake River watershed [28]. Occupancy models cannot use sites with missing site-specific covariates, and water and invertebrates were collected for only 14 of the 20 sites; thus, values for the missing sites were interpolated. Table S1 lists all the candidate covariates used in model selection.

The R packages ‘unmarked’ and ‘AICcmodavg’ were used to build occupancy models for each bird species using an information theoretic approach [49]. A pre-defined list of candidate models was ranked using the Akaike’s Information Criterion (AIC) to account for small sample sizes, such that:

$$AI{C}_c=-2\log \left(L\left(\widehat{\theta }\right)\right)+2K\left(\frac{n}{n-K-1}\right)$$

where L refers to the maximum likelihood for the parameter estimates $\widehat{\theta }$, given the data and model structure; K is the number of estimated parameters; and n is the sample size [49]. Lower AICc values indicate a more parsimonious model. For a given species, candidate models were ranked according to AICc and the difference was calculated between the lowest AICc value and other models’ AICc values, known as the Δ AICc or “AICc difference”. Models with a Δ AICc of less than 2 are considered to have substantial support for being the actual best model, while models with a Δ AICc between 4 and 7 have considerably less support [49]. The weights based on the Δ AICc were calculated such that all the individual weights summed to 1 for a given dataset:

$$w_i=\raisebox{1ex}{$\exp (-\frac{1}{2} {\mathsf{\Delta }}_i)$}\!\left/ \!\raisebox{-1ex}{${{\displaystyle \sum }}_{r=1}^R\exp (-\frac{1}{2} {\mathsf{\Delta }}_i)$}\right.$$

where R is the number of candidate models. Higher weights indicate a higher “weight of evidence” for a model being the actual best model, given the set of candidate models [49].

A hierarchical method was used in which only the survey-specific covariates were modeled before proceeding to the overall best model selection. For the forest and shrub covariates, the best scale was determined, before proceeding to overall model selection. Site-specific covariates that measured similar characteristics (e.g., PC1 of dissolved metals and CCU, or PC1 of invertebrate metals and Pb in invertebrates) were not included in the same model to avoid multicollinearity. Although these covariates measure similar characteristics, it was necessary to test all covariates since, for example, the CCU does not include information on REEs. The “best” models were defined as those with a Δ AICc of less than 3. These “best” models may have unstable parameters or a low predictive skill; therefore, the model-averaged coefficient estimates (betas) and their 95% confidence intervals for parameters found in the “best” models were determined. A variable was considered to be a stable predictor of occupancy if the 95% confidence interval for the beta estimate did not overlap zero. The predictive skill of the best models was also evaluated by conducting a leave-one-out cross-validation (LOOCV) on both the best model and the null model (the model with no survey- or site-specific parameters). The best model was trained on all data points except one; the occupancy of the left-out point was then predicted based on the observed survey- and site-specific covariates and compared to the observed occupancy. This step was repeated by dropping each point in turn to calculate the root-mean-square error (RMSE) for each model.

3. Results

3.1. Trace-Metal and REE Concentrations in Water and Benthic Invertebrates

All results are presented in downstream order, beginning with the Upper Snake, Deer Creek, the Lower Snake, and finally Peru Creek. The complete results for pH, specific conductance, hardness, and dissolved trace-metal concentrations are presented in Table S1. All dissolved REE concentrations can be found in Table S2. The total mass of the combined benthic invertebrate samples and the sampling effort are presented in Table S3. Trace-metal and REE concentrations in invertebrates can be found in Tables S4 and S5, respectively. All of the Moran I values for dissolved trace metals and REEs and for trace metals and REEs in benthic invertebrates were determined to have an absolute value of less than 0.3, and for most of the trace metals and REEs in invertebrates the p values exceeded 0.05 (Tables S6 and S7). Similar to Pearson coefficients, Moran I values can range from 1 to −1 and the low Moran I values for these analytes indicate a random spatial distribution.

Stoneflies in the genus Zapata dominated the Upper Snake and Peru Creek, whereas the most abundant benthic macroinvertebrates in Deer Creek were Baetis mayflies, consistent with previous observations in this watershed (Table S8). Invertebrate biomass in the Lower Snake was extremely sparse. As a result of the difference in the functional feeding group of the dominant species in acidic versus neutral pH streams, these streams were analyzed separately.

In the Upper Snake, the pH values were very acidic: between 3.05 and 3.23. The pH values were slightly basic throughout Deer Creek, ranging from 7.66 to 8.05 (Figure 2). In the Lower Snake, below the confluence of Deer Creek, the pH values were less acidic, returning to circum-neutral below the inflow of Saints John Creek (see Figure 1). The pH values along Peru Creek ranged from 4.31 to 5.27, with higher values below tributary inflows.

The Upper Snake had the highest specific conductance, ranging between 189.5 and 286.5 μS (Figure 2). Specific conductance was much lower in Deer Creek. The Lower Snake had intermediate values between 126.9 and 147.5 μS. Specific conductance in Peru Creek ranged between 167.3 and 230.4 μS.

The CCU values based on the trace-metal concentrations are shown in Figure 3. Dissolved trace-metal and REE concentrations followed a similar pattern throughout the watershed. Figure 4 and Figure 5 present the dissolved concentrations of REEs with high and medium maximum concentrations, respectively. All the REEs in the high category were “light-group” REEs, defined as having an atomic number between 57 and 64.

The concentrations of trace metals in invertebrates are presented in Figure 6, and the concentrations of the more abundant REEs in invertebrates are presented in Figure 7 and Figure 8. In the invertebrates, the concentrations of trace metals and REEs did not follow patterns as clear as were observed in the dissolved concentrations. In general, the REEs that had low dissolved concentrations were also low in invertebrates.

Concentrations of trace metals and REEs in the Upper Snake increased downstream (Figure 4 and Figure 5). A large increase—usually a doubling or tripling in concentration—occurred for Al, Cd, Mn, Ni, Zn, and most REEs between US3 and US4 at the inflow of the “tributary of interest”. The concentrations of many elements were higher at US4 than at any other site. Iron concentrations overall were extremely high and decreased from a maximum of 2508 μg/L at US1 to 1696 μg/L at US4. CCU values averaged 58.5 at US1 and US3, and leapt to 132.8 at US4 (Figure 3).

In the Upper Snake, the Mn, Zn, and most REE concentrations in invertebrates increased downstream, but there was no marked increase at US4, below the tributary of interest (Figure 6, Figure 7 and Figure 8). Meanwhile, Al, Cu, Ni, Pb, Eu, Lu, and Sc peaked at US3 and had lower concentrations at US4. Iron was distinct in that it decreased downstream, with extremely high invertebrate Fe concentrations at US1 of 34,900 μg/g dry mass. Other elements that had maximum concentrations in invertebrates in the Upper Snake were Pb, Ce, Nd, Gd, Sm, and Eu.

In Deer Creek, trace-metal and REE concentrations were relatively low (Figure 3, Figure 4 and Figure 5). Cadmium, Ni, Pb and all REEs were below 1 μg/L; Al, Cu, Mn, and Zn were all below 11 μg/L; and the highest Fe concentration was 124 μg/L. Metal concentrations were relatively constant going downstream. CCU values at all Deer Creek sites were low (Figure 3).

Despite the relatively low dissolved metal and REE concentrations, the highest concentrations of a few trace metals and REEs in invertebrates in the whole watershed occurred in Deer Creek (Figure 6, Figure 7 and Figure 8): Cd peaked at 4.12 μg/g dry mass at D5; Pr peaked at 1.93 μg/g at D1; and La peaked at 8.43 μg/g at D1. Manganese and Zn concentrations in invertebrates were also higher in Deer Creek than in the Upper Snake and REE concentrations were within the same range.

In the Lower Snake, trace-metal and REE concentrations were less than in the Upper Snake (Figure 3, Figure 4 and Figure 5). White aluminum oxide precipitate was observed on the streambed at LS1. CCU values ranged between 17.2 and 39.7 in the Lower Snake.

In the Lower Snake, invertebrate trace-metal and REE concentrations were intermediate or high compared to those for the headwater streams (Figure 6, Figure 7 and Figure 8). Concentrations of several metals (Al, Mn, Cu, and Ni) and REEs (Sc, Y, Dy, Er, Yb, Tb, Lu, Ho, and Tm) were higher in the Lower Snake than at any other study site. Most metals and REEs decreased in concentration between LS1 and LS3. No invertebrates were found at LS5, which was the study site below the inflow of Saints John Creek.

In Peru Creek, the dissolved concentrations of many trace metals and REEs were higher at P1 than at any other site in the watershed (Figure 3, Figure 4 and Figure 5), with Pb concentrations being particularly high. The concentrations of all trace metals and REEs decreased between P1 and P2, where a major tributary joins Peru Creek. P1 had the highest CCU value (95.1). CCU values at P2 and P4 were similar and averaged 46.4.

In Peru Creek, trace-metal and REE concentrations in invertebrates (Figure 6, Figure 7 and Figure 8) were consistently low at P1 and increased at P4, which is the opposite pattern of the dissolved concentrations. Although dissolved Pb concentrations were the highest in Peru Creek, Pb concentrations in invertebrates were not. Zinc was the only metal that had higher invertebrate concentrations in Peru Creek than in other streams.

3.2. Relationships among Trace-Metal and REE Concentrations in Water and Invertebrates

The PCA of dissolved metal and REE concentrations indicated that each of the four basins had distinct characteristics—sites were neatly clustered by basin (Figure 9A). The first principal component explained 82% of the variation in dissolved metal concentrations and separated out the two sites with the highest metal and REE concentrations (P1 and US4) from the other sites. The loadings on PC1 (Figure 9B) were all positive, indicating that the overall magnitude of metal and REE concentrations drove the variation along PC1. The second PC explained 15% of the variation and separated the Upper Snake from the other basins, especially Peru Creek. The metal loadings on PC2 indicated that differences in Fe and Pb influenced the clusters in Figure 9A. Metal and REE loadings are presented in Table S6.

The PCA of trace-metal and REE concentrations in benthic macroinvertebrates did not show a distinct clustering by basin (Figure 10A). The first PC explained 60% of the variation. The loadings on PC1 were all negative, indicating that the overall magnitude of the concentrations drove the variation along PC1. Thus, because the concentrations in Deer Creek and the Upper Snake were similar, the PCA did not separate these two streams. The second PC explained 21% of the variation. The arrows for Fe, Zn, and Mn were the shortest, indicating that they had less of an influence than other elements (Figure 10B). Pb was the trace metal most strongly aligned with the REEs. Metal and REE loadings in invertebrates are presented in Table S7.

The clustering of sites in the PCA for tissue concentrations in benthic invertebrates was not similar to the clustering in the PCA for dissolved concentrations. A notable example is the comparison for Cd and Pb, which are known toxicants. Dissolved lead concentrations were lower than those for Cd, but in invertebrates, Pb concentrations were much higher. For REEs, however, high dissolved concentrations (Figure 4 and Figure 5) generally corresponded to high concentrations in invertebrates (Figure 6 and Figure 7). The exceptions were Gd and Sc, which had medium dissolved concentrations (Figure 5) but were high in invertebrates. In fact, concentrations of Ce, Gd, La, Nd, Sc, and Y in invertebrates were as high or higher than concentrations of Cd and Pb.

Regressions for acidic sites in the Upper Snake, Peru Creek, and the Lower Snake showed that dissolved concentrations were generally poor predictors of those metals and REEs in benthic macroinvertebrates (primarily the shredder Zapata). One exception was Fe (

Table 1. r² and p-values for significant univariate regressions of log-transformed concentrations of metals in water and invertebrates for acidic sites on the Snake River and Peru Creek (n = 7 except for cadmium, where n = 5). Water concentrations were averaged across sites for those sites where invertebrate samples were combined; * indicates a significant p-value (alpha = 0.05).

Metal	Direction of Relationship	Adjusted r2	p-Value
Relationships for acidic sites
Cadmium	Negative	0.829	0.020 *
Iron	Positive	0.717	0.010 *
Lead	Negative	0.516	0.042 *

), with iron concentrations in invertebrates at acidic sites positively correlated with dissolved iron concentrations (p = 0.010). In contrast, the concentrations of Cd and Pb in invertebrates were negatively correlated with the dissolved concentrations of these two metals. This negative relationship may reflect the protective effect of the low pH and high Fe concentrations associated with higher concentrations of Cd and Pb. For all other elements, including REEs, there were no positive or negative significant relationships between the dissolved and invertebrate concentrations (results not shown).

In Deer Creek, where dissolved metal and REE concentrations were low and the dominant invertebrate is a collector–gatherer/scraper (Baetis), there were no elements with significant relationships between concentrations in the water and invertebrates.

3.3. Riparian Bird Distribution and Occupancy Models

Twelve bird species were detected during point counts.

Table 2. Birds detected in the Snake River watershed. Naïve occupancy refers to the ratio of sites where the species was detected in at least one survey. Occupancy models were not built for the greyed-out species.

Common Name	Scientific Name	Naïve Occupancy
American Crow	Corvus brachyrhynchos	0.20
American Robin	Turdus migratorius	0.80
Dark-eyed Junco	Junco hyemalis	0.80
Lincoln’s Sparrow	Melospiza lincolnii	0.55
Mountain Chickadee	Poecile gambeli	0.80
Northern Flicker	Colaptes auratus	0.05
Pine Grosbeak	Pinicola enucleator	0.05
Ruby-crowned Kinglet	Regulus calendula	0.75
Steller’s Jay	Cyanocitta stelleri	0.10
White-crowned Sparrow	Zonotrichia leucophrys	0.65
Wilson’s Warbler	Cardellina pusilla	0.65
Yellow-rumped Warbler	Setophaga coronata	0.65

provides scientific and common names for these species, as well as the naïve occupancy (i.e., the number of sites where the species was detected in at least one survey divided by the number of sites). Eight species were detected at more than 20% of the study sites.

Complete abundance data for all twelve species is provided in Tables S11–S13. The Moran I values for the bird species were all less than 0.065, indicating that the data did not exhibit spatial autocorrelation (Table S14). The American robin and dark-eyed junco were detected at all Deer Creek sites, four sites each in the Upper Snake and Peru Creek, and three sites in the Lower Snake. Lincoln’s sparrow was found at four of five sites along Deer Creek, three in the Lower Snake, and only two each in the Upper Snake and Peru Creek. The mountain chickadee was found at all sites in both Deer Creek and the Upper Snake, but only three sites each along the Lower Snake and Peru Creek. The ruby-crowned kinglet was detected at all Deer Creek sites, four Upper Snake sites, and only three sites each along the Lower Snake and Peru Creek. The white-crowned sparrow and Wilson’s warbler were detected at four of five sites in all basins except Peru Creek, where they were only detected at one site. The yellow-rumped warbler was found at all Upper Snake sites but only three Deer Creek sites, three Peru Creek sites, and two Lower Snake sites.

Values for site-specific covariates and survey-specific covariates are reported in Tables S15 and S16, respectively. Sites along the Snake River and Deer Creek were a mix of forest and wetland shrub habitat, with one high-elevation site (US1) classified as grass. All Peru Creek sites were dominated by forest habitat. The percent of shrub habitat tended to decrease as the scale increased from 100 m to 300 m, whereas forest habitat tended to increase as the radius increased.

Using these habitat parameters, plus the water and invertebrate parameters, occupancy models were built for the eight species detected at more than 20% of the study sites. One clear relationship emerged from these models. Three wetland-shrub habitat specialists were detected in the watershed: Wilson’s warbler, Lincoln’s sparrow, and the white-crowned sparrow. None of these three species were detected at US5, LS5, P2, P4, or P5 (among other sites), which were the only five sites with zero percent shrub within 100 m. Otherwise, these species were commonly detected at Deer Creek, Upper Snake, and Lower Snake sites with shrub habitat, indicating that habitat, but not water quality or invertebrate prey quality, is an important factor for these species during the breeding season. There were no obvious patterns in the occurrence of the other five species. The American robin is a habitat generalist while the dark-eyed junco, mountain chickadee, ruby-crowned kinglet, and yellow-rumped warbler all prefer forest habitats. However, all these species were found at sites with a variety of values for percent forest and shrub. The best models for each species (models with a Δ AICc of less than 3) are listed in

Table 3. Occupancy models with a Δ AICc of less than 3 for each species. Models are named based on their site-specific (Ψ) and survey-specific covariates (p), and Ψ(.), p(.) indicates models with no covariates.

Species	Model Names	AICc	Δ AICc	Weight	Cumulative Weight
American Robin	Ψ(Elevation), p(Stemp)	77.715	0.000	0.825	0.825
Dark-eyed Junco	Ψ(.), p(Stemp)	85.609	0.000	0.318	0.318
	Ψ(.), p(.)	86.950	1.341	0.163	0.481
	Ψ(Elevation), p(Stemp)	87.948	2.339	0.099	0.579
	Ψ(InvertPb), p(Stemp)	88.294	2.685	0.083	0.662
	Ψ(Forest300), p(Stemp)	88.590	2.981	0.072	0.734
	Ψ(Shrub100), p(Stemp)	88.598	2.989	0.071	0.805
Lincoln’s Sparrow	Ψ(Shrub300), p(Julian)	54.974	0.000	0.462	0.462
	Ψ(Forest100), p(Julian)	56.801	1.827	0.185	0.647
Mountain Chickadee	Ψ(.), p(Julian)	67.352	0.000	0.253	0.253
	Ψ(PC1inverts), p(Julian)	68.214	0.861	0.165	0.418
	Ψ(Shrub100), p(Julian)	68.459	1.106	0.146	0.564
	Ψ(Elevation), p(Julian)	68.805	1.453	0.123	0.686
Ruby-crowned Kinglet	Ψ(.), p(Julian + Noise)	81.668	0.000	0.376	0.376
	Ψ(PC1inverts), p(Julian + Noise)	84.364	2.696	0.098	0.473
White-crowned Sparrow	Ψ(Forest100 + InvertPb), p(Julian)	55.342	0.000	0.463	0.463
	Ψ(Shrub100 + InvertPb), p(Julian)	55.348	0.006	0.462	0.925
Wilson’s Warbler	Ψ(Shrub300), p(.)	58.425	0.000	0.417	0.417
	Ψ(Forest100), p(.)	58.453	0.028	0.411	0.828
Yellow-rumped Warbler	Ψ(.), p(Noise)	80.543	0.000	0.182	0.182
	Ψ(.), p(.)	80.747	0.204	0.164	0.346
	Ψ(Forest300 + CCU), p(Noise)	81.047	0.504	0.141	0.487
	Ψ(Forest300), p(Noise)	81.159	0.617	0.133	0.620
	Ψ(CCU), p(Noise)	81.619	1.077	0.106	0.726
	Ψ(PC1water), p(Noise)	82.484	1.941	0.069	0.795
	Ψ(PC1inverts), p(Noise)	82.659	2.117	0.063	0.858
	Ψ(Elevation), p(Noise)	83.096	2.554	0.051	0.909
	Ψ(Shrub300), p(Noise)	83.374	2.831	0.044	0.953

, and model-averaged estimates of the coefficients (betas) for parameters in those models are found in

Table 4. Modeled-averaged parameter coefficient estimates from occupancy models. Both site-specific (Ψ) and survey-specific (p) covariates are included. Only estimates of parameters in occupancy models with a Δ AICc of less than 3 are included. An asterisk (*) indicates a stable estimate, i.e., the confidence interval does not overlap zero.

Species	Parameter	Model-Averaged Beta Estimate	95% Confidence Interval
American Robin	Elevation (Ψ)	34.37	−128.79–197.54
	Stemp (p)	0.07	−0.01–0.14
Dark-eyed Junco	Elevation (Ψ)	1.22	−2.19–4.64
	InvertPb (Ψ)	−0.09	−0.42–0.24
	Forest300 (Ψ)	−0.02	−0.17–0.12
	Shrub100 (Ψ)	−0.02	−0.11–0.07
	Stemp (p)	−0.07	−0.15–0
Lincoln’s Sparrow	Shrub300 (Ψ)	0.15	−0.01–0.3
	Forest100 (Ψ)	−0.03	−0.08–0.03
	Julian (p)	0.33 *	0.06–0.59
Mountain Chickadee	PC1inverts (Ψ)	0.78	−0.92–2.48
	Shrub100 (Ψ)	−0.06	−0.18–0.06
	Elevation (Ψ)	1.03	−0.85–2.91
	Julian (p)	−0.23 *	−0.35–−0.11
Ruby-crowned Kinglet	PC1inverts (Ψ)	−0.25	−0.81–0.32
	Julian (p)	0.1	−0.01–0.21
	Noise1 (p)	−2.21 *	−3.85–−0.58
White-crowned Sparrow	InvertPb (Ψ)	3.72	−7.38–14.82
	Shrub100 (Ψ)	3.4	−6.25–13.01
	Forest100 (Ψ)	−1.52	−6.42–3.38
	Julian (p)	0.2 *	0.06–0.34
Wilson’s Warbler	Shrub300 (Ψ)	2.38	−50.72–55.48
	Forest100 (Ψ)	−0.8	−17.84–16.24
Yellow-rumped Warbler	Forest300 (Ψ)	−0.08	−0.22–0.06
	CCU (Ψ)	0.05	−0.06–0.15
	PC1water (Ψ)	0.33	−0.68–1.34
	PC1inverts (Ψ)	0.50	−0.97–1.97
	Elevation (Ψ)	0.60	−0.82–2.03
	Shrub300 (Ψ)	0.05	−0.12–0.22
	Noise1 (p)	−1.35	−2.96–0.25

. For certain species, one or two models were clearly better than the others. For others, the AICc weights offered support for a suite of models. The site-specific covariates that appeared in the set of best models varied by species and included both water quality and habitat variables. The 95% confidence intervals for all site-specific covariates, for all species, included zero, indicating that none of these variables were stable predictors of occupancy. However, this outcome may be due to the small sample size. Despite having wide confidence intervals, many of the parameter estimates make sense based on the known habitat preferences for particular species. Finally, the leave-one-out cross-validation (LOOCV) results showed that the best models for four of the eight species did a better job at predicting occupancy than the null models for those species (

Table 5. Leave-one-out cross-validation model skill. Lower root-mean-square error (RMSE) indicates that the model did a better job of predicting the occupancy of sites not included in the model’s training data. An asterisk (*) indicates that the best model was better than the null model, which has no survey- or site-specific covariates.

Species	RMSE of Null Model	Best Model	RMSE of Best Model
American Robin	0.443	Ψ(Elevation), p(Stemp)	0.229 *
Dark-eyed Junco	0.436	Ψ(.), p(Stemp)	0.442
Lincoln’s Sparrow	0.524	Ψ(Shrub300), p(Julian)	0.447 *
Mountain Chickadee	0.423	Ψ(.), p(Julian)	0.424
Ruby-crowned Kinglet	0.465	Ψ(.), p(Julian + Noise)	0.469
White-crowned Sparrow	0.510	Ψ(Forest100 + InvertPb), p(Julian)	0.358 *
Wilson’s Warbler	0.505	Ψ(Shrub300), p(.)	0.004 *
Yellow-rumped Warbler	0.507	Ψ(.), p(Noise)	0.533

).

4. Discussion

Overall, the water quality in the study area was similar to many prior studies of this watershed [4,12,28,31,36,37]: acidic sub-watersheds (the Upper Snake and Peru Creek) contained high levels of dissolved trace metals and REEs, while Deer Creek contained low levels of dissolved metals and REEs. As in previous studies, the dominant benthic macroinvertebrate species varied by basin [12,36], which may influence the differences across the sub-watersheds because trace-metal concentrations in freshwater benthic invertebrates have been found to vary with functional feeding guilds [50]. Stoneflies, which dominate the Upper Snake and Peru Creek basins, are shredder–detritivores that consume leaves that have fallen into the stream and have been colonized by bacteria and fungi [51]. As such, the stoneflies may avoid consuming elevated metal concentrations, as demonstrated by the finding that trace-metal and REE concentrations in invertebrates were not well correlated with concentrations in the water. Mayflies, which dominate the relatively pristine Deer Creek basin, are collector–gatherer/scrapers that consume periphyton biomass that may be enriched in trace metals. These mayflies are rarely found in streams with elevated metals, so this taxon is usually considered to be sensitive to poor water quality. However, there were no significant relationships between metal and REE concentrations in the water and invertebrates at Deer Creek sites.

The comparison of invertebrate trace-metal concentrations in this study conducted in 2017 with data for five of the sites studied in 2004 by Yang [36] does show that the concentrations of Al, Cu, Fe, and Ni increased at all five sites. Al concentrations were over an order of magnitude higher in both Deer Creek and the Upper Snake in 2017. Increases in Deer Creek over the thirteen years were smaller, and Cd, Mn, Pb, and Zn concentrations increased at some sites but not others. One of the main patterns that Yang [36] found was a large increase in invertebrate trace-metal concentrations at the Middle Deer site compared to the surrounding sites in Deer Creek. In 2017, Al, Fe, and Mn also peaked at this site (D3) compared to D1 and D4, but other trace metals had decreased. Despite these increases, which may be associated with the increasing concentrations of trace metals over time in the watershed [3,4], it is important to consider that the benthic invertebrate taxa currently present in the watershed may be adapted to these high concentrations of trace metals and REEs that may have occurred during warm periods in the past.

In the context of potential adaptation, it is important to note that the concentrations of REEs in benthic invertebrates found in this ARD-impacted watershed are much higher than other reported values. In comparison to the REE concentrations found by Amyot et al. [18] for invertebrates in pristine temperate lakes in Canada, the REE concentrations for all the stream reaches in the Snake River watershed are 3–4 orders of magnitude higher. Furthermore, the REE concentrations found in this study are 1–2 orders of magnitude higher than the concentrations reported by Pastorino et al. [17] for invertebrates in six stream sites in Italy receiving inputs from nearby agricultural and industrial activities. However, similar to the findings of Pastorino et al. [17], the light REEs generally exhibited higher concentrations than the heavy REEs.

There may be a potential for the bioconcentration of REEs in the aquatic ecosystems of the Snake River watershed. Although there are no native fish populations present in the Snake River, Deer Creek does support native fish and the Lower Snake is occasionally stocked with rainbow trout as a recreational fishery [10]. The concentrations of some REEs in invertebrates were as high or higher than those for Cd or Pb, which are known toxicants; thus, REEs are an emerging concern beyond the widely recognized impacts of ARD.

Of the eight avian species detected at more than 20% of the study sites, the mountain chickadee is the only resident species, although it may move to lower elevations during the winter. The other seven species all migrate south in the fall and arrive back in the Rocky Mountains between mid-March and late-May [52]. All species are reported to eat at least some arthropods, including insects, during the breeding season. The home ranges for these eight species are generally much smaller than the 600-m buffer between the sites. For example, mountain chickadees and white-crowned sparrows have territories of 6–7 ha, which is equivalent to a circle with radius ~150 m [53,54]. Sabo [55] found that territories for the ruby-crowned kinglet, yellow-rumped warbler, and dark-eyed junco ranged from 0.6 to 1.1 ha.

Occupancy models for eight insectivorous avian species provided little evidence that bird presence in the Snake River watershed was influenced by water quality or metal and REE concentrations in invertebrates, which represent potential prey during the breeding season. This is despite large variation in habitat and water-quality parameters across the study area. For three species (Lincoln’s sparrow, Wilson’s warbler, and the white-crowned sparrow), the best models did predict occupancy better than the null model. All three of these species are habitat specialists—they only breed in willow (Salix spp.) wetlands such as those found in parts of the Snake River watershed [56,57,58]. Therefore, it makes sense that the best models included percent shrub or percent forest landcover. For four species (dark-eyed junco, mountain chickadee, ruby-crowned kinglet, yellow-rumped warbler), the best models did not predict occupancy any better than the null models. This suggests that the factors that are influencing the occupancy of these species were not included in the candidate models, potentially due to the widespread distribution of these species in the study area. The best model for the American robin also predicted occupancy better than the null model, but this may have been an artefact of a small sample size; the model suggested that robin occupancy increases with increasing elevation, but robins are known to thrive at all elevations up to those found in the study area [42].

The only potential indicator of a relationship between avian occupancy and metal concentrations was the best model for the white-crowned sparrow. In addition to percent forest cover within 100 m, the best model for this species included lead concentrations in invertebrates. Given that Custer et al. [28] found that Pb concentrations in the livers of tree swallows had been high enough to inhibit some enzyme activity, the presence of invertebrate Pb concentrations in the best model supports the further study of Pb bioavailability and bioaccumulation.

Besides the white-crowned sparrow, other species did not show any relationship between occupancy and metal and REE concentrations in aquatic invertebrates. This finding, suggesting a limited impact on riparian birds, is consistent with the results of the study by Kraus et al. [29] which included 35 streams and showed that concentrations of Cd and Pb were generally lower in adult aquatic insects compared to the aquatic insect larvae—thereby limiting the trophic transfer to predatory riparian spiders. Furthermore, the results of this study provide a broader context for interpreting the results of a previous study in the Snake River watershed, which examined birds that used nest boxes set up for the study [28]. Ninety percent of the birds using the nest boxes were tree swallows (Tachycineta bicolor) or violet–green swallows (T. thalassina) [28]. The authors of the previous study found that benthic macroinvertebrates contained elevated concentrations of metals at one site on Deer Creek and that tree swallows did not use nest boxes at that site, indicating a possible link between water quality and avian habitat selection [28]. However, no swallows were observed anywhere in the study area in 2017. Tree swallows are rarely found above 9000 feet in Colorado [42], but they have been known to breed in nest boxes outside their typical range [59]. Thus, the study area may be at too high an elevation for swallows to breed without the incentive of a nest box. This is relevant because swallows primarily obtain food by capturing insects while flying. All of the songbirds considered in the occupancy models of the current study primarily obtain food by gleaning insects from vegetation or the ground [52]. Due to these behavioral differences, it is likely that most of the biomass consumed by swallows is of aquatic origin, i.e., insects that have emerged directly from the stream [60]. The songbirds detected in 2017 have a more diverse diet that includes spiders, caterpillars (Lepidopetra larvae), aphids, and some vegetable matter (e.g., [55,61]).

The metal concentrations of non-aquatic avian prey items in the study area are not known but it is certainly possible that they contain lower metal concentrations than aquatic invertebrates. Therefore, by eating a variety of invertebrate prey, as well as seeds and berries, songbirds may be exposed to lower metal concentrations than aerial insectivores such as swallows. On the other hand, studies have shown that predatory invertebrates such as spiders are important vectors for the trophic transfer of metal contaminants. For example, Kraus et al. [62] showed that spiders had higher concentrations of zinc, copper, and cadmium than their aquatic prey. However, because insect emergence decreases as metal contamination in streams increases, especially for mayflies, the total metal exported by insects is the lowest at the most contaminated sites [62]. These studies suggest that understanding avian exposure to contaminants requires a complete picture of the food chain.

There are many other possible reasons beyond a broad diet that water quality, specifically elevated metal and REE concentrations in streams and benthic macroinvertebrates, would have no effect on bird occupancy. First, not all metals are harmful. Of the elements reported in this study, only Pb and Cd are known to be harmful to birds [28], although data on the toxicity of rare earth elements are scarce. Metals such as Cu, Fe, Mn, and Zn are essential or beneficial for birds and are easily metabolized. Additionally, birds may not be exposed to elevated metals for enough time to feel the effects. Except for the mountain chickadee, all of the modeled species are migrants that spend at most half a year in the study area. The winter season could be enough time for birds to metabolize the metals and remove them from their systems. However, philopatry is common in birds [63], so it is probable that most individuals return to the same breeding grounds year after year and are potentially repeatedly exposed to contaminants. Certain species are also arriving to breeding grounds earlier each year due to changing climate cues at low elevations [64], meaning that birds could be exposed to contaminants for a longer time, increasing the chances of a chronic toxicity.

There are also many reasons that birds breeding adjacent to ARD streams could, in fact, be experiencing negative effects that would not be detectable from occupancy modeling. For example, Mulvihill et al. [27] found that Louisiana waterthrushes were still present and breeding along ARD streams, but that their territories had to be larger to account for the unfavorable habitat quality. In addition, while fecundity was the same at ARD sites and control sites, birds at ARD sites produced smaller clutches and nestlings grew more slowly [27]. The point count data obtained in the current study would not have captured all these effects.

5. Conclusions

In the Colorado Mineral Belt and other mountainous regions worldwide, acid rock drainage (ARD) has a pervasive and deleterious influence on water quality. In this study, the uptake of both trace metals and REEs by benthic invertebrates in ARD-impacted streams was considered as a potential factor influencing the presence of breeding birds in the Snake River watershed. One important and novel finding is that the concentrations of some REEs in invertebrates were as high or higher than those for Cd or Pb, which are known toxicants. Furthermore, the PCA indicated that REE concentrations in invertebrates were aligned with Pb concentrations in invertebrates. Overall, given the dearth of information about REEs in aquatic organisms, further study of the availability and transfer of REEs through ARD-impacted aquatic ecosystems is warranted.

Another important finding was that dissolved metal and REE concentrations were not good predictors of concentrations in benthic macroinvertebrates, even when acidic and pristine sites were considered separately. This lack of a relationship is likely due to differing exposure and uptake by the resident species of benthic macroinvertebrates. Furthermore, while the benthic invertebrates present in the ARD-impacted streams did have high concentrations of metals and REEs in their tissues, this study did not find evidence of an impact on the presence of riparian songbirds in the watershed, although samples sizes were small. The broad diet of the songbirds and potentially lower metal concentrations in the adult invertebrates compared to the aquatic juvenile insects may account for this result. The monitoring of ARD sites during and after remediation often focuses on the stream environment [65]. The finding of a limited impact of ARD on riparian birds may be useful in assessing whether or not monitoring of these sites during remediation of mining impacts may need to expand to include an evaluation of the metal and REE contents of riparian birds and other terrestrial wildlife.