Why invasive Patagonian beavers thrive in unlikely habitats: a demographic perspective

Abstract

Understanding the demography of an invasive species is crucial to better guide managers seeking to slow the spread of the invader. Habitat differences can affect demographic rates, which may in turn impact the speed of the invasion, but this has been rarely addressed. We studied the demography of invasive North American beavers (Castor canadensis) in 2 contrasting habitat types of the island, forest and steppe, on Tierra del Fuego in southern Patagonia. We used repeated observations, mark-resight methods, telemetry, and camera traps to estimate colony size and demographic rates of beavers in the 2 habitats. Colony size and the number of offspring (“kits”) produced per colony per year were higher in the steppe, contrary to the belief that forest is better habitat. This may be the result of the longer time since invasion in the forests of Tierra del Fuego and that the forest subpopulation is showing density-dependent regulation. Survival of beavers was high in all age classes and was higher than survival rates recorded in North America. Our work shows that plasticity of habitat use and predator release have likely facilitated beaver invasion in Patagonia. The higher productivity and detectability of beavers in the steppe call for active management in a habitat previously assumed to be subprime.

Entender la demografia de las especies invasoras es crucial para guiar a los gestores que intentan frenar la expansión de los invasores. Las diferencias en el tipo de hábitat pueden afectar las tasas demográficas, lo que a su vez puede impactar en la velocidad de la invasión, pero esto ha sido raramente investigado. Estudiamos la demografía de castores norteamericanos invasores (Castor canadensis) en dos hábitats contrastantes, bosque y estepa, en la isla de Tierra del Fuego en el sur de Patagonia. Utilizamos observaciones repetidas, métodos de marcado-reavistaje, telemetría y cámaras trampa para estimar el tamaño de colonia y las tasas demográficas de castores en los dos tipos de hábitat. El tamaño de colonia y el número de crías producido por colonia por año fue mayor en la estepa, contrario a la creencia de que el bosque es un mejor hábitat. Creemos que este puede ser el resultado de un mayor tiempo de invasión en los bosques de Tierra del Fuego y que la subpoblación del bosque está mostrando regulación denso-dependiente. La supervivencia de castores fue alta en todas las clases de edad y fue más alta que las tasas de supervivencia registradas en Norteamérica. Nuestro trabajo muestra que la plasticidad en el uso de hábitat y la ausencia de depredadores han probablemente facilitado la invasion de castores en Patagonia. La alta productividad y detectabilidad de castores en la estepa llaman a un control activo en un hábitat que antes había sido considerado subóptimo.

biological invasions, Castor canadensis, density dependence, habitat heterogeneity, Patagonia, predator release

All habitats are heterogeneous at some scale, and organisms respond to that heterogeneity. Habitat differences can affect demographic rates, which will in turn determine distribution and abundance of a species across a heterogeneous landscape. Although studies on environment-specific demography are widespread for plants, such studies focused on animals, and mammals in particular, have been comparatively scarce, perhaps due to difficulties posed by detection and animal movement that are only now being addressed (Royle et al. 2013). For biological invasions, habitat heterogeneity can affect population growth rates and dispersal at the invasion front, impacting the rate of spread (Shigesada et al. 1986; Dewhirst and Lutscher 2009). Thus, understanding the effect of habitat heterogeneity on demography is particularly important in predicting the spread of biological invasions. In this paper, we quantify differences in colony sizes and demographic rates of introduced North American beavers (Castor canadensis) in 2 very different habitats in Patagonia: forest and steppe.

In addition to relatively fixed differences between habitats, changes to a habitat after it has been invaded may impact demographic rates, affecting both the growth of local populations and the rate at which an invasion proceeds. In models with spatially homogeneous dispersal (Skellam 1951; Kot et al. 1996), negative density dependence behind the invasion front (which may result from habitat degradation) does not influence the speed of an invasion. However, other models have shown that density dependence can generate fluctuating rates of spread by increasing dispersal distance in areas where resources have been depleted (Dwyer and Morris 2006). When assessing effects of different habitats on invading species, we must acknowledge that differences between habitats in the time since invasion may be confounded with inherent differences in the quality of those habitats before they were invaded.

