STUDY AREA

The study area was located in northern Labrador and Québec in the eastern Canadian Low Arctic. Sampling took place at 29 sites within either the Koroc River Basin or smaller catchments draining to the Saglek and Nachvak fjords in the vicinity of Torngat Mountains National Park Reserve (58°24′N to 59°6′N and 62°47′W to 65°26′W; Fig. 1). The Koroc River flows west to Ungava Bay, whereas the remaining streams largely flow in an easterly direction to the Labrador Sea. However, barriers to dispersal (in the form of mountains) were found throughout the Koroc, Saglek, and Nachvak regions, and were not limited to the provincial border that separates the Koroc region from the coastal areas. Sites within the Koroc River Basin were on larger rivers (catchment area ranged between 152 and 3407 km², with a median catchment area of 542 km²), while eastern sites were generally located along smaller streams (catchment area ranged between 13 and 845 km², with a median catchment area of 56 km²). Elevations of sites ranged from 15 to 366 m a.s.l. The climate of this region is Arctic, such that streams and rivers are ice-covered from late October until June. A significant spring freshet follows ice break-up, and flows generally decline through the summer except during periods of rainfall (Prowse et al., 2006c). All sites were above the treeline except for 5 sampling locations in the Koroc Basin that were at or just downstream of the treeline. Water source varied between sites, but streams were generally nival or lake-fed.

DATA COLLECTION

