Beaver Colony Density Trends on the Chequamegon-Nicolet National Forest, 1987 – 2013

The North American beaver (Castor canadensis) is a managed species in the United States. In northern Wisconsin, as part of the state-wide beaver management program, the Chequamegon-Nicolet National Forest removes beavers from targeted trout streams on U.S. Forest Service lands. However, the success of this management program has not been evaluated. Targeted removals comprise only 3% of the annual beaver harvest, a level of effort that may not affect the beaver population. We used colony location data along Forest streams from 1987–2013 (Nicolet, northeast Wisconsin) and 1997–2013 (Chequamegon, northwest Wisconsin) to assess trends in beaver colony density on targeted trout streams compared to non-targeted streams. On the Chequamegon, colony density on non-targeted trout and non-trout streams did not change over time, while colony density on targeted trout streams declined and then stabilized. On the Nicolet, beaver colony density decreased on both non-targeted streams and targeted trout streams. However, colony density on targeted trout streams declined faster. The impact of targeted trapping was similar across the two sides of the Forest (60% reduction relative to non-targeted trout streams). Exploratory analyses of weather influences found that very dry conditions and severe winters were associated with transient reductions in beaver colony density on non-targeted streams on both sides of the Forest. Our findings may help land management agencies weigh more finely calibrated beaver control measures against continued large-scale removal programs.