Beavers were introduced to the island of Tierra del Fuego in 1946 to “enrich” the native fauna and foster a fur trade and have since spread throughout the archipelago. More recently, they began spreading in continental Patagonia (Graells et al. 2015). In North America, beavers inhabit a wide range of habitats from the subarctic to the Rio Grande, and by removing trees and damming streams they create new habitat and increase species richness at the landscape scale (Wright et al. 2002; Rosell et al. 2005; Cunningham et al. 2007). But in southern Patagonia, beavers negatively impact recruitment of Nothofagus spp., the dominant trees in riparian forests (Anderson et al. 2006; Martínez Pastur et al. 2006); create entirely new habitats that allow establishment of other invasive species (Anderson et al. 2006; Henn et al. 2014); and affect food webs of streams (Anderson and Rosemond 2007, 2010). Beavers first established in the forests in which they were initially introduced but by the 1990s began to establish in the adjacent steppe (Skewes et al. 2006).

Impacts of beavers on subantarctic forests have been well investigated, but no detailed studies on demography of invasive beavers in any habitats in Patagonia have been conducted to date. Demographic information is essential to better understand and manage the ongoing beaver invasion. For instance, predator release has been suggested as a major driver of the invasion (Wallem et al. 2007); here, we test this hypothesis by comparing survival rates of beavers in Patagonia with those of beavers in North America. If predation regulates beaver populations in their native range, this comparison may give us information on the level of culling needed to control further spread of the invasion. Beavers also are known to adjust fecundity in response to habitat quality and increased density. Gunson (1970) found that beavers in low-quality habitats show high frequency of resorption of embryos resulting in fewer offspring (“kits”) per female beaver. Bergerud and Miller (1977) and Payne (1984b) showed that mean litter size decreased as density increased in Canada, but high densities also can cause delayed reproduction of juveniles due to decreasing availability of territories (Bergerud and Miller 1977; Busher and Lyons 1999). Differences in fecundity between habitats will likely affect the spread of the invasion in heterogeneous landscapes.

In Patagonia, successful establishment of beavers in the semiarid steppe has challenged previous studies assuming that steppe was unsuitable habitat. Lizarralde et al. (2004) and Wallem et al. (2007) suggested that low density of beaver colonies in the steppe could reflect poor habitat quality, although animal density can be a misleading indicator of habitat quality (Van Horne 1983). Skewes et al. (2006) reported lower densities of colonies and assumed smaller colony sizes in the steppe compared to the forest when estimating total abundance of beavers in Tierra del Fuego. If true, low population growth rates in the steppe could slow the spread of the invasion, but the few data reported on the rate of spread suggest that it is at least as fast and perhaps faster in the steppe (Skewes et al. 2006), which could be due to differences in demography, movement, or both.

Heterogeneity within major habitat types also can affect demographic rates. At the scale of individual beaver colonies, research in northern hemisphere forests has shown an association between food availability and colony size and recruitment. For example, Fryxell (2001) found larger colony sizes and more kits per colony in areas with high woody cover and aquatic vegetation, but the relationship between vegetation and beaver demography has never been assessed in Patagonia. Within-habitat differences in vegetation variables can improve our understanding of what drives demographic differences between forest and steppe. Conversely, absence of differences within or between habitats can be a sign of invader plasticity, which has been cited as a major predictor of invasion success in mammals (González-Suárez et al. 2015).

The goal of this paper is to compare demography of invasive beavers within and between habitats in Patagonia, and between Patagonia and North America. To do this, we estimated age-specific survival, using mark-resight methods and telemetry; the proportion of breeding pairs producing kits, using camera traps; the mean number of kits produced per breeding pair; and the size and trends in size of beaver colonies in the forest and steppe using repeated observations. We used these data to address 4 questions:

1. Do colony sizes and demographic rates currently differ between forest and steppe habitats?

2. Can we identify features within these habitats that are associated with demographic variation?

3. How might the longer history of beaver presence in forests influence the current colony sizes and demographic rates?

4. Do demographic rates of beavers differ between Tierra del Fuego and North America?

Materials and Methods

Study sites.

We defined 4 study sites on the main island of Tierra del Fuego, 2 in the forest and 2 in the steppe (Fig. 1). Sites in the forest were separated by at least 150 km from sites in the steppe. Forest is dominated by 3 tree species in the genus Nothofagus and has an understory with low species richness. Annual precipitation in the forest sites ranges from 500 to 600 mm, and snow covers the ground typically between May and September. High cloud cover and proximity to the sea generate conditions of high relative humidity throughout the year and temperatures that vary between 4°C to 15°C in summer and −3°C to 6°C in winter (Martínez Pastur et al. 2006). Beavers have been present in our forest sites for at least 50 years and have impacted between 20% and 40% of the stream length (Anderson et al. 2014) with densities among the highest reported in the literature (> 4 colonies/km²—Lizarralde 1993; 1.91 colonies/km of stream length—Skewes et al. 2006).

Fig. 1.

Location of the study sites in the forest (dark gray) and the steppe (light gray) in Tierra del Fuego, Argentina.

