Complementary responses of stream fish and benthic macroinvertebrate assemblages to environmental drivers in a shale-gas development area

Sampling was conducted in the Kennebecasis and Pollett watersheds in southeastern New Brunswick, Canada (Fig. 1). The area is primarily forested, with some agriculture and a small percentage of land cover composed of settlements, infrastructure, and industry, including shale-gas exploration and extraction. Headwater regions of the Kennebecasis and Pollett rivers comprise morainal and colluvial surficial geology deposits, whereas the mid to lower reaches are dominated by morainal, glaciofluvial, and alluvial deposits (Rampton 1984). Bedrock geology in the western part of the study area is characterized by Early Carboniferous sedimentary bedrock (Mississippian era) underlain by the Lower Carboniferous Albert formation, which has potential for conventional natural gas and unconventional (shale gas) extraction (Macauley et al. 1984, 1985; Leblanc et al. 2011). The eastern part of the study area is dominated by Late Carboniferous sedimentary bedrock (Pennsylvanian era), which differs from oil shale bedrock with respect to mineral content and organic composition (Dyni 2003). Older classes of bedrock, including sedimentary Devonian-Carboniferous bedrock and volcanic and intrusive bedrock from the Neoproterozoic era, are found along the southern boundary of the study area (see Lento et al. 2019 for a full description of bedrock geology in the study area).

An extensive spatial survey was completed in late August and early September 2015 to characterize fish and BMI assemblage structure, water chemistry, and physical habitat in the Kennebecasis and Pollett watersheds (n = 18 and 6 stations, respectively; Fig. 1). Selection of sampling areas utilized existing geographic information system (GIS) data for topography, geology, and land use information, complemented by local knowledge. Site selection followed a stratified random sampling design to ensure coverage over a variety of catchment sizes across a natural habitat gradient. Mainstem and tributary stations were selected to cover a range of sub-catchment areas from 7.25 to 292.73 km² (Table 1). Stations were classified by underlying geological age as Early Carboniferous, Late Carboniferous, or Older Classes (Devonian-Carboniferous and Neoproterozoic) based on the dominant age of geology in a 1-km upstream catchment buffer (see Lento et al. 2019). Several of the Early Carboniferous stations were located in the vicinity of oil and gas wells (Fig. 1), though active fracking activities were halted in the area in 2014 after the provincial government imposed a moratorium, and many wells were out of operation prior to the moratorium. Groundwater methane levels measured in private wells in the area were elevated compared to surrounding areas, but there was no correlation with proximity to oil and gas wells, which indicated that these were naturally high methane levels (Loomer et al. 2016). Furthermore, spills associated with unconventional shale-gas production have not been reported for the area, and chloride levels in surface waters were not correlated with proximity to oil and gas wells; thus, stations were considered to be least disturbed condition sites (sensu Stoddard et al. 2006).

Water grab samples were collected following provincial guidelines and were analyzed at the New Brunswick Department of Environment and Local Government’s analytical lab. Samples were analyzed for a standard suite of metals, nutrients, and ions following standardized analytical procedures (Table S1). A habitat survey was conducted at each station following protocols of the national Canadian Aquatic Biomonitoring Network (CABIN; Environment Canada 2012), including measurement of stream width, depth, velocity, and bed slope; description of in-stream and riparian vegetation and periphyton coverage; and a rock walk (Table S2). During the rock walk, 100 particles were chosen at random from within the sampled habitat, and the length of the intermediate axis length was recorded in mm (after Wolman 1954). Substrate measurements were grouped by grain size (using the Wentworth scale; see Environment Canada 2012) to quantify the percent composition of different size classes (e.g., sand, gravel, cobble) and median particle size (D₅₀).