Benthic macroinvertebrates were sampled at each stream site in August 2007 by using a time-limited, traveling kick-and-sweep technique, following the standard CABIN (Canadian Aquatic Biomonitoring Network) approach ( http://cabin.cciw.ca). A 400-µm-mesh kicknet was held downstream while the collector disturbed the substrate by kicking and rubbing rocks. The collector moved through the riffle in a zigzag pattern with the D-net, disturbing the substrate for a period of 3 min (Reynoldson et al., 2007). Invertebrate samples were removed from the net and preserved in 95% ethanol for transport.

In the laboratory, macroinvertebrate samples were sorted in their entirety and identified by a Society for Freshwater Science-certified taxonomist. Quality assurance procedures were in place to ensure 95% sorting efficiency. Organisms were identified to the lowest practical taxonomic level (genus when possible). Because of their high abundance, Chironomidae were subsampled until 150 individuals had been obtained from each sample. These individuals were identified to the lowest practical level, and the proportion of Chironomidae taxa in the subsample was used to estimate the total abundance for each chironomid taxon in the total sample.

Environmental data for each sample location were collected at landscape, meso-, and local scales. Due to financial and logistical constraints associated with sampling in such remote regions, field collection of environmental data was necessarily limited. As a result, variables at the larger scales (landscape and meso-scales) were calculated using geospatial data. At the landscape scale, geospatial data were summarized over the entire catchment area upstream of the sample location and classified as catchment-scale variables (Fig. 2, part a). Following Steel et al. (2010), meso-scale data were summarized for a 1-km-radius catchment buffer upstream of each sampling location and classified as reach-scale variables (Fig. 2, part b). At the local scale, physical and chemical data were collected at the sampling location and were classified as site-scale variables (Fig. 2, part c).

FIGURE 1.

Location of study area and the 29 sample sites in the Koroc, Nachvak, and Saglek watershed regions of Québec and Labrador, Canada. The inset panel shows the location of the study area in Canada.

FIGURE 2.

Examples of the three spatial scales over which environmental data were summarized as indicated by the gray-shaded area for (a) catchment, (b) 1-km-radius buffer reach, and (c) site scales.

TABLE 1

Site-scale variables measured for each sample reach, along with their description and measurement units. Variables chosen for inclusion in the RDAs are indicated with *.

To gather site-scale data, water samples were collected at each of the 29 sampling locations and analyzed for nutrients, anions, cations, and metals at Environment Canada's National Laboratory for Environmental Testing in Burlington, Ontario, Canada, following standardized methods (Environment Canada, 2008; Table 1). Chlorophyll a was measured by collecting a composite sample of periphyton from 10 rocks selected at random at each site. A scalpel was used to remove periphyton from within a 12.9 cm² template placed on the upper surface of each rock, and all periphyton scrapes were combined and preserved in 95% ethanol for transport to the laboratory. In the laboratory, the combined sample was homogenized, filtered through a GF/C filter, and the chlorophyll a concentration was determined by extracting the filter and retained material in an 80 °C bath of 90% ethanol for 7 min, then measuring fluorescence with a Turner Designs, model 10 series fluorometer (Sartory, 1982). Chlorophyll a concentrations (g m^-2) were calculated for the combined sample to obtain an average site value.

Substrate composition at each site was characterized by measuring the intermediate axis (b-axis) of approximately 200 particles at the surface of the streambed to the nearest 0.5 cm (0.1 cm for particles less than 1 cm). Particles were chosen at random throughout each study site by walking along transects in the stream reach and measuring substrate particles at the tip of every second boot step. The proportion of particles in the size range of sand (<0.2 cm), gravel (0.2 to <6.4 cm), cobble (6.4 to <25.6 cm), and boulder (2≥25.6 cm) was used to summarize substrate composition. All measured or estimated chemical and physical variables were classified as site-scale variables for our analysis (Table 1).

ArcGIS (Version 9.3, ESRI, St. Paul, Minnesota, U.S.A.), a Geographical Information System (GIS), was used to extract a series of geospatial variables for each sample location. Geospatial variables have increasingly been used to quantify large-scale environmental drivers in studies of hierarchical relationships (Johnson and Host, 2010). A continuous ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer;  http://asterweb.jpl.nasa.gov/gdem.asp) 30-m-resolution Digital Elevation Model (DEM) was used for watershed delineation. The DEM was processed to remove all depressions through a combination of filling and breaching. The upstream catchment for each of the sampling sites was delineated using the ESRI Watershed Delineation toolbox (Noman, 2007). Geospatial variables were extracted for each of the delineated catchments using tools from the Spatial Analyst toolbox and extraction tools from the Geospatial Modelling Environment program (Baier and Neuwirth, 2007; Beyer, 2010; R Development Core Team, 2010). The extracted geospatial variables included bedrock and surficial geology, long-term average annual precipitation, land cover, and landscape measurements (see Table 2 for details). Field-collected measurements were not possible for these larger scales because the remoteness of these locations precluded the opportunity for collecting in situ measurements of additional variables. Each geospatial scale was comprised of a unique set of variables.

DATA PREPARATION

Many early instar insects could not be identified past the level of family or subfamily. To avoid taxonomic redundancy due to mixed-level identification, we evaluated whether the family- or subfamily-level (which included organisms identified to subfamily for Chironomidae, and to family for remaining taxa) could be used to adequately describe the macroinvertebrate assemblages. The effect of classifying individuals at a coarser taxonomic level was evaluated by comparing the spatial arrangement of sample sites in ordination space when different taxonomic levels were used. Results (not shown) indicated that there were no significant differences (at α = 0.05) when either subfamily- or family-level taxonomy was used instead of the genus-level. Due to the strong similarity that was evident in the results based on genus-level and subfamily-level taxonomy, we chose to use the subfamily-level for subsequent analyses. Taxa found in the stream sites are listed in Table 3.

Relative abundance was calculated for each taxonomic level to reduce the influence of extreme values of total abundance across sites and focus on differences in taxonomic composition. Relative abundance values were log₁₀ (x + 1) transformed in order to downweight the influence of both dominant and rare taxa. Invertebrate data were also used to calculate 35 biological metrics in order to compare the association of environmental variables with relative abundance data and summary measures of community structure. Biological metrics included measures of richness, diversity, and composition at multiple taxonomic levels (e.g., order, family, and chironomid subfamily; Table 4). Metric values were either 1og₁₀ (x + 1) or arcsin (√x) transformed depending on which transformation was found to be most appropriate for the data (Table 4). Pearson correlation coefficients were used to remove redundant metrics from the analysis. Strongly correlated |r| > 0.8) metrics were deemed to be redundant, and metrics that shared stronger correlations with other metrics were removed. This resulted in a subset of 24 structural metrics to be used in further analyses to quantify and summarize the biological variability in the benthic community (Table 4).