Sites in the steppe vary in cover of the dominant woody shrub mata negra (Chilliotrichum diffusum; range: 0–54%, mean: 28.7), used by beavers for food and dam construction. Annual precipitation ranges from 200 to 300 mm and snow rarely accumulates in winter. Temperatures are similar to those of the forest sites, but strong winds and lower precipitation increase evapotranspiration and favor the growth of vegetation adapted to xeric conditions. First reports of beaver colonies in our study sites in the steppe were in the early 1990s (Lizarralde 1993) and densities reported for the steppe are very low (< 0.15 colonies/km of stream length—Skewes et al. 2006).

Field sampling.

From 2011 to 2015, we studied 25 beaver colonies in each habitat type. At each site, we observed a minimum of 10 colonies, except for 1 of our sites in the steppe where we only had 5 colonies. We incorporated this last site in the steppe at it had slightly different features than the other (lower precipitation), but repeated counts (see next section) did not significantly differ between sites. A beaver colony usually consists of a breeding pair, newborn kits, and juveniles born the previous year (Bradt 1938). Beaver kits are born around October in the southern hemisphere and juveniles may disperse after they are a year old to establish new colonies.

We conducted 2–3 repeated observations of each colony between mid-January and mid-April of each year between 2012 and 2014. Observations were performed with binoculars or monocular by a trained observer on nonrainy days for 2 h before dusk, when beavers are active. Observation spots were selected to have good visibility of the focal lodge. In our study sites, hunting is rare and beavers are not wary and thus are easy to observe. We recorded the total number of individuals at each colony and the number of kits. Kits are easy to distinguish by size and behavior from juveniles and adults.

Statistical models to estimate colony size and number of kits.

We used N-mixture models to assess differences in colony size and number of kits per colony per year between and within habitats (Royle 2004; Kery et al. 2005). N-mixture models allow one to estimate abundance in spatially replicated populations using temporally repeated counts.

N-mixture models are hierarchical Poisson regression models where abundance is modeled as a Poisson variable, but they attach a logistic regression model to account for imperfect detection (Kéry and Schaub 2012). Generally, the model consists of 2 parts, one to model the ecological state (abundance) and another to model the observation process (counts). In particular, we assume that

$$\begin{array}{cc}{N}_i~Poisson(\lambda )& \text{Ecological process}\end{array}$$

1

$$\begin{array}{cc}{y}_{ij}|N_i~Binomial(N_i,p)& \text{Observation process}\end{array}$$

2

where N_i is the true local abundance at colony i, drawn from a Poisson distribution with mean λ, y_ij is the number of individuals observed at site i at time j (counts), and p is the probability of detection for each individual. In our case, we used colony size as a measure of local abundance.

To account for temporal trends in colony size for each habitat, we simplified the model formulated by Kery et al. (2009). This is an open N-mixture model, in which changes in colony size may occur between years but we assume no changes within years so that counts within years provide information on the detection probability. We modeled λ and p as a function of covariates such that:

$$\log ({\lambda }_{ik})=a_{habita{t}_i}+r_{habita{t}_i}\times (k-1)$$

3

$$logit(p_i)=c_{habita{t}_i}$$

4

Equation 3 indicates that colony size can vary by habitat at site i and year k. $a_{habita{t}_i}$ is an intercept that represents the mean log colony size at year 1 in the type of habitat at site i, habitat i, and $r_{habita{t}_i}$ represents the constant annual population growth rate in that habitat type. Equation 4 assumes detection probability p_i at colony i depends only on the habitat at colony i, as we did not have enough data to estimate a different detectability for each year. Specifically, we assumed that $p_i=c_{habita{t}_i},$ a constant probability of detection given the type of habitat present at site i. We applied this log-linear Poisson model to the colony size data to estimate habitat-specific average colony sizes, growth rates, and detectabilities.

To study how vegetation variables affected colony size within habitats, we ran separate single-season N-mixture models for each year and habitat. For our models, we chose to incorporate vegetation variables to predict colony size and we assumed a constant probability of detection p for each habitat. Assuming constant detectability was preferred as these models may be hard to fit with small sample sizes, particularly if the number of predictors is high, and we were more interested in determinants of colony size than of detectability. Models for the forest had the understory cover, its height, and the median tree diameter as potential predictors of colony size. Models for the steppe included shrub cover and height of the shrub cover as potential predictors.

To model the number of kits produced per colony per year, we again used N-mixture models but assumed that the true number of kits varied among colonies according to a zero-inflated Poisson distribution. Under this model, colonies may have no kits as a consequence of both a Bernoulli process and a Poisson process. The hierarchical model has the structure