BMI samples were collected following standard CABIN protocols (Environment Canada 2012). The operator held a 400-μm mesh kick-net downstream while travelling through the stream channel in a zig-zag fashion and disturbing the substrate for a period of 3 min (Environment Canada 2012). BMI samples were randomly subsampled in the laboratory using the Marchant box method (Environment Canada 2014) until a minimum of 300 organisms was counted. Non-Chironomidae (midge) taxa were sent to an independent Society for Freshwater Science (SFS)-certified taxonomist for identification to genus level where possible, and Chironomidae taxa were sent to EcoAnalysts, Inc. to be identified to species by SFS-certified taxonomists (see Lento et al. 2019 for more details). Samples contained large numbers of taxa that were too small (early instars) to identify to the level of genus and to avoid mixed-level taxonomy within an order (which may lead to misrepresentation of trends), data were rolled up to the level of subfamily for Chironomidae and the level of family or higher for all other organisms. Bowman and Bailey (1997) demonstrated that analysis of BMI communities at the level of genus did not improve the ability to detect patterns in the data.

Fish community surveys were conducted by using a backpack electrofishing unit (Smith-Root LR-24, Vancouver, WA, USA) to fish an extent of the stream site containing at least one riffle and one pool or run in a single upstream pass (stream area ranged from 50 to 830 m²). The effort was standardized by maintaining power within a range of 65–75 W (voltage varied depending on conductivity of the water) and sampling time ranged from 419 to 3009 s. All fish were identified to species level.

Temperature loggers (Onset HOBO UA-001-08, UA-002-64, and U24-001, Bourne, MA, USA) were deployed in stream sites in May 2015 to capture longer-term temperature trends. Loggers were anchored by rebar driven into the substrate and were housed in light gray PVC pipe to avoid effects from solar radiation and shield them from damage caused by shifting substrate. The loggers recorded data every 30 min from May to November. Logger data were examined for records that indicated the logger was out of water (e.g., spikes in temperature or light levels, or decrease in conductivity to near 0 for the combination temperature/conductivity loggers). These records were removed, and remaining data were used to calculate summary temperature metrics for the summer period (June–August), including the number of degree days above 19 °C and the coefficient of variation (CV) of temperature. Catchment area upstream of each station was delineated by using the Hydrology tools in the Spatial Analyst toolbox to analyze a 30 m resolution Digital Elevation Model (DEM; Canada Digital Elevation Data; open.canada.ca/data) in ArcGIS (version 10.3, ESRI, St. Paul, Minnesota, USA). Additional geospatial variables were considered for analysis, including catchment vegetation and slope, but in-stream measures were retained to reduce the number of variables for the analysis. Surficial geology was also characterized from available geospatial layers but was correlated with the bedrock geology classes and therefore removed from the analysis to avoid masking differences in water chemistry and physical habitat variables among geologic age classes.

The PCA of water chemistry and habitat variables identified clear gradients in the data that were driven by a contrast between stations with high levels of ions, nutrients, and suspended sediment particles and stations that were more strongly associated with habitat descriptors. The first axis of the PCA, which explained 31.9% of the variance among stations, separated most Early Carboniferous stations from those underlain by Late Carboniferous or older classes of bedrock, as Early Carboniferous stations were generally positively correlated with conductivity, alkalinity, hardness, and ions (including Cl, Na, K, Ca, and SO₄) as well as turbidity and the habitat variable bankfull depth (Fig. 2). On the opposite end of the first axis gradient, stations underlain by Late Carboniferous or older classes of bedrock were associated with a number of habitat variables, including higher D₅₀ (large median particle size), larger bankfull width, presence of coniferous trees, more variable temperature (higher CV), increased dominance of deciduous trees, and higher cover of periphyton and boulders (Fig. 2).