Environmental variables were log₁₀ (x) or arcsin (√x) transformed as appropriate. Pearson correlations between environmental variables were examined at each scale (catchment, reach, site) in order to select a subset of explanatory variables with minimal redundancy. A pair of variables with |r| > 0.7 was considered to be strongly correlated (0.7 was used instead of 0.8 to ensure that the number of environmental variables did not approach the number of dependent variables), and the variable that shared more strong correlations with other variables was removed to avoid redundancy. The number of selected variables ranged from 7 to 13 depending on the scale that was being considered (Tables 1 and 2).

TABLE 2

Summary of all extracted geospatial variables categorized as climate, landscape measurements, bedrock geology, land cover, and surficial geology descriptors. Variable abbreviations, descriptions, and measurement units are indicated along with the RDA scale for which each variable was used (if any). Data sources are indicated for each descriptor category.

ANALYTICAL APPROACH

Indirect gradient multivariate analysis was used in order to examine the spatial separation of sites based on macroinvertebrate relative abundance or biological metrics. Detrended correspondence analysis (DCA) was run separately for relative abundance data and biological metrics to determine the proper response model for each set of dependent variables. The gradient lengths within our study sites were less than 3.5 standard deviations (SDs) for both relative abundance data and biological metrics (gradient length = 3.02 and 1.29 SDs, respectively), which suggested that a linear response model would be most appropriate (Legendre and Legendre, 1998). We used principal components analysis (PCA) to examine associations among the communities in the stream sites based on relative abundance or metrics.

Direct gradient analysis was used to examine the association between environmental variables at each of the three spatial scales and community composition as described by relative abundance data or biological metrics. Redundancy analysis (RDA) was run separately for each spatial scale and each set of dependent variables. Initial analysis revealed that catchment area was a strong driving factor in the structuring of invertebrate communities in our study area, and that its effect overshadowed the effects of other variables. Catchment area was strongly correlated with longitude (r = - 0.70) and reflected a landscape gradient from the small coastal streams of the Saglek region in the east to the large inland Koroc River in the west. In order to remove the effect of catchment area and focus on the effects of the other environmental variables, catchment area was used as a covariable in all RDAs. The variance inflation factor (VIF) for each variable in the RDA was used to identify redundant variables for each spatial scale; a variable with a VIF > 8 was considered to be redundant and was removed from the analysis. Mean elevation was the only variable found to have a high VIF and was removed from the catchment-scale analysis. At each scale, RDA axis eigenvalues were compared with PCA eigenvalues to determine the proportion of unconstrained community variance explained by each RDA axis. RDA biplots were constructed to visually examine the associations of environmental variables at each scale with the invertebrate community and stream sites.

TABLE 3

Taxa collected at the sample sites, their abbreviations for PCA and RDA biplots, and the number of sites at which each taxon was found.

Variance partitioning (see Borcard et al., 1992, for details) was used to compare the proportion of variance explained by the three scales of environmental variables and determine whether variables at small or large scales contributed more to the structuring of benthic communities at our study sites. One variable was removed from the catchment scale in order to have the same number of variables as in the site and reach scales (7), as suggested in Borcard et al. (1992). In order to choose which variable should be removed from the analysis, we examined eigenvalues and the RDA biplots to determine which variables appeared to explain the least amount of community variance. For the analysis using relative abundance, aspect was removed, while minimum elevation was removed for the analysis using biological metrics. Variance partitioning was accomplished by running RDAs using different combinations of covariables (see Borcard et al., 1992, for details) in order to calculate the proportion of variance explained by each scale alone, as well as the proportion of variance that was explained by the combined effect of the small and large scales. Catchment area was used as a covariable in each of the variance partitioning RDAs.

TABLE 4

Biological metrics considered for the analysis, their abbreviations in PCA and RDA biplots, and data transformations. Metrics chosen for inclusion in the PCAs and RDAs are indicted with *. Note that strongly correlated (|r| > 0.8) metrics were deemed to be redundant, and metrics that shared stronger correlations with other metrics were removed. Abbreviations used include RA (relative abundance) and P/A (presence/ absence).

COMMUNITY STRUCTURE

