Near-Term Effects of Repeated-Thinning with Riparian Buffers on Headwater Stream Vertebrates and Habitats in Oregon, USA

Abstract

We examined the effects of a second-thinning harvest with alternative riparian buffer management approaches on headwater stream habitats and associated vertebrates in western Oregon, USA. Our analyses showed that stream reaches were generally distinguished primarily by average width and depth, along with the percentage of the dry reach length, and secondarily, by the volume of down wood. In the first year post-harvest, we observed no effects of buffer treatment on stream habitat attributes after moderate levels of thinning. One of two “thin-through” riparian treatments showed stronger trends for enlarged stream channels, likely due to harvest disturbances. The effects of buffer treatments on salamanders varied among species and with habitat structure. Densities of Plethodon dunni and Rhyacotriton species increased post-harvest in the moderate-density thinning with no-entry buffers in wider streams with more pools and narrower streams with more down wood, respectively. However, Rhyacotriton densities decreased along streams with the narrowest buffer, 6 m, and P. dunni and Dicamptodon tenebrosus densities decreased in thin-through buffers. Our study supports the use of a 15-m or wider buffer to retain sensitive headwater stream amphibians.

1. Introduction

Forest management practices are a leading concern for world biodiversity. The International Union for the Conservation of Nature (IUCN) has been tracking species conservation for 50 years, with their Red List established as the global authority for species status assessments. Amphibians and fishes are two highly vulnerable chordate taxa (32% of amphibian species are globally threatened, and 43% are experiencing declines [1]; 41% of amphibian species are threatened by other criteria [2]; >40% of freshwater fish species are threatened [3]). In particular, both aquatic-dependent taxa are of great concern in forest ecosystems, with many traditional forestry practices, such as clearcut logging without riparian protective areas, listed as major threats to freshwater fishes at the global scale [4]. For amphibians, habitat loss and degradation is the leading threat [5], with 68% of world amphibians associated with forests: ~5000 species occur in forests [6] of the 7307 species known today [7]. Furthermore, amphibians and fishes are considered “biodiversity indicators” of ecosystem health for the larger aquatic and riparian communities with which they are associated [8]. The societal value of amphibians is growing in this regard, as they have been considered “canaries in the coal mine” or environmental sentinels, signaling system health for other species or ecosystem services (e.g., in North American forest ecosystems [9,10,11,12]). Additionally, forest amphibians have recently been shown to play an important role in forest carbon cycling, facilitating greater carbon storage by reducing leaf litter predators [13].

There are numerous challenges to managing forests for both commodity production and ecological integrity and ecosystem services, including species and their habitats [14,15,16,17,18,19,20]. Retention of ecological services provided by freshwater aquatic systems in forests epitomizes these challenges, because aquatic-dependent species, functions and processes can be directly and indirectly affected by traditional timber-harvest practices (e.g., [21,22,23,24,25]). Concerns for forest aquatic resources include maintenance or restoration of water quality (temperature [26,27,28,29,30]; sedimentation [31,32,33,34]) and suitable habitats for biodiversity (e.g., leaf litter inputs [35,36,37]; down wood [38,39,40,41]).

Aquatic-riparian management zones and, specifically, no-entry buffers around streams and lakes are mitigations applied in diverse contexts worldwide, but especially in North America, termed as “best management practices” to address water quality metrics and aquatic biodiversity conservation [42]. In the Pacific Northwest region of North America, where forests occur on relatively steep topographies with small headwater streams dominating the aquatic system [43], streamside riparian reserves are applied to retain habitat conditions for a diverse set of species and habitat values [14,23,39,44,45,46,47]. New science emerging from empirical studies of amphibians and fishes and their habitat elements is leading the charge to understand riparian management efficacy. In particular, amphibian studies are greatly advancing our knowledge of the effects of riparian buffer management, with: (1) forest stand-scale experiments addressing the effects of alternative buffer widths with alternative upslope timber harvest activities on resident fauna and their habitats; (2) populations being monitored sometimes for many years following treatments; and (3) the spatial setting of studies spanning regions and different forest contexts (review: [47]; also [48,49,50,51,52,53]). Although the inference of each study is often narrow, the extended value of this growing combined dataset is being realized, especially as the role of amphibians as potential environmental sentinels is considered. Hence, findings on how amphibians respond to forested buffers with upland timber harvest are useful to consider as other aquatic-dependent ecosystem services come under scrutiny. Several questions remain, however, including the long-term effects of riparian buffers of different widths, with various upland management regimes, on stream-associated vertebrates and their critical habitat attributes.

The Density Management and Riparian Buffer Study in western Oregon, USA (DMS) has contributed to the growing knowledge base for amphibian responses to riparian buffers in managed forests by tracking several salamander species, in addition to frogs and fish. Elements of this before-after-control study examined the responses of aquatic vertebrates to upland thinning with alternative no-entry streamside buffer widths, including: (1) instream and stream bank salamander counts pre-treatment [54] and after Year 1, 2, 5 and 10 post-thinning [53,55]; and (2) upland salamanders after Years 1 and 2, and 5 and 6 post-thinning [55,56,57,58].

At DMS study sites, a second thinning has now occurred, providing the opportunity to continue to track the efficacy of alternative riparian buffer widths relative to multiple upland forest management entries and tracking the potential effects on aquatic vertebrates and their habitat conditions. DMS has had dual goals of advancing the knowledge of both riparian management alternatives and upland forest restoration approaches, forcing an intersection of land and water management [59]. The original intent of the upland forest density management portion of DMS was to examine ways to accelerate the development of late-successional forest conditions by thinning young managed stands [60], including the development of stand heterogeneity and multi-layered canopies [61]. To achieve these objectives, the DMS study plan called for two forest thinning phases, an initial thinning and a second-entry re-thinning [60]. This experimental design stemmed from prior analyses of tree-growth trajectories and early-successional tree densities in regional old-forest stands. When those old-forest stands were young, sustained tree growth beyond age 50 years in stands with low live tree densities may have been important to attain large tree diameters in stands aged 100–300 years [62,63].

Herein, we report the analyses of the effects of three alternative streamside riparian buffers with a second upland thinning harvest on headwater stream habitat attributes, animal-habitat associations and animal densities. Eight study sites were re-thinned with this second-entry harvest occurring ~12 years after the initial thinning. We analyzed habitat differences and animal-habitat associations pre- and post-second thinning to examine if changes to stream habitat attributes, such as the pool:riffle ratio, stream width, depth, substrate type and down wood volume, account for potential shifts in animal density patterns. We were particularly interested in whether the animal count patterns by species and species groups from Year 10 after the first thinning ([53]; i.e., decreases in stream-bank salamanders and Plethodon dunni, in particular) were retained after the second thinning or whether a new signature emerged from the effects of the different buffers on stream-dependent vertebrate densities. In addition, we examined the potential effects of a new thinned riparian buffer treatment. These “thin-through” buffers may provide insights for forest restoration approaches across the entire riparian zone, which may have long-term ecological benefits, such as accelerating the development of large, stream-bank conifers and future large down wood recruitment into streams.

2. Methods

2.1. Study Area

Our experimental riparian buffer study is replicated at eight sites along the Coast and Cascade Ranges in western Oregon, USA (Figure 1). Located in the western hemlock (Tsuga heterophylla) zone, the climate is characterized as Mediterranean with mild, wet winters and warm, dry summers [64]. Soils are characterized as well- to poorly-drained Ultisols and Inceptisols.

Figure 1. Locations of eight Density Management and Riparian Buffer Study sites in western Oregon, USA, where we examined the effects of a second-thinning harvest with different riparian buffer management designs on instream and stream bank aquatic-dependent vertebrates.

Figure 1. Locations of eight Density Management and Riparian Buffer Study sites in western Oregon, USA, where we examined the effects of a second-thinning harvest with different riparian buffer management designs on instream and stream bank aquatic-dependent vertebrates.