$$\begin{array}{l}\begin{array}{cc}\text{Level 1} (\text{kits produced}, \text{colony} i, \text{year} k)& \end{array}\\ z_{i,k}~Bernoulli(\Omega )\end{array}$$

5

$$\begin{array}{l}\begin{array}{cc}\text{Level 2 (true number of kits}, \text{colony} i, \text{year} k\text{)}& \end{array}\\ K_{i,k} ~ Poisson(z_{i,k}{\kappa }_k)\end{array}$$

6

$$\begin{array}{l}\text{Level 3} (\text{observed number of kits}, \text{colony} i, \text{year} k)\\ x_{i,k} ~ Binomial(K_{i,k}, p_i)\end{array}$$

7

where Ω is the zero inflation parameter (assumed to be constant for all colonies and years in a given habitat), $z_{i,k}$ is a random variable that can be 0 or 1 for colony i in year k (if $z_{i,k}=0,$ no kits are produced, but if $z_{i,k}=1,$, kits still may not be produced, based on the outcome of the Poisson process), ${\kappa }_k$ is the Poisson mean number of kits per colony in year k, $K_{i,k}$ is the true number of kits in colony i in year k (and includes the probability that a colony produces 0 kits), $p_i$ is the probability of detection of each kit in colony i (which only depends on the habitat type in colony i), and $x_{i,k}$ is the observed number of kits in colony i in year k. According to this zero-inflated Poisson model, the probability that colony i produces no kits in year k is $(\text{1}-\Omega )+\Omega \exp \text{(}-{\kappa }_k).$ As before, we estimated each of the parameters using the appropriate link function. To model variation in the number of kits within habitats, we took the same approach that we followed with colony sizes, so that only number of kits—not detectability—was modeled as a function of vegetation covariates.

To independently estimate the proportion of colonies producing kits, we used camera traps in the last year of our study. To do this, we selected 22 colonies (11 in each habitat) and placed camera traps for 3 consecutive nights in each colony. The number of colonies was limited by the presence of good sites to place camera traps, which were selected after watching where beaver families spent more time in our repeated observations. We considered each night to be an observation, and we assigned a 1 to nights when kits were observed and 0 otherwise. We used occupancy models to estimate the proportion of colonies with kits in each habitat. Details are provided in Supplementary Data SD1, and details on parameter estimation of all our models are provided in Supplementary Data SD2.

At each colony, we also quantified abundance of food and dam-building material by placing between 4 and 6 equally spaced 20-m transects perpendicular to the shoreline on the pond where the main lodge was located. Aquatic vegetation is scarce in streams in Tierra del Fuego and was therefore not considered in our sampling. We measured woody understory cover (in forest) and shrub cover (in steppe) by dividing the total length of the vegetation intercepted by the line transect by the transect length. We also estimated the mean height of the vegetation intercepted by the transect.

We measured forest structure at our forest sites as follows. Centered at 7.5, 17.5, 27.5, 37.5, and 47.5 m along a 50-m transect, we located 5 × 10 m quadrats (with the 5-m sides parallel to the transect). Beavers mostly forage within a 50-m distance from the pond (Fryxell 2001). At each quadrat, we recorded the number of live trees higher than 150 cm, diameter of these trees at 20 cm above the ground, and number and diameter of dead trees cut by beavers.

Estimation of survival using mark-resight and telemetry.

Every year between December and April, we livetrapped beavers in the same beaver colonies we observed, using snares with a relaxing lock (McKinstry and Anderson 1998). After trapping, we waited at least 3 days before resuming repeated observations and camera trapping to allow the colony to resume normal activity. Snares were set the day before and were checked early in the morning to minimize the time beavers were held in the traps.

Captured beavers were manually restrained by 2 people without using anesthesia, and we covered their eyes to relieve stress. We placed uniquely colored ear tags (Dalton ID) on captured beavers that could be seen with binoculars when observing colonies. We wiped each beaver’s ears with a diluted chlorhexidine solution prior to piercing to prevent infections. We used the repeated observations we detailed in the previous section to gather information on survival of beavers through mark-resighting.

Both to estimate survival and to track their movements (not reported here), we attached beaver-tail transmitters with a mortality sensor (Advanced Telemetry Systems, Isanti, Minnesota) to juvenile and adult beavers (Arjo et al. 2008). Transmitters weighed 35 g and had a battery life of 500 days. Transmitters were fitted quickly by drilling a 5-mm hole with a cordless drill in the beaver’s tail, offset from the center line where the tail is composed mainly of fatty tissue. We applied a solution of sulfadiazine to the hole to avoid potential infections and facilitate healing. We determined sex of each beaver by palpation of the baculum (Osborn 1955), which was later corroborated using molecular methods with hair samples (Goldberg et al. 2011). We weighed beavers with a spring scale and released them at the site where capture occurred. All capture and handling protocols were approved by the Duke University Animal and Care Committee, the wildlife office of the province of Tierra del Fuego, Argentina and followed ASM guidelines for the use of wild mammals in research (Sikes et al. 2016).