The first axes of the PCA ordinations of relative abundance data and biological metrics accounted for a large proportion of community variance, although biological metrics appeared to provide a stronger statistical summary of macroinvertebrate community structure. The first three axes of the PCA ordination of relative abundance accounted for 58.4% of macroinvertebrate community variance, with the first axis alone accounting for 26.6% of the community variance. In contrast, the first three axes in combination accounted for 75.4% of the community variance in the PCA ordination of biological metrics, while the first axis accounted for 45.3% of the variance in the macroinvertebrate community.

Stream sites were highly spatially separated in ordination space and covered the entire gradient of axes I and II in the ordination based on relative abundance data (Fig. 3, part a), whereas in the ordination based on biological metrics, very few sites were spatially separated and the majority of sites were clustered around the origin (Fig. 3, part b). Despite using catchment area as a covariable, the Koroc River (KR) sites were generally clustered together. There was a stronger distinction among Nachvak Fjord (NF) and Saglek Fjord (SF) sites, particularly in Figure 3, part a.

Few taxa contributed to the strong gradients in community structure; in the PCA biplot based on relative abundance data, the majority of taxa appeared centrally clustered, with only a small number of taxa, primarily Ephemeroptera and Diptera families, determining the spatial separation of sites in ordination space (Fig. 3, part a). The chironomid subfamily Diamesinae appeared to be an important taxon in the ordination plots of both relative abundance data and biological metrics. Although the abundance of Diamesinae showed a weak, positive association with abundances of Baetidae, it was spatially separated from most taxa in ordination space and was negatively associated with several invertebrate families (Fig. 3, part a). Similarly, the metrics that described the presence/absence or abundance of Diamesinae were negatively associated with metrics describing the richness and relative abundance of several taxonomic groups including Ephemeroptera and Plecoptera (Fig. 3, part b). Moreover, Diamesinae metrics were negatively associated with taxonomic diversity and distinctness, which suggested that this subfamily generally dominated the sites in which it was found.

FIGURE 3.

Principal components analysis biplots of (a) relative abundance and (b) biological metrics data from 29 sites located in the Koroc (KR), Nachvak Fjord (NF), and Saglek Fjord (SF) watershed regions of Québec and Labrador, Canada. Taxon abbreviations in (a) are listed in Table 3. Biological metric abbreviations in (b) are listed in Table 4. Data points for taxa/metrics that were only weakly correlated with axes I and II were not labeled in order to improve readability of the plot.

ASSOCIATION WITH ENVIRONMENTAL VARIABLES

Community structure was highly correlated with environmental variables at all three spatial scales. The first three axes of the relative abundance RDAs together explained 74.8–81.1% of the constrained community variance (Table 5). Higher proportions of variance were explained when biological metric data were used; the first three axes together explained 88.5–91.9% of the constrained variance in the community (Table 5). When the RDA eigenvalues were adjusted by the total PCA eigenvalues, the first three relative abundance RDA axes were found to explain 21.0–29.2% of the unconstrained variance in the macroinvertebrate community, whereas for biological metrics, the percent explained unconstrained variance ranged from 29.3% to 40.6% (Table 5). Relative abundance data were most highly correlated with site- and catchment-scale variables, and biological metrics were most highly correlated with site-scale variables.

At each spatial scale, environmental variables in the RDA ordinations had strong associations with the first and second axes (Figs. 4 and 5). At the catchment scale, the variables positively associated with the first axis for both relative abundance and metrics included average slope, long-term annual precipitation, relative area of glaciers, and bedrock geology (a gradient from sedimentary to volcanic bedrock) (Figs. 4, part a, and 5, part a). The second axis appeared to primarily represent a gradient in catchment aspect for biological metrics (Fig. 5, part a). At the reach scale, land cover descriptors were most important. In particular, bare soil dominated the first axis for biological metrics (Fig. 5, part b) and both axes for relative abundance (Fig. 4, part b). At the site scale, the first RDA axis for both ordinations primarily represented a gradient in substrate size, ranging from sand to boulders, whereas the second axis displayed a gradient in nutrients and primary productivity (Figs. 4, part c, and 5, part c).

TABLE 5