At the start of the study, forests were dominated by second-growth 30–70-year-old Douglas-fir (Pseudotsuga menziesii) trees with varying abundances of western hemlock. Other conifer species, such as western redcedar (Thuja plicata), and hardwood species, including bigleaf maple (Acer macrophyllum), red alder (Alnus rubra), Pacific dogwood (Cornus nuttallii), Pacific madrone (Arbutus menziesii) and golden chinquapin (Chrysolepis chrysophylla), were minor components of the overstory. Forest understory communities were characterized by dominant species, including sword fern (Polystichum munitum) and salal (Gaultheria shallon), that vary in relative abundance in association with Oregon Oxalis (Oxalis oregonum) at mesic sites and ocean spray (Holodiscus discolor), common whipplea (Whipplea modesta) and white hawkweed (Hieracium albiflorum) at the driest sites [64,65,66].

Stands regenerated naturally following clearcut and seed tree harvests lacking riparian buffers, with the exception of portions of two sites (Keel Mountain and OM Hubbard), in which Douglas-fir seedlings were planted [60]. Slash management following harvest operations varied among sites with some evidence of burning. There is evidence of historical grazing prior to canopy closure at one site (Green Peak), and three sites were precommercially thinned (study site descriptions: [60]).

2.2. Experimental Design

In conjunction with two consecutive upland forest thinning treatments, three no-harvest riparian buffer treatments included a one site-potential tree-height buffer (“1Tree”, ~70 m on each side of streams), a variable-width buffer (“VAR”, 15-m minimum width on each side of streams) and streamside retention buffer (“STRM”, 6 m on each side of streams; see [53,55]). The second thinning occurred 12 years after the first thinning. With the exception of Perkins Creek, an older stand, which had been thinned 20 years earlier to 250 trees/ha (tph), the first thinning treatment in the moderate-density (MD) upland areas, implemented between 1997 and 2000, retained ~200 tph from a pre-treatment density of ~430–600 tph [60]. At that time, Perkins Creek was thinned to 100–150 tph in its second entry. The second thinning (third for Perkins Creek) occurred between 2009 and 2011 and further reduced overstory tree densities to ~85 tph. At each site, all riparian buffer treatments were paired with an untreated control (i.e., no upland thinning; Table 1).

Table 1. The number of replicate stream reaches with alternative riparian buffer treatments (1Tree (one site-potential tree-height buffer), variable, streamside, thin-through) and untreated controls at eight study sites with upland forests thinned twice to moderate residual tree densities and one site thinned twice to high residual tree densities, in western Oregon, USA; parentheses indicate that these control reaches are represented above in the controls for the moderate-density analysis.

Site	Control	1Tree	Variable	Streamside	Thin-Through
Moderate-density upland retention
Callahan Creek	2	1	2	2	2
Delph Creek	2	0	1	1	0
Green Peak	2	1	1	1	0
O.M. Hubbard Creek	3	2	2	0	0
Keel Mountain	2	1	1	1	3
Perkins Creek	2	0	2	2	0
N. Soup Creek	2	1	1	1	0
Ten High	2	1	3	2	0
Total	17	7	13	10	5
High-density upland retention
Keel Mountain	(2)	0	3	0	4

Table 1. The number of replicate stream reaches with alternative riparian buffer treatments (1Tree (one site-potential tree-height buffer), variable, streamside, thin-through) and untreated controls at eight study sites with upland forests thinned twice to moderate residual tree densities and one site thinned twice to high residual tree densities, in western Oregon, USA; parentheses indicate that these control reaches are represented above in the controls for the moderate-density analysis.

Site	Control	1Tree	Variable	Streamside	Thin-Through
Moderate-density upland retention
Callahan Creek	2	1	2	2	2
Delph Creek	2	0	1	1	0
Green Peak	2	1	1	1	0
O.M. Hubbard Creek	3	2	2	0	0
Keel Mountain	2	1	1	1	3
Perkins Creek	2	0	2	2	0
N. Soup Creek	2	1	1	1	0
Ten High	2	1	3	2	0
Total	17	7	13	10	5
High-density upland retention
Keel Mountain	(2)	0	3	0	4

During the second thinning of MD overstory, density in upland areas was reduced to a new residual density of ~85 tph. An additional “thin-through” (TT) riparian buffer treatment (i.e., without a no-entry riparian buffer zone) was added to five stream reaches at two sites. These five reaches converted previous no-entry riparian buffer widths of two site-potential tree heights (“2Tree” riparian buffers; ~145 m on each side of streams: Callahan Creek and Keel Mountain study sites; Figure 1) to a thinned riparian buffer. Within the 2Tree buffer width; the overstory density was reduced from ~430–600 tph to ~150 tph during the second thinning treatment. We chose a 150-tph riparian thinning density based on lessons-learned from our initial upland density of 200 tph, where the restoration benefits of opening the canopy appeared to be subsiding within about a decade of the thinning due to canopy closure occurring. During the implementation of the thin-through buffer treatment, timber harvesting operators were asked to yard trees away-from or suspended-over streams; avoiding ground disturbance to the stream prism. For analyses, these TT riparian treatments were paired with the other riparian buffer treatments within the MD upland thinning treatment (Table 1).

Finally, to further assess the effects of the interactions between residual upland tree density and riparian buffer width, we replicated the variable-width and thin-through buffer treatments within a second upland thinning treatment of a higher residual overstory tree density (high-density upland treatment: HD). During the first thinning, this upland treatment area was thinned from 430–600 tph to 300 tph, and during the second thinning, it was thinned to 150 tph; hence, residual densities are higher than in our moderate-density (MD) treatment described earlier. At the outset of our project in 1994, this HD treatment was similar to a typical commercial thinning prescription in our region. The HD treatment can thus serve as an additional reference for the contrasting effects of MD treatments aimed at accelerating the development of old-growth forest conditions with HD treatments aimed at wood production [67]. While HD buffer treatments are not as well replicated among sites due to poor stream representation, the Keel Mountain study site contained numerous streams, and 7 stream reaches were treated with variable-width buffers in the HD unit during the first thinning. During the second thinning, the variable-width buffer treatment was retained as a no-entry buffer zone along 3 of 7 stream reaches, and a thin-through (TT) was implemented at the remaining 4 reaches. These TT riparian areas in the HD uplands were thinned to a residual overstory tree density of 150 tph from a pre-thinning density of ~430–600 tph, consistent with the overstory tree density achieved in the associated HD upland thinning after the second thinning and the TT treatment in the MD upland thinning. We compared the HD buffer treatments (VAR and TT) to the same unthinned control reaches at Keel Mountain that were used for comparison with buffer treatments in the MD upland thinning (Table 1).

2.3. Field Methods

Habitat and animal surveys were conducted each year during rainy spring conditions (March–June) along 59 reaches distributed across eight sites (Table 1). Habitat surveys were conducted along entire reach lengths using a modified Hankin and Reeves [68] approach and preceded animal sampling [54,55]. Within reaches, streams were classified into units by surface stream water flow category: pool, riffle, cascade, dry or other. Unit dimensions were measured (i.e., length, width, average depth). Reach habitats were characterized by calculating average dimensions (e.g., width, depth). Reaches were classified into 5 “hydrotypes”: (1) perennial reaches with surface flow in both spring and summer; (2) summer intermittent, with continuous flow in spring and discontinuous flow in summer; (3) perennial-ephemeral, having continuous flow in spring, but being dry in the summer; (4) intermittent, with discontinuous flows in both spring and summer; and (5) intermittent-ephemeral, with discontinuous flow in the spring and dry in the summer. Drainage basin areas were estimated by delineating basins using coordinates from Global Positioning Systems and topography maps in ArcMAP Geographic Information System software v.10 [69].