We checked beaver survival using telemetry every 2 weeks between November and April in each habitat. When we detected a mortality signal, we located the radiotransmitter to confirm death or tag loss. In addition to this, we flew over the study area in a small airplane equipped with a telemetry antenna at the end of each field season to relocate animals that may have dispersed.

Estimation of survival of kits.

Because only individuals heavier than 7 kg can be radiotagged, we separately estimated survival of kits and survival of juveniles and adults. We used weight to define age classes and followed the criteria established by Feldman (2015), such that individuals that weighed more than 5 kg at the time of capture were considered juveniles older than 1 year. Survival of kits was therefore assessed only through mark-resight data, using a Cormack-Jolly-Seber model with a robust design (Gimenez et al. 2007). As the total number of kits captured was small (n = 18), we pooled kits from forest (n = 5) and steppe (n = 13) and estimated yearly survival and detection assuming they did not differ by habitat.

Estimation of survival of juveniles and adults.

We captured 91 adults and juveniles (forest: n = 30, steppe: n = 61) and radiotagged 51 of them. As both juveniles and adults could have had either ear tags only or both ear tags and radio tags, we fitted a multistate capture-recapture model to estimate monthly survival (Lebreton et al. 2009). We transformed monthly survival to yearly survival and approximated the variance using the delta approximation (Powell 2007). Our data set was built by compiling individual encounter histories during 52 months of study. We defined 5 observation states for individual beavers: alive with radio tag, alive without radio tag (ear tags only), dead with radio tag, dead without radio tag, and not observed. Details on the construction of our model are given in Supplementary Data SD2. More generally, the model included estimates of survival using telemetry, survival using mark-resight data, the probability of transitioning from the telemetry group to the mark-resight group (i.e., transmitter stopped working or was lost), and habitat specific probabilities of detection. Because of limited sample size, we assumed survival was the same in forest and steppe, a reasonable assumption as adults and juveniles do not have known predators in either of these habitats in Tierra del Fuego.

Results

Colony size was larger in the steppe than in the forest in all 3 years that repeated observations were conducted (mean colony size ± SD: forest, 4.89 ± 0.49; steppe, 7.06 ± 0.91; Fig. 2). Colony size remained relatively constant over years in the forest (posterior mean growth rate r_forest, 95% credible interval [CRI]: −0.04, −0.22–0.14), while there was a slight increase over time in colony size in the steppe (posterior mean growth rate r_steppe, 95% CRI: 0.12, −0.01–0.27). Probability of detection was higher in the steppe compared to the forest (posterior mean, 95% CRI: forest, 0.49, 0.29–0.64; steppe, 0.57, 0.40–0.70).

Fig. 2.

Mean yearly colony size of beavers (± posterior SD; Castor canadensis) for the forest (solid line) and the steppe (dashed line) sites in Tierra del Fuego, 2012–2014.

Within-habitat heterogeneity impacted colony sizes differently in the forest and the steppe. In the forest, only the median tree diameter was consistently and negatively related with colony size across years, although the effect was close to 0 the second year (Table 1, Supplementary Data SD1). Consistent with this result, cut trees were smaller in diameter than uncut trees. In the steppe, none of the vegetation variables showed a significant and consistent effect in our study (Table 2, Supplementary Data SD1). Interestingly, we recorded some of the largest colony sizes in the steppe in areas with little or no shrub cover, indicating that beavers are able to feed on, and even construct dams from, herbaceous vegetation alone.

We also found a trend toward a higher number of kits per colony in the steppe compared to the forest (Fig. 3). The probability of detection of kits in the steppe was much higher than in the forest (posterior mean, 95% CRI: steppe, 0.63, 0.29–1; forest, 0.19, 0.04–0.49). Zero-inflated Poisson models of the number of kits fitted for each year resulted in very imprecise coefficient estimates, probably because of the combination of small sample size and a high number of covariates (Tables 3 and 4, Supplementary Data SD1). Thus, we were unable to assess the effect of within-habitat heterogeneity on the number of kits. Overall, the mean litter size was similar to values reported by studies in the native range, with the number of kits in the forest falling in the lower end of the range (Table 1).

Fig. 3.

