Restricted cross-scale habitat selection by American beavers.

Animal habitat selection, among other ecological phenomena, is spatially scale dependent. Habitat selection by American beavers Castor canadensis (hereafter, beaver) has been studied at singular spatial scales, but to date no research addresses multi-scale selection. Our objectives were to determine if beaver habitat selection was specialized to semiaquatic habitats and if variables explaining habitat selection are consistent between landscape and fine spatial scales. We built maximum entropy (MaxEnt) models to relate landscape-scale presence-only data to landscape variables, and used generalized linear mixed models to evaluate fine spatial scale habitat selection using global positioning system (GPS) relocation data. Explanatory variables between the landscape and fine spatial scale were compared for consistency. Our findings suggested that beaver habitat selection at coarse (study area) and fine (within home range) scales was congruent, and was influenced by increasing amounts of woody wetland edge density and shrub edge density, and decreasing amounts of open water edge density. Habitat suitability at the landscape scale also increased with decreasing amounts of grass frequency. As territorial, central-place foragers, beavers likely trade-off open water edge density (i.e., smaller non-forested wetlands or lodges closer to banks) for defense and shorter distances to forage and obtain construction material. Woody plants along edges and expanses of open water for predator avoidance may limit beaver fitness and subsequently determine beaver habitat selection.

Habitat selection by animals often exhibits scale-dependent patterns like many other ecological phenomena (Johnson 1980; Wiens 1989; Levin 1992). Advances in spatial technologies and the relative ease of fitting resource selection functions have led to a proliferation of recent habitat selection studies (Northrup et al. 2013). However, many of these studies do not include multi-scale analyses (McGarigal et al. 2016). While some multi-scale studies found inconsistent habitat selection patterns between coarse and fine scales (Corriale and Herrera 2014; Peters et al. 2015), others have shown congruence between spatial scales (Crampton and Sedinger 2011; Prokopenko et al. 2017). However, congruent selection patterns between spatial scales may not translate to habitat types or temporal scales (Crampton and Sedinger 2011). Detecting habitat selection relies on scales of measurement and analysis, and one scale is often insufficient to predict habitat selection at another scale (Mayor et al. 2009).

Variables explaining fine-scale habitat selection can be influenced by behavioral decisions made at a coarse spatial scale (McGreer et al. 2015), which we termed the down-scaling effect. Conversely, behavioral decisions made at fine scales may dictate patterns at larger spatial scales (Jedlikowski et al. 2016), which we refer to as the up-scaling effect. Key resources such as food, cover, and water may affect habitat selection processes at different spatial scales according to life history traits (Coreau and Martin 2007; Jedlikowski et al. 2016) and species-specific ecological requirements (Perez-Garcia et al. 2014). While animals may select habitat structures at different spatial scales (Boyce 2006), the hierarchy of habitat selection should reflect factors affecting individuals’ fitness (Rettie and Messier 2000). For instance, predator avoidance at large spatial scales has been shown to shape habitat selection at fine spatial scales in woodland caribou Rangifer tarandus caribou, which may ultimately limit fitness (Rettie and Messier 2000). Similarly, elk Cervus elaphus avoided areas of high wolf predation risk at landscape scales, whereas fine-scale habitat selection focused on food resources adjacent to human residence (Hebblewhite and Merrill 2009). Scale-dependent or hierarchical habitat selection has a utility in understanding how animals respond to spatiotemporal variation in resource available and predicting animal spatial distributions. However, hierarchical habitat selection has been challenged because the underlying theory dictates that habitat selection at one level may constrain selection at other levels (Mayor et al. 2009). Unconstrained multi-scale analysis provides an opportunity to identify scale independent effects of landscape structure on habitat selection.

The American beaver Castor canadensis (hereafter, beaver) is a large semiaquatic nocturnal rodent found throughout much of North America (Baker and Hill 2003). Although classified as a “choosey” generalist with their presence in many types of ecosystems, beavers are closely associated with water and wetlands (Baker and Hill 2003; Müller-Schwarze and Sun 2003). Beavers are herbivores, feeding on a variety of woody, non-woody, terrestrial, and aquatic vegetation (Baker and Hill 2003). As a wetland ecosystem engineer, beavers removed a large biomass of woody plants at a rate of 1.4 mg per individual per ha per year in Minnesota, USA (Johnston and Naiman 1990), and reduced aquatic plant biomass and plant litter by 60% and 75%, respectively, in Georgia, USA (Parker et al. 2007). Beaver herbivory may alter edges in riparian or bottomland forests, and beavers often forage within a 60-m distance from water (Donkor and Fryxell 1999; Haarberg and Rosell 2006; Steyaert et al. 2015). Likewise, fine-scale space use (e.g., within home ranges) of semiaquatic mammals may depend on the distributions of food resources in proximity to water (Campbell et al. 2013; Corriale and Herrera 2014).