The second axis of the PCA explained 15.8% of the variance among stations and contributed to separation within each bedrock geology age class (Fig. 2). The positive end of the gradient was associated with degree days above 19 °C, average depth, catchment area, nutrients (TN, TOC, nitrate), colour, and aluminum (Fig. 2). Several Early Carboniferous stations (ST2, KB7, KB6, and PL4) followed a gradient in water chemistry along this axis. On the negative end of the gradient, stations were associated with high slope, dominance of coniferous vegetation, and % silt/clay (Fig. 2). Stations such as TC1, TC2, and KB4 (Older Class, Early Carboniferous, and Late Carboniferous, respectively) were negatively correlated with parameters such as TN, TOC, TP, colour, and turbidity and associated total metals such as Al and Fe (Fig. 2). In contrast, some stations on the Pollett River (PL1, PL2, and PL6) were positively associated with TN, TOC, and colour along the second axis while still reflecting a negative association with ions along the first axis. The third axis (not shown), which explained 11.2% of the variance among stations, was positively correlated with TP (associated with Early Carboniferous stations) and periphyton cover (associated with Late Carboniferous stations).

Multivariate assessment of BMI community structure indicated a strong separation of stations along the first and second axes that was similar to patterns in the water quality and habitat PCA. The first axis, which explained 31% of the variance in BMI samples, separated the Older Class stations and a small number of Early and Late Carboniferous stations with higher bankfull width and larger substrate size (stations PL1, PL2, PL6, TC1, and TC2) from the majority of other stations (Fig. 3A). Stations that were positively correlated with the first axis were associated with Odonata and nine families of Trichoptera (Apataniidae, Helicopsychidae, Hydropsychidae, Lepidostomatidae, Leptoceridae, Odontoceridae, Philopotamidae, Polycentropidae, and Psychodidae) whereas only four Trichoptera families were strongly associated with the negative end of the gradient (Brachycentridae, Glossosomatidae, Hydroptilidae, and Rhyacophilidae; Fig. 3A). Plecoptera families appeared evenly distributed along the length of the first axis (including Cloroperlidae, Perlidae, Perlodidae, and Pteronarcyidae on the positive end and Capniidae, Leuctridae, and Nemouridae on the negative end; Fig. 3A), suggesting a wide range of preferences to the underlying environmental gradient. The second PCA axis, which described 16.9% of community variance, was positively correlated with a group of predominantly Early Carboniferous stations that had high conductivity and high concentrations of ions, and was negatively correlated with most Late Carboniferous stations (Fig. 3A). Along the second axis, a large number of taxa were positively associated with the Early Carboniferous stations, including Coleoptera (Elmidae), several families of Ephemeroptera (including Ephemerellidae, Leptohyphidae, and Leptophlebiidae) and Trichoptera (Goeridae, Glossosomatidae, Hydropsychidae, Psychomiidae, and Uenoidae), and most noninsect groups, such as Oligochaeta, Bivalvia, Hirudinea, and Acari (Fig. 3A).

The PCA based on fish relative abundance explained a great deal of variation in the community data, but did not appear to distinguish strongly among bedrock geology age classes (Fig. 3B). The first axis, which explained 65.6% of the variance in the fish community, was positively associated with Pollett River stations PL1, PL2, and PL6, and more weakly correlated with the remaining stations (Fig. 3B). This axis was primarily driven by a negative association between Slimy Sculpin (Cottus cognatus; on the negative end of the gradient) and Blacknose Dace (Rhinichthys atratulus; on the positive end of the gradient), though other stations on the positive end of the gradient were also positively correlated with taxa such as American Eel (Anguilla rostrata), White Sucker (Catostomus commersonii), and Common Shiner (Luxilus cornutus; Fig. 3B). The second axis explained much less community variance (16.7%) but represented a gradient in Brook Trout (Salvelinus fontinalis) abundance along which most stations differed. Brook Trout were positively correlated with Threespine Stickleback (Gasterosteus aculeatus) but were weakly or uncorrelated with other fish species (Fig. 3B).