Mean yearly number of kits (± posterior SD; Castor canadensis) for the forest (solid line) and the steppe (dashed line) sites in Tierra del Fuego, 2012–2014.

Table 1.

Mean litter sizes of North American beavers (Castor canadensis) reported in scientific literature.

Mean litter size	Location	Reference
1.7–2.7a	Newfoundland (Canada)	Bergerud and Miller (1977)
2.7	Ohio (United States)	Svendsen (1980)
1.5–2.7a	Newfoundland (Canada)	Payne (1984b)
3.5	Illinois (United States)	Bloomquist and Nielsen (2010)
2.5	Illinois (United States)	Havens et al. (2013)
1.7 (forest)–2.7 (steppe)	Tierra del Fuego (Argentina)	This study

^aRange reported over 3 years of study.

Camera traps helped us to assess with greater confidence whether the proportion of colonies with kits differed between habitats. Given that we had a small data set and could not fit covariates for both the proportion of colonies with kits and the probability of detection, we determined which model to fit as follows. We compared the fraction of colonies where kits were recorded at least once with the raw probability of detection, which was calculated as the number of nights a kit was recorded by our camera traps divided by the total number of camera trap nights. We found that the raw proportion of colonies producing kits did not vary as much between habitats (forest, 0.72; steppe, 0.63) as the raw probability of detection (forest, 0.5; steppe, 0.76). Thus, we decided to fit an occupancy model where detection, but not the proportion of colonies with kits, varied by habitat. Under this model, the corrected proportion of colonies with kits was 0.76 (95% CRI: 0.68–1), whereas estimation of the proportion of colonies producing kits using estimates of the zero-inflated Poisson model was close to 0.72 . Detection (now of kits, not of all individuals as above) was again higher in the steppe (posterior mean, 95% CRI: forest, 0.43, 0.24–0.65; steppe, 0.71, 0.42–0.9).

Estimates of survival of kits were imprecise given the small sample size (mean survival, 95% CRI: 0.66, 0.39–0.94). An interesting detail of our data is that mortality sensors in the radio tags allowed us to separate mortality from permanent emigration. Of our radiotagged beavers, only 5 died (2 in the forest and 3 in the steppe). The yearly survival estimate derived from telemetry was very high (posterior mean, SD: 0.97, 0.02) and apparent survival estimated from mark-resight data (which includes permanent emigration) was much lower (posterior mean, SD: 0.68, 0.06), suggesting permanent emigration was common, as one would expect for juveniles and young adults. Both survival of kits and the pooled survival rate of adults and juveniles were much higher than survival rates estimated from radiotelemetry data in North America (Table 2).

Table 2.

Mean survival rates of juvenile and adult beavers (Castor canadensis) reported in scientific literature, from studies using telemetry.

Juvenile survival	Adult survival	Location	Reference
0.29–0.59a,b	0.29–0.59	Illinois (United States)	Havens et al. (2013)
0.36	0.66	Wyoming (United States)	McKinstry and Anderson (2002)
0.43		Illinois (United States)	McNew and Woolf (2005)
0.67		Wyoming (United States)	VanDeelen and Pletscher (1996)
0.55–0.76c	0.76–0.87	Illinois (United States)	Bloomquist and Nielsen (2010)
0.82	0.88	Massachusetts (United States)	DeStefano et al. (2006)
0.89		Illinois (United States)	Cleere (2005)
0.97a	0.97	Tierra del Fuego (Argentina)	This study

^aJuveniles and adults were pooled to obtain a single survival estimate.

^bEstimates for the 1st and 2nd year of the study.

^cSurvival was estimated for females and males separately.

Discussion

In our study, we used several robust statistical methods to best estimate colony sizes and demographic rates of invasive beavers in different habitats in Patagonia. Colony size and the number of kits produced per colony were higher in the steppe compared to the forest, but we were unable to detect differences in all other demographic rates between habitats. Survival of juvenile and adult beavers was very high relative to the native range, suggesting predator release may be a major driver of the invasion. In the following paragraphs, we discuss the main reasons for the success of the invasion, examine the evidence for density dependence, and make recommendations based on our results for how the beaver invasion in Patagonia should be managed.

Several factors explain why beavers have been so successful in Patagonia. First, beavers are very plastic in the type of habitats they can use. The larger colony size and greater number of kits produced, along with a trend toward increasing colony sizes in the steppe, highlights not only that beaver populations can thrive in the steppe but that steppe can play an important role in the spread of the invasion. Trends in the number of kits are similar to those found for colony sizes and we believe that larger colony sizes in the steppe are at least partially related to greater kit production there. Although Fig. 2 is not conclusive regarding differences in the number of kits between forest and steppe, 2 additional pieces of evidence confirm the trend. First, the proportion of kits captured in the steppe was almost 30% higher than in the forest. Second, using camera traps, detection of kits along highly conspicuous beaver trails was also higher in the steppe, and variation in detection probability is often caused at least in part by variation in abundance (Royle and Nichols 2003).