Down wood was sampled primarily during the fall season (October–November) due to spring-season time constraints. The diameter and length of down wood pieces ≥10 cm in diameter were estimated visually, and each piece was categorized into one of five decay classes [70]. For a subsample of pieces during each survey (i.e., at each reach and year), width and length were measured accounting for a total of 18% of the total number of pieces recorded. To correct for potential error in estimates and differences among observers, we adjusted estimates using no-intercept linear relationships between measured and estimated values (R² = 0.94 and 0.95 for length and diameter, respectively). The percentage of down wood volume present within three stream zones was recorded: Zone 1, the current wetted channel; Zone 2, outside the wetted channel and within the estimated channel at bankfull stages of water flow; and Zone 3, down wood suspended above the stream prism at the estimated bankfull stage and outside Zones 1 and 2. For our analyses, we calculated total volume and the volume of stream wood in early (Classes 1–2) and late (Classes 3–5) stages of decay, assuming a cylindrical geometry and summed across Zones 1–3 (i.e., within and above the bankfull stage water level).

After habitat surveys were completed, stream-associated vertebrate species were surveyed within streams and along stream banks within a subsample of 10 units per reach [53,54,55], with the exception of two reaches at the Perkins Creek site, where >10 units were sampled due to extended reach lengths, because a larger sample was considered potentially necessary to capture patchy animal occurrences. Reach segments ~15 m from the bottom and top-most portions of the reach were excluded from sampling. Dry units and units with obstructed access (due to vegetation or down wood) were not sampled. Selected units reflected the relative proportions of slow- and fast-water flow classes (i.e., pools, riffles/cascades) in the reach, determined from the prior habitat surveys. The first unit was selected randomly from the first 5 units included at the downstream end of the reach. Moving upstream, subsequent units were selected such that every n-th unit of a type was sampled, where n = the number of units by type (i.e., pool or riffle) in the remaining upstream reach to be sampled, divided by the remaining units [68]. Stream and stream bank sampling were conducted separately at each unit.

For stream animal sampling, units were either sampled by hand (43 reaches measured twice) or by electrofishing (12 reaches measured twice). Only one reach was sampled once by electrofishing and once by hand [53]. Electrofishing was restricted to fish-bearing streams. For hand sampling, animals were captured using dip nets placed downstream from cover objects and herded into nets by hand as objects were moved and replaced in a “soft-touch” approach. For both types of surveys, samplers moved upstream, placing captures in buckets as they proceeded through the unit. Once captures were identified to species, they were returned to the unit. Unit length and width was determined and animal densities per square meter calculated.

Stream bank surveys were conducted for 5 min on each side of the sampled unit within 2 m of the water’s edge. Surveyors searched in, on and under substrates, moving up- and down-stream to use the full 5 min effectively. Animals were placed in plastic bags upon capture, identified and released to the capture location. Handling times were subtracted from the 5-minute search times. Similar to stream surveys, search areas were estimated and species densities were calculated on a per-square-meter basis.

2.4. Statistical Analysis

Pre-treatment data from 10 years after the first thinning and post-treatment data one year after the second thinning were used to characterize the effects of riparian buffer treatments associated with upland forest thinning on stream-associated vertebrates and key instream habitat structures. Because habitat variables may co-vary non-linearly in headwater streams, we identified the main dimensions of habitat variability among streams, riparian buffer and thinning treatments, as well as years, using non-metric multidimensional scaling (NMS) ordinations [71]. The distance matrix was calculated using the Sørenson index of (dis)similarity based on differences in the following quantitative geomorphic habitat characteristics: stream width and depth, average pool depth, average riffle depth, patterns of surface stream flow (i.e., %riffle and %pool by length and area, pool:riffle ratio by length, %dry length), substrate (i.e., percentage fine, small gravel, large gravel, cobble, boulder and bedrock), and total volume of down wood and volumes in early stages (1–2) and late stages (3–5) of decay. The final two-dimension solution was based on the rotation that maximized the variation explained by ordination axes. Axes were interpreted by examining the relationships with the characteristics in the ordination main matrix and secondary exogenous environmental variables (e.g., gradient). Changes in habitat structure associated with the second thinning treatment were evaluated by assessing individual reach and buffer treatment by upland thinning treatment (i.e., MD, HD or untreated control) trajectories and using pairwise comparisons from multi-response permutation tests (MRPP). MRPP tests the null hypothesis that there is no difference between two or more groups by comparing observed differences within groups, with expected differences representing the mean difference for all possible partitions of the data [71]. Finally, associations between individual species and species groups with reach-scale habitat characteristics were assessed by examining the correlations of individual species densities with ordination axes. The NMS ordination and MRPP tests were completed using PC-ORD version 6 [72].

Then, to examine how the second upland thinning treatment interacted with riparian buffers and habitat structure to affect the population densities of instream and stream bank vertebrate species, we used generalized linear mixed model analysis of covariance, including fixed effects of buffer treatment (TT, STRM, VAR, 1Tree or CON), year (pre- or post-second thinning treatment, with the “pre” year being 10 years post-first thinning and the “post” year being one year post-second thinning). Because the effects may vary with the differences in the initial habitat characteristics and stream geomorphology, ordination axis scores, which captured such synthetic habitat variation, were used as a covariate “habitat.” Additionally, the method (hand search vs. electrofishing for stream surveys) was included as a fixed effect [53]. We tested for all main effects: method, habitat, buffer and year; and for two- and three-way interactions among buffer treatment, year and habitat. We examined effects on: (1) all species combined; (2) all amphibian species; (3) all fish species—cutthroat trout (Oncorhynchus clarkii), trout species (Oncorhynchus), sculpin species (Cottidae); (4) stream-breeding amphibian species—coastal tailed frog (Ascaphus truei), coastal giant salamander (Dicamptodon tenebrosus) and torrent salamanders (Rhyacotriton spp.), (5) terrestrial-breeding amphibian species—clouded salamander (Aneides ferreus), Oregon slender salamander (Batrachoseps wrighti), Ensatina (salamander) (Ensatina eschscholtzii), Dunn’s salamander (Plethodon dunni) and western red-backed salamander (P. vehiculum); (6) pond-breeding species—rough-skinned newt (Taricha granulosa), northwestern salamander (Ambystoma gracile) and northern red-legged frog (Rana aurora); and (7) dominant individual species—D. tenebrosus, Rhyacotriton spp., P. dunni and P. vehiculum. To account for the effects of differences in fluvial geomorphology among reaches on population densities, we standardized animal counts to a per-square-meter basis for analyses by dividing counts by stream and stream bank unit area. We included random effects for site and stream (nested in site).

In this model, the effects of the second thinning are supported by evidence for a buffer × year interaction, whereas the main effects of buffer treatment support differences among riparian buffer treatments that are independent of the second thinning (e.g., due to effects of the first thinning or pre-existing variation among treatments). The main effects of year (with no interaction) suggest temporal variation that is not related to treatments. We interpret the evidence as follows: p < 0.01 = strong evidence; p < 0.05 = moderate evidence; and p < 0.15 = weak evidence. Due to our small sample size, as well as inherent difficulties in sampling aquatic vertebrate species, we interpret trends with weak evidence. To account for the zeros, we added a constant (half the smallest non-zero observation) to all density values prior to analyses. We used generalized linear mixed-models assuming a log-normal distribution and identity link function, using the GLIMMIX procedure in SAS version 9.3 [73]. We were constrained by the use of emerging analysis techniques of occupancy and detection probability; due to limited funding, we opted to retain all study sites and reaches, but were only able to conduct single-pass animal surveys of each stream reach. However, we previously reported [53] the results of a case study of occupancy and detection probabilities of three amphibian species from our survey design: instream Dicamptodon tenebrosus and Rhyacotriton spp. occupancies were 0.83 and 0.9, respectively, and detection probabilities ranged 0.35 to 0.53 for Dicamptodon and was 0.65 for Rhyacotriton; whereas bank-dwelling Plethodon dunni had occupancy ranging from 0.45 to 0.85 and detection probabilities ranging from 0.3 to 0.45.