Procrustes analysis indicated that the fit of the BMI PCA (rotational ordination) to the fish PCA (target ordination) was more similar than could be obtained by chance ( $m_{12}^2$  = 0.56; p = 0.001 after 999 permutations). The largest residual difference was found for station PL2 (residual vector length = 0.313), and two other stations had moderately long residual vectors (e.g., PL1 and PL5; Fig. 4), suggesting that the similarity of these stations to other stations differed depending on whether fish or BMI data were considered. However, residual vectors were short for most stations, indicating significant community concordance that was evident as a strong overall similarity between the fish ordination and the BMI ordination.

Baseline criteria were established through analysis of confidence intervals in water quality and biotic metrics, with the differences predominantly found between Early Carboniferous and Older Class stations. Conductivity, which was strongly associated with Early Carboniferous stations in the abiotic PCA, was higher on average in Early Carboniferous stations (mean = 167.9 μS/cm) than in Older Class stations (mean = 41.62 μS/cm), and 95% confidence intervals for these geologic age classes did not overlap (Table 4). Similarly, TP was higher on average in Early Carboniferous stations (mean = 0.013 mg/L) than in Older Class stations (mean = 0.003 mg/L) and confidence intervals did not overlap (Table 4). Nitrate differed significantly between Early and Late Carboniferous stations (results not shown), but the 95% confidence intervals of Older Class stations overlapped with both other classes. Estimates of the normal range for other water quality parameters overlapped between Early Carboniferous, Late Carboniferous, and Older Class stations, which indicated that there was insufficient distinction among geologic classes to allow their use for establishing baseline criteria from the normal range.

There were a number of BMI metrics for which it was possible to develop baseline criteria from the normal range, with all differences found between Early Carboniferous and Older Class stations. Total abundance of BMI was highest on average in Early Carboniferous stations (mean = 25 046.2) and the 95% confidence interval for this geologic class did not overlap with Older Class stations (mean = 8534; Table 4). Taxonomic richness was also higher on average in Early Carboniferous stations (mean = 37.9) than in other geologic classes, and the 95% confidence interval did not overlap with that for Older Class stations (mean = 30.2; Table 4). The remaining metrics that were used to establish baseline criteria reflected the stronger association of Early Carboniferous stations with noninsect taxa and weaker association with EPT taxa in the BMI PCA. The % EPT was higher on average in Older Class stations (mean = 73.5%) than other geologic age classes, particularly Early Carboniferous (mean = 47.4%), and confidence intervals for Early Carboniferous and Older Class stations did not overlap (Table 4). In contrast, the % Chironomidae (the other dominant taxonomic group in these systems) was significantly higher on average in Early Carboniferous stations (mean = 34.5%) than Older Class stations (mean = 18.6%), as was the % Oligochaeta (Early Carboniferous mean = 1.25% and Older Class mean = 0.19%; Table 4).

In contrast to the BMI, there were no fish metrics that differed significantly among geologic age classes, which supported the generally weaker differentiation of Early Carboniferous stations that was evident in the fish PCA. Although species richness was higher in Early Carboniferous stations than Older Class stations, the 95% confidence intervals for these two geologic age classes indicated no significant differences (Early Carboniferous confidence interval range: 4.1–6.8; Older Class confidence interval range: 0.7–5.3). It was therefore not possible to establish baseline criteria for fish metrics to distinguish among geologic age classes and characterize a distinct normal range for Early Carboniferous stations.

Forward selection of variables for Redundancy Analysis resulted in similar sets of environmental variables for analysis of BMI and fish (8 of the 9 final variables were the same in both RDAs), though the relative importance of those variables differed. Notably, the set of variables that was significant in the RDA model differed, as BMI were significantly (p < 0.05) associated with (in order of importance) temperature (degree days over 19 °C), conductivity, and catchment area, whereas fish were significantly associated with catchment area, turbidity, and TN (Table 5). Conductivity, which was strongly associated with Early Carboniferous stations and significant in the BMI RDA, was included in the final fish model, but was one of the weakest variables in the model (Table 5). Average velocity and TP were not included in the final model for either BMI or fish, and only the fish RDA included a substrate size variable (% pebble), whereas slope was only selected in the BMI RDA (Table 5).