Variance in relative abundance and biological metric data explained by environmental variables at each spatial scale (catchment, reach, and site). Variance is reported in terms of % constrained and % unconstrained variance explained by the first three RDA axes.

The dominant gradient at each spatial scale was not associated with the majority of taxa and metrics. At the catchment scale, this was characterized by a clustering of taxa and metrics at the negative end of the primary gradient (Figs. 4, part a, and 5, part a). At smaller scales, most taxa and metrics were distributed along secondary environmental gradients in vegetation and productivity (Figs. 4, parts b and c, and 5, parts b and c). In contrast, the dominant environmental variables along the first axis generally had a strong association with the chironomid subfamily Diamesinae at each spatial scale (Figs. 4 and 5).

Associations between sites and environmental variables were similar for both relative abundance and biological metric data (results not shown). A small number of sites was consistently associated with the primary gradient along the first axis at all scales; SF05, NF01, and NF03 were positively correlated with variables such as the relative area of glaciers, bare soil, and the proportion of sand. At the catchment scale, KR sites were tightly clustered, whereas NF and SF sites were more separated along the first and second axes. At smaller scales, there was a stronger separation of KR sites along gradients of vegetation (reach scale) and productivity (site scale). Productivity (primarily Ch1 a and TN) appeared to strongly contribute to the separation of all sites at the smallest scale.

VARIANCE PARTITIONING

The catchment, reach, and site scales each explained greater amounts of variance in relative abundance data individually (21.3%, 15.7%, and 24.4%, respectively; Fig. 6, part a) than did their combined effects (8.3%; Fig. 6, part a). The reach scale was most weakly correlated with macroinvertebrate relative abundance; however, the strong individual effects of each scale indicated that all three scales explained unique gradients in community structure. A larger percentage of variance was explained by combined effects of the three scales when biological metric data were used (17.6%; Fig. 6, part b). The site scale stood out as the only scale that explained much more variance individually (27.6%; Fig. 6, part b) and thus the only scale to explain unique aspects of community structure when biological metrics were used.

HIERARCHICAL FRAMEWORK OF ENVIRONMENTAL DRIVERS

In our study systems, biological composition was strongly influenced at multiple spatial scales by habitat variables within a hierarchical framework. The strength of the unique and joint contributions differed for each spatial scale, which indicated the biological importance of variables at multiple scales in these stream systems. One of the strongest gradients evident in the RDAs was related to environmental drivers that were associated with habitat structure and instability; these included the relative area of glaciers, bare soil, and the proportion of sand. This gradient was evident across all scales, demonstrating the consistency of habitat-associated influences and the strength of the hierarchical nesting. At the catchment and reach scales, slope, precipitation, the relative area of glaciers, and the relative area of bare soil all have implications for habitat stability, as they may affect overland flow and runoff rates within the system (Knighton, 1998). In addition, the site-scale gradient in substrate size has implications for both habitat structure and instability as it can be directly related to stream power and the potential movement of particles of different sizes (Knighton, 1998). The range of variables contributing to this habitat structure and instability gradient across the hierarchical framework is indicative of the gradient's potential importance to benthic community structure in Low Arctic systems. Further evidence of this habitat structure and instability gradient could have been obtained through analysis of flow regime and additional habitat structural variables, but these variables could not be measured at our study sites due to their remote, northern location. However, other studies in Arctic systems (Friberg et al., 2001) and temperate systems (Brosse et al., 2003; Weigel et al., 2003; Sandin and Johnson, 2004) have indicated the importance of habitat variables, including those related to instability, across multiple spatial scales. For example, Weigel et al. (2003) evaluated the macroinvertebrate communities in a set of minimally disturbed temperate streams and found a primary gradient in habitat instability from largely coarse substrates to fine particles. The strength of this gradient was evident at both reach and catchment scales, similar to the gradient in habitat instability that appeared to structure the macroinvertebrate communities in our Low Arctic streams.

FIGURE 4.