3. Results and Discussion

3.1. NMS Ordination

In analyses of instream habitat attributes and animal densities, the final NMS ordination solution had a final stress of 12.4, with a low proportion of randomized runs, with stress equal to or lower than this value (minimized stress in randomized data = 28.9, p = 0.020). The two-dimensional ordination explained 94% of the variability on the distance matrix (84% on Axis 1, 10% on Axis 2). Axis 1 contrasted perennial stream reaches with greater average depth (r = 0.89) in perennial Hydrotypes 1 and 2 against intermittent stream reaches with a greater %dry length (r = −0.84; Table 2, Figure 2). Axis 2 contrasted stream reaches with a greater percentage of length and area in pools (r = 0.73 and 0.79, respectively) against reaches dominated by riffles (r = −0.53 and −0.74, for length and area, respectively) with greater volumes of down wood (r = −0.40) and, particularly, down wood in early stages of decay (r = −0.54).

Figure 2. Non-metric multi-dimensional scaling ordination of 59 stream reaches described by 24 habitat characteristics; arrows show habitat and animal density patterns with axes, in positive and negative directions of abundance; species acronyms: CT, cutthroat trout, Oncorhynchus clarkii; DITEa, aquatic Dicamptodon tenebrosus, found in streams; DITEt, terrestrial D. tenebrosus, found on stream banks; ALLSP, all species; STRM, all instream species; AMPH, all amphibians; RHsp., Rhyacotriton spp.; ENES, Ensatina eschscholtzii.

Figure 2. Non-metric multi-dimensional scaling ordination of 59 stream reaches described by 24 habitat characteristics; arrows show habitat and animal density patterns with axes, in positive and negative directions of abundance; species acronyms: CT, cutthroat trout, Oncorhynchus clarkii; DITEa, aquatic Dicamptodon tenebrosus, found in streams; DITEt, terrestrial D. tenebrosus, found on stream banks; ALLSP, all species; STRM, all instream species; AMPH, all amphibians; RHsp., Rhyacotriton spp.; ENES, Ensatina eschscholtzii.

Table 2. Correlations of variables in the ordination of reach characteristics (i.e., main matrix) with non-metric multidimensional scaling ordination axes (Figure 2).

Reach Characteristics	Correlation with Ordination Axes
	Axis 1	Axis 2
Dimensions of sampled habitat units
Average width (m)	0.59 ***	0.04
Average depth (m)	0.89 ***	0.50 ***
Pool depth (m)	0.88 ***	0.46 ***
Riffle depth (m)	0.81 ***	0.32 ***
Percentage of sampled habitat units
%Riffle length	−0.33 ***	−0.13
%Riffle area	−0.33 ***	−0.28 ***
%Pool length	0.64 ***	0.32 ***
%Pool area	0.61 ***	0.39 ***
Percentage cover of substrates (reach scale)
%Fine	0.01	−0.16 *
%Small gravel	−0.52 ***	−0.05
%Large gravel	0.56 ***	0.22 **
%Cobble	0.63 ***	0.11
%Boulders	0.69 ***	0.02
%Bedrock	0.58 ***	−0.02
Reach-scale characteristics
%Riffle length	0.42 ***	−0.53 ***
%Riffle area	−0.58 ***	−0.74 ***
%Pool length	0.76 ***	0.73 ***
%Pool area	0.70 ***	0.79 ***
%Dry length	−0.84 ***	0.04
Pool:riffle (length)	0.52 ***	0.75 ***
Average length of dry units	−0.80 ***	0.05
Stream wood (reach-scale)
Volume in early stages of decay (Classes 1, 2)	−0.51 ***	−0.54 ***
Volume in late stages of decay (Classes 3, 4, 5)	−0.36 ***	−0.30 ***
Total volume of wood	−0.43 ***	−0.40 ***

Asterisks indicate the strength of evidence for non-zero correlation; * p < 0.05; ** p < 0.01; *** p < 0.001.

Table 2. Correlations of variables in the ordination of reach characteristics (i.e., main matrix) with non-metric multidimensional scaling ordination axes (Figure 2).

Reach Characteristics	Correlation with Ordination Axes
	Axis 1	Axis 2
Dimensions of sampled habitat units
Average width (m)	0.59 ***	0.04
Average depth (m)	0.89 ***	0.50 ***
Pool depth (m)	0.88 ***	0.46 ***
Riffle depth (m)	0.81 ***	0.32 ***
Percentage of sampled habitat units
%Riffle length	−0.33 ***	−0.13
%Riffle area	−0.33 ***	−0.28 ***
%Pool length	0.64 ***	0.32 ***
%Pool area	0.61 ***	0.39 ***
Percentage cover of substrates (reach scale)
%Fine	0.01	−0.16 *
%Small gravel	−0.52 ***	−0.05
%Large gravel	0.56 ***	0.22 **
%Cobble	0.63 ***	0.11
%Boulders	0.69 ***	0.02
%Bedrock	0.58 ***	−0.02
Reach-scale characteristics
%Riffle length	0.42 ***	−0.53 ***
%Riffle area	−0.58 ***	−0.74 ***
%Pool length	0.76 ***	0.73 ***
%Pool area	0.70 ***	0.79 ***
%Dry length	−0.84 ***	0.04
Pool:riffle (length)	0.52 ***	0.75 ***
Average length of dry units	−0.80 ***	0.05
Stream wood (reach-scale)
Volume in early stages of decay (Classes 1, 2)	−0.51 ***	−0.54 ***
Volume in late stages of decay (Classes 3, 4, 5)	−0.36 ***	−0.30 ***
Total volume of wood	−0.43 ***	−0.40 ***

Asterisks indicate the strength of evidence for non-zero correlation; * p < 0.05; ** p < 0.01; *** p < 0.001.

Species correlations with ordination axes (Table 3) revealed that aquatic D. tenebrosus (r = 0.42, p < 0.001), O. clarkii (r = 0.40, p < 0.001) and A. truei (r = 0.31, p < 0.001) were most strongly associated with larger perennial reaches (Hydrotypes 1 and 2) and positively related to Axis 1 (Figure 2). Rhyacotriton spp. (r = −0.34, p < 0.001) and the terrestrial-occurring D. tenebrosus (r = −0.22, p < 0.01), however, were negatively related to Axis 1 and associated with the intermittent stream reaches (Hydrotypes 3, 4 and 5). Rhyacotriton spp. (r = −0.35, p < 0.001) and E. eschscholtzii (r = −0.29, p < 0.001) were negatively related to Axis 2 and associated with reaches characterized by a greater proportion of riffles and volume of down wood. These findings are highly consistent with the relationships detected in the mid-1990s at 13 DMS sites analyzed before the first thinning treatment was implemented [54]. The terrestrial Dicamptodon association with intermittent streams is a new finding at DMS sites, supporting the contention that terrestrial adults migrate upstream, since juvenile aquatic forms are associated with larger perennial reaches downstream. The E. eschscholtzii relationship with down wood is not surprising, as this species is known for its affinity for down wood at the DMS sites [56,57,74].

With the exception of the thin-through reaches, successional vectors of individual reaches and changes in centroid locations suggested that the second thinning did not lead to dramatic differences among buffer treatments in habitat structure (Figure 3). Visual patterns were confirmed by the absence of significant differences in habitat structure between pre- and post-treatment years in the MRPP analysis (all p > 0.1). However, based on visible differences in centroid locations, there appears to be an overall tendency among all four riparian buffer treatments in the MD uplands and untreated controls for reaches to shift slightly positively along NMS Axis 1 and negatively along NMS Axis 2 (Figure 3). Larger positive shifts along Axis 1 for the MD thin-through reaches suggest that a thin-through buffer strategy without a no-entry zone along streams leads to more considerable broadening and deepening of headwater reaches.

