UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Stream invertebrates in northwest British Columbia : an assessment of the relative importance of forest… Bennett, Shauna Ann 2010

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2011_spring_bennett_shauna.pdf [ 1.62MB ]
Metadata
JSON: 24-1.0071507.json
JSON-LD: 24-1.0071507-ld.json
RDF/XML (Pretty): 24-1.0071507-rdf.xml
RDF/JSON: 24-1.0071507-rdf.json
Turtle: 24-1.0071507-turtle.txt
N-Triples: 24-1.0071507-rdf-ntriples.txt
Original Record: 24-1.0071507-source.json
Full Text
24-1.0071507-fulltext.txt
Citation
24-1.0071507.ris

Full Text

Stream invertebrates in northwest British Columbia: an assessment of the relative importance of forest harvesting and environment factors at local and landscape scales  by  SHAUNA ANN BENNETT B.Sc., University of British Columbia, 1993  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF SCIENCE in The Faculty of Graduate Studies (Forestry)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  December 2010  © Shauna Ann Bennett, 2010  Abstract Streams are strongly linked to their terrestrial environments through processes which influence habitat structure and food availability at several, nested scales. While forest practices can affect the hydrological, sedimentation and disturbance regimes, the relationship between community structure and forest practices across broad scales of natural environments is not well understood. In this study, the relative influence of environment and forest harvest attributes on stream community composition was examined in 143 unique 1st to 6th order streams, across a broad geographic scale in northwest B.C. Environment and forest harvest variables explained unique aspects of the community composition using canonical correspondence analysis (CCA), although landscape level environmental variables representing catchment topography and hydrology were selected first using a stepwise procedure. Using partial CCA, 20% of the variation in community structure between sites was explained by environment variables, 12% by forest harvest variables, and 4% was shared between the two sets of variables. The low proportion of shared variance suggested low redundancy between the two sets of variables. Variables describing recent forest harvest (15 years prior to sampling or less) explained unique aspects of the community composition compared to variables describing older forest harvest (more than 15 years prior to sampling), perhaps suggesting partial recovery of the stream community after 15 years. The 143 stream sites were sorted into three LAND USE groups based on the proportion of forest harvest within the catchment basin. Using catchment area as a covariate, relative abundance and rarified taxa richness increased by 30% and 20% respectively, while the ii  proportion of EPT individuals decreased, in catchments with > 15% harvest compared to undisturbed catchments. Distance-based ordination scores (nonmetric multidimensional scaling) changed significantly between the three LAND USE classes of forest harvest on both axes. The effects of forest harvest on streams were confounded with natural environmental gradients. Generally, the catchment basins with >15% harvest were larger, with less topographic relief, more lakes and wetlands, and less rapid drainage characteristics when compared to undisturbed catchments.  iii  Table of Contents Abstract ..................................................................................................................................... ii Table of Contents ..................................................................................................................... iv List of Tables ........................................................................................................................... vi List of Figures ........................................................................................................................ viii Acknowledgements ................................................................................................................... x 1.  Introduction ................................................................................................................ 1  1.1  Community ecology and disturbance.................................................................. 1  1.2  Stream ecology.................................................................................................... 3  1.3  Geomorphology, forest harvest and response of stream communities ............... 4 1.3.1  Intensity and scale of forest harvesting influence ........................................... 6  1.3.2  Biomonitoring and bioassessment ................................................................... 7  1.4 2  Thesis objectives and approaches ....................................................................... 8 Methods .................................................................................................................... 11  2.1  Study area.......................................................................................................... 11  2.2  Site descriptions ................................................................................................ 12  2.3  Field surveys ..................................................................................................... 13  2.4  Data analyses .................................................................................................... 21  3 3.1  2.4.1  Canonical correspondence analysis ............................................................... 22  2.4.2  Ordination of samples by nonmetric multi-dimensional scaling .................. 24  2.4.3  Analysis of covariance (ANCOVA) ............................................................. 25  Results ...................................................................................................................... 29 Principal components analysis of environmental variables .............................. 29 iv  3.2  Canonical correspondence analysis .................................................................. 30  3.3  Effects of land use on biological metrics .......................................................... 42  3.4  Variation of environmental variables with land use ......................................... 44  4  Discussion ................................................................................................................ 48  4.1  Influence of environment on community structure ........................................... 48  4.2  Comparison of two scales of forest harvest influence ...................................... 50  4.3  The effects of forest harvesting on invertebrate communities .......................... 51  4.4  Covariation of forest harvest and environmental gradients .............................. 55  4.5  Implications for bioassessment ......................................................................... 56  4.6  Future directions ............................................................................................... 57  5  Conclusions .............................................................................................................. 59  References ............................................................................................................................... 60 Appendix A Habitat variables and descriptions..................................................................... 73 Appendix B GIS variables and descriptions .......................................................................... 77 Appendix C Macroinvertebrate taxa list. ............................................................................... 82 Appendix D Component loadings on significant factors, Principal Components Analyses .. 84  v  List of Tables Table 1 The number of sites sampled in each of 6 years categorized by LAND USE class where LAND USE 1 = no harvest, LAND USE 2 <15% harvest, LAND USE 3 ≥ 15% harvest in the catchment basin (regardless of year of harvest). .............................................. 14 Table 2 Reach-scale variables selected as potential predictor variables................................ 15 Table 3 Landscape-scale environment variables selected as potential predictor variables. .. 17 Table 4 Percentage of catchment basin harvested by area in Windfall Creek. The total catchment basin area was 8.5 km2. The cutblock area was inverse distance weighted (d+1)-1 using either the estimated flow distance from the block centroid or the Euclidean distance from the block centroid (Equation 1). Unweighted or weighted cutblock areas were divided by the total catchment area and multiplied by 100 to give a percentage value, which was summed for the catchment to give % of catchment harvested................................................ 20 Table 5. Spearman‘s rank correlation coefficients for each of the catchment metrics for forest harvest land use. The metrics include total unweighted harvest (PCT_HARVEST), inverse-distance-weighted flow distance (FLO_PCT_HARV), and inverse-Euclidean distance-weighted (EUC_PCT_HARV) and road density (ROAD_DENS). (n=145). All values are significant ( critical = r(1,140) = 0.139 for α = 0.05) ................................................... 20 Table 6 Forest harvest disturbance variables selected as potential predictor variables. ........ 21 Table 7 Summary of untransformed mean biological metric values (standard deviation in brackets) for 3 different land use groups where group 1 has no harvest, group 2 is <15% harvest and group 3 is ≥15% harvest in the catchment basin exclusive of cutblock age. ...... 26 Table 8 Mean values (standard deviation in brackets) and range of size and forest harvest disturbance descriptors for 143 sites in three LAND USE groups: group 1, no forest disturbance; group 2, <15% forest harvest; and group 3, >15% forest harvest within the catchment area. ....................................................................................................................... 28 Table 9 Summary of principal components analysis on three sets of environmental variables. Full variable names can be found in Table 2 and Table 3. Component loadings on each of the interpreted factors can be found in Appendix D. ............................................................... 30 vi  Table 10 Description of environment and forest harvest variables used for investigating relationships with invertebrate community composition. The first analyses using CCA and BIO-ENV started with all environment variables (LAND, LOC, and GEO variables) and common forest harvest variables (marked with B or x respectively) and either the forest harvest by age class (FH) or inverse-distance-weighted (Euclidean distance) by age class (IDW) variables. Variables marked under Stepwise CCA or Stepwise BIO-ENV columns at right were selected using the stepwise procedures in CANOCO or PRIMER. ...................... 32 Table 11 Results (p values) of Analyses of Covariance (ANCOVAs) to compare treatment means with land use group as a factor and catchment area (log transformed) as a covariate. (α=0.05, DF=2 for land use groups, DF=1 for covariate, DF=139 for experimental error). If the land use effect was significant (p < 0.05, shown in bold font) then the adjusted means and SEM are shown for each LAND USE group. For variables that did not meet the assumption of normality, the results (p value) of a Kruskal-Wallis non-parametric one-way ANOVA are shown. ..................................................................................................................................... 42 Table 12 Summary of results comparing principal components across three LAND USE groups (α=0.05, DF=2 for land use groups, DF=140 for experimental error). For principal components that did not meet the assumptions of normality and equal variance, a nonparametric Kruskal-Wallis test of significance was run and the values shown are the KruskalWallace test statistics (marked with an asterix). Significant p values are shown in bold, where the null hypothesis of no difference in mean values across treatments was rejected. . 46 Table 13 Adjusted means and standard errors (SEM) for Analyses of Variance (ANOVAs) to compare principal components across the three LAND USE groups (α=0.05, DF=2 for land use groups, DF=140 for experimental error) in cases where the land use effect was significant (p < 0.05, shown in Table 12). Post-hoc Tukey tests were used to compare means between land use groups, and where two values were not significantly different they share the same letter. ........................................................................................................................ 47  vii  List of Figures Figure 1 Location of the 143 sampling sites (black circles) in northwest British Columbia. 12 Figure 2 Forest harvest within the Windfall Creek watershed. For each cutblock, the centroid is shown as a black dot. The flow distance from the centroid to the sampling site was measured in two parts; dT, the distance from the centroid to the most downstream stream network point, and dI, the distance from that point to the sampling site. In addition, the Euclidean distance (direct distance as the crow flies) from each cutblock centroid to the sampling site was measured (not shown)................................................................................ 19 Figure 3 Ordination biplot of families (blue triangles) and abiotic variables (red vectors) resulting from CCA. The log transformed family level data was constrained by select environmental variables and all unweighted forest harvest variables. Axes 1 and 2 explain 8.4% and 3.9% of the variation in the abiotic variables, respectively. Taxa codes are listed in Appendix C. Abiotic variable codes are listed in Table 10. ................................................... 33 Figure 4 Percentages of variation of the family data matrix explained by environment and by forest harvest through variance partitioning with pCCA. The column on the left represents total forest harvest in the catchment basin. The column on the right represents inversedistance weighted (Euclidean) forest harvest in the catchment basin..................................... 35 Figure 5 Partial ordination of family data (triangle symbols) constrained by environment variables (quantitative = red arrows) after removing the effect of forest harvest variables (covariables). Axes 1 and 2 explain 5.2% and 2.9% of the family data. Taxa codes are listed in Appendix C. Abiotic variable codes are listed in Table 10. ............................................... 36 Figure 6 Partial ordination of family data (triangle symbols) constrained by forest harvest variables (quantitative = red arrows, nominal (presence/absence or dummy variable = red triangle) after removing the effect of environmental variables (covariables). Axes 1 and 2 explain 2.2% and 1.7% of the family data. Taxa codes are listed in Appendix C. Abiotic variable codes are listed in Table 10. ...................................................................................... 37 Figure 7 Nonmetric multidimensional scaling biplot for stream macroinvertebrate communities based on Bray-Curtis similarities between pairs of sites with family-level logtransformed abundance data. The number of taxa was reduced to exclude rare species. ...... 38 viii  Figure 8 Bubble plot overlay showing the relative abundance of the families identified in BIO-ENV on the NMDS ordination of family-level data. Bubble size indicates taxon abundance relative to maximum abundance (largest bubble size). Sites have been marked by LAND USE group (1, 2 and 3). NMDS1 is along the horizontal axis (forest harvest disturbance decreases along axis) and NMDS2 is along the vertical axis (forest harvest disturbance increases along axis). ........................................................................................... 41 Figure 9 Least square means (and standard error of the means) of rarified richness (ES512) adjusted by catchment area for each LAND USE category where 1 is no harvest (n = 29), 2 is <15% harvest (n = 82) and 2 is ≥15% harvest (n = 32). ......................................................... 43  ix  Acknowledgements Foremost, I would like to thank my supervisor, Dr. John Richardson. I am extremely grateful to him for sharing his remarkable knowledge of ecology and enthusiasm for life-long learning and for his easy-going, approachable nature. I greatly appreciate the considerable efforts of my thesis supervisor, Dr. John Richardson and my committee, Drs. Dan Moore and Allan Carroll, in reviewing and editing this document. Thank you to BC Environment for providing the data for this project. Financial support for the collection of data for this project was provided by the Forest Investment Account Forest Sciences Program to Ian Sharpe (BC Environment, Smithers, BC), and by the BC Ministry of Environment, BC Timber Sales, Westfraser Mills Ltd., Coast Tsimshian Resources Limited Partnership and Canada Resurgence Ltd. Benthic invertebrate enumerations were completed by Danusia Dolecki (Vancouver, BC) and Sue Salter (Summerland, BC). Thank you to Simon Norris of Hillcrest Geographics, Victoria, BC, for completing all GIS analyses and creating the related figures. A special thanks to Megan Hoole, GIS Analyst at BC Timber Sales, Jon Schulz, Forest Planner at Coast Tsimshian Resources Limited Partnership, and Lance Loggin, Silviculture Forester at Westfraser Mills Ltd. for providing access to their spatial datasets for forest operations in Tree Farm Licence (TFL) areas. Thank you to Drs. Robert Bailey and Adam Yates for discussions of anthropogenic stressor gradients as I developed my thesis proposal. A special thanks to Ian Sharpe, for starting me along the path which has become my career, and to Chris Perrin, for his mentorship and guidance along the path, and also to Margaret Branton, for her friendship and support. Thank you to Jason and my family, for your love, patience and encouragement.  x  1. Introduction 1.1 Community ecology and disturbance Stream communities are shaped by habitat availability, food availability and disturbance regimes, and also by the regional species pool and past history (Resh et al. 1988, Giller and Malmqvist 1998). There is a significant portion of the literature focused on describing stream morphology (Frissell et al. 1986) and disturbance regimes (Resh et al. 1988, Hobbs and Huenneke 1992) and how these influence stream communities (Death and Winterbourn 1995, Lake 2000, Allan 2004, Ward 1998, Splinter et al. 2010). While there are many different hypotheses, it is clear that the processes linking stream communities to the surrounding forests are influenced by the physical environment at different scales (Manel et al. 2000, Dauwalter et al. 2008, Richardson 2008, Beechie et al. 2010). Stream habitats, from interstitial spaces, to pools and riffles, to reach, to segment and finally to stream network, are inherently nested (Frissell et al. 1986, Allan 2004) and the physical environment influences processes that shape the biotic community at all levels (Beechie et al. 2010). Habitat availability is strongly linked to channel morphology, which is a function of the hydrological regime, the sediment regime and riparian vegetation (Elosegi et al. 2010). Riparian vegetation enhances bank stability and reduces erosion (Larsen et al. 2009) while also providing in-stream structure and sediment traps through large wood inputs (Hassan et al. 2005, Richardson 2008). Channel morphology is influenced indirectly by climate, catchment size and topography and soil and bedrock properties which directly influence the hydrological regime, the sediment regime and riparian vegetation. For example, Dauwalter et al. (2008) related fish community structure to reach-scale geomorphology, but only after the effects of climate and geography were removed. Stream networks have well described abiotic and biotic longitudinal patterns due to changes in channel morphology and discharge, from the stream headwaters to the mouth (Lowe et al. 2006). For example, water chemistry (e.g. dissolved salts, nutrient levels and pH) generally changes along a longitudinal gradient (from headwaters to mouth), reflecting changes in vegetation, geology, soils, climate and with anthropogenic influences (Giller and Malmqvist 1998).  1  Stream communities are strongly dependent on their terrestrial surroundings for nutrient supply. Organic matter resources in streams come from several sources including groundwater, soils, and riparian vegetation. Riparian vegetation affects food availability by modifying light inputs and temperatures, which can limit primary production, and by supplying leaf litter, a coarse particulate organic matter that serves as a source of biological energy (Kiffney et al. 2003, Richardson et al. 2005). The hydrologic regime plays a role in food availability through groundwater transport of dissolved nutrients. Topography affects hydrological connectivity which influences the concentrations of nutrients exported to surface water (Kreutzweiser et al. 2008a). Steeper watersheds have higher hydrologic connectivities and lower groundwater residence times, which leads to lower concentrations of transported dissolved nutrients. Groundwater residence time affects dissolved nutrient leaching and water chemistry, both of which are also related to surface geology since rates of weathering minerals vary with rock type (Kreutzweiser et al. 2008a). Dissolved organic carbon is also contributed from groundwater throughout the catchment. The dynamic nature of climate and landscape produce natural disturbance regimes in streams. Disturbance can be defined as ‗any relatively discrete event in time that is characterized by a frequency, intensity, and severity outside a predictable range, and that disrupts ecosystem, community or population structure and changes resources or the physical environment’ (Resh et al. 1988). Both natural (e.g. flooding, fires, avalanches) and anthropogenic factors (e.g. acid mine drainage, hydroelectric power generating dams, forest harvesting) can cause disturbances in streams which influence stream communities through changes to the physical environment (Allan 2004). The effects of specific disturbances on streams are generally predictable (e.g. species richness declines with chemical pollution). However, the intensity and longevity of impacts on stream communities can vary widely among regions since they vary with geomorphology and hydrologic regime (Resh et al. 1988). For example, the effects of nutrient enrichment on streams in West Virginia varied by ecoregion with differences in surface geology (Zheng et al. 2008). Anthropogenic disturbance tends to covary with natural landscape gradients (Allan 2004). An example of this is agricultural land use, which tends to concentrate in valley bottoms where the topography, soils and hydrology are more conducive to crop production.  2  Stream communities are also influenced by the pool of species available for recolonization, a concept known as the regional species pool (Giller and Malmqvist 1998). A combination of physical factors and the regional species pool are thought to determine community structure in streams. Poff (1997) suggests that local processes select the species present in any given location through a series of ‗landscape filters‘ at hierarchical scales which are constrained by the regional species pool.  1.2 Stream ecology The benthic (bottom-dwelling) stream community includes arthropods (insects, mites, scuds and crayfish), molluscs (snails, limpets, mussels and clams), annelids (segmented worms), nematodes (round worms), and Platyhelminthes (flatworms) (Hauer and Resh 1996). Aquatic invertebrates are an important trophic link in stream ecosystems (Hauer and Resh 1996). Invertebrates play an important role in transforming autochthonous and allochthonous organic matter into forms usable by animals higher up the food chain (Merritt and Cummins 1996). Aquatic invertebrates known as ―shredders‖ break down coarse particulate organic matter (e.g. leaf litter) into fine organic matter for consumption by collector-gatherers and collector-filterers and colonization by microbes (Heard and Richardson 1995). Other invertebrates known as ―scrapers‖ or grazers feed on the biofilm, a mix of algae, bacteria and fungi in a gelatinous matrix that covers most surfaces in a stream (Giller and Malmqvist 1998). Shredders, scrapers and collectors are prey for higher level animals (e.g. the plecopteran family Perlidae and fish). Stream invertebrates can be organized by taxonomic group or by functional traits; relationships with habitat (shape, mode of existence), feeding (e.g. mode of food acquisition), or by life cycle attributes (e.g. length of life cycle, number of reproductive cycles) (Cummins and Merritt 1996). Compiling stream community data by traits or functional groups allows greater comparability between sites where the species pool might differ (e.g. over large spatial scales) (Statzner and Bêche 2010). In addition, trait analysis can provide insight into the causal links between changes in environment and community, which is especially effective if the mechanistic effects and expected response have been predicted a priori (Statzner and Bêche 2010). For example, chironomids are predominantly associated with  3  fine sediments, and have been shown to increase in abundance with increases in fine sediments and organic matter (Thompson et al. 2009).  1.3 Geomorphology, forest harvest and response of stream communities The effects of forest harvesting on stream communities occur at several different spatial and temporal scales and can be direct (physical, chemical or biological) or indirect (e.g. via food webs) (Richardson 2008). The type and intensity of effects vary with the intensity of harvesting (e.g. clearcut logging versus shelterwood logging) and also with the climate and landform. For example, partial harvest in the riparian buffer has had no detectable effects on leaf litter decomposition (Kreutzweiser et al. 2010, but see Lecerf and Richardson 2010) or fish and macroinvertebrate assemblages (Chizinski et al. 2010) but there is strong support for effects of clearcut harvest in the riparian buffer on the stream biota (Kiffney et al. 2003, Danehy et al. 2007, Lecerf and Richardson 2010). At the catchment scale, forest harvesting can cause erosion and mass wasting, and lead to changes in the hydrologic regime. Forestry roads and cutblocks disturb and expose soils, which can accelerate geomorphic processes, promoting erosion and mass wasting. Sedimentation in streams can degrade habitat and water quality. The intensity of erosion and mass wasting events are positively associated with increases in precipitation, elevation and hillslope gradient, and with south-facing cutblocks (Litschert and MacDonald 2009). Exposed bedrock and soils weather more quickly, which can increase nutrient transport to the stream via groundwater. Forest canopies may dampen the effect of rainfall on stream flows through interception and help regulate water table levels and soil saturation through evapotranspiration (Moore and Wondzell 2005). It follows that canopy interception and evapotranspiration rates are decreased with forest harvest, which results in wetter soils and increased leaching of dissolved nutrients to groundwater which support increases in primary production in a stream (Richardson 2008, Kreutzweiser et al. 2008a). Snowfall accumulation is greater in cutblocks, where melt rates are also increased due to increased solar radiation, potentially increasing winter storm and freshet peak flows and timing (Moore and Wondzell 2005). Roads can affect the stream community by diverting rainfall to the stream more  4  quickly via ditches, which promotes erosion and sedimentation (Moore and Wondzell 2005). It can be difficult to separate the effects of roads and harvest on stream communities, as the land use activities often occur nearly simultaneously. At the reach scale, removal or reductions in riparian vegetation can lead to changes in geomorphology, nutrient dynamics and temperature regime of a stream. Large woody debris recruited from the riparian area provides structure instream that can trap suspended sediments and organic matter, and create diverse habitats (Hoover et al. 2010). Riparian trees can limit primary productivity and maintain cooler stream temperatures through shading and screening the amount of light radiation reaching the stream (Kiffney et al. 2003). Groundwater from adjacent soils transport dissolved organic matter and nutrients to the stream (Richardson 2008). Leaf litter from riparian vegetation (detritus) provides course particulate organic matter (CPOM) which is decomposed to provide energy for food-webs (Richardson 1992). The effects of riparian harvest on stream function likely depend on the size of the stream. The riparian vegetation is more strongly coupled with small streams due to an increased riparian to stream ratio (Richardson et al. 2010). In small, closed-canopy streams, removal of riparian trees can result in higher instream temperatures due to increased exposure to light (Gomi et al. 2006). Larger size and earlier emergence timing in the mayfly Baetis spp. were linked to increased temperatures in an exposed reach compared to a forested reach (Imholt et al. 2010). Over the long term, recruitment of large wood to streams is reduced. Deciduous trees often replace harvested evergreen conifers along the stream edges after harvest, changing the type, quantity and timing of the leaf litter inputs (Richardson et al. 2005, Kominoski et al. 2010) and the riparian zone hydrology (Moore and Wondzell 2005). A change in leaf litter inputs from coniferous to deciduous (red-alder) in Alaska was associated with increases in invertebrate richness, density and biomass compared with old and young riparian coniferous forests (Hernandez et al. 2005) confirming the strong link between streams and riparian vegetation in headwater streams. Much of the research on biological responses to forest harvest activities has focused on small, headwater streams and the largest changes are observed in catchments that have been harvested with small stream buffers, or that have been clearcut to the stream edges. In small  5  streams (1st and 2nd order), macroinvertebrate community responses to forest practices have generally included increased densities and biomass, and shifts in functional groups (Stone and Wallace 1998, Kiffney et al. 2003, Hernandez et al. 2005, Jackson et al. 2007, Nislow and Lowe 2006, Banks et al. 2007, Danehy et al. 2007, Moldenke and Ver Linden 2007, Thompson et al. 2009, Wilkerson et al. 2010). However, in some cases, there was no change in assemblage structure (Herlihy et al. 2005) or a decline in abundance (Kreutzweiser et al. 2008b) or both (Davies et al. 2005). In cases of increased densities and biomass, there were sometimes increases in Chironomidae (Diptera) densities (Kiffney et al. 2003, Nislow and Lowe 2006). The increased biomass and abundance response has been attributed to a shift from heavily shaded stream to an open canopy stream, which leads to increased solar inputs and primary production (Hernandez et al. 2005, Thompson et al. 2009). Reserves of standing timber on either side of a stream in logged areas serve to filter nutrients from groundwater while maintaining essential habitat for aquatic vertebrates and provide shade to the stream to help regulate temperatures (Kiffney et al. 2003, Gomi et al. 2006). However, biological response to forest harvest has been observed in streams even with intact riparian buffers (Kiffney et al. 2003, Martel et al. 2007, Kreutzweiser et al. 2008b, Smith et al. 2009, Lecerf and Richardson 2010). Fewer studies have investigated the biological response to forest harvest activities in larger streams with intact riparian buffers. However, Martel et al. (2007) observed shifts in community structure, increased densities and decreased taxa richness in third order and larger streams, particularly for caddisflies (Trichoptera).  1.3.1 Intensity and scale of forest harvesting influence One challenge with describing forest harvesting within a catchment is that it can vary in proximity to the stream, intensity (e.g. ranging from clear-cuts to single tree selection helilogging), and time since harvest. There are several types of silvicultural systems with varying intensities of harvest used in B.C. including clearcut, seed tree, patch cut, shelterwood and selection. Clearcut and seed tree are conventional systems that remove >90% of the trees from an area. The remaining types typically retain more trees over a harvest area. A group of small openings (less than 1 ha) define a patch cut system. Shelterwood systems harvest trees in stages so that new cohort of trees grows up under the shelter of larger, older trees. Removal of single trees or small groups of trees (opening size  6  less than 2 tree lengths) is known as a selection system. Within these systems, there are ground, cable and aerial techniques for yarding trees and within those categories there are several options for which type of equipment to use ranging from standard heavy equipment to horses. The scale or area of influence (e.g. whole watershed, riparian corridor or partial watershed based on concentric zones or groundwater travel zones, etc.) must be defined before calculating landscape-level variables. While land use within a catchment area influences stream condition, development nearer to the sampling station may have a greater effect on community assemblages (King et al. 2005) supporting the utility of a distance weighting approach to quantifying landscape activity (Van Sickle and Burch Johnson 2008). Martel et al. (2007) investigated the effect of spatial scale for detecting timber harvest effects on stream invertebrates and found that upstream catchment basin was more effective than varying immediate upstream concentric zones.  1.3.2 Biomonitoring and bioassessment Stream invertebrates are a central part of stream ecology because they are a link between organic matter resources and fish (Hauer and Resh 1996). The diversity of species, their individual habitat preferences and tolerances for disturbance, and their sedentary nature make this community ideal for bioassessment. For well over a century, aquatic invertebrates have been recognized as important indicators of disturbance, beginning in the late 1800s as indicators of water pollution (Davis 1994). Benthic invertebrates are commonly the basis for the predictive models as they are ubiquitous, easy and inexpensive to sample, and generally include species sensitive to human disturbance (Karr and Chu 1999). Increasingly sophisticated methods for predicting invertebrate community structure, species presence, and various community metrics (e.g. Simpson‘s diversity index) have been developed to aid traditional chemical and physical impact assessment tools. In British Columbia, the federal and provincial governments have been working on a reference condition approach (RCA) to stream bioassessment (Reynoldson et al. 2001, Mazor et al. 2006, Perrin et al. 2007). The reference condition approach (RCA) has been developed to overcome impact assessment design difficulties typical of stream bioassessment (Bailey et al. 1998, Bailey et al. 2004). For example, impact assessment studies strive to employ either, a before-after (BA), a 7  control-impact (CI) or, a combined before-after-control-impact (BACI) design, which is often not possible in stream systems due to changes along a longitudinal stream gradient. The RCA allows prediction of the expected invertebrate community composition at a site using natural environmental descriptors to pair it with a cluster of biologically similar reference sites (Bailey et al. 2007). While bioassessment has been an effective tool for several types of land use (placer mining, agriculture, urbanization), this has not been shown for forest harvesting and there is contradicting evidence of changes in stream invertebrate community composition associated with forest harvesting (no change, Herlihy et al. 2005; change, Davies et al. 2005, Kreutzweiser et al. 2008b). A better understanding of the effects of forest harvesting activities on the benthic invertebrate community composition, and the environmental variables confounded with those effects, would benefit forest managers by allowing them to identify and select thresholds for specific stressors or combinations of stressors to guide sustainable forest management or indicate system recovery (Richardson 2008).  1.4 Thesis objectives and approaches Empirical studies have helped us learn about the responses of invertebrate communities and indicator groups to certain stressors (e.g. sedimentation). However, the relationship between community structure and forest practices across broad scales of natural environments is not well understood. Observational studies allow us to draw inferences from patterns in broad spatial and temporal scales (Underwood et al. 2000) but they are not without challenges. Often, one challenge is the lack of experimental controls. The predisturbance condition of test sites has not been measured and must be inferred from the condition at undisturbed sites. Another challenge with observational studies over broad geographic scales is the high degree of natural variation, often co-varying with land use, which can mask disturbance effects. Statistical techniques that allow the effects of covariates to be held constant are useful in this type of study design (Allan 2004, Martel et al. 2007). A third challenge is that multiple environmental stressors are often interacting and it can be difficult to assign causation (Allan 2004). This study examined variation in stream invertebrate communities across a large geographic area with environmental and forest disturbance variables. Based on the 8  prediction that both environmental (e.g. size, surface geology, topography of catchment) and forest harvesting variables would influence invertebrate communities, I used partial canonical correspondence analysis (pCCA) to partition the variance explained in the familylevel community structure between the two sets of variables (Borcard et al. 1992). The objective was to determine which factors influence community structure over a broad geographic scale and to describe how the community structure is related to these variables. A second ordination technique, nonmetric multidimensional scaling (NMDS), was used to investigate relationships between community structure and critical environmental and forest harvesting variables using BIO-ENV in PRIMER software (Clarke and Ainsworth 1993, Clarke et al. 2008). While pCCA allows selection of the most importance variables though a step-wise procedure, the NMDS ordinations are the basis of predictive modelling for bioassessment in British Columbia and Canada. A second objective was to compare unweighted and inverse-distance weighted (IDW) measures of forest harvesting influence on community composition. I predicted that the influence of forest harvest in the watershed would decline as a function of distance and that the IDW measure of forest harvest would be more strongly related to the invertebrate community composition than an unweighted forest harvest calculation. Variance partitioning in canonical correspondence analysis was used to compare the percent variation attributable to forest harvest under the two scales of influence. Many of the predicted effects of forest harvest on stream communities are expected to be size dependent. I tested biological endpoints (e.g. taxa richness) across three forest harvest land use classes using catchment area as a covariate. The comparison among land use classes was predicted to show increases in relative abundance of individuals and proportion of chironomids due to increased of food resources, and decreases in taxa richness due to sedimentation. This study focuses on three main objectives. First, to determine which of the environment and forest harvest disturbance variables were most strongly associated with variation in the benthic invertebrate community structure and partition the variation explained with each group of variables. Secondly, to determine whether the proportion of  9  variation in community structure explained by the analyses was different between unweighted and distance-weighted scales of forest harvest influence. Thirdly, to determine whether biological metrics (e.g. taxa richness, relative abundance) changed predictably with forest harvesting using watershed size as a covariate.  10  2 Methods Approximately 300 streams were sampled as part of a BC Environment project to develop a multivariate predictive bioassessment model for the Skeena Region (Skeena BEAST; Perrin et al. 2007). Potential sites were selected through a mapping exercise based on accessibility (most sites were accessed by vehicle with roughly 10% of sites reached by helicopter or boat) and level of disturbance. The goal was to capture a broad range of anthropogenic disturbance types and intensities, and an equal number of undisturbed sites (see Perrin et al. 2007 for details). The original set of 300 sites was reduced to 143 sites that were selected for inclusion in this study, based on the following 5 criteria. 1. The site must be below treeline. 2. Forest harvest was the only known land use disturbance (sites with urban, agriculture or mining land use were removed). 3. Invertebrate samples were sub-sampled using a size-fractionated method (see section 2.3 for details). 4. Sites that were sampled immediately (within 500 m) downstream of lakes, wetlands or beaver impoundments were removed. 5. A complete set of habitat data must be available for the site.  2.1 Study area The study was located in northwest British Columbia (study centroid 127.42° W, 54.86° N) and covered an area of approximately 139,000 km2 (Figure 1). The 143 sites extended from roughly 100 km east of Prince George (121.52°W) to approximately 18 km westsouthwest of the Bell II crossing (second crossing of the Bell River travelling north on the Stewart-Cassiar highway where a heli-ski lodge is located) and 85 km north of Stewart (130.06°W). Latitude of the sites ranged from 53.11°N to 57.09°N.  11  Figure 1 Location of the 143 sampling sites (black circles) in northwest British Columbia.  2.2 Site descriptions Sites were sampled from two terrestrial ecozones, including 79 sites in the Montane Cordillera ecozone and 64 sites in the Pacific Maritime ecozone. The descriptions of the areas were from the Ecological Stratification Working Group (1995) document. In the Pacific Maritime ecozone, sites were sampled in the Coastal Gap, Nass Basin and Nass Ranges ecoregions. The three ecoregions cover the Kitimat Ranges and coast areas of the Coast Mountains. Forests in the three ecoregions are primarily western red cedar and western hemlock in the west, transitioning to mountain hemlock, lodgepole pine, Engelmann spruce and subalpine fir in the east. Maximum annual precipitation ranges from 4500 mm in the Coastal Gap ecoregion, to 2500 mm in the Nass Basin and then to 1500 mm in the Nass Ranges. The Nass Basin and Nass Ranges represent a transitional zone between coastal and interior climates.  12  In the Montane Cordillera ecozone, samples were collected from four main ecoregions including the Skeena Mountains, Bulkley Ranges, Fraser Plateau and Fraser Basin. The Skeena Mountains lie between the rugged coastal ranges to the west and the Omineca Mountains to the east. Interior western red cedar and western hemlock are still common but there is a transition to lodgepole pine, Engelmann spruce and subalpine fir in the east. Similarly, the Bulkley Ranges is a high elevation, mountainous ecoregion that serves as a transition zone to interior areas. In both the Skeena Mountains and Bulkley Ranges, mean annual precipitation ranges from 1500 mm in the west to 600 mm in the east. The Fraser Plateau and Fraser Basin ecoregions represent more interior conditions where the mean annual precipitation ranges from 600 mm in the west to 250 mm in the east. White spruce, lodgepole pine, trembling aspen and Douglas-fir forests are most common in the Fraser Plateau, while trembling aspen, paper birch, lodgepole pine and spruce are more common in the gently rolling Fraser Basin. Overall, the study area covered a broad range of conditions, from the rugged coastal mountains to the interior rolling plateaus and plains of north-central British Columbia.  2.3 Field surveys Benthic invertebrate samples were collected from streams in late summer (range of sampling dates from August 12 to September 16); with 7 streams in 2009 sampled in early fall (September 30 and October 1). Sampling was generally conducted upstream of road crossings. Riffle habitat in each stream was sampled for 3 minutes using a 400 m mesh kick-net (30 cm width). Sampling in larger rivers was restricted to the margins for safety reasons. Samples were preserved and shipped to a Vancouver-based laboratory for identification and enumeration of the invertebrates. Size-fractionated (2 mm sieve, 250 µm sieve) sub-sampling was used to reduce the amount of time needed to sort each sample. Either the whole sample (both fractions) or a minimum of 200 animals from each size fraction were identified to the genus level where possible, except for chironomids which were counted at the family level. The analyses in this study were run with family-level taxonomic resolution of invertebrate data because it matched the taxonomic level of predictive bioassessment models. Reece et al. (2001) found that family and genus level data were equally sensitive to seasonal variation in predictive models for streams in the Fraser  13  River basin. Similarly, a comparison of assemblage-environment relations across various taxonomic levels and transformations conducted by Heino (2008) found that the family-level log- transformed data had the strongest (and significant) correlation with a subset of environmental variables, and that the family and genus level resemblance matrices were highly correlated. The biological data were compiled into a matrix of family abundance by sample. Field surveys were conducted between 2004 and 2009. Sampling year was not included as a potential predictor variable but sites from each LAND USE class were sampled in each year (Table 1). Due to the nature of sampling remote locations, stream sampling in a given year was not geographically independent (e.g. in 2009 all samples were collected from tributaries to the Bell-Irving or Lower Nass rivers near Bell II). Table 1 The number of sites sampled in each of 6 years categorized by LAND USE class where LAND USE 1 = no harvest, LAND USE 2 <15% harvest, LAND USE 3 ≥ 15% harvest in the catchment basin (regardless of year of harvest).  LAND USE1  2004 8  2005 7  2006 4  2007 1  2008 7  2009 2  Total 29  LAND USE2  18  23  10  10  19  2  82  LAND USE3 Total  4 30  13 43  6 20  4 15  2 28  3 7  32 143  Local habitat variables collected during the site visit have been listed in Appendix A (from Perrin et al. 2007). Unfortunately, not all variables were available for all sites. In particular, there were many missing data for flow, water chemistry and substrate variables and therefore, these were not considered in this study. Mean values of three measurements were recorded for wetted and bankfull widths. Stream gradient (%) was estimated with a hand-held clinometer. The local reach was defined as roughly six times the bankfull width, and the proportion of pools, riffles, cascades and glides within the reach area were visually estimated. From these measurements, several variables were included as potential predictor variables representing channel morphology, as shown in Table 2. Variables included average bankfull width, channel ratio, the proportion of pools and the proportion of cascades. Wetted width was not included as it varies with  14  many factors (e.g. recent precipitation, time of day), and is somewhat redundant with bankfull width which was considered to be a more stable model variable. Instream macrophyte and periphyton coverage were considered important indicators of nutrient availability and trophic status of the stream (e.g. oligotrophic, eutrophic). Macrophyte coverage was left as a categorical variable and periphyton coverage, which was recorded as one of 5 categories ranging from 1 (rocks not slippery, no colour) to 5 (rocks mostly obscured by algal mat, algal mass may be in long strands), was converted to a presence-absence variable. Riparian vegetation plays an important role in providing litter, channel stability and shade in streams. Variables representing the dominant structural stage and composition of the riparian vegetation were included as potential predictor variables. Table 2 Reach-scale variables selected as potential predictor variables. Variable code Gradient %  Description Stream gradient (%)  % pools  Area instream as pools (%)  % cascades  BWAve  Area instream as cascades (%) Macrophyte coverage, 5 categories 0%, 1 to 25%, 26 to 50%, 51 to 75% or 76 to 100%. Periphyton coverage noticeably slippery (category 3), very slippery (category 4) or rocks mostly obscured (category 5), yes = 1, no = 0 Canopy coverage, 5 categories 0%, 1 to 25%, 26 to 50%, 51 to 75% or 76 to 100% Structural stage of dominant vegetation, tree cover greater than 10% yes = 1, no = 0 Presence / absence of dominant riparian vegetation as conifers (present = 1, absent = 0) Mean bankfull width (m)  Stream Order  Stream order  CHANRATIO  Channel ratio (wetted width / channel width)  MC PERI_35 CC SSCOV10 RipVeg_Con  Landscape level variables were measured using geographic information systems (GIS) analyses. The catchment basin, the area of the watershed upstream of the sampling point, was delineated using the Fresh Water Atlas (formerly Corporate Watershed Base) layer (first order watersheds) and the TRIM 25 m Digital elevation model (25 m DEM), a  15  digital model of the ground surface topography. It was a three step process that involved extracting all first order polygons upstream of the sample point, using the 25 m DEM to refine the boundary to terminate at the sample point, and then combining the first order watersheds into a single polygon (pers. comm., Simon Norris, Hillcrest Geographics, Victoria, BC, March 20, 2009) which was then given a unique identifier. The watershed polygons were then used to run a GIS overlay. Data layers were compiled for each of the variables listed in Appendix B. Variables were compiled using ArcGIS 9 geographic information systems (GIS) software developed by ESRI© (http://esri.com/). Spatial datasets were accessed through the Province of BC spatial data directories known as the Land and Resource Data Warehouse (LRDW) and the Integrated Land Management Bureau (ILMB) ARC warehouse, and the federal spatial database known as the Canadian Soil Information System (CANSIS). All GIS analyses were performed by Simon Norris (Hillcrest Geographics, Victoria, BC). Landscape-scale variables included as potential predictor variables in this study are summarized in Table 3. Several variables were selected to represent regional topography, geology and land cover which have a strong influence on stream geomorphology (Montgomery 1999). Surface geology was represented by rock types, local surface form and dominant drainage characteristics. Sedimentary and volcanic rocks were expected to be important variables as these types of rock tend to weather more quickly than others (e.g. metamorphic) and contribute phosphorus to the stream. The average aspect (degrees) of the catchment basin was used to calculate the heat load index where warmer slopes have values close to 1, and cooler slopes have values close to 0 (McCune et al. 2002).  16  Table 3 Landscape-scale environment variables selected as potential predictor variables. Variable code ECOZONE  Description Dummy variable, 2 ecozones (0=pacific maritime, 1=montane cordillera)  LATITUDE  Decimal degrees  LONGITUDE  Decimal degrees  % SEDIM  Area of catchment at sedimentary rock (%)  % INTRU  Area of catchment as intrusive rock (%)  % VOLCA  Area of catchment as volcanic rock (%)  % METAM  Area of catchment as metamorphic rock (%)  % ULTRA  Area of catchment as ultramafic rock (%)  Dom_DRAIN  % WETLANDS  Dominant drainage class (1=imperfect, 2=moderate, 3=rapid, 4=well drained) Dominant local surface form class (1=dissected, 2=hummocky, 3=inclined, 4=level, 5=rolling, 6=ridged, 7=steep, 8=terraced, 9=undulating). Area of catchment as wetlands (%)  %LAKES  Area of catchment as lakes (%)  % ALPINE  Area of catchment as alpine (%)  ELEVATION  Elevation at sampling site (masl)  RELIEF  Maximum catchment elevation - elevation at sampling site (m)  DRN_DNSTY  Drainage density (km/km2)  % hillslope > 60%  Proportion of watershed with hillslopes greater than 60%  BARREN_PA  Barren land greater than 1% of catchment basin? Yes = 1, no = 0  HEAT LOAD INDEX  Aspect converted to scale of 0 to 1, Heat Load Index  AREA  Area of catchment basin (km2)  ICE_PA  Area of catchment as ice (%) converted to presence absence variable  dom_LOCSF  For each unique cutblock in a catchment basin, the type of logging (e.g. clearcut), age and area of cutblock, flow-distance and Euclidean-distance from the block centroid to the sample site were compiled from three different sources: the Vegetation Resources Inventory (VRI), Tree Farm Licence (TFL) and Reporting Silviculture Updates and Land Status Tracking System (RESULTS). The TFL spatial dataset was provided by forest licence holders to the GIS analyst, while the VRI and RESULTS databases were available through the LRDW. In this study, the proximity of harvest activities to the stream was highly variable but the intensity of logging did not vary at the cutblock level and was not further considered as a potential predictor variable. Clear-cut logging accounted for 99.3% of the total harvest in all catchments (591 km2), with the remaining 0.7% being attributable to selection and shelterwood logging (4 km2).  17  The flow distance was modified from Van Sickle and Burch Johnson (2008) and consisted of two parts, the distance from the cutblock centroid to the nearest downstream waterbody connected to the sampled stream (DT), and the distance from the point of entry to the nearest downstream waterbody to the sample site (DI) (Figure 2). The total percent harvest in each catchment basin was calculated three ways: without distance weighting (%FH), and using inverse-distance-weighting (Equation 1) of each cutblock area based on flow distance (%FH_FLO), and Euclidean distance (%FH_EUC). An example of the distance weighted calculations for one watershed has been prepared and is shown in Table 4.  18  Figure 2 Forest harvest within the Windfall Creek watershed. For each cutblock, the centroid is shown as a black dot. The flow distance from the centroid to the sampling site was measured in two parts; dT, the distance from the centroid to the most downstream stream network point, and dI, the distance from that point to the sampling site. In addition, the Euclidean distance (direct distance as the crow flies) from each cutblock centroid to the sampling site was measured (not shown).  19  Equation 1  Inverse Distance Weighted Area of cutblock = A∙(d + 1)-1 Where A= cutblock area (km2), and d = distance (km)  Table 4 Percentage of catchment basin harvested by area in Windfall Creek. The total catchment basin area was 8.5 km2. The cutblock area was inverse distance weighted (d+1)-1 using either the estimated flow distance from the block centroid or the Euclidean distance from the block centroid (Equation 1). Unweighted or weighted cutblock areas were divided by the total catchment area and multiplied by 100 to give a percentage value, which was summed for the catchment to give % of catchment harvested.  Cutblock ID 48358 28661 28667 28660 28648 28649 28621  Cutblock Area (km2) 0.136 0.103 0.147 0.330 0.931 0.514 0.799  Flow Distance (dT + dI) (km) 5.40 5.12 5.30 4.41 3.55 2.26 2.72  Euclidean Distance (km) 4.58 2.88 2.84 3.70 1.89 0.82 2.39  Total % of Watershed Harvested  % of Catchment Harvested 1.6 1.2 0.2 3.9 10.9 6.0 9.4  % of Catchment Harvested (area weighted by Inverse Flow Distance) 0.2 0.2 0.0 0.7 2.4 1.9 2.5  % of Catchment Harvested (area weighted by Inverse Euclidean Distance) 0.3 0.3 0.0 0.8 3.8 3.3 2.8  33.2  8.0  11.3  Number of Years from Harvest to Sample Collection 13 15 15 16 17 17 22  Peterson et al. (2010) found that the two models for inverse-distance weighting were highly correlated. Similarly, in this study the two distance weighted measures were highly correlated (Spearman rank correlation, r = 0.996, Table 5). Since the flow based distance could not be calculated exactly as planned (DT was a direct distance from the centroid rather than a flow path based distance), further comparisons proceeded with only the Euclidean inverse-distance weighted and unweighted data. Table 5. Spearman’s rank correlation coefficients for each of the catchment metrics for forest harvest land use. The metrics include total unweighted harvest (PCT_HARVEST), inverse-distance-weighted flow distance (FLO_PCT_HARV), and inverse-Euclidean distance-weighted (EUC_PCT_HARV) and road density (ROAD_DENS). (n=145). All values are significant (rcritical = r(1,140) = 0.139 for α = 0.05)  PCT_HARVEST FLO_PCT_HARV EUC_PCT_HARV ROAD_DENS  PCT_HARVEST  FLO_PCT_HARV  EUC_PCT_HARV  0.941 0.957 0.873  0.996 0.797  0.814  20  Using GIS, the length of unbuffered stream in each cutblock was measured. Only waterbodies connected to the stream network were included in the analyses (i.e. lakes, wetlands or streams that did not have a surface water connection to the stream network were excluded). The length was recorded as the same value whether logging was on one side or two of the unbuffered waterbody. The total length of roads in each catchment basin was compiled from the VRI layer. All variables considered as potential predictor variables have been summarized in Table 6. Table 6 Forest harvest disturbance variables selected as potential predictor variables. Variable code RH_10PA  Age 1-20y  Description Presence of riparian harvest in 10 years prior to sampling Proportion of riparian harvest in 5 years prior to sampling (length of riparian harvest / stream length) Proportion of age class 1 forest (1-20years) in watershed  Age 21-40y  Proportion of age class 2 forest (21-40years) in watershed  Age >250y  Road_Dens  Proportion of age class 9 forest (>250years) in watershed Mean cutblock age for all blocks in the catchment basin. In catchments with no harvest, a value of 99 years was used to represent an older forest. Kilometers of roads per square kilometers watershed area (km/km2)  FH % harvest  Proportion of catchment as forest harvest (%)  IDW % Harvest  Proportion of euclidean-flow-weighted harvest in catchment (%)  %FH <15y  Proportion of total % harvest by age, less than 15 years  %FH 15-29y  Proportion of total % harvest by age, 15 - 29 years  %FH >30y  Proportion of total % harvest by age, greater than 30 years  % IDW <15y  Proportion of inverse-direct-distance weighted % harvest by age, less than 15 years  % IDW 15-29y  Proportion of inverse-direct-distance weighted % harvest by age, 15 - 29 years  %IDW >30y  Proportion of inverse-direct-distance weighted % harvest by age, greater than 30 years  RH_5p  BLKage  2.4 Data analyses Before relating environmental and disturbance variables to the biota, the number of abiotic variables had to be reduced. Principal components analyses were used to reduce groups of correlated variables into fewer principal component variables. Principal components analysis (using SYSTAT 11) was run separately on each of three groups of environmental variables: landscape scale variables, local stream habitat variables and geological variables. Separating the variables into groups prior to PCA, rather than running a  21  single PCA, was done to permit variance partitioning of environmental data groups. All principal components with an eigenvalue greater than 1 were interpreted as being significant.  2.4.1 Canonical correspondence analysis Canonical correspondence analyses (CCA) were performed to relate taxonomic community structure at the family level (software: CANOCO version 4.55, Biometris – Plant Research International, Wageningen, The Netherlands). Prior to analyses, taxa that occurred at fewer than 7 sites were removed from the dataset. This reduced the overall number of taxa from 73 to 44, while removing only 0.02% of individuals from the total abundance (Appendix C). The data were log10(x + 1) transformed to increase the weighting of taxa with low or moderate abundance by downweighting the highly abundant families. The additional option to ‗downweight rare species‘ in CANOCO was not selected after preliminary results showed that less total variation was explained with that option selected. A preliminary detrended correspondence analysis showed that the gradient lengths in the data were relatively short (1.299 to 1.513 standard deviation units of species turnover, or SD). The gradient length ―is a measure of how unimodal the species response are along an ordination axis‖ (ter Braak and Smilauer 2002). This is a measure of beta diversity, the amount of change in family composition along a directly measured gradient of environment (or time). Where gradient lengths are less than 4 SD, the risk of an arch effect in the CCA ordination is low. Direct gradient analysis with biplot scaling (site scores are rescaled such that the mean is zero and the variance is one) was selected to promote an ordination of taxa with an optimal environmental basis. Environmental and forest harvest variables selected through the stepwise procedure were plotted along with species on biplots. The distances between species and distances between vectors and species on the biplots approximates their Chi-square distances which allows direct spatial interpretation of any relationships between species and environment points from the biplots (McCune et al. 2002). The probability of type I error (α), was set at 0.1 for stepwise selection of the significant environment variables using Monte Carlo permutation tests.  22  Partial Canonical Correspondence Analysis (pCCA) was used to partition the total explained variance in community structure into 3 components, environmental, forest disturbance, and shared variation (Borcard et al. 1992). Partial, constrained ordinations are run by removing the effects of selected variable(s) as covariates through multiple linear regression (Borcard et al. 1992). The fractions of variance were calculated with 4 separate CCA runs where the sum of all canonical eigenvalues using the given sub-set of abiotic variables is divided by the sum of all eigenvalues (total inertia). The total inertia is the variation in the family level data that can be explained by the abiotic variables. Partitioning the variance was a four step procedure. First, the CCA was run constraining the family-level abundance data by the abiotic PC variables (landscape, local and geological). Second, the CCA was run constraining the family-level abundance data by the forest harvest variables. Third, a partial CCA was run constraining the ordination by the abiotic PC variables selected using forward-selection after removing the effects of the forest harvest variables. Fourth, a partial CCA was run constraining the ordination by the forest harvest variables selected using forward-selection after removing the effects of the abiotic PC variables. The total variance explained by the abiotic variables (environment and forest harvest) was calculated as the variance explained from step 1 plus the variance explained from step 4 (or step 2 plus step 3). Steps 3 and 4 provide the variation explained by the environment and forest harvest variables respectively. Shared variation between the environment and forest harvest variables was the difference between steps 1 and 3, or steps 2 and 4 (should be similar). The significance of the pCCA outputs were tested using a Monte Carlo randomization test of the first axis and all axes together (499 unrestricted permutations under the reduced model, α = 0.05). All analyses were performed twice, once with unweighted forest harvest variables and a second time with inverse-direct-distance weighted forest harvest variables. The purpose was to compare the two measures of forest harvest to see if one or the other explained a greater proportion of variation in the family data. Principal components were used to represent environmental variables (landscape, local and geology) in order to retain the maximum amount of information in the fewest possible variables. Forest disturbance variables were not reduced using PCA since they were few in number and interpretation of the results was simplified by using individual variables. 23  In each step of the calculations, unrestricted Monte Carlo permutation tests (499 permutations) on the overall trace statistics were performed to test if one or several predictor variables could explain a significant part of the variance in community structure.  2.4.2 Ordination of samples by nonmetric multi-dimensional scaling Non-metric multidimensional scaling (NMDS) was used to further examine the invertebrate community structure at the sites. NMDS was used in addition to CCA, to investigate relationships between community structure and critical environmental and forest harvesting variables using BIO-ENV in PRIMER software (Clarke and Ainsworth 1993, Clarke et al. 2008). While CCA allowed selection of the most important environmental variables through a step-wise procedure with tests of significance and variance partitioning, the NMDS ordinations added value as they underlie the predictive modelling for bioassessment in British Columbia. The samples were ordinated using non-metric Multidimensional Scaling (NMDS) based on the Bray-Curtis similarity matrix of the log10 (x + 1) transformed data. The NMDS ordination is a graphic representation of the rank order and relative distance between the sites. The corresponding stress value indicates how well the graphic represents the similarity matrix and, ultimately, whether or not an interpretation based on the ordination will be reliable. For 2-dimensional ordinations, stress values less than 0.3 with a large number of sites (>50) suggest that interpretation of the ordination is reasonable, although stress values less than 0.2 are preferred. Scores along each NMDS axis (NMDS1 and NMDS2) were saved for each site. The relationship of environmental and forest disturbance variables with the NMDS ordination was explored with a routine called BIO-ENV in the statistical software PRIMER (version 6, Primer-E, Plymouth Marine Laboratory, Plymouth, UK). BIO-ENV uses a stepwise procedure (BVSTEP) to select a sub-set of the abiotic variables that optimizes a match with the biotic similarity matrix underlying the MDS ordination. This procedure can also be run with the families as variables to select a sub-set of taxa that are more influential in driving the ordination pattern (Clarke and Warwick 2001).  24  The BIO-ENV procedure differs from CCA in several ways. CCA partitions the variance in the abiotic variables to see which variables have the strongest relationship with the biota and to express the proportion of the total variance in the original data represented by the ordination. NMDS is a more flexible ordination approach that has fewer assumptions than other ordination algorithms about the data and the response to the abiotic variables (McCune et al. 2002). The NMDS ordination is a graphical representation of all the variation in the biotic data while the CCA represents only the variation in community structure that is explained by the selected abiotic variables. NMDS ‗preserves the rank order of amongsample dissimilarities in the rank order of distances‘ (Clarke 1993).  2.4.3 Analysis of covariance (ANCOVA) Invertebrate community descriptors included family-level taxa richness (rarified to the lowest relative abundance of 512 individuals, ES512), relative abundance, Simpson‘s diversity index (Equation 2), Simpson‘s evenness (Equation 3), Ephemeroptera- PlecopteraTrichoptera (EPT) richness, proportion of EPT individuals, proportion of Chironomidae, proportion of trichopterans, and trichopteran taxa richness (Table 7). Rarefaction was used to standardize taxa richness to a common effort across all samples, since richness is known to increase with the area sampled and the number of individuals in the sample, called the collector‘s curve phenomenon (Colwell and Coddington 1995). These community metrics are commonly used in biomonitoring programs and often have predictable responses to specific types of disturbance (Lenat 1988, Karr and Chu 1999). For example, EPT taxa richness decreases with increasing urban and agricultural land use disturbance (Lenat 1988), and both richness and abundance of individuals in these insect orders have been shown to vary with forest harvesting (Hernandez et al. 2005, Martel et al. 2007). Simpson‘s Diversity index is a measure of the abundance and taxonomic richness of a community while Simpson‘s Evenness index is a measure of the distribution and relative contribution of each taxon present to the overall community (McCune et al. 2002). Both Simpson‘s Diversity index and evenness range from 0 to 1, with values close to 1 signifying a more diverse and even community. All community metrics met the assumptions of normality and equal variance.  25  Equation 2 s  Simpson’s Diversity Index =  1 – Σ (pi)2 i–1  where pi is the proportion of the ith taxon at the station, and S is the total number of taxa at the station  Equation 3 s  Simpson’s where Evenness Index =  1 / Σ (pi)2 / S i–1  Table 7 Summary of untransformed mean biological metric values (standard deviation in brackets) for 3 different land use groups where group 1 has no harvest, group 2 is <15% harvest and group 3 is ≥15% harvest in the catchment basin exclusive of cutblock age. Biological Metric  LAND USE 1 n = 29  LAND USE 2 n = 81  LAND USE 3 n = 33  Relative Abundance  3547 (1968)  3560 (2673)  5214 (4064)  Rarified Family Richness (ES512)  16.6 (3.3)  17.0 (3.8)  19.0 (3.1)  Simpson's Diversity  0.77 (0.11)  0.77 (0.13)  0.80 (0.09)  Simpson‘s Evenness  0.28 (0.09)  0.27 (0.09)  0.27 (0.09)  EPT Richness (ES512)  12.1 (1.9)  11.6 (2.1)  12.2 (1.95)  Proportion of EPT Individuals  0.78 (0.22)  0.78 (0.17)  0.68 (0.16)  Proportion of Chironomidae  0.16 (0.19)  0.15 (0.16)  0.20 (0.18)  Trichoptera Richness  3.3 (1.4)  3.5 (1.3)  3.9 (1.3)  Proportion of Trichopterans  0.09 (0.09)  0.07 (0.07)  0.08 (0.07)  NMDS1  0.300 (0.487)  0.139 (0.743)  -0.613 (0.761)  NMDS2  -0.298 (0.508)  0.018 (0.593)  0.217 (0.740)  Considering both the skewness and kurtosis of the sampling distribution of each variable along with p values from a Shapiro-Wilk normality test, the variables family taxa richness, EPT richness, NMDS1 and NMDS2 met the assumption of normality. One variable, relative abundance, was log10 transformed to reduce skew. Two variables, proportion of trichopterans and proportion of chironomids, were square root-transformed to meet the 26  assumption of normality. An arcsine (√x) transformation was also tried for proportional data, but resulted in less improvement. Since the variable, proportion of EPT individuals, did not meet the assumption of normality and normality was not improved with several transformations, a non-parametric test was used. Analyses of covariance (ANCOVAs) were conducted to determine the effect of forest harvest on streams after accounting for potential differences due to stream size, with post-hoc Tukey tests to compare means between pairs. Sites were arbitrarily categorized into three land use (LAND USE) classes: 1 as shown in Table 8. While age of forest harvest was not included in the land use categorization, there was no significant difference in average block age between LAND USE 2 and LAND USE 3 (2 sample t test, DF = 112, p = 0.702). Taxa richness and abundance have been shown to vary predictably with stream width (Paller et al. 2006). Catchment area was used as a covariate to account for possible differences due to stream size in separate analyses of covariance (ANCOVAs) to investigate the relationship of invertebrate community descriptors with land use category in SYSTAT 11. Across all sites, the linear relationship between catchment area and bankfull width was highly significant (R2 = 0.51, F1, 141 = 146.7, p < 0.001). Bankfull width is potentially influenced by land use and confounded with stream gradient, so catchment area was selected as the covariate to account for potential changes in invertebrate community descriptors with size. The GLM procedure was used to estimate the model with the interaction term between treatment (LAND USE) and the covariate (L10_area) in the model, to ensure that there was no significant interaction between the treatment and covariate. There was no significant interaction (e.g. p = 0.565 for the interaction term with rarified richness as the dependent variable) indicating that the assumption of homogeneity of slopes is reasonable.  27  Table 8 Mean values (standard deviation in brackets) and range of size and forest harvest disturbance descriptors for 143 sites in three LAND USE groups: group 1, no forest disturbance; group 2, <15% forest harvest; and group 3, >15% forest harvest within the catchment area. LAND USE 1 Catchment Descriptor Catchment Area (km2)  LAND USE 2  LAND USE 3  n=29 30.7 (42.7) 0.9 - 167.0 8.6 (5.8) 1.5 - 26.7  n=82 n=32 52.8 (51.8) 38.7 (38.3) 1.1 - 198.1 2.1 - 147.1 Bankfull Width (m) 13.5 (14.3) 6.5 (4.0) 1.5 - 78.3 1.5 - 20.2 Forest Harvest (%) 5.9 (4.3) 31.8 (14.1) 0 0.3 - 14.9 15.8 - 76.3 Average Cutblock Age1 99.4 (21.0) 15.3 (9.5) 16.2 (7.3) 16 - 99 1 - 40 0 - 33 Road Density (km/km2) 0.02 (0.04) 0.34 (0.31) 1.25 (0.59) 0 - 0.2 0 - 1.43 0.2 - 3.1 1 In catchment basins without any harvest, an average cutblock age of 99 years was assigned to approximate a mature forest in these watersheds.  28  3 Results The relative abundance of invertebrate individuals ranged from 512 to 17,532, while the number of unique family-level taxa ranged from 11 to 29. The number of unique trichopteran families ranged from 0 to 8, with a mean value of 1.3 while the unique number of EPT families ranged from 8 to 18, with a mean value of 13 taxa. On average, individuals in the EPT orders accounted for 76% of all individuals, although the proportion ranged from 11.2 to 99.8%. The proportion of chironomids ranged from 0 to 88%, although on average they represented 17% of individuals. The average Simpson‘s diversity value was relatively high (0.78) and ranged from 0.22 to 0.91. Simpson‘s evenness ranged from 0.08 to 0.52, with a mean value of 0.27 at all sites.  3.1 Principal components analysis of environmental variables The 14 landscape variables were reduced to 4 interpretable factors which explained 67% of the total variation in the dataset (Table 9). The 11 local scale variables were reduced to 5 interpretable factors which explained 66% of the total variation in the original variables. The 7 geological variables were reduced to 3 interpretable factors which explained 65% of the total variation in the original variables.  29  Table 9 Summary of principal components analysis on three sets of environmental variables. Full variable names can be found in Table 2 and Table 3. Component loadings on each of the interpreted factors can be found in Appendix D.  Variable Set  Principal Component Name  % Total Variance Explained  Variables with high positive (+) or negative (-) loadings  Interpretation  LANDSCAPE  LANDpc1  33.7  % hillslope >60% (+), relief (+), longitude (-),% alpine(+)  catchment topography  LANDSCAPE  LANDpc2  13.4  Area (+), % lakes (+),% barren (+), % wetlands (+)  catchment basin size and hydrology  LANDSCAPE  LANDpc3  10.9  ecozone (+), elevation (+)  geographical location  LANDSCAPE  LANDpc4  9.0  Latitude  LOCAL  LOCpc1  22.1  LOCAL  LOCpc2  14.0  Latitude (+) Bwave (-), stream order (-), canopy cover (+) %cascades (+), gradient % (+)  LOCAL  LOCpc3  10.7  %Pools (-), Macrophyte Cover (+)  Habitat diversity  LOCAL  LOCpc4  9.8  cover >10% (+), rip veg CON (+)  mature, coniferous cover  LOCAL  LOCpc5  9.3  chanratio (-), peri35 (+)  Nutrients, trophic status  GEOLOGY  GEOpc1  28.1  % sedimentary (-),% volcanic (+)  Gradient of nutrient weathering rock  GEOLOGY  GEOpc2  20.0  % intrusive (+), DOMdrain (+), locSF (+)  catchment drainage characteristics  GEOLOGY  GEOpc3  16.6  % metam (-), % ultra (-)  non-weathering rock  stream size and cover Stream power  3.2 Canonical correspondence analysis Roughly 29% of the variation in the family data matrix was explained by all the environment (landscape, local and geology) and unweighted forest harvest variables (Table 10) using CCA. The first two axes of the 4 canonical eigenvalues were interpreted. Significant variables (Monte Carlo permutation n = 499, α= 0.1) with the strongest relationship to the community (represented by the arrow length) included catchment basin size and hydrology, steepness of the catchment basin, stream size (LANDpc2, LANDpc1, LOCpc1 in Table 9) and % forest harvest less than 15 years prior to sampling (Figure 3). Road density and forest age class 1 – 20 years (Age 1 – 20 y) had long arrows but were not statistically significant during the forward step-wise procedure. The riparian and forest harvest variables grouped together on the biplot and were positively correlated with species 30  axis 1. The principal component variable representing local stream size (LOCpc1) was positively correlated with species axis 2 (weighted correlation coefficient, 0.45). Vectors representing the percentage of forest harvest 15 to 29 years prior to sampling, and more than 30 years prior to sampling, were correlated. The principal component variable representing channel ratio (ratio of wetted to bankfull width) and periphyton abundance (LOCpc5) grouped with the forest harvest variables along with the pc variable representing the proportion of sedimentary and volcanic rocks (GEOpc1). The landscape variable representing the steepness of the catchment (LANDpc1) was correlated with the pc variable representing local stream gradient (LOCpc2). The step-wise procedure reduced the number of environment variables to 9 and the number of forest harvest variables to 7 (Table 10) in two separate CCA analyses. Significant variables common to both CCA analyses were marked with a ‗B‘ in Table 10. Selected variables specific to either the unweighted forest harvest or the inverse-distance weighted forest harvest CCA runs were marked as ‗FH‘ or ‗IDW‘ respectively.  31  Table 10 Description of environment and forest harvest variables used for investigating relationships with invertebrate community composition. The first analyses using CCA and BIO-ENV started with all environment variables (LAND, LOC, and GEO variables) and common forest harvest variables (marked with B or x respectively) and either the forest harvest by age class (FH) or inverse-distance-weighted (Euclidean distance) by age class (IDW) variables. Variables marked under Stepwise CCA or Stepwise BIO-ENV columns at right were selected using the stepwise procedures in CANOCO or PRIMER. Abiotic Description Variables Label Environment Variables LANDpc1 catchment basin topography LANDpc2 catchment basin size and land cover (% lakes, %wetlands, %barren) LANDpc3 geographical location (ecozone and elevation) LANDpc4 Latitude LOCpc1 LOCpc2 LOCpc3 LOCpc4 LOCpc5  stream size and cover steepness / high gradient (% cascades, gradient) % pools and macrophyte coverage mature, coniferous cover UV exposure / stream cover  GEOpc1 GEOpc2 GEOpc3  Gradient of nutrient weathering rock catchment drainage characteristics Metamorphic, ultra  Forest Harvest Variables RH_5p Proportion of riparian harvest 5 years prior to sampling (length of stream harvest / total stream length) RH_10PA Presence/absence of riparian harvest 10 years prior to sampling rd_dens Age 1-20y Age 21-40y Age >240y  Road density in the catchment basin (km/km2) Proportion of watershed with forest aged 1-20 years Proportion of watershed with forest aged 21-40 years Proportion of watershed with forest aged greater than 240 years  BLKage  Average cutblock age in the catchment  Proportion of forest harvest by age class % FH < 15y Proportion of catchment harvested 14 years prior to sampling % FH 15-29y Proportion of catchment harvested 15 to 29 years prior to sampling % FH >30y Proportion of catchment harvested more than 30 years prior to sampling Inverse-distance-weighted Proportion of forest harvest by age class % IDW <15y Proportion of catchment harvested 14 years prior to sampling % IDW 15-29y Proportion of catchment harvested 15 to 29 years prior to sampling % IDW >30y Proportion of catchment harvested more than 30 years prior to sampling % Variation explained by all abiotic variable in CCA (FH) Correlation with NMDS (FH)  Stepwise CCA  Stepwise BIO-ENV  B B  x x  B B  x  B B B  x  B  B B  X  B x B B FH FH FH  IDW IDW IDW 29% 0.374  32  Figure 3 Ordination biplot of families (blue triangles) and abiotic variables (red vectors) resulting from CCA. The log transformed family level data was constrained by select environmental variables and all unweighted forest harvest variables. Axes 1 and 2 explain 8.4% and 3.9% of the variation in the abiotic variables, respectively. Taxa codes are listed in Appendix C. Abiotic variable codes are listed in Table 10.  While Euclidean distance between the cutblock centroid and the sample site ranged from 70 m to 24 km, using an inverse-distance-weighted approach to calculating the percent forest harvest in the catchment made no difference in the total percentage of variation in the family data matrix explained by forest harvest variables (Figure 4). After removing the effect of environmental variables, the effect of forest harvesting explained 12% of the variation in the invertebrate community. After removing the effect of forest harvesting, the effect of the environmental variables explained 20% of the variation in the invertebrate community. 33  Within each of those variable groups, 3 to 4% of the explained variation was shared between the two sets of abiotic variables. A large portion of the community variation could not be explained by the abiotic variables (73%). In all pCCAs, the variation in the family data matrix explained by the abiotic variables was significant (p = 0.002) as evaluated by an unrestricted Monte Carlo permutation test on the reduced model (n = 499). After removing the effect of the significant forest harvest variables, the landscape, local and geology variables were shown in Figure 5. The LANDpc1 variable was negatively correlated with species axis 1 (weighted correlation coefficient, -0.57) while the LANDpc2 and LANDpc3 variables were weakly correlated with species axis 2 and species axis 3 respectively (weighted correlation coefficients, 0.47 and 0.40). The remaining significant variables were weakly correlated with the species axes. Overlap of the LANDpc1 and LOCpc2 arrows indicated strong covariation of the two variables. After removing the effect of the significant landscape, local and geological variables, the forest harvest variables had a greater range (Figure 6). The variable representing the presence of riparian harvest 10 years prior to sampling (RH_10PA) was a nominal variable (presence-absence) and was represented as a point, rather than an arrow like the continuous variables. The proportion of forest 21 to 40 years old (Age 21-40y) was weakly correlated with species axis 1 (correlation coefficient, 0.36). All other variables were weakly correlated to the species axes. Generally, variables representing recent harvest were positively correlated with species axis 2 and were opposite to variables representing harvest more than 15 years old which were negatively correlated with species axis 2. Family centroids located close to an arrow are associated with that variable. Psychodidae and Ceratopogonidae appear to be associated with riparian harvest within 10 years and forest harvest greater than 30 years old, respectively. The dipteran family Psychodidae appears to be associated with the maximum block age, which represents mature forests.  34  100%  90%  80%  70%  Percent of Variation  72  73  Unexplained  60% Forest Harvest 50%  Environment + Forest Harvest  40%  Environment  30%  20%  10%  8  8  4  3  16  16  Unweighted Forest Harvest  Inverse-distance-weighted harvest  0%  Figure 4 Percentages of variation of the family data matrix explained by environment and by forest harvest through variance partitioning with pCCA. The column on the left represents total forest harvest in the catchment basin. The column on the right represents inverse-distance weighted (Euclidean) forest harvest in the catchment basin.  35  Figure 5 Partial ordination of family data (triangle symbols) constrained by environment variables (quantitative = red arrows) after removing the effect of forest harvest variables (covariables). Axes 1 and 2 explain 5.2% and 2.9% of the family data. Taxa codes are listed in Appendix C. Abiotic variable codes are listed in Table 10.  36  Figure 6 Partial ordination of family data (triangle symbols) constrained by forest harvest variables (quantitative = red arrows, nominal (presence/absence or dummy variable = red triangle) after removing the effect of environmental variables (covariables). Axes 1 and 2 explain 2.2% and 1.7% of the family data. Taxa codes are listed in Appendix C. Abiotic variable codes are listed in Table 10.  NMDS ordination showed a wide spread of sites in space (Figure 7). The stress value for the ordination was relatively high (0.26), although this is partially a function of the large sample size (n = 143). Using the BVSTEP procedure in the BIO-ENV routine, the community structure was most strongly correlated with the proportion of riparian harvest five years prior to sample collection (RH_5p) (BVSTEP correlation coefficient = 0.245). A subset of 6 variables including proportion of riparian harvest 5 years prior to sample collection, LANDpc1, LANDpc3, Age 1-20y, LANDpc2 and LOCpc2 had a combined BVSTEP 37  correlation coefficient of 0.363 (significance level = 0.2%). There was considerable overlap in the significant variables selected by the stepwise procedure in CCA and the combination of variables with the best correlation to the NMDS ordination using stepwise BIO-ENV procedure (Table 10).  Figure 7 Nonmetric multidimensional scaling biplot for stream macroinvertebrate communities based on Bray-Curtis similarities between pairs of sites with family-level log-transformed abundance data. The number of taxa was reduced to exclude rare species.  A separate run of the BVSTEP routine relating the family abundance matrix to the BrayCurtis similarity matrix indicated that 26 of the 44 taxa were highly influential in the ordination (BVSTEP combined correlation coefficient = 0.953, p < 0.002). Taeniopterigydae was the most correlated with the ordination (BVSTEP stepwise correlation coefficient = 0.409, p < 0.01). A correlation coefficient of nearly 0.8 could be achieved with a combination of 11 taxa including Elmidae, Chironomidae, Ephemerellidae, Leptophlebiidae, Capniidae, Chloroperlidae, Perlidae, Taeniopterygidae, Hydropsychidae, Rhyacophilidae and Torrenticolidae (not in order of importance). Only three of the influential taxa did not belong to the EPT orders. Bubble-plot overlays are shown in Figure 8 for 5 of the first 8 taxa 38  selected using the BVSTEP procedure. The relative abundance of the family Taeniopterygidae tended to be higher in undisturbed and moderately disturbed sites (LAND USE 1 and 2) along the NMDS1 axis as shown in Figure 8. In contrast, the stonefly predator Perlidae tended to be more abundant in sites with forest harvest disturbance (LAND USE 2 and 3) along both axis 1 and 2, although not all sites with forest harvest housed this family. The relative abundance of the mayfly family Leptophlebiidae increased with forest harvest along the axis NMDS1. Elmidae beetles appeared to increase with increasing forest harvest along both NMDS1 and NMDS2 axes. The caddisfly family Rhyacophilidae, which occurred in 140 of 143 sites, appeared to increase along NMDS2, although not necessarily with LAND USE groups. The stonefly families Capniidae and Chloroperlidae and the caddisfly family Hydropsychidae are not shown, but patterns for these taxa were most similar to that of the family Perlidae.  39  40  Figure 8 Bubble plot overlay showing the relative abundance of the families identified in BIO-ENV on the NMDS ordination of family-level data. Bubble size indicates taxon abundance relative to maximum abundance (largest bubble size). Sites have been marked by LAND USE group (1, 2 and 3). NMDS1 is along the horizontal axis (forest harvest disturbance decreases along axis) and NMDS2 is along the vertical axis (forest harvest disturbance increases along axis).  41  3.3 Effects of land use on biological metrics Invertebrate abundance (log10 transformed to improve normality) differed among land use groups (ANCOVA, F2,139 = 3.991, R2 = 0.065, p = 0.021) after adjustment for stream size (Table 11). Least squares means (adjusted for covariate) and standard error (SE) of the mean were calculated by back-transforming the log10-transformed outputs. Using backtransformations creates asymmetrical errors around the adjusted mean, so the range of SE is shown. After adjusting for catchment area, relative abundance was not significantly different at undisturbed sites (LAND USE 1, mean 3177 individuals per sample, SE 2793 to 3614, n = 29) and sites with < 15% harvesting (LAND USE 2, 2742, SEM 2541 to 2958, n = 82; Tukey group 1 versus 3, p = 0.310, groups 1, 2 p = 0.603). Relative abundance was significantly greater at sites with > 15% harvesting (LAND USE 3, mean 4112 individuals per sample, SE 3639 to 4645, n = 32; Tukey group 2 versus 3, p = 0.013). Table 11 Results (p values) of Analyses of Covariance (ANCOVAs) to compare treatment means with land use group as a factor and catchment area (log transformed) as a covariate. (α=0.05, DF=2 for land use groups, DF=1 for covariate, DF=139 for experimental error). If the land use effect was significant (p < 0.05, shown in bold font) then the adjusted means and SEM are shown for each LAND USE group. For variables that did not meet the assumption of normality, the results (p value) of a Kruskal-Wallis nonparametric one-way ANOVA are shown.  Biological Metric Relative Abundance Rarified Family Taxa Richness (ES512) Simpson's Diversity  Land Use Effect P value  Covariate Effect P value  LAND USE 1 n = 29  LAND USE 2 n = 82  LAND USE 3 n = 32  0.021  0.121  3177 a (2793 – 3614)  2742 a (2541 – 2958)  4112 b (3639 – 4645)  0.005  0.009  16.2 ± 0.7 a  17.1 ± 0.4 a  19.1 ± 0.6 b  0.369  0.050  0.322 ± 0.13 a  0.134 ± 0.08 a  -0.634 ± 0.12 b  Simpson‘s Evenness  0.873  0.035  EPT Richness (ES512) Proportion of EPT Individuals1 Proportion of chironomids2 Trichoptera Richness1 Proportion of Trichopterans NMDS1  0.578  <0.001  0.001 0.338  0.570  0.320 0.954  <0.001  < 0.001  0.194  NMDS2 <0.001 -0.210 ± 0.11 a -0.007 ± 0.06 ab 0.209 ± 0.10 b 0.018 1 2 Did not meet the assumption of normality, Did not meet assumption of equality of slopes.  42  Family taxa richness of invertebrates (rarified to the lowest relative abundance of 512 individuals) differed among land use groups after adjustment for stream size (ANCOVA, F2,139 = 5.516, R2 = 0.111, p = 0.005) as shown in Figure 9. The covariate effect was also significant, suggesting that the use of ANCOVA was justified. Undisturbed streams had the lowest size-adjusted mean taxa richness (LAND USE1, 16.2 ± 0.7 SE, n = 29) with similar values in streams with <15% harvesting (LAND USE 2, 17.1 ± 0.4 SE, n = 82). Sites with >15% harvest in the catchment basin had significantly higher size-adjusted mean taxa richness (LAND USE 3, 19.1 ± 0.6 SE, n = 32; Tukey group 1 versus 2, p = 0.454; group 1 versus 3, p = 0.004; groups 2, 3 p = 0.023).  Least Squares Means 22  ES512  20 18 16 14 12  1  2 LANDUSE  3  Figure 9 Least square means (and standard error of the means) of rarified richness (ES512) adjusted by catchment area for each LAND USE category where 1 is no harvest (n = 29), 2 is <15% harvest (n = 82) and 2 is ≥15% harvest (n = 32).  Simpson‘s diversity and evenness were not significantly different between the three LAND USE groups (ANCOVAs, F2, 138 = 1.005, R2 = 0.041, p = 0.369 and F2,138 = 0.136, R2 = 0.032, p = 0.873). A non-parametric one-way ANOVA was significant for the proportion of EPT individuals with LAND USE (Kruskal-Wallis test statistic 13.481, p = 0.004). While the 43  overall proportion of EPT individuals was significantly lower at the LAND USE 3 sites, there was no change in the unique number of EPT taxa (EPT richness) between any of the three LAND USE groups (ANCOVA, F2, 138 = 0.551, R2 = 0.120, p = 0.578). The proportion of chironomids (square-root transformed to improve normality) did not differ among land use groups (ANCOVA, F2,139 = 1.092, R2 = 0.018, p = 0.338). The covariate effect was not significant, suggesting that the ANCOVA was not needed as there was no effect of stream size on the results. The proportion of trichopterans (square-root transformed to improve normality) was not significantly different between the three LAND USE groups (ANCOVA, F2, 138 = 0.047, R2 = 0.138, p = 0.954). Trichoptera richness did not meet the assumption of normality. The normality of the variable was not improved with several common transformations, so a nonparametric test was used to test the effect of LAND USE on the variable. A non-parametric one-way ANOVA was not significant for Trichoptera richness with LAND USE (KruskalWallis test statistic 2.278, p = 0.320). Variables consisting of NMDS ordination scores from a 2-dimensional plot (NMDS1 and NMDS2, Figure 7) were both significantly different between the LAND USE groups (NMDS1 ANCOVA, F2, 138 = 17.967, R2 = 0.212, p < 0.001; NMDS2 ANCOVA, F2, 138 = 4.125, R2 = 0.228, p = 0.018). The covariate was only significant along the second axis (NMDS2, p < 0.001). Along NMDS1, the adjusted mean scores for LAND USE 1 (0.322 ± 0.132 SEM) and LAND USE 2 (0.134 ± 0.077) were not significantly different (Tukey, p = 0.442), but both were significantly greater than LAND USE 3 (-0.634 ± 0.123) (Tukey, p < 0.001). For NMDS2, LAND USE 1 (-0.210 ± 0.107) was significantly less than LAND USE 3 (0.209 ± 0.100) (Tukey, p = 0.012), but there was no significant difference in LAND USE 2 (-0.007 ± 0.063) from either group (Tukey group 1 versus 2, p = 0.240; group 2 versus 3, p = 0.163).  3.4 Variation of environmental variables with land use Seven of the twelve variables were significantly different between treatment groups (Table 12 and Table 13). LANDpc1 differed among land use groups (ANOVA, F2,140 = 44  13.681, R2 = 0.163, p < 0.001). Undisturbed catchments were the steepest with similar values in streams with <15% harvesting (Tukey group 1 versus 2, p = 0.996). Catchments with >15% harvest in the catchment basin were significantly less steep (Table 13, Tukey group 1 versus 3, p < 0.001, group 2 versus 3, p < 0.001). LANDpc2 differed among land use groups (ANOVA, F2,140 = 3.456, R2 = 0.047, p = 0.047). Post-hoc comparison of least squares means showed that undisturbed sites (LAND USE 1) were significantly smaller catchment size with less lakes and wetlands cover, from sites with <15% forest harvest (Tukey group 1 versus 2, p = 0.038). Catchments with >15% forest harvest had the highest LANDpc2 values, but were not significantly greater than values from undisturbed catchments (Tukey group 1 versus 3, p = 0.059) or catchments with <15% forest harvest (Tukey group 2 versus 3, p = 0.963). There were significant differences in LANDpc3 between all three LAND USE groups (ANOVA, F2,140 = 16.797, R2 < 0.001, p = 0.194) with a decreasing trend from higher elevations in the Pacific Maritime ecozone to lower elevations in the Montane Cordillera ecozone. Post-hoc comparison of least squares means showed that undisturbed sites (LAND USE 1) were significantly greater than sites with forest harvest (Tukey group 1 versus 2, p < 0.001; group 1 versus 3, p < 0.001) and sites with < 15% forest harvest were significantly greater than sites with > 15% forest harvest (Tukey group 2 versus 3, p = 0.007). The LANDpc4 variable represented latitudinal change. LANDpc4 differed among land use groups (ANOVA, F2,140 = 3.778, R2 = 0.051, p = 0.025) but the change was not linear across the LANDUSE groups. Post-hoc comparison of least squares means showed that LANDUSE 1 sites had the highest latitudes, and LANDUSE 2 (Tukey group 1 versus 2, p = 0.025) had the lowest latitudes, with latitudes for LANDUSE 3 sites in between (Tukey group 1 versus 3, p = 0.261; groups 2, 3 p = 0.626). At the local scale, two variables were significantly different across the LAND USE groups; LOCpc2 represented an increasing gradient of channel slope (% cascades and stream gradient), and LOCpc3 represented an increasing gradient of the proportion of pools and macrophyte cover. LOCpc2 differed among land use groups (ANOVA, F2,140 = 14.13, R2 = 0.199, p < 0.001). Post-hoc comparison of least squares means showed that LAND USE 1 and 2 were higher gradient streams (Tukey group 1 versus 2, p = 0.199), and both differed from LAND USE 3 which had lower gradient streams (Tukey group 1 versus 3, p < 0.001, 45  2,3 p < 0.001). A non-parametric one-way ANOVA was significant for the LOCpc3 with LAND USE (Kruskal-Wallis test statistic 8.179, p = 0.017). GEOpc 2 represented catchment drainage characteristics from well drained to poorly drained. The GEOpc2 variable differed among land use groups (ANOVA, F2,140 = 3.836, R2 = 0.052, p = 0.024) but the change was not linear across the LANDUSE groups. Post-hoc comparison of least squares means showed that LANDUSE 1 was not significantly different from LANDUSE 2 or LANDUSE 3 (Tukey group 1 versus 2, p = 0.968; 1,3 p = 0.111) but they were significantly different from each other (Tukey group 1 versus 3, p = 0.018). Table 12 Summary of results comparing principal components across three LAND USE groups (α=0.05, DF=2 for land use groups, DF=140 for experimental error). For principal components that did not meet the assumptions of normality and equal variance, a non-parametric Kruskal-Wallis test of significance was run and the values shown are the Kruskal-Wallace test statistics (marked with an asterix). Significant p values are shown in bold, where the null hypothesis of no difference in mean values across treatments was rejected. Principal Component  Interpretation  F value  p value  R2  LANDpc1  Topography/steepness of catchment  13.681  <0.001  0.163  LANDpc2  Catchment size and natural land cover  LANDpc3  Ecozone, elevation  LANDpc4  3.456  0.034  0.047  16.797  < 0.001  0.194  Latitude  3.788  0.051  LOCpc1  Stream size (width, order) and cover  2.428  0.025 0.092  LOCpc2  Channel slope (% cascades, gradient)  14.13  < 0.001  0.199  LOCpc3  Habitat diversity (pools, macrophyte cover)  *8.179  LOCpc4  Mature, coniferous cover  0.075  0.017 0.928  0.001  LOCpc5  UV exposure / stream cover  0.694  0.501  0.01  GEOpc1  Sedimentary to volcanic gradient  *1.373  0.503  GEOpc2  Catchment drainage characteristics  3.836  GEOpc3  Metamorphic / ultramafic geology  *0.575  0.024 0.75  0.034  0.052  46  Table 13 Adjusted means and standard errors (SEM) for Analyses of Variance (ANOVAs) to compare principal components across the three LAND USE groups (α=0.05, DF=2 for land use groups, DF=140 for experimental error) in cases where the land use effect was significant (p < 0.05, shown in Table 12). Posthoc Tukey tests were used to compare means between land use groups, and where two values were not significantly different they share the same letter. Principal Component Variable  LAND USE 1  LAND USE 2  LAND USE 3  0.229 ± 0.171 a  0.212 ± 0.102 a  -0.750 ± 0.163 b  LANDpc2  Description Topography/steepness of catchment Catchment size and natural land cover  -0.426 ± 0.183 a  0.093 ± 0.109 b  0.147 ± 0.174 b  LANDpc3  Ecozone, elevation  0.742 ± 0.168 a  -0.029 ± 0.100 b  -0.598 ± 0.160 c  LANDpc4  Latitude Channel slope (% cascades, gradient) Catchment drainage characteristics  0.368 ± 0.182 a  -0.183 ± 0.108 b  0.136 ± 0.173 ab  0.475 ± 0.167 a  0.141 ± 0.100 a  -0.722 ± 0.159 b  0.096 ± 0.181 ab  0.147 ± 0.108 a  -0.406 ± 0.172 b  LANDpc1  LOCpc2 GEOpc2  47  4 Discussion 4.1 Influence of environment on community structure The general goal of this type of large-scale, observational study is to define and describe the abiotic variables that affect the processes that shape community assemblages (Beisner et al. 2006). Understanding the relative influence of environmental factors on community structure at various spatial scales (e.g. reach, catchment, landscape) is a critical first step in evaluating the effects of disturbance on the biota (Resh et al. 1988, Beechie et al. 2010). Previous large-scale aquatic studies have found that abiotic variables (geomorphic, geographic, stream chemistry and land use) have accounted for 19 to 50% of the variation explained in macroinvertebrate assemblages and indicators in streams (Johnson et al. 2004, Sandin and Johnson 2004, Kratzer et al. 2006, Johnson et al. 2007) and wetlands (Brazner et al. 2007). In terrestrial systems, the results have been generally similar. In floral assemblages, abiotic variables have accounted for 35 to 45 % of the variation explained in upland grass assemblages (Vandvik and Birks 2002), vegetation assemblages in sphagnum bogs (Lachance and Lavoie 2004), and in ground flora in plantation forests (French et al. 2008). Similarly, in a study of terrestrial spiders, abiotic variables have accounted for 49% of the variation explained in communities (Bowden and Buddle 2010). In this study, environmental and forest harvest variables explained 29% of the variation in the family level assemblage data, which is consistent with previous studies. The portion of unexplained variance was high, but this is not uncommon for assemblage data (Borcard et al. 1992). The large amount of variation unexplained by the environmental and forest harvest variables suggests that there might have been relevant factors not taken into account or that there is a large degree of stochastic variation in the dataset (McCune et al. 2002). Catchment-scale variables had a strong influence on invertebrate communities, followed by reach-scale variables in two different analyses (pCCA and BIO-ENV). Stream invertebrate communities varied with the topography or ‗steepness‘ and size of the catchment, although the proportion of area as lakes, wetlands and barren lands were also important. This was in agreement with a study of streams in New Zealand that found catchment scale variables including relief ratio, basin shape and catchment area were the best 48  predictors of invertebrate assemblages (Townsend et al. 2003). In addition, Brosse et al. (2003) found that relief ratio was the most important catchment-scale variable for predicting richness and evenness. In this study, longitude covaried with topography, likely a reflection of the range in geographic location of sites from the Interior Plateau and across the Coast Mountains. Stream invertebrates varied with terrestrial ecozone, which was not unexpected since these classifications represent climate, geology and vegetation. At the reach scale, stream gradient and the proportion of habitat as cascades were highly influential. Mature coniferous riparian cover and in-stream periphyton cover were significant but less influential variables. Stream size (e.g. bankfull width, stream order), proportion of pools and macrophyte coverage were not important environmental drivers of community structure in this study. Overall, reach-scale variables were less important than landscapescale variables in this study. There is evidence that reach-scale variables play a significant role in driving community composition (e.g. Brosse et al. 2003, Johnson et al. 2007, Beechie et al. 2010). However, local- and reach-scale variables were not well represented in this study. In particular, bedform particle size and water chemistry variables (e.g. pH, dissolved organic carbon, dissolved nutrients) were not included in the analysis, and may have increased the proportion of variation explained in community structure. Another possibility is that reach-scale variables are coupled tightly with habitat structure and could be more important drivers of species traits (e.g. body shape, functional feeding group) rather than community structure. Richards et al. (1997) found that reach-scale variables (% shallow areas, stream size, canopy cover, large wood, and % fines) were important predictors of species traits. Surface geology, including area of specific rock types (e.g. sedimentary), local surface form and drainage characteristics contributed very little to the variation explained in the biota. Sedimentary and volcanic rocks were expected to be important variables as these types of rock tend to weather more quickly than others (e.g. metamorphic) and contribute phosphorus to the stream. In other studies, geology has been a strong environmental driver of species variation, explaining up to 17% of variation in the species assemblage (Kratzer et al. 2006). It is not clear why surface geology variables were not more important in this study. Perhaps, the streams are nitrogen limited so that the addition of phosphorus does not 49  result in changes to the invertebrate community. It is also possible that the surface geology variables were redundant with other environment variables and therefore, were not selected in the canonical correspondence analysis.  4.2 Comparison of two scales of forest harvest influence In the field of landscape ecology, there have been a number of studies focused on defining the scale of influence for particular land uses to improve the understanding of the relationship between land use and the biota (Johnson and Host 2010). These studies address the question of defining the portion of the watershed (e.g. catchment, riparian zone) that best reflects the observed changes in the biota (Burcher 2009). The premise underlying distanceweighted land use is that there should be a gradual reduction of influence of land use on a stream site with increasing distance upslope and upstream away from the site. Initial work focused on comparing the catchment-, riparian buffer- and reach-scale of land use (see Allan 2004). More recent work has incorporated distance-weighted calculations of land use using flow-pathways to try and derive the most realistic representation of the true scale of influence (Van Sickle and Burch Johnson 2008, Peterson et al. 2010). Strong associations between the stream assemblages (e.g. fish, invertebrates, diatoms) and distance-weighted land use variables have been shown for urban and agricultural land use (King et al. 2005, Van Sickle and Burch Johnson 2008). In this study, the unweighted proportion of forest harvest in the catchment was compared to an inverse-distance-weighted proportion of forest harvest in the catchment. The inversedistance-weighted forest harvest variable did not explain a higher proportion of variation in the community structure of family-level data. The result suggests that there are not areas of greater influence in the catchment and that invertebrate communities are responding to forest harvesting in the whole catchment rather than an ―effective catchment area‖ specific to forest harvest land use (Johnson and Host 2010). The result may also suggest that a Euclideandistance based weighting scheme was not an effective representation of the potential influence of forest harvesting. Nutrient leaching and sediment transport are linked to geology, precipitation and catchment morphology. It seems reasonable that the scale of influence would vary with geomorphology (e.g. a cutblock 10 km from the stream may have more influence in a steep watershed than a flat one) and that a distance-weighting scheme 50  should account for these factors. Perhaps a flow-path based or groundwater travel time based representation would be more strongly linked to community structure. Van Sickle and Burch Johnson (2008) suggest a two part weighting scheme that separates the groundwater flow component from the in-stream component. Burcher (2009) suggests using groundwater travel time zones (e.g. 30 min zone, 90 minute zone) which reflect the distance travelled by groundwater over a set time period, which reflects both the lateral and longitudinal nature of streams. In a similar approach, one recent study used topography and hydrological connectivity to define the lateral extent of the expected reach scale area of influence which they successfully linked to invertebrate taxonomic and feeding guild composition (LeCraw and Mackereth 2010). One drawback to these methods is that they require intensive quality assurance checking of GIS analyses to ensure flow distances are accurately represented, especially in catchments with low relief (Burcher 2009, Van Sickle and Burch Johnson 2008). Both the reach-scale approach by LeCraw and Mackereth (2010) and the inversedistance-weighted approach are likely to be more effective in small watersheds (Johnson and Host 2010). The large range of catchment basin sizes in this study or the lack of consideration for catchment topography and hydrological connectivity may be reasons why there was no improvement in linking the forest harvest land use with the biota using an inverse-distance-weighted approach.  4.3 The effects of forest harvesting on invertebrate communities The effects of forest harvesting act at different scales within the catchment basin. The quality of aquatic habitat is determined by the stability and complexity of the stream channel, which are influenced by changes at the catchment (hydrological and sediment regimes) and reach (large wood inputs) scales (Hassan et al. 2005). The quality and quantity of organic matter inputs are also controlled by both catchment- and reach-scale processes through surface and subsurface flows, sedimentation processes and riparian dynamics. Beechie et al. (2010) suggest that the catchment- and reach-level processes are nested so that watershedscale processes affect channel morphology which, in turn, determine habitat structure, water quality and the biotic assemblage. For example, forest harvest roads upslope of the stream act at the catchment scale through changes to surface and subsurface flows, which can lead to changes in both the flow regime and sedimentation regimes. Changes in these regimes can 51  lead to changes in stream morphology, and, ultimately, affect the quality of habitat and the biota. Results from pCCA after removing the effect of environmental variables showed that forest harvesting in the catchment and within the riparian area (reach-scale) both influenced community composition. Harvest within the riparian zone was identified as an important influence on community structure in both ordination methods. The proportion of riparian zones harvested 5 years prior to sampling was the most highly correlated of all variables to the NMDS ordination of sites (using BIO-ENV). The proportion of riparian zones harvested less than 5 years prior to sampling and the presence of riparian zones harvested less than 10 years prior to sampling did not affect the community in the same way (Figure 6). This change in harvesting effect in riparian zones over time possibly represents recovery of riparian function perhaps associated with growth of the deciduous vegetation that is likely to replace the original coniferous vegetation. Organic matter inputs can be substantially reduced following clearcut in the riparian zone, although litter inputs have been shown to recover by 10 years post harvest (Wipfli et al. 2007). Forest harvesting at the catchment scale was expected to influence the biota through changes to the hydrologic and sediment regimes. The proportion of forest harvesting in the first 15 years prior to sampling had a different effect on the community assemblage than the proportion of forest harvesting 16 to 30 years prior to sampling and more than 30 years prior to sampling. Similarly, in a long-term paired-watershed study, large differences in community structure and function still existed 16 years post-harvest, although by 26 years post-harvest, the reference and disturbed communities were much more similar (Ely and Wallace 2010). In this study, the differences in effects on the community with time since harvest might indicate partial hydrologic recovery associated with forest regrowth 15 years post-harvest. The average projected height for cutblocks 0 to 15 years old was 1.4 m compared to an 8.2 m projected height for cutblocks 15 to 29 years old, and a 13.6 m projected height for cutblocks greater than 30 years of age (TFL data only). A review by Moore and Wondzell (2005) suggested that there may be large gains in hydrologic recovery in the first 10 to 15 years (60-75%) due to initial forest regrowth, with full recovery requiring a much longer time interval (e.g. 30 to 80 years). 52  Both family-level taxa richness and relative abundance metrics were significantly greater with increased forest harvest ( ≥ 15%). Increases in invertebrate abundance and richness suggest an increase in local productivity perhaps due to nutrient enrichment through subsurface flows, increased solar radiation promoting algal growth or changes to quality or quantity of litter inputs or possibly increased heterogeneity of the stream (e.g. increased habitat complexity). As an example, an increase in taxa richness was related to a change in the bioavailability of litter (from coniferous to deciduous) in Pacific Northwest streams (Kominoski et al. 2010). However, there are many other reported cases of decreases in taxa richness with forest harvest (Martel et al. 2007, Kreutzweiser et al. 2008b, Zhang et al. 2009) perhaps due to sedimentation effects. Forest harvesting can have direct and indirect effects on the stream habitat and food availability and, depending on the specific balance of changes in some processes and the geomorphology of the stream, we might expect different outcomes in terms of the effect on the biota. For example, Davies et al. (2009) found that forest harvest led to changes in channel morphology that lead to decreased FPOM storage potential and a deficit in FPOM, which ultimately decreased invertebrate abundance. However, Danehy et al. (2007) found that forest harvest led to increased total nitrogen instream which lead to increased invertebrate abundance. It follows that forest harvest effects on invertebrate abundances in the literature are variable, with some reports of increased abundances (Stone and Wallace 1998, Nislow and Lowe 2006, Danehy et al. 2007, Martel et al. 2007) and biomass (Fuchs et al. 2003), and some studies reporting decreased abundances (Davies et al. 2005, Zhang et al. 2009). Taxonomic richness has been shown to decrease with logging while numerical density increased (Martel et al. 2007). In particular, Martel et al. (2007) found a marked decrease in the abundance of Trichopterans and an increase in the abundance of chironomids with forestry activities. Catchments with higher forest harvest intensity had lower proportions of EPT individuals but there was no change in EPT taxa richness, Trichoptera richness or proportion of Trichopterans. Stone and Wallace (1998) found no change in EPT richness with logging but Martel et al. (2007) reported a significant decrease in Trichoptera richness. In general, there seems to be support for a lack of change in richness (Stone and Wallace 1998) and decrease 53  in proportions of total EPT (Danehy et al. 2007) or proportions of individual taxa from the EPT orders (Kreutzweiser et al. 2008b) with forest harvesting. Community structure (mean NMDS ordination axis scores) was significantly different between undisturbed sites and sites in catchments with ≥ 15% forest harvest along both axes. Changes in distance-based measures of community structure are often reported with logging (Davies et al. 2005, Kreutzweiser et al. 2008b, Reid et al. 2010) although in a study of 167 first-order streams there was no association between NMDS scores and logging activity (Herlihy et al. 2005). In one study, the change in community structure was accounted for by increases in Diptera, Mollusca and Oligochaetes, and decreases in Ephemeroptera, that were linked to increases in epilithon biomass, temperature and proportion of fines in the streambed (Reid et al. 2010). The results are positive for the reference condition approach to predictive bioassessment suggesting that this type of bioassessment tool is sensitive to changes in community structure with forest harvest land use. Results from this study and one other (Danehy et al. 2007) suggest that Simpson‘s diversity and evenness are not effective indicators of forest harvesting. A lack of change in diversity and evenness measures coupled with a shift in the community assemblage suggests that as taxa are lost from the stream invertebrate community, they are replaced by others. Diversity and evenness might be better indicators in systems where there is no replacement of taxa, or a shift in dominance by a few tolerant taxa as might be expected with other types of land use (e.g. in the presence of toxic stressors). There was no change in the proportion of chironomids with increased proportion of forest harvest in this study. Chironomids are generally thought to be insensitive to sedimentation, likely because of compensation by tolerant species replacing more sensitive species within the family (Shaw and Richardson 2001, Kreutzweiser et al. 2005). Large increases in the proportion of Chironomidae have been reported within the first couple of years of disturbance and with only low to moderate intensities of forest harvesting (Kiffney et al. 2003, Nislow and Lowe 2006, Martel et al. 2007). Martel et al. (2007) suggest that chironomids thrive in a changing environment because they are small, have short life cycles and are generalist feeders. It is possible that the coarse age classes used in this study (e.g. % 54  forest harvest 0 to 14 years prior to sampling) masked short term increases in chironomids (less than 15 years), or that the family-level taxonomic resolution masked changes at the genus-level that might be important. Several families of invertebrates were correlated with the gradient of forest harvest from undisturbed sites, to sites with > 15% harvest in the catchment. Taeniopterygidae, Rhyacophilidae, Leptophlebiidae, Elmidae, Capniidae, Chloroperlidae, Hydropsychidae and Perlidae families have potential as indicators of forest harvest. Taeniopterygidae is more abundant in undisturbed streams while the other families increased with forest harvest disturbance.  4.4 Covariation of forest harvest and environmental gradients Stream networks have well described abiotic and biotic longitudinal patterns due to changes in channel morphology and discharge, from the stream headwaters to the mouth (Lowe et al. 2006). For example, water chemistry (e.g. dissolved salts, nutrient levels and pH) generally changes along a longitudinal gradient (from headwaters to mouth), reflecting changes in vegetation, geology, soils, climate and with anthropogenic influences (Giller and Malmqvist 1998). Covariation of anthropogenic gradients with natural gradients is relatively common (e.g. Richards et al. 1997, Yates and Bailey 2010) and can make interpretation of disturbance effects more difficult. In this study, tests of environmental gradients across forest harvest LAND USE groups illustrated that the natural and anthropogenic gradients were highly confounded. Relative to unharvested areas, forest harvesting activities were concentrated in larger, lower elevation Pacific Maritime ecozone catchment basins with gentler hillslopes, less relief, and more lakes and wetlands, and where the stream channels were lower gradient with a low proportion of habitat as cascades. While there were significant differences in latitude and surface geology between forest harvesting LAND USE groups, the trends were not linear although they tend to show that forest harvest was concentrated in lower latitude, less rapidly drained catchments with a lower proportion of intrusive rock. All of these environmental gradients also influenced invertebrate community composition (as shown through pCCA) except for the catchment drainage characteristics. 55  The gradient representing the proportion of intrusive rock, the dominant drainage type and the dominant local surface form (GEOpc2) varied significantly between forest harvesting LAND USE groups, but did not explain a significant portion of the variation in the biota. This suggests that these catchment properties do not shape the stream communities. Similarly, LOCpc3, representing an instream habitat gradient from fewer pools and abundant macrophytes to more pools and fewer macrophytes, varied significantly between forest harvesting LAND USE groups but was not significantly related to variation in the stream biota. Two other variables, LOCpc4 and GEOpc3, explained some variation in the stream biota, but were not significantly different across forest harvest LAND USE groups. LOCpc4 represented mature coniferous riparian forest with mature tree cover greater than 10%. It was not surprising that this variable did not vary with forest harvest land use since it was measured over a length of stream 6 times the bankfull width at each sampling site, and most streams were not sampled directly in cutblocks. GEOpc3 represented a gradient of catchment area dominated by ultramafic to metamorphic surface geology. As a next step, it would be useful to know how the intensity of forest harvesting disturbance and recovery vary along the natural gradients. While the results from the previous section showed significant changes in taxa richness and abundance metrics, these results may be reflecting changes in environmental gradients rather than forest harvest disturbance.  4.5 Implications for bioassessment Significant differences in mean NMDS ordination axis scores with forest harvest ≥ 15% in a catchment provides support for the use of reference condition approach based predictive models (e.g. BEnthic Assessment of SedimenT, Reynoldson et al. 1995) for assessing forest harvest impacts. The reference condition approach clusters undisturbed sites based on community assemblages (Reynoldson et al. 1995) or taxa occurrence (Parsons and Norris 1996) and then uses a discriminant function analysis to find the best subset of abiotic variables (called predictor variables) to discriminate between the biologically based clusters and predict potentially disturbed ―test‖ sites to a proper reference group. Each individual test 56  site is plotted in NMDS ordination space with the sites in the reference group in which it was predicted to belong, and the community assemblage is expected to fall inside a 90% probability ellipse drawn around the reference sites for the predicted group (Reynoldson et al. 1997). Overall, the results of this study suggest a significant link between stream invertebrate community structure and forest harvesting. However, these differences could be masked by natural variability of environmental drivers. Understanding the environmental drivers of community structure in streams is important for using the assemblage as an indicator of environmental condition in bioassessment. It would be useful to consider the catchment and reach-level variables that drive natural variation in community composition when choosing reference sites. If the bioassessment model covers a broad spatial area, it is important that reference sites cover the same range of physical and environmental conditions as the potentially disturbed sites. In addition, including a broad suite of abiotic variables in the model building (discriminant function analysis) would be important. Feist et al. (2010) found that a single set of abiotic predictors for their study area was not achievable because landscape predictors of salmon abundance varied by basin. This emphasizes the importance of understanding the environmental drivers specific to an area of interest. In this study, the important landscape, reach, and geology variables covaried with forest harvesting. In a study of streams in Ontario, Yates and Bailey (2010) found a similar covariation of natural and anthropogenic gradients and provided an argument for the use of predictive bioassessment models which classify reference sites by the biotic assemblage, removing the need for landscape-level stratification of the sites (e.g. by ecoregion). However, it is not clear whether the biologically-based reference groups would eliminate the potential of the key abiotic variables (e.g. topographic relief) to mask variability in community structure related to land uses. This is likely a function of the geographic scope of the predictive model, the reference site coverage within the project area, and the abiotic variables selected for the predictive modelling step of the bioassessment model building stage.  4.6 Future directions There are many questions still unanswered, including whether or not the effects of forest harvesting are greater for some types of streams than others? In addition, it would be useful 57  to assess the effectiveness of biological traits in this type of study as it seems likely that the factors that shape stream communities (e.g. habitat complexity, food availability, disturbance regimes) would vary along the natural gradient. Answering these questions would help refine methods for detecting forest harvest effects using bioassessment approaches.  58  5 Conclusions Family-level richness and relative abundance of stream invertebrates were significantly greater in catchments with > 15% forest harvest, compared to unharvested catchments. The results suggest that stream communities may benefit through increased food availability or habitat diversity under the observed forest harvest disturbance regimes. Both natural abiotic and forest harvesting variables were shown to be important in shaping stream invertebrate communities. Forest harvesting at the catchment and riparian scales influenced stream community structure. Inverse-distance-weighting the forest harvest in the catchment was not linked more strongly to the aquatic community as predicted. For environment variables, landscape-scale variables were more important than reach-scale variables in this study, although this may reflect the poor representation of reach scale variables in the analyses. Topography, basin area and land cover (% lakes, % wetlands, % barren), elevation and ecozone were important at the catchment scale. At the reach scale, % cascades and stream gradient were important. Forest harvesting land use covaried with this environmental gradient, and the effects of forest harvest are potentially masked by the natural variation in communities along the environmental gradient. Relative to unharvested areas, forest harvesting activities were concentrated in larger, lower elevation Pacific Maritime ecozone catchment basins with gentler hillslopes, less relief, and more lakes and wetlands; and where the stream channels were lower gradient with a low proportion of habitat as cascades. The results stress the importance of regional calibration of bioassessment tools so that land use effects are not masked by regional gradients.  59  References Allan, J.D. 2004. Influence of land use and landscape setting on the ecological status of rivers. Limnetica 23(3-4):187-198. Arscott, D.B., J.K. Jackson, and E.B. Kratzer. 2006. Role of rarity and taxonomic resolution in a regional and spatial analysis of stream macroinvertebrates. J. N. Am. Benthol. Soc., 25(4):977-997. Bailey, R.C., M.G. Kennedy, M.Z. Dervish and R.M. Taylor. 1998. Biological assessment of freshwater ecosystems using a reference condition approach: comparing predicted and actual benthic invertebrate communities in Yukon streams. Freshwater Biology 39:765-774. Bailey, R.C., R.H. Norris and T.B. Reynoldson. 2004. Bioassessment of Freshwater Ecosystems – Using the Reference Condition Approach. Kluwer Academic Publishers, Boston. 170 pp. Bailey, R.C., Reynoldson, T.B., Yates, A.G., Bailey, J. and S. Linke. 2007. Integrating stream bioassessment and landscape ecology as a tool for landuse planning. Freshwater Biology 52:908-917. Banks, J.L., J. Li and A. Herlihy. 2007. Influence of clearcut logging, flow duration, and season on emergent aquatic insects in headwater streams of the central Oregon Coast Range. J. N. Am. Benthol. Soc. 26(4):620-632. Beechie, T.J., D.A. Sear, J.D. Olden, G.R. Pess, J.M. Buffington, H.Moir, P.Roni, and M.M. Pollock. 2010. Process-based principles for restoring river ecosystems. Bioscience 60(3):209-222. Beisner, B.E., P.R. Peres-Neto, E.S. Lindstrom, A. Barnett, and M.L. Longhi. 2006. The role of environmental and spatial processes in structuring lake communities from bacteria to fish. Ecology 87(12):2985-2991. 60  Borcard, D., P. Legendre and P. Drapeau. 1992. Partialling out the spatial component of ecological variation. Ecology 73(3):1045-1055. Bouden, J.J., and C.M. Buddle. 2010. Spider assemblages across elevational and latitudinal gradients in the Yukon Territory, Canada. Arctic 63(3): 261-272. Brazner, J.C., N.P. Danz, G.J. Niemi, R.R. Regal, A.S. Trebitz, R.W. Howe, J.M. Hanowski, L.B. Johnson, J.J.H. Ciborowski, C.A. Johnston, E.D. Reavie, V.J. Brady, and G.V. Sgro. 2007. Evaluation of geographic, geomorphic and human influences on Great lakes wetland indicators: a multi-assemblage approach. Ecological Indictors 7:610635. Brosse, S., C.J. Arbuckle and C.R. Townsend. 2003. Habitat scale and biodiversity: influence of catchment, stream reach and bedform scales on local invertebrate diversity. Biodiversity and Conservation 12:2057-2075. Burcher, C.L. 2009. Using simplified watershed hydrology to define spatially explicit ‗zones of influence‘. Hydrobiologia 618:149-160. Chizinski, C.J., B. Vondracek, C.R. Blinn, R.M. Newman, D.M. Atuke, K.Fredricks, N.A. Hemstad, E. Merten and N. Schlesser. 2010. The influence of partial timber harvesting in riparian buffers on macroinvertebrate and fish communities in small streams in Minnesota, USA. Forest Ecology and Management 259:1946-1958. Clarke, K.R., and M. Ainsworth. 1993. A method of linking multivariate community structure to environmental variables. Mar. Ecol. Prog. Ser. 92:205-219. Clarke, K.R., and R.M. Warwick. 2001. Change in marine communities: an approach to statistical analysis and interpretation, 2nd edition. PRIMER-E: Plymouth. Clarke, K.R., P.J. Somerfield, and R.N. Gorley. 2008. Testing of null hypothesis in exploratory community analyses: similarity profiles and biota-environment linkage. Journal of Experimental Marine Biology and Ecology 366:56-69.  61  Cummins, K.W., and R.W. Merritt. 1996. Ecology and distribution of aquatic insects. In An Introduction to the Aquatic Insects of North America, Third Edition. R.W. Merritt and K.W. Cummins (editors). Kendall/Hunt Publishing Company, Dubuque, Iowa. 862 pp. Danehy, R.J., S.S. Chan, G.T. Lester, R.B. Langshaw, and T.R. Turner. 2007. Periphyton and macroinvertebrate assemblage structure in headwaters bordered by mature, thinned and clearcut douglas-fir stands. Forest Science 53(2): 294-307. Dauwalter, D.C., D.K. Splinter, W.L. Fisher, and R.A. Marston. 2008. Biogeography, ecoregions, and geomorphology affect fish species composition in streams of eastern Oklahoma, USA. Environ. Biol. Fish 82: 237-249. Davies, P.E., L.S.J. Cook, P.D. McIntosh and S.A. Munks. 2005. Changes in stream biota along a gradient of logging disturbance, 15 years after logging at Ben Nevis, Tasmania. Forest Ecology and Management 219: 132-148. Davis, W.S. 1994. Biological Assessment and Criteria: Building on the Past. In Biological Assessment and Criteria: Tools for Water Resource Planning and Decision Making. W.S. Davis and T.P. Simon (editors). CRC Press, Inc., Boca Raton, Florida, USA. Death, R.G., and M.J. Winterbourn. 1995. Diversity patterns in stream benthic invertebrate communities: the influence of habitat stability. Ecology 76(5):1446-1460. Ecological Stratification Working Group 1995. A National Ecological Framework for Canada. Agriculture and Agri-Food Canada, Research Branch, Centre for Land and Biological Resources Research and Environment Canada, State of the Environment Directorate, Ecozone Analysis Branch, Ottawa/Hull. Report and national map at 1: 7 500 000 scale. Elosegi, A., J. Diez and M. Mutz. 2010. Effects of hydromorphological integrity on biodiversity and functioning of river ecosystems. Hydrobiologia 657:199-215.  62  Ely, D.T., and J.B. Wallace. 2010. Long-term functional group recovery of lotic macroinvertebrates from logging disturbance. Can. J. Fish. Aquat. Sci. 67:11261134. Feist, B.E., E.A. Steel, D.W. Jensen, and D.N.D. Sather. 2010. Does the scale of our observational window affect our conclusions about correlations between endangered salmon populations and their habitat? Landscape Ecol. 25:727-743. Feld, C.K., and D. Hering. 2007. Community structure or function: effects of environmental stress on benthic macroinvertebrates at different spatial scales. Freshwater Biology 52:1380-1399. French, L.J., G.F. Smith, D.L. Kelly, F.J.G. Mitchell, S. O‘Donoghue, S.F. Iremonger, and A. McKee. 2008. Ground flora communities in temperature oceanic plantation forests and the influence of silvicultural, geographic and edaphic factors. Forest Ecology and Management 255:476-494. Frissell, C.A., W. J. Liss, C.E. Warren, and M.D. Hurley. 1986. A hierarchical framework for stream habitat classification: viewing streams in a watershed context. Environmental Management 10(2):199-214. Fuchs, S.A., S.G. Hinch, and E. Mellina. 2003. Effects of streamside logging on stream macroinvertebrate communities and habitat in the sub-boreal forests of British Columbia, Canada. Can. J. For. Res. 33:1408-1415. Giller, P.S., and B. Malmqvist. 1998. The biology of streams and rivers. Oxford University Press Inc., Oxford, New York. Gomi, T., R.D. Moore, and A.S. Dhakal. 2006. Headwater stream temperature response to clear-cut harvesting with different riparian treatments, coastal British Columbia, Canada. Water Resour. Res. 42, W08437, doi:10.1029/2005WR004162. Goodman, K.J., M.A. Baker and W.A. Wurtsbaugh. 2010. Mountain lakes increase organic matter decomposition rates in streams. J. N. Am. Benthol. Soc. 29(2):521-529. 63  Hassan, M.A., M. Church, T.E. Lisle, F. Brardinoni, L. Benda, and G.E. Grant. 2005. Sediment transport and channel morphology of small, forested streams. Journal of the American Water Resources Association (JAWRA) 41(4):853-876. Hauer, F. R., and V.H. Resh. 1996. Benthic macroinvertebrates. In Methods in Stream Ecology F.R. Hauer and G.A. Lamberti (editors). Academic Press Ltd., San Diego, California, USA. Heino, J. 2008. Influence of taxonomic resolution and data transformation on biotic matrix concordance and assemblage-environment relationships in stream macroinvertebrates. Boreal Environment Research 13:359-369. Herlihy, A.T., W. J. Gerth, J. Li and J.L. Banks. 2005. Macroinvertebrate community response to natural and forest harvest gradients in western Oregon headwater streams. Freshwater Biology 50:905-919. Hernandez, O., R.W. Merritt, and M.S. Wipfli. 2005. Benthic invertebrate community structure is influenced by forest succession after clearcut logging in southeastern Alaska. Hydrobiologia 533:45-59. Hobbs, R.J., and L.F. Huenneke. 1992. Disturbance, diversity, and invasion: implications for conservation. Conservation Biology 6(3):324-337. Hoover, T.M., L.B. Marczak, J.S. Richardson and N. Yonemitsu. 2010. Transport and settlement of organic matter in small streams. Freshwater Biology 55:436-449. Imholt, C., C.N. Gibbons, I.A. Malcolm, S. Langan and C. Soulsby. 2010. Influence of riparian cover on stream temperatures and the growth of the mayfly Baetis rhodani in an upland stream. Aquatic Ecology 44:669-678. Jackson, C.R., D.P. Batzer, S.S. Cross, S.M. Haggerty, and C.A. Sturm. 2007. Headwater streams and timber harvest: channel, macroinvertebrate, and amphibian response and recovery. Forest Science 53(2):356-370.  64  Johnson, L.B., and G.E. Host. 2010. Recent developments in landscape approaches for the study of aquatic ecosystems. J. N. Am. Benthol. Soc. 29(1):41-66. Johnson, R.K., M.T. Furse, D. Hering and L. Sandin. 2007. Ecological relationships between stream communities and spatial scale: implications for designing catchmentlevel monitoring programmes. Freshwater Biology 52:939-958. Johnson, R.K., W. Goedkoop and L. Sandin. 2004. Spatial scale and ecological relationships between the macroinvertebrate communities of stony habitats of streams and lakes. Freshwater Biology 49:1179-1194. Karr, J.R., and E.W. Chu. 1999. Restoring Life in Running Waters: Better Biological Monitoring. Island Press, Covelo, California. 206 pp. Kiffney, P.M., J.S. Richardson, and J.P. Bull. 2003. Responses of periphyton and insects to experimental manipulation of riparian buffer width along forest streams. Journal of Applied Ecology 40:1060-1076. King, R.S., M.E. Baker, D.F. Whigham, D.E. Weller, T.E. Jordan, P.F. Kazyak, and M.K. Hurd. 2005. Spatial considerations for linking watershed land cover to ecological indicators in streams. Ecological Applications 15(1):137-153. Kominoski, J.S., L.B. Marczak, and J.S. Richardson. In press. Riparian forest composition affects stream litter decomposition despite similar microbial and invertebrate communities. Ecology. [doi:10.1890/10-0028.1]. Kratzer, E.B., J.K. Jackson, D.B. Arscott, A.K. Aufdenkampe, C.L. Dow, L.A. Kaplan, J.D. Newbold and B.W. Sweeney. 2006. Macroinvertebrate distribution in relation to land use and water chemistry in New York City drinking-water-supply watersheds. J. N. Am. Benthol. Soc. 25(4):954-976. Kreutzweiser, D., E. Muto, S. Holmes, and J. Gunn. 2010. Effects of upland clearcutting and riparian partial harvesting on leaf pack breakdown and aquatic invertebrates in boreal forest streams. Freshwater Biology 55:2238-2252. 65  Kreutzweiser, D.P., K.P. Good, S.S. Capell, and S.B. Holmes. 2008b. Leaf-litter decomposition and macroinvertebrate communities in boreal forest streams linked to upland logging disturbance. J. N. Am. Benthol. Soc. 27(1):1-15. Kreutzweiser, D.P., P.W. Hazlett and J.M. Gunn. 2008a. Logging impacts on the biogeochemistry of boreal forest soils and nutrient export to aquatic systems: a review. Environ. Res. 16:157-179. Kreutzweiser, D.P., S.S. Capell and K.P. Good. 2005. Macroinvertebrate community responses to selection logging in riparian and upland areas of headwater catchments in a northern hardwood forest. J. N. Am. Benthol. Soc. 24(1):208-222. Lachance, D., and C. Lavoie. 2004. Vegetation of Sphagnum bogs in highly disturbed landscapes: the relative influence of abiotic and anthropogenic factors. Applied Vegetation Science 7:183-192. Lake, P.S. 2000. Disturbance, patchiness and diversity in streams. J. N. Am. Benthol. Soc. 19(4):573-592. Larsen, S., I.P., Vaughan and S.J. Ormerod. 2009. Scale-dependent effects of fine sediments on temperate headwater invertebrates. Freshwater Biology 54:203-219. Lecerf, A., and J.S. Richardson. 2010. Litter decomposition can detect effects of high and moderate levels of forest disturbance on stream condition. Forest Ecology and Management 259(12):2433-2443. Lecraw, R., and R. Mackereth. 2010. Sources of small-scale variation in the invertebrate communities of headwater streams. Freshwater Biology 55:1219-1233. Lenat, D.R. 1988. Water quality assessment of streams using a qualitative collection method for benthic macroinvertebrates. J. N. Am. Benthol. Soc. 7(3):222-233.  66  Litschert, S.E., and L.H. MacDonald 2009. Frequency and characteristics of sediment delivery pathways from forest harvest units to streams. Forest Ecology and Management 259:143-150. Lowe, W.H., G.E. Likens and M.E. Power. 2006. Linking scales in stream ecology. Bioscience 56(7):591-597. Manel, S., S.T. Buckton and S.J. Ormerod. 2000. Testing large-scale hypotheses using surveys: the effects of land use on the habitats, invertebrates and birds of Himalayan rivers. Journal of Applied Ecology 37:756-770. Martel, N., M. A. Rodriguez and P. Berube. 2007. Multi-scale analysis of responses of stream macrobenthos to forestry activities and environmental context. Freshwater Biology 52:85-97. Mazor, R. D., T. B. Reynoldson, D. M. Rosenberg, and V. H. Resh. 2006. Effects of biotic assemblage, classification and assessment method on bioassessment performance. Can. J. Fish. Aquat. Sci. 63:394-411. McCune, B., J.B. Grace and D.L. Urban. 2002. Analysis of Ecological Communities. MJM Software Design, Gleneden Beach, Oregon, USA. Melody, K. J., and J. S. Richardson. 2007. Riparian forest harvesting and its influence on benthic communities of small streams of sub-boreal British Columbia. Can. J. For. Res. 37:907-918. Merritt, R.W., and K.W. Cummins. 1996. An introduction to the aquatic insects of North America. Third Edition. Kendall/Hunt Publishing Company, Dubuque, Iowa. Moldenke A.R., and C. Ver Linden. 2007. Effects of clearcutting and riparian buffers on the yield of adult aquatic macroinvertebrates from headwater streams. Forest Science 53(2):308-319.  67  Montgomery, D.R. 1999. Process domains and the river continuum. Journal of the American Water Resources Association 35(2):397-410. Moore, R.D., and S. M. Wondzell. 2005. Physical hydrology and the effects of forest harvesting in the Pacific Northwest: a review. Journal of the American Water Resources Association (JAWRA) 41(4):763-784. Nislow, K.H., and W.H. Lowe. 2006. Influences of logging history and riparian forest characteristics on macroinvertebrates and brook trout (Salvelinus fontinalis) in headwater streams (New Hampshire, U.S.A.). Freshwater Biology 51:388-397. Paller, M.H., W.L. Specht and S.A. Dyer. 2006. Effects of stream size on taxa richness and other commonly used benthic bioassessment metrics. Hydrobiologia 568:309-316. Parsons, M., and R.H. Norris. 1996. The effect of habitat-specific sampling on biological assessment of water quality using a predictive model. Freshwater Biology 36:419434. Perrin, C.J., S. Bennett, S. Linke, A.J. Downie, G. Tamblyn, B. Ells, I. Sharpe, and R.C. Bailey. 2007. Bioassessment of streams in north-central British Columbia using the reference condition approach. Report prepared by Limnotek Research and Development Inc. and BC Ministry of Environment for the BC Forest Science Program. 121p. Peterson, E.E., F. Sheldon, R. Darnell, S.E. Bunn, and B. D. Harch. In press. A comparison of spatially explicit landscape representation methods and their relationship to stream condition. Freshwater Biology, no. doi: 10.1111/j.1365-2427.2010.02507.x. Poff, N.L. 1997. Landscape filters and species traits: towards mechanistic understanding and prediction in stream ecology. J. N. Am. Benthol. Soc. 16(2):391-409. Reece, P.F., T.B. Reynoldson, J.S. Richardson and D.M. Rosenberg. 2001. Implications of seasonal variation for biomonitoring with predictive models in the Fraser River catchment, British Columbia. Can. J. Fish. Aquat. Sci. 58:1411-1418. 68  Reid, D.J., J.M. Quinn and A.E. Wright-Stow. 2010. Responses of stream macroinvertebrate communities to progressive forest harvesting: influences of harvest intensity, stream size and riparian buffers. Forest Ecology and Management 260:1804-1815. Resh, V.H., A.V. Brown, A.P. Covich, M.E. Gurtz, H.W. Li, W.Minshall, S.R. Reice, A. L. Sheldon, J.B. Wallace and R.C. Wissmar. 1988. The role of disturbance in stream ecology. J. N. Am. Benthol. Soc. 7(4):433-455. Reynoldson, T.B., D.M. Rosenberg, and V.H. Resh. 2001. Comparison of models predicting invertebrate assemblages for bio-monitoring in the Fraser River catchment, British Columbia. Can. J. Fish. Aquat. Sci. 58:1395-1409. Reynoldson, T.B., R.C. Bailey, K.E. Day and R.H. Norris. 1995. Biological guidelines for freshwater sediment based on BEnthic Assessment of SedimenT (the BEAST) using a multivariate approach for predicting biological state. Australian Journal of Ecology 20:198-219. Reynoldson, T.B., R.H. Norris, R.H. Resh, K.E. Day, and D.M. Rosenberg. 1997. The reference condition approach: a comparison of multimetric and multivariate approaches to assess water-quality impairment using benthic macroinvertebrates. J. N. Am. Benthol. Soc. 16(4):833-852. Richards, C., R.J. Haro, L.B Johnson, and G.E. Host. 1997. Catchment and reach-scale properties as indicators of macroinvertebrate species traits. Freshwater Biology 37:219-230. Richardson, J.S. 1992. Food, microhabitat, or both? Macroinvertebrate use of leaf accumulations in a montane stream. Freshwater Biology 27:169-176. Richardson, J.S. 2008. Aquatic arthropods and forestry: effects of large-scale land use on aquatic systems in Nearctic temperate regions. Can. Entomol. 140:495-509.  69  Richardson, J.S., R.E. Bilby, and C.A. Bondar. 2005. Organic matter dynamics in small streams of the Pacific Northwest. Journal of the American Water Resources Association (JAWRA) 41(4):921-934. Richardson, J.S., Y. Zhang, and L.B. Marczak. 2010. Resource subsidies across the landfreshwater interface and responses in recipient communities. River Research and Applications 26:55-66. Sandin, L., and R.K. Johnson. 2004. Local, landscape and regional factors structuring benthic macroinvertebrate assemblages in Swedish streams. Landscape Ecology 19: 501-514. Shaw, E.A., and J.S. Richardson. 2001. Direct and indirect effects of sediment pulse duration on stream invertebrate assemblages and rainbow trout (Oncorhynchus mykiss) growth and survival. Can. J. Fish. Aquat. Sci. 58:2213-2221. Smith, B.J., P.E. Davies, and S.A. Munks. 2009. Changes in benthic invertebrate communities in upper catchment streams across a gradient of catchment forest operation history. Forest Ecology and Management 257: 2168-2174. Splinter, D.K., D.C. Dauwalter, R.A. Marston, and W.L. Fisher. 2010. Ecoregions and stream morphology in eastern Oklahoma. Geomorphology 122:117-128. Statzner, B., and L.A. Bêche. 2010. Can biological invertebrate traits resolve effects of multiple stressors on running water ecosystems? Freshwater Biology 55(Suppl. 1):80-119. Stone, M.K., and J.B. Wallace. 1998. Long-term recovery of a mountain stream from clearcut logging: the effects of forest succession on benthic invertebrate community structure. Freshwater Biology 39:151-169. Systat Software, Inc. 2004. Systat for Windows, version 11. SYSTAT Software, Inc., Chicago, IL, U.S.A.  70  ter Braak, C.J.F., and P. Smilauer. 2002. CANOCO Reference manual and CanoDraw for Windows User‘s guide: Software for Canonical Community Ordination (version 4.5). Microcomputer Power (Ithaca, NY, USA), 500pp. Thompson, R.M., N.R. Phillips, and C.R. Townsend. 2009. Biological consequences of clear-cut logging around streams – Moderating effects of management. Forest Ecology and Management 257: 931-940. Townsend, C.R., S. Doledec, R. Norris, K. Peacock and C. Arbuckle. 2003. The influence of scale and geography on relationships between stream community composition and landscape variables: description and prediction. Freshwater Biology 48:768-785. Underwood, A.J., M.G. Chapman, and S.D. Connell. 2000. Observations in ecology: you can‘t make progress on processes without understanding the patterns. Journal of Experimental Marine Biology and Ecology 250:97-115. Van Sickle, J., and C. Burch Johnson. 2008. Parametric distance weighting of landscape influence on streams. Landscape Ecology 23:427-438. Vandvik, V., and H.J.B. Birks. 2002. Partitioning floristic variance in Norwegian upland grasslands into within-site and between-site components: are the patterns determined by environment or by land-use? Plant Ecology 162:233-245. Vinson, M.R., and C.P. Hawkins. 1996. Effects of sampling area and subsampling procedure on comparisons of taxa richness among streams. J. N. Am. Benthol. Soc. 15(3):392-399. Vinson, M.R., and C.P. Hawkins. 2003. Broad-scale geographical patterns in local stream insect genera richness. Ecography 26:751-767. Ward, J.V. 1998. Riverine landscapes: biodiversity patterns, disturbance regimes and aquatic conservation. Biological Conservation 83(3):269-278.  71  Wilkerson, E., J.M. Hagan, and A. A. Whitman. 2010. The effectiveness of different buffer widths for protecting water quality and macroinvertebrate and periphyton assemblages of headwater streams in Maine, USA. Can. J. Fish. Aquat. Sci. 67:177190. Wipfli, M.S., J.S. Richardson and R.J. Naiman. 2007. Ecological linkages between headwaters and downstream ecosystems: transport of organic matter, invertebrates, and wood down headwater channels. Journal of the American Water Resources Association (JAWRA) 43(1):72-85. Yates, A.G., and R.C. Bailey. 2010. Covarying patterns of macroinvertebrate and fish assemblages along natural and human activity gradients: implications for bioassessment. Hydrobiologia 637:87-100. Zhang, Y., J.S. Richardson and X. Pinto. 2009. Catchment-scale effects of forestry practices on benthic invertebrate communities in Pacific coastal streams. Journal of Applied Ecology 46:1292-1303. Zheng, L., J. Gerritsen, J. Beckman, J. Ludwig and S. Wilkes. 2008. Land use, geology, enrichment, and stream biota in the eastern ridge and valley Ecoregion: implications for nutrient enrichment criteria development. Journal of the American Water Resources Association (JAWRA) 44(6):1521-1536.  72  Appendix A Habitat variables and descriptions The following table of habitat variables and descriptions is from Perrin et al. (2007). Variable  Method or Standard Used  Weather Current and recent weather  Documentation based on personal observation; choices include: storm, rain, showers, overcast, and sunny  General Site Information  Location  A hand-held global positioning system was used to record latitude and longitude (in decimal degrees) and elevation. A site description was noted and site diagram completed.  Photographs  Taken looking upstream, downstream, across, and at the substrate (substrate photo included a 50cm quadrat for scale)  Water Quality Air Temperature  Recorded from thermometer or handheld meter  Water Temperature  Recorded from handheld meter (YSI model 63)  pH  Recorded from handheld meter (YSI model 63)  Specific Conductance  Recorded from handheld meter (YSI model 63)  Dissolved Oxygen  Recorded from handheld meter (Handy MK II) Collected according to sampling procedures outlined in the B.C. Ambient Freshwater and Effluent Sampling Manual (RIC, 1997). Sample bottles included:  Water Samples  1 L for general ions (alkalinity, chloride, true colour, sulphate, total dissolved solids, total suspended solids and turbidity), 250 mL for nutrients (ammonia, nitrite, nitrate, organic and total nitrogen; orthophosphorus, total phosphorus), 250 mL for total organic carbon (H2SO4 preservative added in the field), and 250 mL for total metals (collected in acid washed bottle, with HNO3 preservative added in the field). All water samples were sent to Maxxam Analytics Inc. by overnight courier, to achieve the recommended 72hr holding time. Metals analysis was done using an ICPMS scan, and all parameters were analyzed according to methods in the B.C. Environmental Laboratory Manual (Horvath, 2005).  Vegetation Cover  73  Variable  Method or Standard Used  Percent Cover  An estimate of the % of the wetted surface area that is covered (within 1 m of the water surface) by each of the following: woody debris, boulders, undercut banks, deep pools, and overhanging vegetation. Modified from the Reconnaissance (1:20 000) Fish and Fish Habitat Inventory: Standards and Procedures - Version 2.0 (RIC, 2001).  Macrophyte Coverage  An estimate of the % of the streambed that is covered by macrophytes (Reynoldson et al. 2004).  Periphyton Coverage  An estimate of periphyton coverage in the running water, using a scale of 1 to 5  Disturbance Indicators  Misc. disturbance indicators  Documentation of the presence of scour, sediment wedges, extensive riffles and limited pools, unvegetated and mid-channel bars, multiple channels, eroding banks, isolated sidechannels, recently formed large woody debris (LWD) jams and LWD parallel to banks. Modified from RIC (2001).  Stream Channel Characteristics Gradient  Measured using a clinometer and reported in %  Habitat Units  An estimate of the % of the channel area within the reach that is occupied by pools, glides, riffles, and cascades.  Stream Widths  Three measurements of the wetted and bankfull width, taken from within the stream reach. Modified from RIC (2001). Values were later averaged to obtain wetted and bankfull width.  Stream Profile (for Discharge Measurement)  A discharge measurement from the kicked area, calculated from velocity and depth measurements made at 5-8 equidistant points across the stream, using a Swoffer 2001 current meter. Field methods were based on procedures in the Manual of Standard Operating Procedures for Hydrometric Surveys in B.C. (RIC, 1998), with velocity measurements taken at a depth of 0.6 x total depth. Discharge calculations followed RIC (1998).  Substrate  Composition  A visual estimate of the % of the stream reach that is covered with various particle sizes (sand, gravel, pebble, cobble, boulder and bedrock), according to the Wentworth Scale. The estimate was made by each of three members of the field crew and average values were recorded.  Embeddedness  An estimate of how embedded the cobbles are in the surrounding fines (measured in the riffle habitat).  Pebble Count  A Wolman Pebble Count (Wolman, 1954), where the intermediate diameter of 100 randomly-selected particles within the stream reach was measured using a ruler. From this data, SYSTAT 11 and MS Excel were used to calculate the median particle size (D50), geometric mean particle size (Dg)  74  Variable  Method or Standard Used and the Fredle Index (FI).  Odours/Oils  Presence documented on field sheet  Riparian Vegetation  Vegetation Types  An estimate of the % of different vegetation types (bare soil, grass/herb, shrub, deciduous, coniferous) present at the site. Modified from RIC (2001). Included documentation of species present.  Structural Stage  An estimate of structural stage, according to RIC (2001).  Canopy Closure  An estimate of canopy closure, according to RIC (2001).  Mountain Pine Beetle Infestation  An estimate of the presence of pine trees and severity of the pine beetle infestation, at the site (riparian) and in the watershed.  Land Use Observations  Observations of land use, erosion, and NPS pollution at or near the site.  Rapid Bioassessment  USEPA Rapid Bioassessment Protocol (RBP)  A scoring of 10 habitat parameters (epifaunal substrate/cover, embeddedness, velocity-depth combinations, sediment deposition, channel flow status, channel alteration, channel sinuosity, bank stability, bank vegetative protection and riparian vegetative zone width), using a scale of 1-20 according to the USEPA RBP Field Sheet (Barbour et. al., 1999).  75  References Barbour, M.T., J. Gerritsen, B.D. Snyder, and J.B. Stribling. 1999. Rapid Bioassessment Protocols for Use in Streams and Wadeable Rivers: Periphyton, Benthic Macroinvertebrates and Fish, Second Edition. EPA 841-B-99-002. U.S. Environmental Protection Agency; Office of Water; Washington, D.C. Horvath, S (editor). 2005. British Columbia Environmental Laboratory Manual — For the Analysis of Water, Wastewater, Sediment, Biological Materials and Discrete Ambient Air Samples. 2005. Water and Air Monitoring and Reporting Section, Ministry of Environment. Resources Inventory Committee (RIC). 1997. BC Ambient Freshwater and Effluent Sampling Manual. Prepared by BC Environment. Resources Inventory Committee (RIC). 1998. Manual of Standard Operating Procedures for Hydrometric Surveys in B.C. Prepared by Ministry of Environment, Lands and Parks, Resources Inventory Branch for the Aquatic Inventory Task Force Resources Inventory Committee. November 2, 1998. Version 1.1. Resources Inventory Committee (RIC). 2001. Reconnaissance (1: 20 000) Fish and Fish Habitat Inventory Standards and Procedures. Prepared by BC Fisheries Information Services Branch for the Resources Inventory Committee April, 2001. Version 2.0. Reynoldson, T.B., C. Logan, T. Pascoe, S. Thompson and S. Sylvestre. 2004. Invertebrate biomonitoring field and laboratory manual for running water habitats. Environment Canada, Canadian Aquatic Biomonitoring Network, Vancouver, BC. 53 p. Wolman. M. 1954. A method of sampling coarse bed material. American Geophysical Union, Transactions, 35:951-956.  76  Appendix B GIS variables and descriptions The following tables list the GIS variables that were available for the project. Not all variables were used. Source Path  ATTRIBUTE  DESCRIPTION  SOURCE DATABASE  ECOREG ECOZONE  Ecoregion Ecozone Area of sample site watershed Length of streams within watershed Stream order at sampling site Area of lakes Area of wetlands Area of ice Distance from sample site to upstream lake outlet (where applicable) Distance from sample site to upstream wetlands outlet (where applicable) Area of forest burned post 1988 Area of forest, burned before 2000 Area classified as agricultural Area classified as alpine Area classified as avalanche chute - below the tree line, devoid of forest growth due primarily to snow avalanches. Usually herb or shrub covered. Area classified as rock barrens, badlands, sand and gravel flats, dunes and beaches where unvegetated surfaces predominate, except where these conditions occur in ALPINE areas. Area classified as mine/extraction site Area classified as Urban Length of roads Number of Minfile mines Elevation (metres) of sampling site  CANSIS CANSIS  1 1  n/a  2  LRDW  3  LRDW LRDW LRDW LRDW  3 4 5 6  LRDW  3 and 4  LRDW  3 and 5  LRDW  7  LRDW  7  LRDW LRDW  8 8  LRDW  8  LRDW  8  LRDW LRDW LRDW LRDW  8 8 9 10  ILMB arcwhse  11  WATERSHED_AREA_HA STREAMS_LENGTH_KM STREAMS_ORDER_MAX LAKES_AREA_HA WETLANDS_AREA_HA ICE_AREA_HA  LAKES_MIN_DIST_METRES  WETLANDS_MIN_DIST_METRES FOR_BURNPAST20YRS_HA FOR_BURNPRE2000_HA AGRICUL_HA ALPINE_HA  AVALAN_HA  BARREN_HA MINING_HA URBAN_HA ROADS_LENGTH_KM MINFILE_NUM ELEV_MAX  77  ATTRIBUTE ELEV_MIN SLOPE_LT30_HA SLOPE_30_50_HA SLOPE_50_60_HA SLOPE_GT60_HA ROCK_SEDIM_HA ROCK_INTRU_HA ROCK_VOLCA_HA ROCK_METAM_HA ROCK_ULTRA_HA  PMDEP_1_HA  PMDEP_A_HA  PMDEP_C_HA  PMDEP_F_HA  PMDEP_L_HA  PMDEP_M_HA  PMDEP_R_HA  PMDEP_W_HA  CFRAG_1_HA CFRAG_A_HA CFRAG_B_HA  DESCRIPTION Maximum elevation of watershed Area of hillslope gradient of 0 to 30% Area of hillslope gradient of 30 to 50% Area of hillslope gradient of 50 to 60% Area of hillslope gradient of 60% and greater Area underlain by sedimentary rock Area underlain by intrusive rock Area underlain by volcanic rock Area underlain by metamorphic rock Area underlain by ultramafic rock Area of parent material mode of deposition = '#' (Not Applicable) Area of parent material mode of deposition = 'A' (Alluvial (Fluvial) Material) Area of parent material mode of deposition = 'C' (Colluvium) Area of parent material mode of deposition = 'F' (Fluvioglacial Glaciofluvial - Material) Area of parent material mode of deposition = 'L' (Lacustrine Material) Area of parent material mode of deposition = 'M' (Morainal Material - Till) Area of parent material mode of deposition = 'R' (Rock) Area of parent material mode of deposition = 'W' (Marine Material) Area of coarse fragment content (%) = '#' (Not Applicable) Area of coarse fragment content (%) = 'A' (<10%) Area of coarse fragment content (%) = 'B' (10-30%)  SOURCE DATABASE  Source Path  ILMB arcwhse  11  ILMB arcwhse  11  ILMB arcwhse  11  ILMB arcwhse  11  ILMB arcwhse  11  LRDW  12  LRDW  12  LRDW  12  LRDW  12  LRDW  12  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  78  ATTRIBUTE CFRAG_C_HA  DRAIN_1_HA  DRAIN_I_HA  DRAIN_M_HA  DRAIN_R_HA  DRAIN_W_HA LOCSF_1_HA  LOCSF_D_HA  LOCSF_H_HA LOCSF_I_HA LOCSF_L_HA LOCSF_M_HA LOCSF_R_HA LOCSF_S_HA LOCSF_T_HA LOCSF_U_HA  various Length of roads (km) Total Area Harvested (ha) AC1 (ha) AC2 (ha)  DESCRIPTION Area of coarse fragment content (%) = 'C' (31-65%) Area of soil drainage attribute = '#' (Not Applicable) Area of soil drainage attribute = 'I' (Imperfectly drained) Area of soil drainage attribute = 'M' (Moderately well drained) Area of soil drainage attribute = 'R' (Rapidly drained) Area of soil drainage attribute = 'W' (Well drained) Area of local surface form = '#' (Not Applicable) Area of local surface form = 'D' (Dissected - a dissected or gullied pattern providing external drainage for an area) Area of local surface form = 'H' (Hummocky or Irregular) Area of local surface form = 'I' (Inclined) Area of local surface form = 'L' (Level) Area of local surface form = 'M' (Rolling) Area of local surface form = 'R' (Ridged) Area of local surface form = 'S' (Steep) Area of local surface form = 'T' (Terraced) Area of local surface form = 'U' (Undulating) Climate variables such as temperature and precipitation for the sample site (not the watershed) Total length of roads within the watershed Total harvested area within the watershed, any date Total area of age class 1 forest (1-20years) Total area of age class 2  SOURCE DATABASE  Source Path  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  13  CANSIS  14  Digital Road Atlas VRI, RESULTS, TFL forest cover  15  VRI, TFL forest cover db VRI, TFL forest cover db  16 16  16  79  ATTRIBUTE  DESCRIPTION forest (21-40years) Total area of age class 3 forest (41-60years) Total area of age class 4 forest (61-80years) Total area of age class 5 forest (81-100years) Total area of age class 6 forest (101-120years) Total area of age class 7 forest (121-140years) Total area of age class 8 forest (14-250years) Total area of age class 9 forest (>250years) Average aspect of watershed  AC3 (ha) AC4 (ha) AC5 (ha) AC6 (ha) AC7 (ha) AC8 (ha) AC9 (ha) Aspect (avg)  Source Path Number 1  SOURCE DATABASE  Source Path  VRI, TFL forest cover db  16  VRI, TFL forest cover db  16  VRI, TFL forest cover db  16  VRI, TFL forest cover db  16  VRI, TFL forest cover db  16  VRI, TFL forest cover db  16  VRI, TFL forest cover db BC TRIM 1:20,000 digital elevation model (25m)  16  SOURCE PATH  2  http://sis.agr.gc.ca/cansis/nsdb/ecostrat/gis_data.html Skeena Bioassessment Watersheds, derived from CWB_WATERSHED_POLY and TRIM DEM  3 4 5 6  WHSE_BASEMAPPING.CWB_STREAM_NETWORKS WHSE_BASEMAPPING.CWB_LAKES_POLY WHSE_BASEMAPPING.CWB_WETLAND_POLY WHSE_BASEMAPPING.CWB_GLACIERS_POLY  7  WHSE_FOREST_VEGETATION.VEG_COMP_LYR_R1_POLY and WHSE_BASEMAPPING.BTM_PRESENT_LAND_USE_V1_SVW  8 9 10 11  WHSE_BASEMAPPING.BTM_PRESENT_LAND_USE_V1_SVW WHSE_BASEMAPPING.DRA_DIGITAL_ROAD_ATLAS_LINE_SP WHSE_MINERAL_TENURE.MINFIL_MINERAL_FILE P:\corp\arcwhse\gdbc\tdem_<250kMAP_TILE>  12 13 14  WHSE_MINERAL_TENURE.GEOL_BEDROCK_UNIT_POLY_SVW http://sis.agr.gc.ca/cansis/nsdb/slc/v3.1.1/intro.html http://sis.agr.gc.ca/cansis/nsdb/ecostrat/district/climate.html https://apps.gov.bc.ca/pub/geometadata/metadataDetail.do?recordUID=45674&recordSet =ISO19115 https://apps.gov.bc.ca/pub/geometadata/metad+ataDetail.do?recordUID=52218&recordS et=ISO19115 https://apps.gov.bc.ca/pub/geometadata/metadataDetail.do?recordUID=52583&recordSet =ISO19115  15  16  80  Cutblock Report Attributes ATTRIBUTE Cutblock ID Area within watershed (ha) Harvest Date Height Class dR1 - Connected streams (m) dR2 - Connected side channel streams (m) dR3 - Connected Lake/Wetland edges (m)  dT  dI  DESCRIPTION Unique identifier for cutblock polygon. (VRI = ), (RESULTS = ), (TFL = ) Total area of given cutblock within watershed (ha) Date of most recent harvest activity for given cutblock Current (projected) height class of layer 1 leading species in harvest polygon Length of streams within the cutblock that drain to outside of the cutblock (m) Length of side channels within the cutblock that drain to outside of the cutblock (m) Length of the edges of lakes and wetlands within the cutblock that drain to outside of the cutblock (m) For cutblocks not containing stream features (dR1 = 0), the distance from centroid of cutblock to the nearest downstream stream. For cutblocks containing water edge features (dR1 > 0), the distance from the centroid of the cutblock to the intersection of the edge of cutblock and the most downstream segment of the stream network within the cutblock. The distance along the stream network from the sampling site to the point on the stream network closest to (and most downstream from) the centroid of the cutblock.  81  Appendix C Macroinvertebrate taxa list. Identifications were completed by Danusia Dolecki, Vancouver, BC and Sue Salter, Summerland, BC. Families that were included in the reduced data set were given a taxon code.  Taxon Code  Order (unless marked by *for class or **for suborder)  Family  Enchyt  Haplotaxida  Total # of Individuals (all samples)  Occurrence (# of sites)  % of all individuals  Enchytraeidae  694  29  0.12%  Haplotaxida  Lumbricidae  118  4  0.02%  Naidid  Haplotaxida  Naididae  1,656  11  0.29%  Tubifi  Haplotaxida  Tubificidae  374  21  0.07%  Lumbri  Lumbriculida  Lumbriculidae  366  28  0.06%  Isotom  Collembola  Isotomidae  104  8  0.02%  Collembola  Poduridae  8  2  <0.01%  Collembola  Sminthuridae  49  9  0.01%  Coleoptera  Amphizoidae  4  2  <0.01%  Coleoptera  Curculionidae  2  3  <0.01%  Coleoptera  Dytiscidae  Coleoptera  Elmidae  Coleoptera  Sminth  Elmida  84  3  0.01%  10,932  36  1.94%  Staphylinidae  21  4  <0.01%  Diptera  Athericidae  1  2  <0.01%  Diptera  Blephariceridae  29  6  0.01%  Cerato  Diptera  Ceratopogonidae  1,602  66  0.28%  Chiron  Diptera  Chironomidae  104,756  143  18.60%  Diptera  Deuterophlebiidae  36  5  0.01%  Diptera  Dixidae  57  4  0.01%  Empidi  Diptera  Empididae  3,163  108  0.56%  Ephydr  Diptera  Ephydridae  180  10  0.03%  Diptera  Oreoleptidae  2  2  <0.01%  Diptera  Pelecorhynchidae  4  2  <0.01%  Diptera  Phoridae  12  3  <0.01%  Psycho  Diptera  Psychodidae  5,541  34  0.98%  Simuli  Diptera  Simuliidae  4,885  81  0.87%  Diptera  Stratiomyidae  9  3  <0.01%  Tipuli  Diptera  Tipulidae  3,146  111  0.56%  Amelet  Ephemeroptera  Ameletidae  2,493  85  0.44%  Baetid  Ephemeroptera  Baetidae  78,392  144  13.92%  Epheme  Ephemeroptera  Ephemerellidae  37,673  141  6.69%  Heptag  Ephemeroptera  Heptageniidae  94,568  143  16.79%  Leptop  Ephemeroptera  Leptophlebiidae  19,788  84  3.51%  Megaloptera  Sialidae  4  2  <0.01%  82  Order (unless marked by *for class or **for suborder)  Family  Odonata  Gomphidae  Capnii  Plecoptera  Chloro  Total # of Individuals (all samples)  Occurrence (# of sites)  % of all individuals  3  3  <0.01%  Capniidae  12,349  92  2.19%  Plecoptera  Chloroperlidae  20,305  140  3.60%  Leuctr  Plecoptera  Leuctridae  4,566  82  0.81%  Nemour  Plecoptera  Nemouridae  54,876  144  9.74%  Plecoptera  Peltoperlidae  128  2  0.02%  Perlid  Plecoptera  Perlidae  2,000  34  0.36%  Perlod  Plecoptera  Perlodidae  4,009  124  0.71%  Plecoptera  Pteronarcyidae  57  6  0.01%  Taenio  Plecoptera  Taeniopterygidae  39,266  94  6.97%  Apatan  Trichoptera  Apataniidae  702  16  0.12%  Brachy  Trichoptera  Brachycentridae  5,509  42  0.98%  Glosso  Trichoptera  Glossosomatidae  6,123  85  1.09%  Hydro  Trichoptera  Hydropsychidae  8,028  114  1.43%  Trichoptera  Hydroptilidae  179  5  0.03%  Trichoptera  Lepidostomatidae  750  14  0.13%  Trichoptera  Leptoceridae  16  2  <0.01%  Limnep  Trichoptera  Limnephilidae  993  52  0.18%  Philop  Trichoptera  Philopotamidae  90  9  0.02%  Taxon Code  Lepido  Trichoptera  Polycentropodidae  8  2  <0.01%  Rhyaco  Trichoptera  Rhyacophilidae  14,259  140  2.53%  Ueniod  Trichoptera  Uenoidae  5,212  39  0.93%  Amphipoda  Crangonyctidae  6  3  <0.01%  Amphipoda  Gammaridae  113  3  0.02%  **Prostigmata  Halacaridae  4  2  <0.01%  **Oribatida  Hydrozetidae  339  20  0.06%  **Oribatida  Oribatidae  55  4  0.01%  Hydryp  **Prostigmata  Hydryphantidae  671  21  0.12%  Hygrob  **Prostigmata  Hygrobatidae  567  22  0.10%  Lebert  **Prostigmata  Lebertiidae  1,382  30  0.25%  Sperch  **Prostigmata  Sperchonidae  2,113  61  0.38%  **Prostigmata  Stygothrombidiidae  **Prostigmata  Torrenticolidae  *Hydrozoa  Hydridae  Veneroida  Sphaeriidae (Pisidiidae)  Basommatophora Basommatophora  Hydroz  Torren Sphaer Planor Planar  32  3  0.01%  3,695  38  0.66%  11  3  <0.01%  1,747  11  0.31%  Physidae  3  3  <0.01%  Planorbidae  55  8  0.01%  Heterostropha  Valvatidae  0  1  <0.01%  Tricladida  Planariidae  2,302  60  0.41%  83  Appendix D Component loadings on significant factors, Principal Components Analyses LANDSCAPE LEVEL Components  1  2  3  4  ECOZONE  -0.691  -0.003  0.562  0.236  % WETLANDS  -0.663  0.500  -0.211  0.129  % LAKES  -0.199  0.559  -0.482  0.183  ICE PA  0.490  0.285  0.318  0.140  % ALPINE  0.713  0.076  0.521  0.119  -0.659  -0.026  0.533  0.399  RELIEF  0.822  0.242  0.155  -0.105  DRN_DENS  0.564  -0.311  -0.015  0.265  0.879  ELEVATION  %HILLSLOPE >60%  0.012  0.183  -0.238  -0.045  0.507  0.244  -0.075  HEAT_INDEX  0.219  0.324  0.089  0.082  LATITUDE  0.353  -0.025  -0.183  0.847  -0.776  -0.004  0.336  -0.301  AREA  0.181  0.844  0.071  -0.074  Eigenvalue  4.716  1.881  1.525  1.262  % Total Variance Explained  33.7  13.4  10.9  9.0  Cumulative % Variance Explained  33.7  47.1  58.0  67.0  1  2  3  4  5  GRADIENT %  0.434  0.520  0.210  -0.424  0.208  % POOLS  0.287  -0.175  -0.682  -0.157  0.018  % CASCADES  0.013  0.736  0.207  -0.351  -0.180  MC  0.165  -0.307  0.643  0.124  0.183  -0.436  -0.181  0.183  -0.031  0.454  CC  0.582  -0.071  0.119  0.204  0.313  SSCOV10  0.051  0.497  -0.061  0.551  0.387  RIPVEG CON  0.324  0.460  -0.248  0.455  -0.010  -0.848  0.311  0.037  -0.005  -0.014  0.373  -0.026  0.318  0.326  -0.661  -0.797  0.091  -0.003  0.277  -0.139  BARREN PA  LONGITUDE  LOCAL LEVEL Components  PERI 35  BWAVE CHANRATIO STREAM ORDER Eigenvalue  2.426  1.538  1.182  1.079  1.019  % Total Variance Explained  22.1  14.0  10.7  9.8  9.3  Cumulative % Variance Explained  22.1  36.0  46.8  56.6  65.9  84  GEOLOGY Components  1  2  3  -0.947  -0.218  0.078  % INTRU  0.242  0.726  -0.284  % VOLCA  0.846  -0.172  0.382  %METAM  0.132  -0.349  -0.736  % ULTRA  0.044  0.108  -0.608  DOM DRAIN  -0.501  0.560  0.126  DOM LOCSF  0.147  0.589  0.046  Eigenvalue  1.964  1.397  1.162  % Total Variance Explained  28.1  20.0  16.6  Cumulative % Variance Explained  28.1  48.1  64.7  % SEDI  85  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0071507/manifest

Comment

Related Items