The spatial arrangement of stations and BMI taxa in constrained ordination space was similar between the RDA biplot and the BMI PCA, which was indicative of the strong contribution of the selected water chemistry and habitat variables to the environmental gradient underlying patterns in BMI community composition. The first axis of the RDA explained 14.6% of the unconstrained variance in BMI assemblage composition, and the second axis explained an additional 8.6% of the unconstrained assemblage variance (Fig. 5A). Several stations in the Pollett River catchment (PL1, PL2, PL4, and PL6) were positively associated with higher temperatures, larger catchment area, higher average depth, and higher TN along the first axis of the RDA (Fig. 5A). A large number of Late Carboniferous stations were negatively correlated with this gradient and with the second axis, and these stations were associated with lower temperatures, smaller catchment area and depth, lower TN, high periphyton coverage and high slope (stations CB1, PL5, ST1, KB4, and SB3A; Fig. 5A). In contrast, six of the Early Carboniferous stations were associated with high conductivity and turbidity along the second axis gradient. Plecoptera and Trichoptera families were distributed along the first axis, as in the BMI PCA, whereas the Early Carboniferous stations on the positive end of the second axis gradient (associated with high conductivity and turbidity) were still positively associated with noninsect groups and a small number of insect families (Fig. 5A).

Separation of stations in the fish RDA showed a weaker grouping of geological age classes, with most classes distributed across the ordination (Fig. 5B). The first axis, which explained 17.4% of the constrained variance, was primarily driven by a gradient in temperature, catchment area, depth, and TN, as in the BMI RDA. However, Pollett River stations were largely uncorrelated with this axis, and the gradient was driven by a negative association of these variables with stations such as ST1, PL5, SB1, MB1, and CB1 (Fig. 5B). Brook Trout, and to a lesser extent Slimy Sculpin, Ninespine Stickleback (Pungitius pungitius), and Blacknose Shiner (Notropis heterolepis) were negatively associated with system size and temperature along the first axis gradient. The second axis indicated a gradient in conductivity and turbidity, both of which were positively associated with a number of Early Carboniferous stations (Fig. 5B). On the other end of that axis, stations in all geological age classes were positively associated with the presence of coniferous vegetation or TN. Brook Trout and Threespine Stickleback were positively correlated with conductivity along the second axis, whereas Burbot (Lota lota), Lake Chub (Couesius plumbeus), and Sea Lamprey (Petromyzon marinus) were positively associated with turbidity (Fig. 5B).

A concordance gradient was created by aligning stations in order of decreasing Procrustes residuals, using the results of the analysis of fish and BMI concordance (i.e., as plotted in Fig. 4). Geology age classes were generally spread across the gradient, although the highest residuals were found in Older Class and Late Carboniferous stations, whereas the lowest residuals were found in Early Carboniferous stations (Fig. 6), indicating stronger similarity in station characterization based on fish and BMI for Early Carboniferous stations. The strongest correlation of environmental variables with the concordance gradient was a negative correlation of residual values with conductivity (r = −0.67) that indicated that concordance between fish and BMI increased as conductivity increased (Fig. 6). Other environmental variables were generally weakly correlated with the concordance gradient (|r| < 0.30; results not shown), but there was evidence of a positive correlation with D₅₀ (r = 0.43) and bankfull width (r = 0.38) that indicated a decline in concordance (increasing Procrustes residuals) with an increase in substrate size and system size.