This trend for a broadening of stream reaches in the thin-through treatment with moderate upland residual tree densities (MD) was not apparent in the thin-through treatment with high residual tree densities (HD), however. Although it is possible that the lower upland retention in the MD thinning treatment may be contributing to this shift, it was noted by our survey crews that physical damage to some MD streams resulted from implementation of thin-through riparian buffers; it appeared that logs were dragged through some MD streams and that harvest machinery entered some MD stream prisms during operations [75]. It is likely that our analyses showing enlarged stream habitats with thin-through buffers in MD uplands, but not HD uplands, could simply be due to differences in very local timber harvest activities. In summary, there appears to be no evidence for a treatment effect of the second thinning on the habitat structures we measured, although patterns of down wood recruitment and stream flow in reaches with no-entry buffers and the pattern for enlarged channels in reaches with thin-through buffers all warrant further consideration.

Table 3. Correlations of species densities with axes from the ordination of reaches described by habitat characteristics.

Species/Group	Correlation with Ordination Axes
	Axis 1	Axis 2
All species combined	0.40 ***	0.15 *
Fish species	0.29 ***	0.18 *
Oncorhynchus spp.	0.22 **	0.21 **
O. clarkii	0.40 ***	0.23 **
Cottidae spp.	0.25 ***	0.15 *
Amphibian species	0.33 ***	0.04
Stream-breeding amphibians	0.32 ***	0.04
A. truei	0.31 ***	0.14
D. tenebrosus (aquatic)	0.42 ***	0.16 *
D. tenebrosus (terrestrial)	−0.26 ***	−0.08
Rhyacotriton spp.	−0.32 ***	−0.35 ***
Terrestrial-breeding amphibians	0.17 *	−0.13
A. ferreus	0.03	0.01
B. wrighti	−0.19 *	−0.15*
E. eschscholtzii	−0.16 *	−0.29 ***
P. dunni	0.22 **	−0.08
P. vehiculum	−0.08	−0.07
Pond-breeding amphibians	0.15 *	0.17 *
A. gracile	0.09	0.08
T. granulosa	0.11	0.13
R. aurora	0.13	0.12

Asterisks indicate the strength of evidence for non-zero correlation; * p < 0.05; ** p < 0.01; *** p < 0.001.

Table 3. Correlations of species densities with axes from the ordination of reaches described by habitat characteristics.

Species/Group	Correlation with Ordination Axes
	Axis 1	Axis 2
All species combined	0.40 ***	0.15 *
Fish species	0.29 ***	0.18 *
Oncorhynchus spp.	0.22 **	0.21 **
O. clarkii	0.40 ***	0.23 **
Cottidae spp.	0.25 ***	0.15 *
Amphibian species	0.33 ***	0.04
Stream-breeding amphibians	0.32 ***	0.04
A. truei	0.31 ***	0.14
D. tenebrosus (aquatic)	0.42 ***	0.16 *
D. tenebrosus (terrestrial)	−0.26 ***	−0.08
Rhyacotriton spp.	−0.32 ***	−0.35 ***
Terrestrial-breeding amphibians	0.17 *	−0.13
A. ferreus	0.03	0.01
B. wrighti	−0.19 *	−0.15*
E. eschscholtzii	−0.16 *	−0.29 ***
P. dunni	0.22 **	−0.08
P. vehiculum	−0.08	−0.07
Pond-breeding amphibians	0.15 *	0.17 *
A. gracile	0.09	0.08
T. granulosa	0.11	0.13
R. aurora	0.13	0.12

Asterisks indicate the strength of evidence for non-zero correlation; * p < 0.05; ** p < 0.01; *** p < 0.001.

Figure 3. Changes in reach habitat structure by buffer treatment supported by non-metric multi-dimensional scaling ordination of 56 stream reaches described by 24 habitat characteristics; Open and solid symbols show reaches prior to and following second thinning treatments, respectively; centroids show mean locations; results for moderate-density thinning (A) are shown separately from high-density thinning (B).

Figure 3. Changes in reach habitat structure by buffer treatment supported by non-metric multi-dimensional scaling ordination of 56 stream reaches described by 24 habitat characteristics; Open and solid symbols show reaches prior to and following second thinning treatments, respectively; centroids show mean locations; results for moderate-density thinning (A) are shown separately from high-density thinning (B).

3.2. Moderate-Density Thinning

For the MD thinning treatment, results of our analysis of co-variance using generalized linear mixed models suggest that the effects of the second thinning were complex and vary among species and with habitat conditions. The effects of the sampling method, buffer treatments (unrelated to second thinning) and the interaction between buffer and habitat structure (as quantified using NMS Axis 1 scores) were observed for the density of all species combined, amphibian species and stream-breeding amphibians (Table 4). In contrast, densities of fish species were primarily related to sampling method and year, and densities of terrestrial-breeding amphibians as a group were related to habitat structure (NMS Axis 1; Table 4). Trends for groups appeared to be related to patterns in the dominant species.

Results provide moderate evidence (e.g., p < 0.05) for the effects of the second thinning treatment on the density of P. dunni occurring along wider stream reaches with deeper pools and higher Axis 1 scores (Table 4: habitat × buffer × year effect). For P. dunni, the dominant terrestrial-breeding and stream-bank dwelling amphibian species, density increased with NMS Axis 1 in all but the 1Tree buffer treatment, which showed a negative relationship with Axis 1 (Figure 4). In all, but the TT treatment, where density declined, density increased following the second thinning. In the TT treatments in the MD uplands, the reduced P. dunni densities may have resulted from the direct effects of the TT harvest on this species or its habitats. Differences between pre- and post-harvest years were apparent only at high axis scores where population densities were higher. The observed increase in P. dunni densities after the second thinning in the other treatments (1Tree, VAR, STRM and the untreated control) suggests that P. dunni populations are not stable at our DMS sites and perhaps are rebounding from the depressed levels reported previously. For instance, at the 10-year mark following the first thinning, counts had decreased in the narrow streamside buffer treatment [53].

Table 4. Results from mixed model analyses of variance testing for the effects of riparian buffer treatments (1Tree, Variable, Streamside and no-buffer/thin-through vs. untreated controls) prior to and following the second thinning treatment (year), while controlling for sampling method (hand search vs. electrofishing) and habitat structure (NMS Axis 1 scores) on stream-associated vertebrate species; DF_Num = degrees of freedom for the numerator; DF_Den = degrees of freedom for the denominator; Pr = probability.

	Moderate-Density Upland Thinning				High-Density Upland Thinning