Therefore, it is plausible to hypothesize that main factors limiting beaver fitness and determining habitat selection at both large and small spatial scales include the availability of woody plants and access to wetland habitats. Nevertheless, beaver habitat selection was rarely examined at multiple spatial scales or with modern probability distribution techniques. Studies of habitat selection by beavers primarily focused on the characterization of dams and lodge site selections using presence-absence or presence-pseudo absence habitat comparisons (Allen et al. 1983; Beier and Barrett 1987).

In this study, our objectives were to determine if beaver habitat selection was specialized to semiaquatic habitats and if variables explaining habitat selection are consistent between landscape and fine spatial scales. Specifically, we predicted that habitat selection by beavers would be positively correlated with edge densities of woody wetland, shrub, open water, and emergent herbaceous wetland at both large and small spatial scales (prediction P1). Additionally, we predicted that beaver habitat selection would be scale dependent, whereby beavers select different habitat variables at large and small spatial scales (prediction P2). We built maximum entropy (MaxEnt) models to relate landscape-scale presence-only data to landscape variables, and used generalized linear mixed models (GLMMs) to evaluate fine spatial scale habitat selection using global positioning system (GPS) relocation data. Explanatory variables between the landscape and fine spatial scale were compared for consistency.

Materials and Methods

Study area