Stations underlain by different geological age classes were characterized by environmental variables that have the potential to drive gradients in biotic assemblages (e.g., Snelder et al. 2004). Early Carboniferous stations, which have high shale-gas potential, were strongly characterized by high levels of several water chemistry variables (e.g., ions, conductivity, turbidity, and metals), which may be driven by underlying geologic composition (Johnson et al. 1997; Reimann et al. 2009). In contrast, stations underlain by Late Carboniferous or older classes of bedrock were more strongly associated with a range of habitat variables that were indicative of gradients in stream width, periphyton cover, and riparian vegetation. These results supported the idea that surface and groundwaters in distinctive geological regions might display characteristic abiotic conditions with unique chemical habitats (Stapinsky et al. 2002; Reimann et al. 2009; Neff and Jackson 2011, 2012). Furthermore, differences in conductivity and TP among geologic age classes in our study allowed the development of baseline criteria for monitoring based on the normal range for stations underlain by Early Carboniferous bedrock.

There was clearly a stronger association of BMI with conductivity, which was a defining feature of stations underlain by Early Carboniferous bedrock geology, whereas fish were more strongly associated with system size (catchment area and depth) and a gradient in nutrients that was not specific to a particular geologic age class. Neff and Jackson (2013) noted similar patterns when comparing stream systems on the Precambrian Shield and off-Shield, with BMI responding more strongly to water chemistry differences across geology types and fish predominantly responding to physical habitat variables. Conductivity and alkalinity played a dominant role in distinguishing between Shield and off-Shield stations in their study, and Neff and Jackson (2011) found an association of several noninsect groups (e.g., Gastropoda, Amphipoda, mites, Hirudinea, and Isopoda) with the higher ion levels in off-Shield stations, similar to our results. Furthermore, Johnson et al. (2015) sampled streams in active shale-gas development areas in the southern United States and found their BMI assemblages to be dominated by Oligochaeta, Chironomidae, and tolerant genera of Ephemeroptera, which corresponds to the BMI patterns seen in our study for Early Carboniferous stations. The predominance of noninsect groups in high conductivity waters underlain by Early Carboniferous bedrock may relate to the dependence of taxa such as gastropods and bivalves on uptake of calcium from the water to maintain shell growth (Havas and Rosseland 1995; Neff and Jackson 2013). Though conductivity appeared to be a less dominant driver for fish, it was positively correlated with both Brook Trout and Threespine Stickleback. Conductivity has been found to correlate with abundance (Scarnecchia and Bergersen 1987) and size (as outlined in Copp 2003) of some fish species including trout, with high conductivity stations equated to more productive systems. The positive association of Brook Trout relative abundance with both conductivity and periphyton in these streams appears to support this idea.

Physical habitat differences among stations were reflected by taxonomic associations of both BMI and fish with gradients in temperature and system size. For example, stations that were cooler were associated with Plecoptera taxa such as Leuctridae, Nemouridae, and Chloroperlidae, which prefer small cool-water streams (McCafferty 1998). Brook Trout were similarly associated with small cool-water streams, consistent with their preferred temperature of 16 °C (Coker et al. 2001). The strongest gradient in the fish PCA was a negative relationship along the first axis between Slimy Sculpin, which have a cool-water preference (11–16 °C; Lyons 1990), and stations PL1, PL2, and PL6, which were large, warm systems, with higher maximum summer temperatures than other systems (ranging from 26 to 27.5 °C). These stations had more degree days above 19 °C, and they were positively correlated with Pteronarcyidae, Perlidae, and Perlodidae, which inhabit larger systems and are more tolerant of warm temperatures (McCafferty 1998). Blacknose Dace and American Eel were the fish species that were most strongly associated with these warm-water stations, and these species have preferred temperatures of 24.6 °C and 19 °C, respectively (Coker et al. 2001). Temperature in our study was not specifically related to particular geological age classes, though O’Sullivan et al. (2019) found that surficial geology, slope, and areas of geologic contact could be used to predict stream temperatures. The Pollett River stations, which had the highest degree days above 19 °C and the most variable temperatures in our study, were large stations that also had low percent cover of riparian vegetation and low percent forest in the upstream catchment. Johnson (2004) found evidence that shading affected maximum temperatures in small stream systems, and the lack of riparian vegetation at the Pollett River stations may have contributed to increased temperatures, driving differences in BMI and fish assemblages.