Effect	DFNum	DFDen	F-Value	Pr > F	DFNum	DFDen	F-Value	Pr > F
All species
Sampling method	1	38	32.5	<0.0001	1	1	6.3	0.241
Habitat	1	38	0.0	0.876	1	1	2.9	0.336
Buffer	4	38	4.5	0.004	2	1	1.2	0.544
Habitat × Buffer	4	38	3.0	0.030	2	1	3.0	0.379
Year	1	38	0.0	0.947	1	1	2.4	0.367
Buffer × Year	4	38	0.2	0.930	2	1	10.3	0.215
Habitat × Year	1	38	2.2	0.149	1	1	61.7	0.081
Habitat × Buffer × Year	4	38	0.8	0.512	2	1	28.1	0.132
Amphibians
Sampling method	1	38	21.6	<0.0001	1	1	6.3	0.241
Habitat	1	38	0.0	0.921	1	1	2.9	0.336
Buffer	4	38	3.9	0.010	2	1	1.2	0.544
Habitat × Buffer	4	38	3.2	0.024	2	1	3.0	0.379
Year	1	38	0.2	0.666	1	1	2.4	0.367
Buffer × Year	4	38	0.3	0.866	2	1	10.3	0.215
Habitat × Year	1	38	2.6	0.117	1	1	61.7	0.081
Habitat × Buffer × Year	4	38	1.0	0.418	2	1	28.1	0.132
Stream-breeding amphibians
Sampling method	1	38	46.9	<0.0001	1	1	15.9	0.156
Habitat	1	38	0.2	0.645	1	1	0.2	0.753
Buffer	4	38	4.8	0.003	2	1	0.7	0.653
Habitat × Buffer	4	38	2.2	0.091	2	1	0.1	0.907
Year	1	38	0.3	0.608	1	1	0.0	0.953
Buffer × Year	4	38	0.7	0.586	2	1	4.9	0.304
Habitat × Year	1	38	2.5	0.126	1	1	1.5	0.437
Habitat × Buffer × Year	4	38	1.8	0.157	2	1	6.9	0.260
Terrestrial-breeding amphibians
Sampling method	1	38	1.46	0.234	1	1	3.42	0.316
Habitat	1	38	4.32	0.045	1	1	0.2	0.732
Buffer	4	38	1.16	0.343	2	1	0.03	0.967
Habitat × Buffer	4	38	1.42	0.246	2	1	3.47	0.355
Year	1	38	0.04	0.843	1	1	1.32	0.456
Buffer × Year	4	38	0.75	0.566	2	1	0.78	0.626
Habitat × Year	1	38	0.78	0.384	1	1	4.27	0.287
Habitat × Buffer × Year	4	38	1.32	0.280	2	1	1.53	0.497
Pond-breeding amphibians
Sampling method	1	38	0.3	0.615	1	1	4.5	0.280
Habitat	1	38	1.5	0.225	1	1	2.6	0.354
Buffer	4	38	0.3	0.910	2	1	29.4	0.129
Habitat × Buffer	4	38	0.5	0.773	2	1	4.1	0.330
Year	1	38	3.6	0.066	1	1	20.8	0.138
Buffer × Year	4	38	0.8	0.523	2	1	22.5	0.147
Habitat × Year	1	38	0.5	0.503	1	1	11.2	0.185
Habitat × Buffer × Year	4	38	0.7	0.633	2	1	2.3	0.423
Fish
Sampling method	1	38	33.19	<0.0001
Habitat	1	38	1.89	0.1777
Buffer	4	38	1.23	0.315
Habitat × Buffer	4	38	1.58	0.1994
Year	1	38	4.17	0.0482
Buffer × Year	4	38	0.82	0.5195
Habitat × Year	1	38	0.19	0.6681
Habitat × Buffer × Year	4	38	0.48	0.7505
Dicamptodon tenebrosus
Sampling method	1	38	17.0	0.000	1	1	3.4	0.317
Habitat	1	38	10.7	0.002	1	1	31.5	0.112
Buffer	4	38	5.8	0.001	2	1	8.6	0.235
Habitat × Buffer	4	38	2.5	0.061	2	1	4.9	0.305
Year	1	38	0.8	0.370	1	1	9.2	0.203
Buffer × Year	4	38	0.5	0.775	2	1	0.7	0.636
Habitat × Year	1	38	2.7	0.109	1	1	7.7	0.220
Habitat × Buffer × Year	4	38	0.3	0.908	2	1	9.9	0.220
Rhyacotriton spp.
Sampling method	1	38	0.4	0.530	1	1	0.0	0.921
Habitat	1	38	8.7	0.006	1	1	0.1	0.820
Buffer	4	38	3.0	0.029	2	1	0.0	0.958
Habitat × Buffer	4	38	1.7	0.177	2	1	0.1	0.956
Year	1	38	1.9	0.176	1	1	5.8	0.251
Buffer × Year	4	38	1.8	0.144	2	1	1.4	0.510
Habitat × Year	1	38	1.8	0.186	1	1	0.0	0.965
Habitat × Buffer × Year	4	38	2.1	0.104	2	1	0.5	0.717
Plethodon dunni
Sampling method	1	38	0.1	0.769	1	1	1.9	0.403
Habitat	1	38	5.6	0.023	1	1	1.0	0.508
Buffer	4	38	1.4	0.251	2	1	0.8	0.626
Habitat × Buffer	4	38	3.1	0.026	2	1	0.7	0.646
Year	1	38	0.0	0.954	1	1	0.0	0.869
Buffer × Year	4	38	2.7	0.043	2	1	0.6	0.674
Habitat × Year	1	38	0.4	0.559	1	1	0.3	0.673
Habitat × Buffer × Year	4	38	3.3	0.021	2	1	0.0	0.996

Table 4. Results from mixed model analyses of variance testing for the effects of riparian buffer treatments (1Tree, Variable, Streamside and no-buffer/thin-through vs. untreated controls) prior to and following the second thinning treatment (year), while controlling for sampling method (hand search vs. electrofishing) and habitat structure (NMS Axis 1 scores) on stream-associated vertebrate species; DF_Num = degrees of freedom for the numerator; DF_Den = degrees of freedom for the denominator; Pr = probability.

	Moderate-Density Upland Thinning				High-Density Upland Thinning
Effect	DFNum	DFDen	F-Value	Pr > F	DFNum	DFDen	F-Value	Pr > F
All species
Sampling method	1	38	32.5	<0.0001	1	1	6.3	0.241
Habitat	1	38	0.0	0.876	1	1	2.9	0.336
Buffer	4	38	4.5	0.004	2	1	1.2	0.544
Habitat × Buffer	4	38	3.0	0.030	2	1	3.0	0.379
Year	1	38	0.0	0.947	1	1	2.4	0.367
Buffer × Year	4	38	0.2	0.930	2	1	10.3	0.215
Habitat × Year	1	38	2.2	0.149	1	1	61.7	0.081
Habitat × Buffer × Year	4	38	0.8	0.512	2	1	28.1	0.132
Amphibians
Sampling method	1	38	21.6	<0.0001	1	1	6.3	0.241
Habitat	1	38	0.0	0.921	1	1	2.9	0.336
Buffer	4	38	3.9	0.010	2	1	1.2	0.544
Habitat × Buffer	4	38	3.2	0.024	2	1	3.0	0.379
Year	1	38	0.2	0.666	1	1	2.4	0.367
Buffer × Year	4	38	0.3	0.866	2	1	10.3	0.215
Habitat × Year	1	38	2.6	0.117	1	1	61.7	0.081
Habitat × Buffer × Year	4	38	1.0	0.418	2	1	28.1	0.132
Stream-breeding amphibians
Sampling method	1	38	46.9	<0.0001	1	1	15.9	0.156
Habitat	1	38	0.2	0.645	1	1	0.2	0.753
Buffer	4	38	4.8	0.003	2	1	0.7	0.653
Habitat × Buffer	4	38	2.2	0.091	2	1	0.1	0.907
Year	1	38	0.3	0.608	1	1	0.0	0.953
Buffer × Year	4	38	0.7	0.586	2	1	4.9	0.304
Habitat × Year	1	38	2.5	0.126	1	1	1.5	0.437
Habitat × Buffer × Year	4	38	1.8	0.157	2	1	6.9	0.260
Terrestrial-breeding amphibians
Sampling method	1	38	1.46	0.234	1	1	3.42	0.316
Habitat	1	38	4.32	0.045	1	1	0.2	0.732
Buffer	4	38	1.16	0.343	2	1	0.03	0.967
Habitat × Buffer	4	38	1.42	0.246	2	1	3.47	0.355
Year	1	38	0.04	0.843	1	1	1.32	0.456
Buffer × Year	4	38	0.75	0.566	2	1	0.78	0.626
Habitat × Year	1	38	0.78	0.384	1	1	4.27	0.287
Habitat × Buffer × Year	4	38	1.32	0.280	2	1	1.53	0.497
Pond-breeding amphibians
Sampling method	1	38	0.3	0.615	1	1	4.5	0.280
Habitat	1	38	1.5	0.225	1	1	2.6	0.354
Buffer	4	38	0.3	0.910	2	1	29.4	0.129
Habitat × Buffer	4	38	0.5	0.773	2	1	4.1	0.330
Year	1	38	3.6	0.066	1	1	20.8	0.138
Buffer × Year	4	38	0.8	0.523	2	1	22.5	0.147
Habitat × Year	1	38	0.5	0.503	1	1	11.2	0.185
Habitat × Buffer × Year	4	38	0.7	0.633	2	1	2.3	0.423
Fish
Sampling method	1	38	33.19	<0.0001
Habitat	1	38	1.89	0.1777
Buffer	4	38	1.23	0.315
Habitat × Buffer	4	38	1.58	0.1994
Year	1	38	4.17	0.0482
Buffer × Year	4	38	0.82	0.5195
Habitat × Year	1	38	0.19	0.6681
Habitat × Buffer × Year	4	38	0.48	0.7505
Dicamptodon tenebrosus
Sampling method	1	38	17.0	0.000	1	1	3.4	0.317
Habitat	1	38	10.7	0.002	1	1	31.5	0.112
Buffer	4	38	5.8	0.001	2	1	8.6	0.235
Habitat × Buffer	4	38	2.5	0.061	2	1	4.9	0.305
Year	1	38	0.8	0.370	1	1	9.2	0.203
Buffer × Year	4	38	0.5	0.775	2	1	0.7	0.636
Habitat × Year	1	38	2.7	0.109	1	1	7.7	0.220
Habitat × Buffer × Year	4	38	0.3	0.908	2	1	9.9	0.220
Rhyacotriton spp.
Sampling method	1	38	0.4	0.530	1	1	0.0	0.921
Habitat	1	38	8.7	0.006	1	1	0.1	0.820
Buffer	4	38	3.0	0.029	2	1	0.0	0.958
Habitat × Buffer	4	38	1.7	0.177	2	1	0.1	0.956
Year	1	38	1.9	0.176	1	1	5.8	0.251
Buffer × Year	4	38	1.8	0.144	2	1	1.4	0.510
Habitat × Year	1	38	1.8	0.186	1	1	0.0	0.965
Habitat × Buffer × Year	4	38	2.1	0.104	2	1	0.5	0.717
Plethodon dunni
Sampling method	1	38	0.1	0.769	1	1	1.9	0.403
Habitat	1	38	5.6	0.023	1	1	1.0	0.508
Buffer	4	38	1.4	0.251	2	1	0.8	0.626
Habitat × Buffer	4	38	3.1	0.026	2	1	0.7	0.646
Year	1	38	0.0	0.954	1	1	0.0	0.869
Buffer × Year	4	38	2.7	0.043	2	1	0.6	0.674
Habitat × Year	1	38	0.4	0.559	1	1	0.3	0.673
Habitat × Buffer × Year	4	38	3.3	0.021	2	1	0.0	0.996

