Are beavers a solution to the freshwater biodiversity crisis?

To determine whether reintroduced beavers, as an example of native herbivorous megafauna, can increase freshwater biodiversity at the landscape scale and to compare effects on two contrasting taxonomic groups.


| INTRODUC TI ON
The decline and extinction of megafauna have been a defining feature of the Anthropocene. In freshwaters, the scale of decline in megafauna (He et al., 2017) is symptomatic of the declines in their biodiversity globally (WWF, 2018). Most megafauna are uniquely affected by top-down pressures such as hunting and persecution, while sharing the bottom-up pressures of habitat loss, fragmentation or deterioration (caused, for example, by land-use change, pollution, invasions and climate change) with the wider biota (Reid et al., 2018). However, there is an argument that modern landscapes, and their attendant extinction crises, are also in some respects a legacy of the loss of megafauna and the impact this has subsequently had on dependent-processes, including energy cycling, ecological dynamics and maintenance of heterogeneity at a range of scales (Doughty et al., 2016). Might resurrecting populations of megafauna therefore go some way to reversing the present biodiversity crisis in freshwaters?
Among freshwater megafauna, beavers are unusual in that their populations are stable (North American beaver Castor canadensis) or expanding (Eurasian beaver Castor fiber) (He et al., 2017). This is largely due to continent-wide protective legislation (e.g. European Habitats Directive) and numerous reintroduction schemes that have together, over the last 50 years or so, reversed a long-term decline in beaver populations (Halley & Rosell, 2002;Naiman, Johnston, & Kelley, 1988). Beavers are also unusual in terms of the extent of their ecosystem engineering activities, although they are by no means unique among the freshwater herbivorous megafauna in playing important ecosystem engineering roles (Bakker, Pagès, Arthur, & Alcoverro, 2016;Bump, 2018;Moss, 2015). Beavers modify existing freshwater habitats (streams, lakes and ponds) by building dams to raise and stabilize water levels, thus maintaining a submerged lodge entrance that reduces exposure to terrestrial predators and assists foraging (Hartman, 1996). Beaver ponds are more dynamic than other permanent wetlands, partly due to a fluctuating hydrological regime caused by leakage from dams and their repair, that results in intermittent exposure of pond margins (Gurnell, 1998). Smaller scale disturbances, including selective herbivory (Law, Jones, & Willby, 2014), fall or wind blow of dead trees, accumulation of woody debris and excavation of channels by beaver all add to the uniqueness of these engineered wetlands (Hood & Larson, 2014). Dams retain sediment and organic matter, thereby modifying nutrient cycling and decomposition dynamics which influences water chemistry and materials transported downstream (Ecke et al., 2017;Naiman et al., 1988;Puttock, Graham, Cunliffe, Elliott, & Brazier, 2017). The physical and biological characteristics of surrounding areas are also altered by inundation, with megafauna playing a significant role in the movement of energy and materials across the aquatic-terrestrial boundary (Johnston & Naiman, 1987;Moss, 2015). Landscape-scale heterogeneity is typically increased through the combination of beaver-engineered and non-engineered habitat, and the coexistence of engineered sites ranging from newly formed to long-abandoned (Willby, Law, Levanoni, Foster, & Ecke, 2018), with potential benefits for multiple taxonomic groups (Rosell, Bozser, Collen, & Parker, 2005).
In view of the societal importance of freshwaters (de Groot, Brander, & Finlayson, 2016), the increasing evidence of how current and emerging pressures are affecting this resource (Dudgeon et al., 2006;Reid et al., 2018), and the unparalleled rate at which biodiversity is now being lost from freshwaters (WWF, 2018), novel approaches and adaptive methods are urgently needed to protect what remains, or restore what has been lost. A promising approach is to reinstate natural processes and allow the ecosystem to self-design (Sayer et al., 2016). The re-establishment of herbivorous megafauna, such as beavers, is one such natural and novel restoration method, that exploits their ecosystem engineering activities to promote biodiversity and heterogeneity, while restoring lost ecosystem services (Law, Gaywood, Jones, Ramsay, & Willby, 2017). Indeed, now that their population status is secure, beaver is increasingly being reintroduced specifically to restore landscape heterogeneity and increase resilience to floods and droughts via direct and indirect habitat transformation (Burchsted, Daniels, Thorson, & Vokoun, 2010;Halley & Rosell, 2002;Nolet & Rosell, 1998).
Evidently, beavers can create hydrologically distinctive features within the landscape (Nummi & Holopainen, 2014;Westbrook, Cooper, & Baker, 2010;Wright, Jones, & Flecker, 2002). As such, they have the potential to increase species richness across several organismal groups (Janiszewski, Hanzal, & Misiukiewicz, 2014;Rosell et al., 2005;Stringer & Gaywood, 2016) and indeed across ecosystems by creating aquatic-terrestrial linkages (Anderson, Paszkowski, & Hood, 2014;Nummi, Kattainen, Ulander, & Hahtola, 2011). Willby et al., (2018, established that beaver ponds supported higher plant species richness at patch and site scales and that turnover in composition was higher between patches in beaver ponds than other wetlands. Moreover, they also found that beetle richness and abundance (but not turnover) were higher in beaver ponds. However, while greater species richness and habitat heterogeneity are positive attributes in conservation, this does not preclude the possibility that such habitats are dominated by common, generalist or non-native taxa at the expense of scarcer, specialist or native ones; we cannot assume that novel habitats (or habitats created by novel means) will necessarily support novel biodiversity.
In this study, we assess the potential implications of beavers, as one of the few remaining widespread representatives of the freshwater herbivorous megafauna, for enhancing regional freshwater biodiversity. Our study was based in Sweden, to which beavers were successfully reintroduced from Norway between 1922 and 39, following their extinction in the 1870s (Hartman, 1996). We use plants and beetles as focal taxa which are ideal study groups since they are taxonomically diverse, indicative of particular environmental conditions and provide a contrast between passive and active dispersers. Beetles also exhibit high functional diversity and are wellknown taxonomically and biogeographically (Bilton, Ribera, & Short, 2019) and can be considered representative of wider macroinvertebrate assemblages (Bilton, Mcabendroth, Bedford, & Ramsay, 2006;Ruhí & Batzer, 2014). Furthermore, beetles are easily live-sorted from sample debris, and the majority of individuals can be identified in the field, thereby reducing destructive sampling. We tested (a) whether the composition, rarity and native status of plants and beetles, and growth strategies of vegetation, differ between beaver ponds and other wetlands co-occurring in the same landscape and (b) whether predictable differences in the physical environment between wetland types related to their origin drive these compositional differences. On this basis, we assess whether reintroducing selected megafauna could aid in the recovery of freshwater biodiversity elsewhere.

| Field sites
The study focused on a 100 × 100 km area between Örebro and Skinnskatteberg, in south-central Sweden (59°30′N, 15°10′W), dominated by managed forests or low-intensity agriculture. Here, valley wetlands formed through stream impoundment by beaver dams (beaver ponds-BP) coexist with other permanent, shallow, standing freshwaters such as small lakes, ponds and river oxbows (other wetlands -OW). A total of 10 BP and 10 OW were sampled.
All BP supported active beaver colonies (indicated by freshly grazed plants or felled trees, canal creation, dam maintenance and lodge construction) and were estimated to have been formed for at least 5 years (from aerial imagery and extent of dead wood). Other wetlands were close (<5 km) to sampled BP, but were not paired with specific sites (see Appendix S1 for a summary of local environmental variables).
The extent of leaf litter, open water, woody debris, bare ground and grazing associated with each plot or sample was scored visually on the 1-5 scale as above, while mean plant height and water depth were determined from replicate measurements. Water conductivity was measured using a multi-range conductivity metre (Hanna instruments HI 9033) calibrated to 25°C. These explanatory variables are comprised of environmental variables that may influence biotic assemblages and were expected to differ between wetland types (Willby et al., 2018).

| Data analysis
To determine whether sufficient waterbodies were surveyed per wetland type for the focal taxonomic groups, the sample coverage was calculated using the 'iNEXT' r package (Hsieh, Ma, & Chao, 2016) based on incidence data per site. Plant cover data were logtransformed (x + 1), and a plot level Bray-Curtis dissimilarity matrix was calculated. Due to low abundance (an average sweep sample contained six individuals (range 1-10), a dissimilarity matrix was calculated using binary transformed data for beetles. Species composition per wetland type was then compared using non-metric multidimensional scaling (NMDS). A permutational multivariate analysis of variance, based on 999 permutations, was used to test for differences in species composition between wetland types for both assemblages. Species characteristic of particular wetlands was identified using multilevel pattern analysis which tests for associations between species patterns and combinations of groups of sites (De Cáceres & Legendre, 2009). Total beta diversity for plants and beetles was partitioned into turnover versus nestedness components based on species abundance and incidence, respectively, using the 'betapart' r package (Baselga & Orme, 2012).
Rarity scores were assigned to plants and beetles using a 1-5 ranking (1 = common and widespread; 2 = common-frequent and fairly widespread; 3 = locally common but scattered; 4 = infrequent; and 5 = rare). Scores for plants were based on descriptions in Mossberg and Stenberg (2018) and assessment of the distribution of records in south-central Sweden and adjacent regions held by the Global Biodiversity Information Facility (GBIF; gbif.org).
Sample-specific rarity scores for plants and beetles were then derived based on the weighted average mean rarity score using weighting by cover (plants) or count of individuals (beetles). Differences in rarity between wetlands were tested using a Kruskal-Wallis test as these data did not meet parametric requirements.
Prior to constrained ordination using redundancy analysis (RDA), all continuous explanatory variables were log-transformed, mean centred and scaled by 1 SD to improve comparability between variables and to reduce the effect of outliers. Correlations between predictor variables were assessed in a correlation matrix and checked for variance inflation. The site × species matrices for plants (logtransformed cover) and beetles were analysed using all continuous predictors, with wetland type added as a categorical variable. In addition, the number of plant species, maximum plant coverage per plot (%) and plant height were added to the beetle RDA to determine whether there were secondary effects of vegetation structure on beetles. An automated, forward selection of predictor variables was conducted on the initial global model with the most parsimonious models being selected based on the significance of each variable (p < .05) using the 'vEgaN' r package (Oksanen et al., 2019).
Plant growth strategies were assigned based on Pierce et al. (2017). This approach uses the major axes of variation in functional leaf traits associated with size and resource economics to represent Grime's (Grime, 1977) opposing plant growth strategies, C (competitor), S (stress tolerator) and R (ruderal), on a continuous scale. Of the 156 species we recorded, 78% could be matched directly with the species documented by Pierce et al. (2017), with a further 13% matched to a closely related documented species with a similar growth form, habitat and life history. For each quadrat, we then calculated the representation of CSR strategies within the vegetation, weighted by the cover scores of the component species, following the approach of Willby, Pulford, and Flowers (2001). Differences in representation of each growth strategy between wetlands at the quadrat scale were then tested using generalized linear mixed effects model with site as a random intercept. Ten quadrats, in which species unassigned to CSR accounted for >30% of the plant cover, were excluded from the analysis. In terms of functional responses, we focused on plants since they are most likely to respond directly to changes induced by beavers, while the pool of species present also span the gradient of variation in growth strategies. Suitable trait data for beetles within this study area were unavailable.

| Compositional differences between wetlands
Estimated sample coverage was generally high (mean = 94%) indicating effective sampling of each taxonomic group per wetland type (Table 1). Within each site, a sufficient number of samples were taken to capture the majority of plant species, and however, sampling more beetles would have resulted in a greater number of species being found (Appendix S2). Half the total species pool of plants (156 species) and 45% of the species pool for beetles (66 species) was shared by both wetland types. For both taxonomic groups, a higher proportion of the total species pool was found only within BP (31% and 30% for plants and beetles respectively), compared to OW (19% and 20% respectively). These general differences can be visualized in the unconstrained ordination (Figure 1a,c). Despite the large overlaps in the hulls for each taxon group between wetland types, the mean species composition (as represented by the centre of each ordispider) differed significantly between wetlands for both plants (p < .001) and beetles (p = .034). Total beta diversity was strongly dependent on turnover for both plants and beetles (96% and 94% respectively), rather than nestedness (4% and 6%).
A total of 37 plant species were significant indicators (p < .05) of a wetland type (Figure 1b

| Environmental basis for differences between wetlands
When both species assemblages were constrained by local environmental variables (see Appendix S1), the separation of the two wetland types was more distinct (Figure 2). In both cases, the overall constrained models were significant (p < .001 (plants); p = .018 (beetles)). For plants, plots from BP were associated with more woody debris, open and bare ground, while those in OW had greater leaf litter, plant height and plant coverage (Figure 2a). Water depth was the only significant environmental variable that explained beetle assemblages, though was driven by one outlying site (Figure 2b). When this outlier was removed, the overall model was not significant (p = .186). Wetland type accounted for a significant proportion of the compositional differences for plants (p < .001), over and above the effect of other variables, but not for beetles (p = .136). However, only 11% of variance in composition was explained in either model.

| Differences in growth strategies between wetlands
No significant differences were found between wetland types in the representation of the competitor growth strategy in the TA B L E 1 Summary of species richness , uniqueness, sampling efficiency and rarity (mean ± SE (range)) per wetland type for each taxon group quadrat-level vegetation (p = .16) (Figure 3a). However, in BP the representation of stress tolerators was significantly lower (p = .01), while ruderals were more common (p = .002) in comparison with OW ( Figure 3b,c). Specifically among the subset of indicator species, the mean representation of growth strategies in BP indicator plant species (23.2 ± 3.7, 29.9 ± 5.3, 46.8 ± 5.4% for CSR respectively) contrasted strongly with the OW indicators (51.5 ± 11.1, 39.3 ± 10.8, 9.1 ± 3.5% for CSR respectively), highlighting a strong characterization of BP vegetation by ruderals and OW vegetation by competitors and stress tolerators.

| D ISCUSS I ON
The loss of megafauna from modern landscapes has contributed to deterioration of ecosystem function and heterogeneity, with cascading negative effects on biodiversity (Doughty et al., 2016). The scale and consequences of this loss often only emerge fully after populations of megafauna have been restored (Bakker & Svenning, 2018). Beavers are increasingly recognized as facilitators of natural ecosystem processes and a keystone species (Ecke et al., 2017;Stoffyn-Egli & Willison, 2011;Stringer & Gaywood, 2016). In this F I G U R E 1 Unconstrained ordination of plants (a and b) and beetles (c and d) for beaver ponds (blue) and other wetlands (red), using non-metric multidimensional scaling (NMDS). Convex hulls enclose each wetland type with 'spider' plots showing spread of samples from the wetland type centroid. NMDS species scores are shown on the right hand plots and coloured for species significantly associated (p < .05) with each wetland; grey = no association, blue = beaver pond and red = other wetland

(d)
study, we found that wetlands engineered by beavers supported a larger and more distinct species pool than other wetlands in the same landscape and that this effect was not simply a product of attracting common, generalist or non-native taxa. Therefore, beaver ponds are not subsets of other wetlands, as indicated by their beta diversity being predominantly explained by turnover, with only a low percentage explained by nestedness. Differences in species and/or growth strategy composition of focal biota between wetlands could be clearly traced to the features of wetlands related to their differing origins and disturbance regime, although effects on plants were more pronounced than those on beetles. Our results suggest that re-establishing beaver populations could be an important mechanism in supporting freshwater biodiversity recovery in degraded landscapes.
In this study, despite negligible differences in some physicochemical variables (see Appendix S1), these wetland types sup-  ** some known from independent studies to be heavily grazed by beavers, for example Menyanthes trifoliata, Nymphaea alba and Schoenoplectus lacustris Milligan & Humphries, 2010;Willby, Perfect, & Law, 2014). Our observation that wetland type remains a significant predictor of vegetation composition after the effects of measured abiotic factors are accounted for is consistent with other evidence that herbivory by beavers (a direct effect), significantly influences wetland vegetation composition (Law et al., 2017 over and above the effects of dam-induced inundation (indirect effects). In contrast, indicator plant species in BP were typically smaller, faster growing ruderals (e.g. Alisma plantago-aquatica, Glyceria fluitans, Rorippa palustris, Callitriche spp.), often associated with shallow water or intermittently exposed or otherwise disturbed margins (Abernethy & Willby, 1999). The differences between wetlands in growth strategy composition of their vegetation clearly reflects the patchiness and ecological dynamism imposed by beavers. The persistence of these compositional differences between engineered versus non-engineered wetlands after beavers disperse to new territories is unclear, although some studies indicate they may last for decades (Bartel, Haddad, & Wright, 2010;Ray, Rebertus, & Ray, 2001).

Conductivity
In contrast to plants, few beetle species were significantly associated with either wetland type; H. heydeni and I. ater, which were associated with BP, are both typical of enriched, well-vegetated standing waters (Foster, Bilton, & Nelson, 2016). The lack of indicator species could be partly explained by the relatively high dispersal ability of beetles (Bilton et al., 2019), as we found no evidence that beaver-related effects on vegetation structure explained differences in beetle assemblages. Bloechl, Koenemann, Philippi, and Melber (2010) and Lundkvist, Landin, and Milberg (2001) also found that variance in water beetle composition across the landscape was low in both artificially created and agricultural ponds. Factors behind habitat selection by beetles are poorly known and are likely to be scale, species, and life stage dependent (Lundkvist et al., 2001;Yee, Taylor, & Vamosi, 2009). One potentially important influence on beetle richness and composition is habitat complexity linked to heterogeneity of vegetation (Gioria, Bacaro, & Feehan, 2011), or features such as beaver-dug canals and dead wood (Hood & Larson, 2014), occurring at finer scales than we measured. Willby et al. (2018) found that BP had 26% more water beetles than OW, thereby suggesting that beetles benefit from habitat heterogeneity at least in terms of their density.
In common with other ecosystem engineers (Romero, Gonçalves-Souza, Vieira, & Koricheva, 2015), the effects of beavers can be taxon specific, but also appear to be context specific.
For example, disturbances may leave some ecosystems more susceptible to invasion by non-native species (Strayer, 2010), as found from beaver studies in North and South America (Lesica & Miles, 2004;Westbrook, Cooper, & Anderson, 2017). By contrast, we found no increase in non-native in beaver wetlands (although non-native plants were intrinsically scarce), in common with other studies of beaver wetlands in the United States (Brzyski & Schulte, 2009;McMaster & McMaster, 2000). Moreover, the evidence that rewilding contributes to biological invasions is mixed, with several factors at play, including a strong influence from local processes (Derham, Duncan, Johnson, & Jones, 2018). This implies that a link between habitat disturbance by ecosystem engineers and biological invasions is not generalizable. The novelty of the megafaunacreated habitats themselves, relative to those that already exist in a landscape, should also be considered context specific (Wright et al., 2002). In the present study, novel wetland habitats (or habitats created by novel means) supported a unique, albeit not necessarily rare, suite of species that were present in the regional species pool and able to colonize this new habitat. By contrast, biota within ponds formed by invasive North American beaver in Chile were largely similar to those found in naturally occurring lentic habitat (Anderson, Vanessa Lencinas, et al., 2014), with ubiquitous, rather than unique species being found in beaver ponds (Anderson & Rosemond, 2007). Evidently, ecosystem responses to invasion by beavers far outside their native range cannot be considered a useful guide to their effects when reintroduced within their native range. Finally, while reviews of the effects of ecosystem engineers generally (Romero et al., 2015) or beavers specifically (Stringer & Gaywood, 2016) tend to reveal positive effects, it is likely that these effects are moderated by factors such as the population size of the engineer and the seasonality of resource supply (Brzyski & Schulte, 2009).

| Implications
Studies of comparative biodiversity across coexisting habitats often infer compositional differences rather than directly quantifying these. Our study demonstrates that, in their natural range, beavers create ponds that, while superficially similar to other shallow wetlands, differ subtly in their physical characteristics and disturbance regime. This results in distinctive species assemblages that are indicative of beaver wetlands, particularly so for plants, rather than being simply a subset of those found in other freshwater habitats in the same landscape. With freshwater biodiversity declining at an unsustainable rate, recognizing the major role that herbivorous megafauna can play in creating novel habitats, increasing spatial heterogeneity and stimulating ecological dynamism is key. However, this role is largely underappreciated (Moss, 2015) and even domesticated livestock now exert only very limited influence on wetlands in most agricultural landscapes due to their routine exclusion by fencing. Habitat engineering by beaver offers a passive, low-tech, megafauna-based ecosystem restoration technique that could be widely implemented under the flagship of rewilding, and which, as this study shows, remains effective even in landscapes affected by agriculture and forestry. Of course beavers alone cannot solve the freshwater biodiversity crisis but wider acceptance of the mounting evidence for the 'healing power' of these and other megafauna will be progress.

ACK N OWLED G EM ENTS
The Swedish Research Council Formas (grant no. 2010-1647) financed the contributions of OL and FE. We thank the Carnegie Trust for a grant to cover travel costs and a University of Stirling Horizon studentship for funding of AL.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data referred to in this article and code used in analyses are deposited in DataSTORRE-the University of Stirling research data repository.