The detection of significant community concordance among BMI and fish has two potential implications: (i) that the current range of environmental drivers in the study systems acts on BMI and fish in a similar manner, resulting in similar groupings of stations regardless of the organism examined and (ii) that there is redundancy in the evaluation of both BMI and fish assemblage structures within the same system, as they both lead to the same conclusion (Jackson and Harvey 1993). In the case of the Kennebecasis and Pollett River watersheds, the evaluation of both BMI and fish assemblages cannot be considered fully redundant as these groups of organisms appeared to respond most strongly to different environmental drivers. Other studies have noted significant concordance among organism groups despite evidence of differing responses to localized drivers (e.g., Jackson and Harvey 1993; Dolph et al. 2011; Bae et al. 2014). In particular, Neff and Jackson (2013) found significant community concordance between fish and BMI assemblages in Shield and off-Shield stream systems despite a stronger association of BMI with water chemistry variables and fish with physical habitat variables. Similar to our study, Neff and Jackson (2013) assessed streams within a localized region where there were distinct geological differences that had a significant effect on the physical and chemical habitat of the stream systems. Where the importance of drivers differs among groups that are concordant (e.g., fish and BMI), each assemblage could be expected to respond differently to future perturbation, which indicates that they provide complementary information about ecosystem condition.

Previous studies that have noted weaker concordance at smaller spatial scales (e.g., within drainage basins or within ecoregions) have attributed the difference in response of organism groups to strong localized differences in small-scale environmental drivers such as water chemistry and physical habitat (Paavola et al. 2003, 2006; Larsen et al. 2012). However, Paavola et al. (2006) suggested that concordance might be evident in small-scale assessments if there is a single, strong environmental gradient driving assemblages in a similar way. In our study, the geological topologies characterized by Early Carboniferous, Late Carboniferous, and older classes of bedrock may have acted in such a manner. The concordance gradient that we identified across stations indicated that the strongest concordance or consistency among fish and BMI was evident in small systems with small substrate size and high conductivity. This association with system size is contrary to the findings of Paavola et al. (2006), who found increased concordance in larger systems in their broadscale regional assessment of streams. However, this may speak to the strength of the geological gradient in the current study and the important distinction of Early Carboniferous (shale-gas producing) areas from those underlain by other geological formations, providing further evidence of the dominant role geology can play in characterizing freshwater habitats (e.g., Neff and Jackson 2011, 2012, 2013).

Though differences in system size, substrate size, and temperature were evident within and among geology age classes, water chemistry was the most important set of variables for characterizing Early Carboniferous stations. In a PCA including only water chemistry data (excluding physical habitat data), the first two axes explained 71.8% of the variance among sites (results not shown) further indicating the importance of water chemistry to distinguish among geology classes. Water chemistry also represents one of the characteristics of surface waters that is most likely to be impacted by potential spills, contamination by untreated flowback fluids (hydraulic fracturing fluids that return to the surface), or the release of treated effluent from shale -gas production (Entrekin et al. 2011; Vengosh et al. 2014). In their review of the potential impacts of hydraulic fracturing activities, Vengosh et al. (2014) noted that surface spills in the United States have caused high levels of organics, metalloids (Se, As), radionuclides (Ra), methane, pH, salts (total dissolved salts, and Cl and Br in particular), and specific conductivity. Wilson et al. (2014) found Cl levels in produced water from shale gas activities to be on average 133 733 mg/L, several orders of magnitude higher than the average Cl levels in the Early Carboniferous stations in the current study (8.38 mg/L), which were elevated compared to Cl levels in Late Carboniferous and Older Class stations (2.34 and 1.63 mg/L, respectively). Lento et al. (2019) found a number of ions had elevated levels in Early Carboniferous stations, consistent with the strong positive association of these stations with conductivity and ions in the present study. Given the already elevated levels of Cl and conductivity in the Early Carboniferous (high shale-gas potential) stations relative to other geological formations in the current study, this suggests the potential for even stronger geological distinction of stations if shale-gas production were to cause impacts to surface waters through wastewater spills. To support monitoring and detection of such impacts, baseline criteria can be developed by identifying the range of values outside which there may be evidence of impairment (normal range; Munkittrick et al. 2009). In our study, confidence intervals for conductivity provided evidence of the distinction among geologic age classes with respect to this chemical parameter, but also indicated the normal range for Early Carboniferous stations within which 95% of the conductivity values fell. With additional monitoring of streams in areas of potential shale gas development, these baseline criteria could be further refined to allow greater accuracy in the detection of impacts.