For D. tenebrosus, there also was weak evidence for an interaction between habitat and buffer treatment (p = 0.061) and for an interaction between habitat and year (p = 0.109; Table 4). Density of the dominant stream-breeding salamander, D. tenebrosus, increased with increases in NMS Axis 1 scores (Figure 4) in all but the TT treatment, suggesting that this species tracks increases in stream reach depth and width (Table 2). Differences among treatments and years were apparent at high NMS Axis 1 scores, where overall densities were higher (Figure 4). In the 1Tree, VAR and STRM buffer treatments, as well as in untreated controls, densities were lower at high axis scores following treatments than prior to treatments. Because salamander densities in controls followed this pattern, it may be reflective of patterns within the larger population at the site rather than a signal of treatment effects, per se. In the TT treatment, the D. tenebrosus results reflected both a shift in habitat structure and a reversed relationship between density and habitat structure: there was a positive relationship between salamander density and NMS Axis 1 prior to treatment and a negative relationship following the second thinning. This result is consistent with our observation that instream habitats of TT reaches were affected during the second-entry timber harvest and suggests that the harvest may have either altered critical habitats for this species or affected the species more directly during operations.

The density of the subdominant stream-breeding salamander species, Rhyacotriton spp., declined with increases in NMS Axis 1 scores, consistent with previous observations that these species are associated with smaller intermittent streams [54]. Results show increases in density at low Axis 1 scores following the second thinning in all but the narrowest buffer treatment (STRM), where density appears to have declined after the second thinning (Figure 4). Hence, a complex pattern emerges suggesting a positive effect of thinning in wider stream buffers and an adverse effect in the narrowest streamside buffers. There was no effect of thinning in the TT treatment on Rhyacotriton spp., although densities were low in these reaches, which may have reduced our ability to examine the effects (Figure 4).

Figure 4. Estimates from generalized mixed model analysis of co-variance testing for the effects of riparian buffer treatments after a second moderate-density thinning treatment: 1Tree, variable (VAR), streamside (STRM), thin-through (TT) and year (pre-treatment = open symbols, solid lines; post-treatment = solid symbols, dashed lines) while controlling for habitat (non-metric multidimensional scaling (NMS) Axis 1) and sampling method (hand sampling vs. electrofishing); the inset for D. tenebrosus shows reaches sampled using electrofishing; where significant effects of the sampling method were not observed, methods were plotted together.

Figure 4. Estimates from generalized mixed model analysis of co-variance testing for the effects of riparian buffer treatments after a second moderate-density thinning treatment: 1Tree, variable (VAR), streamside (STRM), thin-through (TT) and year (pre-treatment = open symbols, solid lines; post-treatment = solid symbols, dashed lines) while controlling for habitat (non-metric multidimensional scaling (NMS) Axis 1) and sampling method (hand sampling vs. electrofishing); the inset for D. tenebrosus shows reaches sampled using electrofishing; where significant effects of the sampling method were not observed, methods were plotted together.

3.3. High-Density Thinning

For the HD upland thinning treatment, there was little evidence for relationships between population densities and fixed effects. There was some evidence for an interaction between habitat and year for all species combined and amphibian species (p = 0.081) and for a three-way interaction between habitat, buffer treatment and year (p = 0.132). This effect appears to be related to the broader range of habitat conditions in the thin-through buffer treatment reaches, resulting in higher densities of animals at high NMS Axis 1 scores, particularly following the thin-through treatment (Figure 5). Because sample sizes were particularly small for this analysis, p-values >0.05 may be reflective of emerging ecological patterns, yet we advise caution in interpreting these results too broadly. Fish were not observed in any HD reaches.

Figure 5. Estimates from generalized mixed model analysis of co-variance testing for the effects of riparian buffer treatments after a second high-density thinning treatment: 1Tree, variable (VAR) and thin-through (TT) and year (pre-treatment = open symbols, solid lines, post-treatment = solid symbols, dashed lines) while controlling for habitat (NMS Axis 1) and sampling method (hand sampling vs. electrofishing). Because significant effects of sampling method (hand sampling vs. electrofishing) were not observed, the methods were plotted together.

Figure 5. Estimates from generalized mixed model analysis of co-variance testing for the effects of riparian buffer treatments after a second high-density thinning treatment: 1Tree, variable (VAR) and thin-through (TT) and year (pre-treatment = open symbols, solid lines, post-treatment = solid symbols, dashed lines) while controlling for habitat (NMS Axis 1) and sampling method (hand sampling vs. electrofishing). Because significant effects of sampling method (hand sampling vs. electrofishing) were not observed, the methods were plotted together.

4. Conclusions