Introduction
In Wisconsin, as in many states, the North American beaver (Castor canadensis) is a managed species. Wisconsin's management plan addresses a wide variety of issues important to citizens a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Chequamegon-Nicolet NF. Both sides of the Forest are located in the northern highland province of Wisconsin [8]. This province is generally flat to gently rolling with little topographic variation [9] due to hundreds of millions of years of erosion combined with the effects of repeated glaciation cycles [8]. Streams in the Chequamegon-Nicolet NF are considered to be headwater streams, with over 80% of streams on both sides of the Forest being stream order 1 or 2, and predominantly less than 15 m wide. The streams typically occur in old glacial drainages and wetlands with broad valleys having low slopes, and sand or gravel as dominant sediment [9]. Stream gradient averages 8.9 m/km on the Chequamegon side and 4.1 m/km on the Nicolet side; because of these low gradients, streams are easily dammed [9]. Water temperatures are regulated by groundwater inflow and degree of vegetative shading [10]. The Forest's vegetation is composed of mesic mixed northern hardwoods [black ash (Fraxinus nigra), American elm (Ulmus americana), red maple (Acer rubrum), sugar maple (Acer saccharum), basswood (Tilia americana), oak (Quercus spp.)], conifer [jack pine (Pinus banksiana), red pine (Pinus resinosa), balsam fir (Abies balsamea), eastern white pine (Pinus The average annual temperature is about 4.4˚C in northern Wisconsin [11]. During more than one-half of the winters, temperatures fall to -40˚C or lower; summer temperatures above 32.2˚C occur on average two to four days in northern Wisconsin [11]. However, there is a significant warming trend in northwestern Wisconsin [12]. The long-term mean annual precipitation is 12.6 cm in northern Wisconsin [11] but summer precipitation has declined while autumn precipitation has increased [12]. Average seasonal snowfall is 164 cm and the average date of first snowfall in northern Wisconsin is in early November [11]. Many of the streams and lakes in northern Wisconsin are ice-covered from late November to early April [11].

Beaver Control and Monitoring Program
The Chequamegon-Nicolet NF classifies its streams into "trout" and "non-trout", based on the Wisconsin Department of Natural Resources (WDNR) categorization of streams into three trout-relevant classes [13]. "Trout" streams are high quality with natural or stocked trout populations (i.e., WDNR Class 1 or 2). "Non-trout" streams correspond to WDNR Class 3 streams (marginal habitat quality with no natural reproduction). As part of the state beaver management plan, the WDNR determined the high-quality trout streams that needed to be protected using beaver removal. The Chequamegon-Nicolet NF determined which of those streams fell within the borders of the Forest and designated trout streams as "targeted" if those were subject to beaver removal and "non-targeted" for streams that were not. If the WDNR-designated trout streams were in wilderness areas (where no management activities are allowed) or ran through private lands in the Forest with no permission to conduct beaver management from the landowner, those trout streams were moved into the "non-targeted" category. Non-trout streams that were tributaries of targeted trout streams were moved into the "targeted" category. Targeted and non-targeted streams are displayed in Fig 1, and the associated lengths for each of the four categories are presented in Table 1. The beaver management program started ten years later on the Chequamegon side (1997) than the Nicolet side (1987). As noted above, the Forest is open for regulated trapping but it does not collect information on how many beavers are removed. The WDNR estimates the state-wide beaver harvest using voluntary responses to a questionnaire mailed to a sample of licensed trappers [1].
The Forest's management program had two components. The first component consisted of contracting the removal of all known beaver colonies and lodges from targeted streams in the spring (April to early May) of each year by USDA-WS [14,15]. Note that beaver control program activities are not subject to IACUC review, as they are land management activities, not research activities. However, control program activities were conducted in accordance with state law and best management practices sponsored by the Association of Fish and Wildlife Agencies as codified in USDA-WS directives [16].
The second component was fall monitoring of beaver locations on targeted and non-targeted streams to assist in locating active beaver colonies on targeted streams for removal the following spring. The Forest contracted with USDA-WS to conduct fall aerial surveys to document the location of active beaver colonies along its land base ( Table 1). The fall flights were land management activities, not research activities, and therefore not subject to IACUC review. Aerial surveys for beaver colonies have been found to be an effective method for examining abundance patterns [6,17].
Aerial fall flights were conducted after leaf-off, usually mid-October-November. Flights were made between 0900 and 1600 (Central Standard Time) at 213-244 m above the surface, flying at 128-161 km/hr. The Nicolet side was flown by the same contract pilot since monitoring began in 1987. The Chequamegon side was flown by a variety of WDNR pilots. Completing both surveys required approximately 40 hours of flight (approximately 24 hours for Nicolet side and 16 hours for Chequamegon side), and had to be completed after leaf-off and before snowfall to adequately observe current beaver activity. Since 2007, flight paths have been tracked using GPS. Because of the short window between leaf-off and first snowfall, flights were conducted on consecutive days if the pilot assessed the weather conditions to be safe for flying. Adjusting the flight schedule for optimal weather conditions to increase observation conditions was not possible because of the short observation window and large land base to cover.
One USDA-WS biologist made all observations on the Chequamegon side since the flights began; on the Nicolet side, there were 4 biologists. On each flight, the biologist documented the location of colonies based on evidence of recent beaver activity. A colony was considered active if there was a lodge with fresh mud or cuttings and/or a food cache with fresh cuttings. Active colonies were recorded along the length of the stream, typically on one pass. However, if there was heavy conifer cover or steep terrain, the pilot would make a second pass to observe the stream from a different vantage point. The pass was conducted such that the stream was on the observer (or right) side of the plane. Active colonies were marked on a map and then digitized from the map into a geographic information system. Accuracy was estimated at ±322 m. Maps of beaver colony locations made from the aerial surveys were made available to the public by request, and the Forest informed the Wisconsin Trappers Association when the maps became available. No information was available about how regulated trappers used those maps to target streams on the Forest (as noted above).
For this study, we used beaver colony count data from 1987-2013 for the Nicolet side and 1997-2013 for the Chequamegon side. We used data for all stream systems flown for !5 years. The response variable was colony density (active beaver colony count/km stream flown) by stream category, year, and side of Forest. Because few non-trout streams were targeted (Table 1), we did not use this category in any further analyses.

Weather Variables
To explore whether weather variability could explain inter-year variation or trends in beaver colony densities on non-targeted streams, we used the beaver literature to determine what variables to consider in the analysis. We were specifically interested in weather variables that affected survival of the 2-year old beavers that annually disperse from their natal dens to establish their own independent colonies. In a population viability analysis, early-age class survival was critical for beaver population persistence, and large variation in early-age class survival rates made the beaver population vulnerable to extinction [18].
A beaver colony consists of a pair of adult beavers and their young. Young stay in the lodge until they are 2 years of age [19], when they disperse in the spring (April-May) [20]. Proximate factors affecting beaver litter size and survival of the young in the den include quantity and quality of available food and severity of winter weather [18,19]. Quantity and quality of food resources is affected by weather conditions [21,22]. For the European beaver (Castor fiber) in Norway, higher rainfall in the fall prior to the spring when females give birth led to smaller litter sizes [22]. The higher rainfall stressed important food plants by waterlogging their roots, resulting in lowered nutrient content for beavers during the reproductive period [22]. Reduced survival of young in the den (leading to fewer 2-year-olds available for dispersal) was associated with warmer spring temperatures and increased spring-summer precipitation, both factors that can reduce the nutritional value of plants [21,22]. However, warmer and wetter conditions during spring dispersal would be favorable for the survival of the 2-year-olds as they move away from the den [23].
In northern areas, beavers rely on winter food caches, generally set aside after the first heavy frost [18]. Severely cold winters can prompt early depletion of food caches, leading to physiological stress and reduced survival [24,25]. This effect can be offset by heavy snowfall, which improves the insulation of the lodge, reducing thermoregulatory energy demands [24,25]; lighter than average snowfall would have the opposite effect.
Thus, we chose four broad weather variables as affecting female production of young and survival of 2-year old beavers: temperature during late spring (when plants start to grow, i.e., time when green-up occurs), precipitation as reflected in soil moisture during the growing season, temperature during the winter, and winter snowfall. All weather variables were lagged by up to three years to reflect environmental conditions and potential food nutritional value when the young were in the lodge (lag of 0 to 2 years) or for the female before breeding (lag of 3 years). We then defined our measures of these variables to reflect conditions in northern Wisconsin (Table 2). In northern Wisconsin, green-up occurs in May, and the growing season for aspen (Populus spp.) extends through September [26]. In northern Wisconsin, the first Table 2. Variables used to assess the ability of weather to explain minor variations in the beaver density index on non-managed streams or function as drivers of trends on those streams. Each variable was computed for the current year and the three preceding years. Variables were standardized against their long-term averages (1983-2013).

Variable description Period of coverage
Mean maximum daily temperature Green-up: May heavy frost occurred, on average, in late August to early September [11] though there was a shift of 3-12 days later in the 2000s [12]. We defined winter to be December-March because these months represent the coldest months in northern Wisconsin, and the months when the heaviest snowfall is expected. To represent severe winter conditions, we used annual cumulative snowfall and number of days below −17.8˚C (0˚F). We standardized all weather variables by their long-term (1983-2013) means. This standard statistical transformation put extreme events on the same scale (i.e., anomalous conditions will have large negative or positive values, regardless of the variable) and also made variables comparable across the sides of the Forest. We used the 30-year time span so that the overall average conditions were based on data encompassing the most variability in weather conditions against which to judge extreme events (i.e., incorporates the effects of climate change [12]).
All temperature data were downloaded from PRISM [27], imported into ArcMap [28] using a Python for ArcGIS (ArcPy) script, and given a NAD1983 projection if the data did not already have a projection. PRISM data are provided at a 30 m x 30 m scale of resolution so there are multiple 30 m x 30 m cells per district. As part of the ArcPy script, the tool "Zonal Stats as Table" was then run using the Forest's management districts as zones. Using the mean function in the tool, data from the 30 m x 30 m cells were averaged over the districts by side of the Forest. For the Chequamegon side, we used a weighted average where the weights were the number of miles flown within a district during the fall flights.
For temperature conditions during the month when green-up occurs, we calculated the mean maximum daily temperature for May of each year. To represent precipitation conditions during the growing season each year, we used the Palmer Drought Severity Index (PDSI), a standardized index for how much water is available to plants [29]. The Nicolet side was covered by the PDSI Northeast Division (4703). On the Chequamegon side, the Washburn/Great Divide district was covered by the PDSI Northwest Division (4701) and the Park Falls and Medford districts were covered by the PDSI North Central Division (4702). An overall Chequamegon PDSI was an average of PDSI for Divisions 4701 and 4702 weighted by the number of miles flown during the fall flights for the districts.
We used cumulative snowfall data from National Weather Service Cooperative Observer Program (COOP) weather stations nearest the Chequamegon-Nicolet NF (Fig 1). Monthly snow data were downloaded from the National Climatic Data Center [30]. Cumulative snowfall for the Chequamegon side was an average weighted by the number of miles flown during the fall flights for the district closest to the COOP stations. Cumulative snowfall for the Nicolet side was the average of the data from the COOP stations.
Number of days with an average daily temperature less than −17.8˚C was calculated from daily temperatures for December-March from the PRISM website [27]. For each year, number of days with an average daily temperature less than −17.8˚C was calculated by month and then summed for December-March.

Statistical Models
We analyzed time trends by side of the Forest (Chequamegon, Nicolet) and by 3 stream categories (targeted trout, non-targeted trout, non-targeted non-trout). For time trend analyses on non-targeted streams, we used the trout or non-trout classification as a variable in the model to determine differences between trout and non-trout streams. To explore the influence of weather, we focused on non-targeted streams to avoid the confounding effect of beaver management. We added each weather variable to either a time trend or a constant model (if there was no time trend).
We used generalized additive models with a gamma of 1.4 to avoid overfitting [31] for all analyses. For the time trend models, we also used segmented linear models to estimate when changes in slope occurred. We used a Gaussian error structure with no autocorrelation in all models. To compare models, we used Akaike Information Criteria corrected for small sample size (AIC c ) [32]. We calculated evidence ratios using AICc weights to determine how good the minimum AIC c model was with respect to the other models [32]. We considered models within 2 AIC c units of the minimum AIC c model to be competitive models [32]. We used the Davies test to determine the significance of the segmented linear models. All analyses were done with R version 3.1.2 [33]. Significance was assessed at a P of 0.05 and trends at 0.10.

Time Trends
On the Chequamegon side, there was no significant linear trend in beaver colony density on non-targeted streams (Year term: t = −1.1, df = 30, P = 0.29), regardless of whether the stream was a trout or non-trout stream (Trout Classification Ã Year interaction term, t = 0.169, df = 30, P = 0.87) (Fig 2a). However, there was a difference in colony density between non-targeted trout and non-targeted non-trout streams (Stream Classification term: t = 4.45, df = 30, P < 0.001). Colony density was 0.226 colonies/km (SE = 0.005, n = 17) on non-targeted trout streams compared to 0.173 (SE = 0.005, n = 17) for non-targeted non-trout streams.

Exploratory Analysis of Beaver Colony-Weather Relationships on Nontargeted Streams
On the Chequamegon side, when we added weather variables to a model controlling for the differences in beaver colony density on the non-targeted streams, there was a significant positive relationship with June-July PDSI in the year prior to the fall flight (s(June-July PDSI lag 1): F = 9.07, df = 5.3, 27.7, P < 0.001; adjusted R 2 = 0.77). However, 2008 beaver colony density was an influential point. Specifically, 2008 had the lowest colony density in the series and was linked to the 2007 June-July PDSI, which was the lowest value in the PDSI series and the only one that fell into the PDSI severe drought category. When the 2008 point was removed, the June-July PDSI in the year prior to the fall flight was no longer significant (June-July PDSI lag1: t = −0.00002, df = 30, P = 0.99). When we added an indicator variable for the 2007 drought to the model adjusting for different densities on the non-managed streams, the influence of the severe drought in 2007 was significant (2007 severe drought: t = −4.87, df = 31, P < 0.0001; adjusted R 2 = 0.627). When we added weather variables to the 2007 drought model, the minimum AIC model linked colony density to May temperature 3 years before the fall flight (i.e., conditions affecting females before pregnancy). A competitive model included severe winter temperatures the year of the fall flight (i.e., winter before dispersal of the 2-yearolds) ( Table 3). For the minimum AIC model, there were higher colony densities in year i when there were cooler May temperatures three years prior to the fall flight (coefficient = −0.006; t = −2.88, df = 30, P = 0.007; adjusted R 2 = 0.698). This model was only 1.6 times better than the competitive model (severe winter temperatures) but 16 times better than the 2007 drought model ( Table 3). The competitive model indicated that there were higher colony densities when the winter before the fall flight was less severe (coefficient = −0.0009; t = −2.7, df = 30, P = 0.011; adjusted R 2 = 0.689).
On the Nicolet side, when we added the weather variables to the time trend model for nontargeted streams, there was a significant non-linear relationship with May PDSI three years prior to the fall flight (s(May PDSI lag 3): F = 3.14, df = 2.3, 45.35, P = 0.037; adjusted R 2 = 0.896) ( Table 4). Specifically, years that were in moderate drought or, conversely, moderately moist or wetter in the years before the 2-year-olds were born were associated with years with lower colony densities when those 2-year-olds dispersed. Using a baseline year of 1990 (to control for year), average predicted colony density when May PDSI values were in the moderate drought category (−2 to −2.99) was 0.196 colonies/km (standard error = 0.0009, n = 7), 0.209 colonies/km (standard error = 0.0015, n = 10) for the moderately moist category (2 to 2.99), and 0.222 colonies/km (standard error = 0.0005, n = 11) for mid-range values (−0.5 to 0.5) There were two competitive models. The first competitive model included August-October PDSI three years prior to the fall flight (i.e., food caching period before the females become pregnant) (s(Aug-Oct PDSI lag 3): F = 2.77, df = 1.9, 42.7, P = 0.064; adjusted R 2 = 0.896) Table 3. Weather models for beaver colony density on streams where no targeted beaver management occurred on the Chequamegon side of the Chequamegon-Nicolet National Forest, 1997-2013. Models listed are the minimum AIC model, competitive models (those models within 2 AC units of the 006Dinimum), and the 2007 drought model. The evidence ratio is the ratio of the Akaike weights for the minimum AIC model and an alternate model.

Evidence ratio
Trout stream type + Indicator of severe drought in 2007 + Deviance of mean maximum May temperature three years before fall flight

Discussion
After 20 years of beaver control, we found a substantial and sustained Forest-wide reduction in beaver colony density on targeted trout streams across the Chequamegon-Nicolet National Forest. The magnitude of the reduction on both sides of the Forest was similar at roughly 60% relative to non-targeted trout streams-the targeted streams were regularly recolonized, but well below pre-control levels. Using the Chequamegon side as the benchmark, it took seven years of effort to achieve the full reduction, which was then sustained over the remainder of the monitoring period. The percentage declines are similar to the declines (but not the numerical harvest levels) documented in beaver pelt harvest records of the mid-1800s when commercial beaver trapping was occurring [34].
The Forest conducted a test in the early 1990s using only regulated trapping on a subset of its trout streams. This approach was not able to meet beaver reduction goals, and the test streams were reassigned to the targeted trapping program. It is an interesting question why regulated trapping was not able to generate a reduction similar to that found in the targeted beaver population. We speculate there was interplay between effort and population dynamics. Absent very high pelt prices, regulated trappers might have skipped pursuing beaver colonies that were located on harder-to-access reaches of the managed streams. This could have resulted in persistent colonies producing dispersing young that could colonize the parts of the streams that were accessible by trappers, keeping the annual fall flight population index data near unmanaged levels. USDA-WS trappers, on the other hand, are paid to get into harder-toaccess areas and so might more effectively reduce the targeted beaver population.
While our analysis was based on a population index rather than a population estimate, our finding of higher colony densities on the Chequamegon side compared to the Nicolet side of the Forest is consistent with beaver population estimates for northern Wisconsin. Historically, after the over-trapping era of the late 1890s, remnant beaver populations remained in the far northwestern part of the state [35]. Currently, the beaver population in northwestern Wisconsin (location of the Chequamegon) is consistently higher than that for northeastern Wisconsin (location of the Nicolet) [36,37]. Our finding of different time trends between the two sides of the Forest for colonies on non-targeted streams is also consistent with the time trends seen in the beaver population in northern Wisconsin. While the beaver population in northern Wisconsin declined overall between 1992 and 2008, the decline was steepest in northeastern Wisconsin and a sharp decline occurred between 1995 and 1998 [36]. The 1998 population estimate in northeastern Wisconsin declined by 47% from the 1995 estimate while the population in northwestern Wisconsin declined only 13% (calculated from [36]). The influence of the 2007 severe drought that we detected on the Chequamegon was likely a contributing factor to the 32% decline in the northwestern Wisconsin beaver population between 2005 and 2008, the steepest decline seen in the northwest beaver population (calculated from [36]). Environmental conditions can be important factors causing mortality in beaver populations [18]. While our evaluation of weather influences was exploratory, the weather signals we found correspond to factors found in other studies that affect the size of the cohort of 2-yearold dispersers, either through a linkage with survival while in the den or to cohort size via female condition [21,22]. In our study, though the specific weather variables associated with beaver colony density on the two sides of the Forest were different, the commonality was that dry conditions during spring and severe winters were associated with reduced beaver colony density in the fall (either the same year or lagged). Drought was implicated in the depression of local beaver populations in western Nebraska [18] and has been suggested as a cause of mortality in beaver in Kansas [38]. Severe winters and starvation have also been noted as principal causes of mortality in northern climates [6].
Beaver populations without extensive trapping can have complex dynamics including declines [39] and nonlinear trends [40]. The pattern we found on the Chequamegon side of the Forest (decline on the targeted trout streams and no pattern on the non-targeted streams) was the expected pattern if regulated trapping was at sustainable levels and randomly distributed across the streamscape, making targeted trapping the differentiating factor affecting colony density on trout streams. We did not find that pattern on the Nicolet side of the Forest. Instead we found substantial declines on the non-targeted streams as well. On the Nicolet, the observed 65% decline in the colony index from 1994-2013 could be attributable to either a set of local factors driving the population down apart from the targeted beaver management program (e.g., local concentrations of regulated trapping, changes in habitat-based carrying capacity) or it could be that targeted streams in the broader landscape had previously provided a source of colonizing beavers for the Forest (i.e., source-sink dynamics [39]). In northeast Wisconsin counties, 24% to 67% of total perennial stream kilometers are under targeted beaver management [1]. In contrast, in northwest Wisconsin, only in one county does targeted beaver management exceed 25% of perennial stream kilometers with most other counties having <10% of stream kilometers under targeted beaver management ( Figure 10 in [1]). While the interaction between streams within the Forest and the broader landscape is not known, we do know that new beaver colonies on the targeted streams in the Forest arise only through new colonization each year and these colonies are removed the following spring, before they can produce dispersing young. Currently, even with declining populations, the Nicolet non-targeted streams are providing colonizers for the targeted streams at the same relative rate as occurs in the Chequamegon system (generating a targeted colony density 60% lower than the non-targeted colony density in both the Chequamegon and Nicolet systems).
Many public land management organizations are faced with optimizing and prioritizing resource allocations to competing management objectives. The question is whether the desired cold-water ecosystem and trout benefits can be maintained if the intensity of targeted beaver management, and thus its cost, is reduced. Our finding that colony densities vary temporally and spatially may help refine management options across the Forest and the state to lower the cost of beaver management. Relying more heavily on targeting densities by region, watersheds, specific locations, or by stream systems may be a means of transitioning to variable trapping pressures over larger regions. For example, stable cold-water groundwater inputs or extensive groundwater inputs combined with vegetative stream shading were able to maintain downstream water temperatures in northern Wisconsin despite the presence of beaver dams [10]. Thus, "judicious removal of a few key beaver impoundments may have a strong effect on downstream temperatures without the loss of diversity associated with large-scale dam removal" [10]. Management programs using variable trapping can also integrate other beaver damage factors beyond fisheries such as infrastructure and bank stabilization. While understanding local impacts of beavers on trout fisheries is an important research topic in the state of Wisconsin [1], further investigation on how local stream characteristics plus surrounding landscapes influence beaver recolonization events and territory settlement patterns will be important to help land management agencies refine large-scale beaver management programs.