The stronger association of BMI with conductivity has important implications for biomonitoring in areas with high shale-gas potential, as it suggests that any changes to ion levels that might result from shale-gas production through spills or overland runoff may be most evident as impacts on BMI assemblages. Noninsect taxa and some families of Ephemeroptera and Trichoptera were highly associated with ion levels and could act as indicators of shifts in conductivity, whereas many Plecoptera in particular were more related to temperature gradients. BMI abundance, taxonomic richness, and the percent composition of Chironomidae, EPT, and Oligochaeta were all significantly different in stations underlain by Early Carboniferous bedrock relative to older classes of bedrock geology, similar to the patterns noted by Lento et al. (2019), and the baseline criteria developed for these metrics can be used to monitor how test stations compare with the expected normal range in the absence of impacts from shale-gas development. Fish assemblages provided complementary information by strongly reflecting other physical and chemical habitat characteristics that might be altered with shale-gas development (e.g., increased turbidity from development activities; Entrekin et al. 2011) and natural gradients (e.g., nutrients) that did not relate as strongly to geologic age. Shifts in fish assemblages may be expected to reflect impacts to the system that are related to more general land use and development changes, such as changes to nutrient and sediment regimes or habitat fragmentation caused by clearing the land for construction of wells and development of supporting infrastructure (Entrekin et al. 2011; Brittingham et al. 2014). Thus, monitoring of both BMI and fish may provide a measure of the relative strength of different anthropogenic impacts, including multiple aspects of shale-gas development, as taxonomic preferences and tolerance levels inform assessment of community shifts.

This study characterized the abiotic habitat and biotic communities of streams in geologic areas with high shale-gas potential, contrasting with streams on neighbouring geologic formations. There was clear evidence of differences in the abiotic habitat across geological age classes, and these differences were associated with variability in BMI and fish assemblage structure, with BMI more strongly reflecting temperature and conductivity and fish relating more to system size and gradients in nutrients and turbidity. Despite the differences in dominant drivers across organism groups, there was evidence of significant community concordance among BMI and fish assemblages, and this concordance increased with increasing conductivity and decreasing system size and substrate size. Geological age, differentiating Early Carboniferous, Late Carboniferous, and older classes of bedrock, appeared to act as a large-scale driver in these systems, contributing to concordance within the biotic community despite differences in localized environmental drivers (i.e., physical and chemical habitat variables), and the results of this study suggest that BMI and fish provide complementary information about ecosystem condition across the range of geology types examined. Effective environmental monitoring in areas of shale-gas development should therefore include assessment of both BMI and fish assemblages to allow for the detection of impacts related to shale-gas extraction as well as general perturbation related to development, with deviations from concordance used to detect assemblage-specific driver–response relationships. Importantly, baseline criteria developed in this paper provide estimates of the normal range for chemical and biotic metrics that might be affected by shale-gas development, thereby supporting the development of diagnostic tools for monitoring. This study highlights the importance of examining biotic response to environmental drivers in multiple assemblages to understand and interpret baseline conditions and ultimately measure meaningful change to adequately, and appropriately, inform management and biomonitoring of freshwater systems.