The lens of both scientists and managers is focusing on the role of forested streamside buffers in the retention of aquatic ecosystem services. Many questions are being posed. How wide do no-entry riparian buffer zones need to be to retain sensitive fish and amphibian species or their critical habitat conditions? Do we need a no-entry riparian buffer if we are employing ecological forestry approaches with green-tree retention, rather than clearcut harvesting? To achieve biodiversity goals in managed forests where wood commodity production is a priority, what is the efficacy of an inner streamside no-harvest zone relative to a thinned outer riparian zone? If such an approach has merit for biodiversity objectives, how wide does each of these areas need to be? In some contexts, riparian thinning has been proposed to accelerate the growth of larger trees for future large down wood recruitment; a habitat benefit to both aquatic-riparian and terrestrial species. However, there is controversy about riparian restoration approaches, including how thinning for restoration purposes might be conducted without risking various potential adverse effects on aquatic-dependent fauna. These questions have been at the forefront of DMS goals to develop and assess integrated riparian and upland forest restoration alternatives [59,76].

In our study, we examined the effects of a second-entry forest thinning harvest to a moderate residual tree density (MD) with four alternative riparian buffer widths relative to untreated controls on headwater stream vertebrates and their habitats. Several key findings emerged. First, with respect to stream habitat attributes: (1) instream habitat did not appear to be strongly affected by the three riparian buffers having a no-entry zone; and (2) stream channels tended to be broader and deeper in thinned buffer zones post-harvest in moderate density thinning treatments, suggesting that thinning affected the stream prism more directly. The pattern we observed for enlarged stream channels within the reaches with a thin-through riparian buffer suggests either direct disturbance to the stream channel during harvest implementation or increased erosion immediately post-harvest.

Second, in the MD treatment, there were patterns of animal densities related to the second thinning, with a complex story developing by species, but not by species groups. The species groups analyses were not informative; a few patterns emerged that were reflective of the dominant species examined here. Both P. dunni and Rhyacotriton species showed patterns of increased densities post-harvest within stream reaches having 70-m and 15-m buffers. The Rhyacotriton pattern is of interest because this genus is of special concern in the Pacific Northwest, with perceived sensitivities to increased water temperatures and sedimentation that may result from timber harvest [77,78]. Rhyacotriton densities increased post-harvest in all stream reaches, but the streamside-retention buffer width, which measures 6 m on each side of streams. This narrow band may be insufficient to retain critical habitat conditions for this animal, such as microclimate [79,80]. However, continued monitoring is required to ascertain if the density decrease we detected is transient or will be sustained over time. Nevertheless, with two thinning treatments, a potential risk of this narrow buffer has been identified for this species.

The dominant stream-bank salamander species, P. dunni, showed the strongest evidence for increases over time in the 70 m and 15 m-buffered reaches, but not in reaches with the thin-through buffer. This finding is of special interest, because P. dunni exhibited a significant drop in numbers in the 6-m buffer reaches previously, 10 years after the first thinning; hence, those depressed numbers appear to have been sustained after the second thinning. Furthermore, the thin-through buffers appear to have adversely affected both P. dunni and D. tenebrosus. With our documented enlargement of thin-through stream reach channels, this suggests that both instream and stream bank salamander habitats were adversely affected by the forestry operations during the thin-through harvest in the MD treatment. We did not quantify stream-bank habitats where P. dunni occur, so this suggestion is not definitively supported by our data. However, our field crew observed tracks from harvest equipment in some streams, places where it appeared that logs had been dragged through streams and habitat disruptions to stream bank areas within thin-through buffers in MD thinning treatment [75]. Because this thin-through effect on animals was not confirmed by data from our HD uplands study component, the results may be due to local reach-specific effects of timber harvest activities [75]. Also, because thin-through treatments were located primarily on smaller streams with steep topography, our design may be prone to haphazard disturbances to small streams. This may be especially the case in areas where timber harvest contractors have not been in the habit of protecting small fishless headwater streams, for example on many private lands [46]. Longer-term (i.e., Year 5 and Year 10) post-treatment monitoring of responses to the second-entry thinning will track these patterns further to examine the longer-term effects on salamander occurrences.

We also report a reduction of D. tenebrosus densities in and along stream reaches across study sites, inclusive of control units. This suggests that there may be a larger population pattern emerging for the managed forest landscapes within which our study sites are placed. Dicamptodon salamanders are not restricted to streams, as aquatic larvae transform to terrestrial adults and occur in forested uplands. They are perceived to range much more broadly than other stream-dependent salamanders, such as Rhyacotriton or P. dunni (P. dunni movements were largely within 5 m of streams and averaged <5 m in a cover-board study at one of our DMS sites [58]). Although home ranges of Dicamptodon salamanders are not precisely known, the stream-breeders have been found several hundred meters from streams [47]. Hence, patterns that we are observing in and along streams may reflect animal responses, not only locally to our stream buffer treatments, but to activities ongoing in a larger area of the forested landscape.

Our study has been implemented coincident with management directions for federal lands from the Northwest Forest Plan [14], on federal lands managed by the USA Bureau of Land Management. These lands are arrayed in a checkerboard land-ownership pattern with private timberlands (e.g., [81]). Since 1994, private lands have moved to a more frequent timber harvest interval, now approaching 30–40 years, whereas federal lands in land use allocations for wood commodity production are on longer intervals for regeneration harvests, with the dominant intervening activity being thinning (e.g., [82]). A spatial pattern has emerged in this checkerboard landscape, such as that depicted by our Keel Mountain study site (Figure 5). The influence of forest cutting on the pattern of structural conditions and age classes is visually apparent, with a signature of buffers along larger streams connecting several patches. Most headwater streams on private lands in the region are not buffered or have considerably less protection during timber harvest [46]. With many of the forest stands on federal lands nearing their rotation ages, increased harvest activities are projected, potentially increasing the extent of cutting within this checkerboard landscape. The importance of our DMS sites is emerging to track the responses of aquatic-dependent vertebrates, not only to our treatments, but perhaps to this larger-scale influence, as well. Dicamptodon density patterns reported in this study may be particularly relevant to track, as their composite stream and upland habitat requirements may be disturbed to a greater extent than other species that we are following.

It is important to weigh both the short- and long-term risks and benefits of forest management actions. For ecological restoration, a short-term risk for a long-term benefit might be an option. However, whereas risks to species may be considered for broader ecological restoration goals, there may be alternative management strategies that provide species, ecological and economic benefits, especially when sensitive species are addressed. For example, P. dunni is a species of concern in Washington State, but not Oregon [83], and Rhyacotriton spp. are of concern in both states [78,83,84]. In western Oregon, the relevance of a 15-m riparian buffer width has been examined from aquatic, riparian and upland perspectives and may provide both short- and long-term benefits to a variety of ecosystem services, including retaining these sensitive species. At DMS sites, riparian forest restoration benefits resulting thinning, including a boost in overstory tree growth within no-entry riparian buffers within ~15 m of the riparian-upland forest edge [85]. In secondary forests, a 15-m stream effect on riparian microclimate has been documented as the distance from streams where gradients in cool, moist near-stream microclimates shifted to more xeric upland conditions [79,80]. These conditions likely contribute to the near-stream zone within 15 m of streams, emerging as important for some plethodontid salamander distributions and life history functions, such as habitat associations and movements [57,58]. Herein, we report an increase in the densities of both P. dunni and Rhyacotriton species after a second thinning treatment within stream reaches with a 15-m buffer. Hence, 15 m appears to be a useful riparian buffer width for the retention or development of riparian biodiversity and ecosystem services, including forest-dependent salamanders with species-of-concern status in thinned forests; yet, caution is still needed, due to uncertainty regarding long-term impacts. Although recent studies have focused largely on the effects of upland thinning regimes with riparian buffers in headwaters, the utility of streamside buffers of varying widths within upland clearcuts and other harvests with lower levels of overstory retention remains largely unresolved for both small headwater streams and larger downstream reaches.

Figure 6. Aerial photograph of the Keel Mountain study site (yellow square), showing forest management activities within the landscape on private and federal lands arrayed in a checkerboard pattern.

Figure 6. Aerial photograph of the Keel Mountain study site (yellow square), showing forest management activities within the landscape on private and federal lands arrayed in a checkerboard pattern.