Redundancy analysis (RDA) biplots of relative abundance data from 29 sites located in the Koroc, Nachvak Fjord, and Saglek Fjord watershed regions of Québec and Labrador, Canada, constrained to environmental variables at three spatial scales: (a) catchment, (b) reach, and (c) site. Environmental variable and taxon abbreviations are as listed in Tables 1 and 3. Data points for taxa that were only weakly correlated with environmental variables were not labeled in order to improve readability of the plot.

FIGURE 5.

RDA biplots of biological metric data from 29 sites located in the Koroc, Nachvak Fjord, and Saglek Fjord watershed regions of Québec and Labrador, Canada, constrained to environmental variables at three spatial scales: (a) catchment, (b) reach, and (c) site. Environmental variable and biological metric abbreviations are listed in Tables 1 and 4. Data points for metrics that were only weakly correlated with environmental variables were not labeled in order to improve readability of the plot.

FIGURE 6.

Results of variance partitioning, showing the percent variance explained by the catchment, reach, and site scales alone, their combined effects, and the unexplained variance for the RDAs using (a) relative abundance and (b) biological metric data. Macroinvertebrate data are from 29 sites located in the Koroc, Nachvak Fjord, and Saglek Fjord watershed regions.

Building on the hierarchical landscape-scale template of stream classification proposed by Frissell et al. (1986), Poff (1997) suggested that hierarchically nested environmental variables may act in succession, with landscape-scale variables acting initially as a filter that determines which organisms can inhabit a particular area. Local-scale variables subsequently act to modify species composition. Although originally described by Poff (1997) to explain the occurrence of species traits, the concept applies to the analysis of community composition. Variables at the landscape scale may initially determine which taxa can inhabit the region (for example, by filtering out taxa that are not adapted to general high-latitude conditions), whereas variables at each successively smaller scale act to further modify the taxonomic composition of the stream (for example, through the availability of suitable habitat patches). While not exclusionary of combined effects of the variables at multiple scales, this concept emphasizes the importance of the individual effect of variables at each scale. Moreover, Townsend et al. (2003) suggested that the combination of the effects of variables at different scales is the most likely determinant of community structure, rather than the effect of one scale alone, and the dominance of each scale as a driver varies among regions (e.g., Brosse et al., 2003; Townsend et al., 2003; Weigel et al., 2003).

In our hierarchical framework, the primary filter acting on community structure was characterized by habitat variables that formed the dominant environmental gradient at each scale. Secondary nested filters were evident at smaller scales and were related to productivity and vegetated land cover. This was particularly evident for sites on the Koroc River. At the catchment scale, there was a strong similarity within KR sites based on the primary filter of large-scale habitat variables. However, at smaller scales, there was a stronger separation among these KR sites along the secondary filter gradients of land cover and productivity. This differentiation of sites would not have been evident if variables at only one scale had been considered. The strength of the unique and joint effects of variables at each scale emphasized the importance of considering a combination of small-, meso-, and large-scale influences on highlatitude stream systems.

INVERTEBRATE RESPONSE TO HIERARCHICAL FILTERS

The chironomid subfamily Diamesinae dominated the ordinations of community structure and provided an example of the strong association of taxa with a gradient in habitat instability. The abundance or presence/absence of Diamesinae was positively correlated with environmental drivers associated with habitat instability, and these associations were evident across all spatial scales. This taxon has been shown to be highly tolerant of extreme habitat conditions that are characteristic of high-latitude and high-elevation streams, namely low water temperatures and high bed instability (Milner et al., 2001). Milner and Petts (1994) developed a conceptual model to describe how macroinvertebrate composition would be expected to change with decreasing distance to a glacier, and therefore decreasing water temperature and habitat stability. In this model, diversity declined with decreasing water temperature and bed stability, and Diamesinae (particularly the genus Diamesa) were suggested to be among the most tolerant organisms. Milner et al. (2001) validated their model with field data, confirming the high tolerance of Diamesinae to low temperatures and to bed instability. Although water temperature data could not be collected at all of our study sites, the available data indicated lower maximum water temperature values and standard deviations during the sampling season in NF and SF area streams (average maximum temperature = 13.55 °C ± 3.19 [n = 4] and 15.37 °C ± 3.64 [n = 3], respectively) than in the KR sites (average maximum temperature = 18.32 °C ± 4.67 [n = 8]). Moreover, the Nachvak and Saglek Fjord regions also contained smaller glaciers not detected because of the GIS data layer resolution, indicating that temperature differences may be due in part to glacial influence. High abundance of Diamesinae was most associated with NF and SF sites, consistent with the available information on water temperature. Although conditions at our study sites were less extreme than those at the glacial sites that were evaluated by Milner et al. (2001), Diamesinae was also found to be positively associated with a habitat instability gradient at all three spatial scales. Diamesinae was particularly dominant at SF05, NF01, and NF03, three sites that were consistently positively correlated with the habitat instability gradient. Furthermore, water temperature data for NF01 indicated a noticeably low maximum temperature with low seasonal variability (average maximum temperature = 10.36 °C ± 2.16).