We conducted the study at Redstone Arsenal (RSA), a 15,342-ha military installation managed by the Department of Defense (DOD), located in Madison County, Alabama, USA (34°39'00.5"N 86°37'52.1"W). RSA was bordered by the city of Huntsville to the North and East, Madison County to the west, and the Tennessee River to the south (McClintic et al. 2014a). We collected data at 11 wetland sites spread across RSA, where beavers were active (Figure 1). The landscape containing the 11 sites was composed of agricultural fields, upland pine forests, mixed forests, bottomland hardwood forests, and developed areas (e.g., DOD infrastructure). Average annual temperatures ranged from 5.2 C° in January to 27 C° in July, with an overall annual average of 16.7 C°. Average annual precipitation was 138 cm, ranging from 9 to 15 cm per month (Huntsville-Decatur International Airport weather station).

Figure 1

Habitat suitability map of American beavers in RSA, Alabama, USA. Stars depict the 11 wetlands where beavers were monitored. Habitat suitability index ranges from 0 to 1 as shown in the legend with values closer to 1 representing more suitable habitat. The second panel is a map of the National Land Cover Data for RSA, Alabama, USA with the same 11 monitoring locations.

Species presence locations

We collected data on American beaver presence locations to build habitat suitability models. We documented beaver presence based on beaver markings, structures, and capture locations in 11 wetlands during two periods, March 2015 to August 2015 and February 2016. Markings and structures included dams, lodges, feeding stations, foraging locations, and castor mounds. We recorded geographic coordinates of all presence locations using a handheld GPS unit (Garmin E-Trex 10, Garmin, Olathe, KS, USA). The geographic locations were projected in the Universal Transverse Mercator Zone 16 North (UTM 16 N) using World Geodetic System 1984 (WGS 84) datum.

We located beaver lodges using a combination of techniques: homing to tagged beaver locations using very high frequency (VHF) telemetry and systematic searches by walking/wading and boat. We also used systematic searches to quantify all other beaver signs in each wetland. We identified foraging locations by searching in quadrants radiating from lodges, whereby we recorded all woody vegetation (≥1 cm) browsed by beavers within 2 m plots. We identified active dams by the presence of fresh construction material (e.g., sticks, mud, and rocks). We identified active castor mounds by sight and smell, as beavers reapply castor to these territorial markings frequently. We categorized feeding stations as locations where beaver processed and consumed food. These were usually logs or small islands just above waterline with accumulations of freshly processed sticks and shavings.

Preparation of environmental variables

We generated raster layers of landscape variables as the covariates of habitat suitability models. We obtained land cover data from the National Land Cover Database (NLCD) 2011 with a resolution of 30 × 30 m per grid cell (Homer et al. 2015), reclassified it into 10 land cover and land use classes, and clipped the data to the extent of RSA in ArcMap 10.3 (ESRI, Redlands, CA, USA). The 10 classifications used in this study were grass (NLCD classes 71 and 72), developed (NLCD classes 21, 22, 23, 24, and 31), deciduous forest (NLCD class 41), evergreen forest (NLCD class 42), mixed forest (NLCD class 43), shrub (NLCD class 51 and 52), cultivated crops (NLCD class 82), woody wetlands (NLCD class 90), emergent herbaceous wetlands (NLCD class 95), and open water (NLCD class 11) (Homer et al. 2015). We used the Booleanisator tool in the Biomapper software package (Hirzel et al. 2002) to create a new raster map by landscape cover type. We used the Circular Analysis program within Biomapper (Hirzel et al. 2002) to calculate edge density and relative frequency of occurrence (0–1). This contextual operator can be viewed as a circular-moving window based on a user-defined buffer. Our buffer size (11.88 ha) was equivalent to annual home range size of American beavers previously reported on RSA (McClintic et al. 2014a). We also used program Distance Analysis in Biomapper (Hirzel et al. 2002) to calculate the exact minimum distances from a grid cell to the nearest grid cell of a land cover type.

We performed a principle component analysis (PCA) on 30 landscape layers (i.e., edge density, frequency of occurrence, and minimum distance for each of 10 land cover types) in R to avoid multicollinearity between land cover variables (Everitt 2004; R Core Team 2016). We selected the number of principal components (PCs) to retain >90% of the total variability of the original data of the 30 landscape variables, and we used landscape PCs to predict suitable beaver habitats across RSA.

Landscape scale habitat suitability modeling

We used presence locations to build habitat suitability models with MaxEnt methods (Elith et al. 2011; Phillips et al. 2006). To avoid pseudoreplication, we parameterized MaxEnt to remove duplicate presence points as well as points that were within 30 m of other presence points (Razgour et al. 2011). We set the number of random pseudo absence locations to 10,000. Training and test locations were randomly chosen by MaxEnt at a ratio of 80% to 20% of the total number of locations. The training data were used to tune MaxEnt parameters, and the testing was used to test the MaxEnt performance.

We used the PCs of landscape variables to build a MaxEnt model for generating a habitat suitability index map of the study area, a black box approach without inferring beaver resource selection. We referred to this model as a predictive model. We also built a MaxEnt model to infer landscape variables selected by beavers through the model selection of 30 original, untransformed landscape variables using a combination of the information-theoretic approach and least absolute shrinkage and selection operator (LASSO). LASSO is a regularization technique for regression, accounting for multicollinearity and avoiding model overfitting (Tibshirani 1996). The second model is referred to as inferential model.

We built the MaxEnt inferential models using the R package MaxentVariableSelection (MVS) (Jueterbock et al. 2016). The function VariableSelection within MVS uses corrected Akaike information criterion (AICc), LASSO, and the relative contribution of variables to model fit to select the most parsimonious model. VariableSelection selects a subset of landscape variables and a regularization multiplier (β) for LASSO, which minimize the AICc value (Jueterbock et al. 2016). LASSO uses the multiplier β to shrink the coefficient of least influential covariates toward zero to avoid overfitting (Merow et al. 2013; Jueterbock et al. 2016). The smaller the value of β, the closer the fit between the projected distribution and the training data set. We parameterized VariableSelection to remove a variable, which explained <5% of the model deviance (Jueterbock et al. 2016). We also parameterized VariableSelection to remove one of the two variables having the absolute Pearson correlation |r| > 0.7. The variable removed had lower contribution to model fit than the variable retained.

Beaver capture and GPS relocation data

We collected GPS location data from free-ranging beavers on RSA to validate habitat suitability maps developed with species presence locations. We captured beavers using Hancock live traps weighing 13.2 kg with dimensions of 71 × 91 × 10 cm in late March (Hancock Trap, Custer, SD, USA). We baited traps with castor scent and/or food lures (Backbreaker or Woodchipper, Dobbins’ Products, Goldsboro, NC, USA). We set traps at 1500 h and checked traps daily by 0900 h. We weighed each captured beaver to the nearest 0.1 kg using an electronic scale (Berkley, Inc., Spirit Lake, IA, USA) and estimated age by body mass (<1 year: ≤ 11 kg; yearling [1 year]: 11–16 kg; sub-adult [1–2 year]: 16–19 kg; and adult [> 2 year]: > 19 kg) (McTaggart 2002).

We attached UltraLITE GPS transmitters (model G10-210, Advanced Telemetry Systems [ATS], Inc., Isanti, MN, USA) integrated with VHF radio telemetry tags (model M3530, ATS, Inc.) to two beavers in each of 5 wetlands dispersed across RSA (n = 10). We located beavers with GPS units weekly during April and May by searching for active VHF beacons. We setup GPS units to record a location every 15 min over a 12-h period from 1800 h to 0600 h, based on an estimated 30-day battery life. We assumed that underwater locations and locations in lodges would not be recorded due to interference with satellites. We used the homing technique to recover transmitters with a VHF mortality signal and to locate transmitters after expected GPS battery failure. We extracted GPS data manually from each transmitter using ATS software Robin Manager version 2.5.14248 (ATS, Isanti, MN, USA). All beaver capture and handling methods were approved by the Institutional Animal Use and Care Committee of the United States Department of Agriculture, Animal and Plant Health Inspection Service, Wildlife Services, National Wildlife Research Center (Protocol # QA-2436).

Validation of MaxEnt habitat suitability models

MaxEnt models randomly separate the complete presence data into the training data for model building and test data for cross-validation. We used the area under the curve (AUC) of receiver operating characteristics (ROCs) with the 20% test data to cross-validate the MaxEnt models (Hilden 1991; Liu et al. 2011). The AUC value ranges from 0 to 1. An AUC value of 0.5 indicates the predictive performance of a random model, and 1.0 for a perfect predictive accuracy of the tested model (Liu et al. 2011).

In the cross-validation, the training and testing data may not be independent of each other. We also used GPS locations of free-ranging beavers to validate the MaxEnt models with the continuous Boyce index (CBI). The CBI is a Spearman correlation between the predicted-to-expected (P/E) ratio of habitat suitability value and mean Habitat Suitability Index (HSI), within a window moving over the HSI range in the predictive MaxEnt model (Boyce et al. 2002; Hirzel et al. 2006). The CBI value ranges from −1 to 1, with 0 being equivalent to predictions by random models and a negative value indicating an incorrect model (Hirzel et al. 2006). The CBI was evaluated with each of 5 sets of GPS locations of free-ranging beavers. Average CBI was computed over the 5 individual evaluations.

Although our presence locations were not randomly sampled as desired for habitat suitability or ecological niche modeling studies (Elith et al. 2011a; Renner et al. 2015), GPS locations from free-ranging beaver can be assumed to generate a spatial stochastic process (Cagnacci et al. 2010). We used the validation test with the GPS locations, independent of the training locations, to test the representativeness of our presence locations for beaver space use.

Fine-scale habitat selection

To evaluate fine-scale habitat selection, we built Resource Selection Functions (RSF) for the 5 GPS-tracked beavers using a use versus available approach. We estimated 100% minimum convex polygon (MCP) home ranges with the R package adehabitatHR (Calenge 2006). The GPS locations within 100% MCP home ranges represented habitat use. We randomly sampled the same number of locations without replacement within MCPs to represent available habitat (Boyce and McDonald 1999; Lele and Keim 2006). We then used GLMMs as RSFs to compare resources or habitat used to resources or habitat available and to determine landscape variables influencing fine-scale habitat selection by beaver with an animal identification number (ID) as random effects following Steyaert et al. (2015). We used landscape variables selected by landscape scale inferential MaxEnt models to build individual RSFs, with one of two highly correlated landscape variables (|r|> 0.7) being included in a RSF (Merow et al. 2013). We used AICc to select most parsimonious RSF for each of the 5 GPS tracked animals, with the lowest AICc representing the most parsimonious model (Burnham and Anderson 2002). A model of <2.0 ΔAICc is considered a competing model. We conducted model selection in a forward stepwise fashion. If the effect directions of landscape variables remained in the final GLMM were consistent with that revealed by the variable response curve produced by MaxEnt, we concluded that the effects of the variables were scale independent. We also fit generalized additive mixed models (GAMMs) to the GPS locations with animal ID as random effects to demonstrate nonlinear effects of the landscape variables selected by GLMMs and compared them to the trends of the MaxEnt response curves.

Results

We detected 334 presence locations (180 foraging locations, 12 feeding stations, 4 dams, 10 lodges, 19 cast mounds, and 109 trapping locations) in the 11 wetlands (Table 1). Of the 10 GPS transmitters attached to captured beaver, 7 were recovered after one month. Of the 7, 5 yielded data with a total of 607 GPS locations. Due to the denning and underwater habits of beavers, GPS relocations were irregular and did not represent 15-min intervals.

Table 1

Type and number of presence locations sampled from wetlands across RSA, Alabama, USA.

Wetland	Total presence per wetland	Main lodge (1= found, 0 = not found)	Number of secondary lodges	Dams	Castor mounds	Foraging locations	Feeding station	Captures
Blueberry (BB)	7	1	3	0	0	0	0	3
Easter Posey Wetland (EPW)	71	1	0	1	0	52	2	15
Thiokol Wetland (TW)	108	1	0	1	19	67	10	10
Martin Road North (MRN)	74	1	1	0	0	60	0	12
Patton Road Pond (PRP)	67	1	1a	3 (2)b	0	55a	0	7
Igloo Pond WetlaND (IPW)	17	0	0	1a	0	0	0	16
Mckinley Range (MK)	5	1	0	1a	0	0	0	3
Hudson Recreation Area (HUD)	55	1a	0	1a	0	44a	5a	4
Hale Road Wetland (HRW)	97	1	2a	3 (1)b	12a	46a	12a	21
Nasa Gun Range (NGR)	4	1	0	0	0	0	0	3
DDT Impoundment Site (DDT)c	39	1	2	2a	1a	16a	0	15

Notes:

location not included in Maxent model.

numbers in parentheses were not used in Maxent model.

2 main lodges.

Our predictive model included the first 15 PCs, which accounted for 90% of total variability in the original landscape variables. The predictive model had excellent accuracy with an AUC value of 0.97. The CBI was 0.97 in the validation test using the 607 GPS locations of 5 beavers. Average CBI was 0.84 (standard deviation [SD] =0.03) over 5 separate tests using the GPS locations of individual beavers. Therefore, validation tests using the GPS locations independent of the training data indicated an excellent predictive power for the model. The suitable habitat of beaver was highly discontinuous and fragmented (Figure 1).

The results of MVS, across 10 MaxEnt models showed that a regularization parameter (β multiplier) of 2.0 had the lowest AICc. Thus, we set the β multiplier to 2.0 for subsequent model selection. Model selection for the best inferential MaxEnt model showed that open water edge density, shrub edge density, woody wetland edge density, and grassland frequency influenced habitat selection by beavers. Variable response curves demonstrated nonlinear effects of the 4 landscape variables on the habitat suitability of beavers (Figure 2).

Figure 2

Response curves of predicted occurrence likelihood of grassland frequency, water edge density, shrub edge density, and woody wetland edge density in the Maxent model for the habitat suitability of American beaver. Each curve represents how the predicted likelihood of habitat suitability changes with increased value of a landscape variables while the other landscape variables are held at averages.

GLMMs showed that water edge density, shrub edge density, and woody wetland edge density affect fine-scale habitat selection of beaver (Table 2). GAMMs demonstrated the nonlinear effects of woody wetland edge density (Figure 3A) and shrub edge density (Figure 3B) but linear effects of open water edge density (Figure 3C) on the habitat suitability of beaver at the fine scale.

Table 2

Model selection using Akaike information criteria for the effects of woody wetland edge density (wwetbd), shrub edge density (shrubbd), water edge density (waterbd), and grassland frequency (grassfq) on the habitat suitability of American beavers in RSA, Northern Alabama, USA. The top 3 of 15 models were presented.

Model	AICc	ΔAICc
y∼wwetbd+shrubbd+waterbd	1422.61	0
y∼waterbd+shrubbd+grassfq+wwetbd	1424.14	1.53
y∼wwetbd+shrubbd	1432.84	10.23

Figure 3

Results of the GAMM examining the effects of (A) woody wetland edge density, (B) shrub edge density, and (C) water edge density on fine-scale habitat selection by American beavers in RSA, Alabama, USA. Woody wetland edge density and shrub edge density exhibit nonlinear responses, whereas water edge density shows an inverse linear response.

Discussion

In this study, we found evidence for scale independent effects of landscape structure on beaver habitat selection by applying a multi-scale approach with no hierarchical constraints. First, the fine-scale spatial distribution of beaver relocations was correlated with mean suitability indices predicted by landscape-scale MaxEnt models (Speakman correlation index or CBI = 0.87). Second, 3 of 4 landscape variables influencing landscape-scale habitat selection had similar effects on the fine-scale habitat selection by beavers (Figures 2 and 3). Of the 3 variables that were congruent between scales, edge densities of woody wetland and shrub were likely critical to forage and lodge construction opportunities for beavers. Scale-dependent habitat selection by beavers dictates that habitat selection at broader scales constrains selection at finer scales (Mayor et al. 2009). By removing hierarchical constraints, habitat selection may be determined by the characteristic scales of resources, such as spatial extents of woody wetland in the case of beavers, or predation risks (Mayor et al. 2009). Although we found congruence in habitat characteristics among spatial scales, our unrestricted approach allowed us to detect and contrast the direction of relationships between scales (Figures 2 and 3).

At the fine scale (i.e., within home ranges), we saw nonlinear relationships with shrub edge densities, suggesting that beavers found intermediate thresholds to meet forage requirements and may have used ecotones between aquatic and terrestrial habitats (Donkor and Fryxell 1999; John and Kostkan 2009). Although beavers move out of water up to 60 m to forage (Donkor and Fryxell 1999), this was unlikely on RSA given the large expanses of wetlands, consistent with negative relationships between habitat selection and distance to open water. With long regeneration time of some woody plants, beaver herbivory may reduce food availability and increase woody wetland edge density, so beaver may include more woody vegetation within their home ranges (Campbell et al. 2005; Brzyski and Schulte 2009; McClintic et al. 2014a). Consequently, beavers may increase the use of the area having high woody wetland edge density (Figure 3A). Although we did not quantify plant vigor in this study, these patterns are consistent with the resource heterogeneity hypothesis previously observed on RSA (McClintic et al. 2014a).

The negative linear relationship with herbaceous wetland edge density at the home range scale suggests that while beavers need open water for protection from predators, they are able to trade-off the size of wetlands, or distance of lodge to bank, for easier access to woody plants (e.g., food and construction resources). This is consistent with central foraging theory and the deliberate movement hypothesis (McClintic et al. 2014b). Beavers also defend their territories from others beavers, often marking their territories along the same aquatic-terrestrial edges where they forage on woody plants. Beaver selection for decreased herbaceous wetland edge density also is consistent with maintaining an “economically defendable” area (Campbell et al. 2005). Thus, minimizing predator risk and travel requirements for food, construction material, and defense all promote beaver fitness. Alternatively, foraging locations used in habitat selection modeling were primarily tree biting marks, which represented winter and early spring habitat use before herbaceous wetland plants remerged.

While our results partially supported our prediction (P1), herbaceous wetlands did not appear to be a major factor influencing beaver habitat selection on RSA. However, this is understandable as aquatic herbaceous plants also occurred in woody wetlands where beavers were found. The amount and timing of aquatic vegetation in beaver diets varies with latitude (Svendsen 1980; Parker et al. 2007; Milligan and Humphries 2010). In the southeastern USA, beavers tend to shift their diet toward non-woody plants in summer months (Parker et al. 2007). Although we did not quantify annual food consumption rates of beaver at our study sites, beavers have been shown to alter landscapes through consumption of large volumes of aquatic (Parker et al. 2007) and woody plants (Johnston and Naiman 1990). We likely missed much of the herbivory on aquatic vegetation during our sampling period. Nevertheless, beavers have been shown to feed on a variety of terrestrial, emergent, and aquatic non-woody vegetation. An exhaustive list of all forage species is unavailable, and would vary between regions (Gallant et al. 2004). About 59% of 180 feeding signs observed in our study were from deciduous trees Acer rubrum, Liquidambar styraciflua, and Nyssa spp. However, only 0.07% was from Ligustrum sinense and Cephalanthus occidentalis, the only shrub species we documented upon which beaver foraged on RSA. Our habitat selection models supported the avoidance of both coniferous and mixed pine-hardwood landscapes, which is consistent with other studies of habitat selection in beaver (Roberts 1981; Gallant et al. 2004).

In summary, multi-scale habitat selection has become a theoretic foundation for understanding animal habitat selection or resource use (Boyce 2006; Johnson 1980). To our knowledge, this is the first known study to evaluate beaver habitat selection using multi-scale modeling techniques. Future research is needed to test these predictions across the range of American beavers and to apply statistical techniques that optimize scale, which is critical to assess scale dependence (McGarigal et al. 2016).