The second reason for the beaver’s invasion success is that, as in North America, they have the ability to utilize the best locations within habitat types, even when they can thrive under a range of conditions. Median tree diameter was related (negatively) with colony size in the forest. Beavers usually selected smaller trees in the forest, as cutting down large trees represents a cost of both greater predation risk and higher processing time (Fryxell and Doucet 1993). Although we did not find any consistent predictors of colony size in the steppe, colony size and the number of kits produced per colony were large even in areas with no woody vegetation cover, which suggests that beavers are extremely adaptable and can rely on herbaceous vegetation only. Other studies showed that the availability of woody vegetation cover does not necessarily correlate positively with site occupancy (Beier and Barrett 1987; Hartman 1996; Suzuki and McComb 1998). Phenotypic plasticity and, more specifically, habitat breadth have been cited as main predictors of invasion success in mammals (González-Suárez et al. 2015). For invasive beavers in Patagonia, the ability to utilize diverse habitats has been an underappreciated contributor to invasion success.

The third reason for the beaver’s invasion success is likely to be release from predators, and perhaps diseases as well. Survival of juveniles and adults was among the highest reported in the literature (Table 2). We argue this is probably a consequence of predator release. Pumas (Puma concolor) could prey upon beavers in continental Patagonia but they are absent on Tierra del Fuego. Another factor that often impacts beavers in their native habitat is the bacterial disease tularemia, but the only study on diseases of beavers in Patagonia did not find tularemia (Skewes et al. 1999). Two out of the 5 deaths we reported in our study were the result of trapping of nuisance beavers or shooting, which is occasional in our study sites and has remained low in most of Tierra del Fuego. In areas of their native range, wolves and coyotes seem to be the main predators of beavers (Muller-Schwarze 2011), but they are absent on Tierra del Fuego.

Particularly for juveniles, dispersal is associated with an increase in mortality (DeStefano et al. 2006; Muller-Schwarze 2011). Mortality can be caused by predation, human-beaver encounters, or intraspecific encounters, among other causes. Human-beaver encounters do occur on Tierra del Fuego, but probably at much lower rates than in populated areas in North America (DeStefano et al. 2006). We only found 1 beaver, at 1 of our steppe sites, which presented injuries that may be ascribed to a fight with a conspecific. We did however find that survival was lower for kits than for juveniles and adults, as other studies report. Payne (1984a) found as few as 48% of beaver kits survive to the first 6 months, and Bloomquist and Nielsen (2010) reported survival rates of 28% for the first 11 months. Although imprecise given the small size of our sample, estimates of kit survival also were higher in our study than these North American estimates. Beaver kits that venture far from the lodge could be preyed upon by foxes (Lycalopex spp.) or raptors in Tierra del Fuego, even though these predators are unlikely to be able to kill the larger juvenile and adult beavers.

Our study adds further evidence in favor of the predator release hypothesis in vertebrates. Although frequently cited as a driver of invasion success, the predator release hypothesis has received only partial support, overwhelmingly biased toward studies of plants (Keane and Crawley 2002; Colautti et al. 2004; Liu and Stiling 2006). In a recent review, Jeschke et al. (2012) revealed that 5 out of 8 studies on vertebrates but 54 out of 106 studies on all species combined show some support for the predator release hypothesis, suggesting that introduced vertebrates are in fact more likely to experience predator release.

Classic invasion theory predicts that successful invaders should exhibit early reproduction and high fecundity, which will lead to high population growth rates (Pimm 1991). However, Sæther et al. (2004) found that longer reproductive lifespans are associated with higher establishment success for birds. In a recent paper, Capellini et al. (2015) showed that long lifespans and large litters are traits that characterize successful introduction and establishment of mammalian invaders. Further, their work showed that spread of mammalian invaders after successful establishment can also be predicted by age at 1st reproduction. All these traits (potentially large litter sizes, long reproductive spans, and early maturation) have been reported for beavers (Boyce 1981) and shed light on other factors that facilitated the spread of beavers in Patagonia.