IMPLICATIONS FOR MONITORING AND MANAGEMENT OF ARCTIC STREAMS

Understanding the spatial scales at which physicochemical variables affect aquatic communities provides important information for management and monitoring of streams. Indeed, the strong relationship between aquatic community structure and environmental drivers at local and landscape scales in our study indicates that future perturbations on such systems may influence biotic communities at multiple spatial scales. Northern aquatic ecosystems may be even more vulnerable to the influence of direct and indirect local landscape disturbances due to the potential for regional and multi-scale effects of climate change. Projected regional increases in temperature will cause potential shifts to the Arctic landscape that may lead to changing land use and landscape disturbances (Prowse et al., 2009). For example, large-scale disturbances associated with resource exploration and development can alter local hydrological patterns, which will directly affect aquatic habitat availability (Prowse and Brown, 2010). Melting of glaciers due to warming may also have significant effects on regional species pools (Jacobsen et al., 2012). In addition, increased nutrients associated with local-scale habitat change, thawing of surface soils, and slumping (e.g., Bowden et al., 2008) caused by large-scale climate warming have the potential to alter productivity and community composition. Locally increased rates of slumping associated with warmer temperatures have already caused greater sedimentation levels in aquatic ecosystems (Lantz and Kokelj, 2008). Moreover, Benstead et al. (2005) demonstrated a dramatic increase in abundance and biomass of macroinvertebrates in response to short-term additions of nutrients within a Low Arctic stream.

Despite the relatively low levels of warming that have occurred in our study area to date (Prowse et al., 2006a), recent models have suggested that the Labrador Sea may be particularly sensitive to warming-related shifts in circulation and salinity, and this may impact inland waters in northern Labrador and Québec (Greene et al., 2008). The results of our study indicate the importance of monitoring environmental variables along the full extent of the hierarchical framework to detect impacts from future warming. Furthermore, climate warming is projected to cause changes to the hierarchical framework that we have described, for example with the northward shifting of the treeline (Culp et al., 2012a), and these changes have the potential to alter the biological importance of drivers at different scales. However, to date few long-term baseline data sets exist for northern regions, particularly for lotic systems, which limits our ability to assess the potential effects of future landscape disturbance and climate change impacts (Prowse et al., 2006c). The development of a strategic long-term monitoring strategy for aquatic ecosystems should build on the limited existing data networks in addition to incorporating large-scale remote sensing information (Culp et al., 2012a).

SUMMARY

This research provides critical baseline data on running water communities in northern Labrador and Québec that may be used for future comparative spatial and temporal analyses. The results of this study suggest that environmental variables at local (site), meso- (reach), and landscape (catchment) scales are important drivers of benthic macroinvertebrate community structure in eastern Canadian Arctic systems. In addition to combined effects of variables at all three scales, the site, reach, and catchment scales each described unique gradients in community structure, emphasizing the importance of the hierarchical framework of environmental variables within this region. There was a primary filter characterized by a clear gradient in habitat structure and instability at large and small scales, indicating the significance of physical structuring variables in the absence of anthropogenic effects. A secondary filter was evident at smaller scales and was related to productivity and land cover. The hierarchical framework that we have described provides a knowledge base for the development of a sustainable long-term monitoring approach for Low Arctic running water ecosystems. By incorporating measures of environmental drivers at multiple spatial scales, future monitoring efforts may more effectively respond to current and future pressures on aquatic biodiversity.