Smaller colony sizes and lower number of kits in the forest may be a sign of density-dependent regulation in beaver populations. The density of colonies at our study sites in forest is certainly at the upper end of the highest densities recorded, and in contrast to the steppe (Pietrek and González-Roglich 2015), most suitable sites in the forest have already been occupied. Our finding that larger tree diameters can negatively affect colony size in the forest combined with the absence of any vegetation effects in the steppe supports the hypothesis that because beavers have occupied the forest for a longer period of time and have already removed smaller trees in many sites, they are now experiencing negative density dependence there but not in the steppe. Smaller litter sizes have been reported at high colony densities (Bergerud and Miller 1977; Payne 1984b). However, we did not observe differences between habitats in the proportion of colonies producing kits in the last year of our study.

In Table 1, we compared litter sizes reported for beavers in North America with those in our study. We have here included only studies in North America that used observations or complete colony removal to estimate the number of newborns. Other studies have looked at the number of embryos or placental scars in females, which overestimate the number of kits. The mean number of beaver kits we found in the steppe is consistent with most estimates from North America (Table 1) and from invasive North American beavers in Europe (Danilov et al. 2011; Parker et al. 2012). The lower number of kits in the forest agrees with the work of Payne (1984b) who reported mean litter sizes as low as 1.5 kits/colony at high densities. Thus, fewer kits in the forest may again signal the onset of density dependence there but not in the steppe.

Juvenile survival and fecundity are thought to be among the 1st demographic rates to decline in response to density dependence in large mammals. Eberhardt (1977) proposed that in marine mammals, increases in density should first affect early-age survival, then fecundity, and finally, adult survival. Gaillard et al. (2000) confirmed that pattern using large ungulates as a model. Although we could not separate survival of kits between forest and steppe, fecundity seems to be declining in forested areas invaded a longer time ago, whereas survival of adults remains high in both habitats. Thus, beavers could be another instance of this general pattern found in larger mammals, but a further study to get a better estimate of survival of kits is needed to support this hypothesis.

Many studies on the spread of invasions assume that negative density dependence does not affect the speed of the invasion at the front. Although, we did not study demographic rates at the front, more recently colonized areas showed higher recruitment than areas invaded a longer time ago, suggesting that the influence of negative density dependence is weaker near the front. While demography near the front may be the major contributor to the speed of the invasion, changes in dispersal as a consequence of changes in density may also affect the speed of the invasion.

Our results have important implications for managing the spread of the invasion of beavers in Patagonia. First, we should leave behind the traditional view of the steppe as suboptimal habitat that cannot support the invasion, and instead encourage active management, particularly in areas of the steppe that can be a source of propagules to mainland Patagonia (Fig. 1). Moreover, as detectability is higher in the steppe, this may be a better place to monitor beavers and to cull the population. Second, high survival is likely a major driver of the invasion. In previous work we have conducted using an integrodifference-equation structured model (Neubert and Caswell 2000), we simulated the beaver invasion using demographic rates of the species in North America. The speed of the invasion was highly sensitive to increase in adult survival and an increase in survival rates had a much higher impact on the invasion speed than did an increase in the mean dispersal distance. Thus, reducing adult survival should be a key aim of management. Third, all else being equal, areas that were invaded longer ago will contribute proportionally less to the spread of the invasion than more recently invaded areas for 2 reasons: lower fecundity and lower proximity to the invasion front. Lastly, the new demographic data on beavers in Patagonia that we have presented here should be useful for parameterizing more realistic population models to better manage the ongoing beaver invasion in continental Patagonia.

Supplementary Data

Supplementary Data SD1.—Description of occupancy models used to estimate the proportion of colonies with kits. Tables including within-habitat variables predicting colony size and number of kits.

Supplementary Data SD2.—Data manipulation and structure of the multistate capture-recapture model.

Supplementary Data SD3.—Beaver captures discriminated by habitat, colony, sex, and age.

Acknowledgments

AGP wants to thank W. F. Morris for mentoring and advice and B. Gardner, O. Gimenez, M. Kery, and M. Schaub for statistical advice. M. Casalinuovo helped with the overflight telemetry surveys. M. Fasanella and C. Fachinetti provided assistance with the molecular sexing. G. Himes Boor provided helpful comments that improved the manuscript. This work was funded by the Center of Latin American and Caribbean studies at Duke, Cleveland Metroparks Zoo, CREOI (Conservation and Research Opportunities International), Duke Graduate School, Idea Wild, a Sue Wijdenes Memorial Grant from Neotropical Grassland Conservancy awarded to AGP, Planete Urgence, Rufford Foundation, Wildlife Conservation Society and CONICET. The authors finally want to thank personnel working at Estancia Cullen, Estancia Sara, and Estancia Harberton for allowing us to conduct research in their properties.

Author notes

* Correspondent:

Associate Editor was Marcus Vieira.