Riparian forests throwback at the Eurasian beaver era: a woody vegetation assessment for Mediterranean regions

Several biotic and abiotic interactions will contribute to riparian ecosystem changes. The impact of Eurasian beaver (Castor fiber) on woody vegetation is still unknown for the Mediterranean biogeographical area. Through a replicable approach applied on a cluster of three rivers, we studied how the tree layer of Mediterranean riparian sites is impacted by the beaver's recent comeback. For each site, we collected data (e.g., stem diameter, species, distance from riverbank) for all standing trees and additional information only for gnawed trees at plot level. Data elaboration allowed to characterise impacts on riparian vegetation. Salix spp. and Populus spp. are the main gnawed species, but sporadically other species can be selected based on their size and spatial distribution (e.g., Alnus glutinosa). Diameter means of gnawed trees are significantly lower than the not gnawed ones. Most of the selected trees have low diameter classes (< 12 cm), even if diameter preferences may vary on the basis of overall stand tree size range and distribution. Over 90% of the gnawed trees are entirely harvested, with stumps as the remaining standing element. Main changes on the overall forest stand occurred in the first ten metres from the riverbank, as beaver gnawing activity is significantly influenced by the interaction among tree distance from the river and diameter size. Our approach can be used as a model system to be implemented in other Mediterranean sites where beaver is expanding, with the aim of predicting mid-term riparian forests vegetation changes.


Introduction
Mediterranean riparian forests are one of most threatened ecosystems (Rodríguez-González et al. 2021).The exposure to land-use changes and to global warming, coupled with new biotic interactions might potentially contribute to structural shifts of these ecosystems in the mid-term.Riparian woody vegetation structure and biodiversity can deeply change following Eurasian beaver (Castor fiber L.) gnawing activities (Piton et al. 2020).Eurasian beaver is a semiaquatic rodent, which spends most of its time in the water, whereas walking on land only to search for food (Campbell-Palmer et al. 2016).The Eurasian beaver is listed in Annex III of the Berne Convention (1979) and in Annexes II and IV of Habitat Directive (92/43/EEC).Both legal instruments require protection and establish the possibility of reintroduction after feasibility studies considering both the acceptability by the general public and potential impacts.Gnawing activities on trees are made essentially for feeding and river dams creation (Puttock et al. 2017;Brazier et al. 2021).Its diet is strictly vegetarian and includes both wooden and herbaceous plants (Nolet et al. 1994;Fustec et al. 2001;Campbell-Palmer et al. 2016).Literature underlines that woody species preferences can vary on the basis of different variables and their interactions (Nolet et al. 1994;Campbell-Palmer et al. 2016).Woody species taxonomy, shoot diameter, and distance from the riverbank seem to play a determinant role (Haarberg and Rosell 2006) at local level (forest stand scale), but their relative influence on beaver activity remains to be clarified.Salicaceae species (Salix spp., Populus spp.) have been reported as the most selected by beaver populations in Hungary (Juhász et al. 2022), Poland (Jackowiak et al. 2020), Czech Republic (Mikulka et al. 2022), and France (Piton et al. 2020).Where Salicaceae species are less available, other species (e.g., dogwoods and alders) are selected (Fryxell and Doucet 1993;Mikulka et al. 2022).In young riparian stands, Nolet et al. (1994) found that mainly shrub willows were selected as being the most abundant plant species locally.Plants belonging to other genera (e.g., Alnus, Corylus, Fraxinus, Populus, and Prunus) were selectively consumed, even if present sporadically within the stand.Trees with diameter up to 30 cm can be occasionally gnawed (O'Connell et al. 2008), but only trees with relatively small diameters (< 15 cm) are significantly selected (Haarberg and Rosell 2006;Janiszewski et al. 2012;Jackowiak et al. 2020).With the increasing distance from the water, the probability of beaver gnawing activity on woody plants decreased (Haarberg and Rosell 2006), with some different trends among the woody plant genera (Mahoney and Stella 2020;Mikulka et al. 2022).Some studies reported that trunk diameter is more important than species composition (Mahoney and Stella 2020), others that the selection of woody plants is more pronounced with increasing distance from water (Fryxell and Doucet 1991).
The Eurasian beaver was once present in a large part of the Palaearctic, ranging from Western Iberian Peninsula to North-western China, and Western Mongolia, throughout all suitable riparian habitat types in forests, tundra, and steppe (Halley et al. 2021).At the start of 1900, the range of the Eurasian beaver was limited to a few, small-sized refugia between France and Mongolia, hosting less than 1200 individuals (Halley and Rosell 2002).Currently, the species shows reproductive populations in most of its original range, mainly due to successful reintroduction programs (Halley et al. 2021;Kodzhabashev et al. 2021;Calderón et al. 2022;Paladi and Cassir 2022).Apart from the refugium of Southern France, in the Mediterranean area, Eurasian beaver has been only very recently recolonising riparian habitats (Halley et al. 2021;Pucci et al. 2021;Calderón et al. 2022).So far, the impact on the European riparian woody vegetation has been described only for the northern and eastern countries (e.g., Hartman 1996;Haarberg and Rosell 2006;Jackowiak et al. 2020), or for a few cases in the France Alps (Piton et al. 2020).No studies have been conducted yet within the Mediterranean biogeographical region (sensu Olson and Dinerstein 1998: see Larsen et al. 2021;Grudzinski et al. 2022).With our study, we aimed to carry out a woody vegetation analysis of Mediterranean riparian forests undergoing beaver recolonization.The specific goals are to understand (i) the impact on woody vegetation in terms of gnawed trees by beaver, (ii) which local variables influence beaver tree gnawing activity and (iii) spatial and dimensional thresholds influencing gnawing activity.Following previous literature, we predicted that willows and poplars would be the preferred species by beavers rather than other genera, and that intermediate diameters would represent the most consumed plants independently by distance from the riverbank.Thus, we tested a novel and replicable approach on a cluster of three rivers located in central Italy, as a model system to study how the tree layer of Mediterranean riparian sites is impacted by beaver recolonization.

Study area
The study was carried out between summer and autumn 2022, within two river basins characterised by the presence of Eurasian beavers, i.e., Ombrone-Merse and Tevere (see map in Fig. 1a).
In the sites it is confirmed the presence of reproductive naturalised populations since the end-2019 (Mori et al. 2021(Mori et al. , 2022a;;Pucci et al. 2021;Viviano et al. 2022).Apart from two single individuals in two areas of North-Eastern Italy (Pontarini et al. 2019;Pucci et al. 2021), the species is actually still confined only in this portion of central Italy and it is in its initial expansion phase.
Vegetation assemblages, flow variability and natural and human disturbances are typical of Mediterranean semi-natural riparian areas (Biondi et al. 2003;Hupp and Rinaldi 2007;Zaimes 2020).As shown by climate diagrams (Fig. 1a), the local climate is typically meso-Mediterranean, with most precipitation occurring during autumn and winter.

Sampling design and data collection
Beaver presence throughout rivers is still very fragmented, thus the maximum length of a riparian forest with continuous beaver activity in the study area was between 60 and 120 m.In the climate diagrams, the red line represents the variations of the mean monthly temperatures, the blue line represents the variations of the monthly precipitation, the blue hatched area the humid period, the red dot area the dry period and the blue filled area represents the wet period.The coloured boxes along the x axis show months where frost is likely (cyan) or definite (blue).Validated and continuous temperature and precipitation data were provided by the Regional Hydrological Sector of Tuscany (SIR) and cover the period 1993-2022 (https:// www.sir.tosca na.it/ consi stenza-rete).For our purposes we selected two SIR meteorological stations located near the study plots and representative of the homogeneous areas.The weather station of Sansepolcro (code T0S11000039, Lat 43.559,Long 12.097) was selected as representative of Tevere river whilst Buonconvento (code T0S11000067, Lat 43.092,Long 11.439) was selected for Merse and Ombrone river.The climate diagrams were created in R using the "climatol" package (Guijarro 2019).b Sampling design exemplification.The picture shows an exemplificative site where two homogeneous stand are surveyed collecting data in two representative plots for each stand.Minimum required distances between plots and stands are also reported 1 3 We selected two riparian stands for each of the three surveyed river sites (for a total of 6 stands (Fig. 1b)).Stands were selected on the basis of the significant beaver activity (i.e., abundant "fresh" signs: Mori et al. 2022a).Within every single river site, each stand was at least 1 km distant from the other, to include the overall heterogeneity of the riparian forest.We collected forest stand data within two rectangular (5 × 20 m) plots with minimum distance among each other of 30 and maximum of 60 m.Plots were located perpendicular to the main river flow with one of the two short sides placed at the edge of the riverbank (Fig. 1b).Mean values of the two plots were assumed as being representative of each riparian forest stand.
Within the plot, we collected data on standing trees (Table 1) adapting previous experiences (e.g., Misiukiewicz et al. 2016;Mikulka et al. 2022) to our study area and study aims.Data were collected in an average time frame of an hour and half (3 operators).The diameter threshold used to select woody elements was ≥ 3 cm at 15 cm above the ground, as in presence of very small diameter beavers often eat or otherwise use the entire plant and thus gnawing activity could be difficult to detect (Crisler and Russell 2010).

Data analysis
Descriptive statistics were conducted at different aggregation levels (river site and stand), through R statistical software (R Core Team 2021).The "ggpubr" (Kassambara and Kassambara 2020) and "ggstatsplot" (Patil 2021) packages were used to test and visualise data distribution, as well as to compare means between groups and to perform correlations tests.To test the influence of different factors on the probability of tree gnawing, a generalised linear mixed model (GLMM) with a binomial distribution was used (Zuur et al. 2009), through the "lme4" package (Bates et al. 2015).We considered the presence/absence of gnawing activity on the woody plant as the response variable, whereas the distance from water and diameter were used as covariates.Genera (for Salix and Populus spp.) or a taxonomic cluster (for all other genera) and their interaction with distance and trunk diameter were included in the model as explanatory variables.Stands (representing woodland variability in each site-river) was used as a random factor.We used the generalised form since data deviates from the normality (Zuur et al. 2009), and Laplace approximation since we have fewer than five random effects (Ju et al. 2020).The Wald statistic test (χ 2 , Zuur et al. 2009, for GLMM) was used to test the effect size of fixed effects.Moreover, the overall explanatory power of the models was estimated by calculating conditional R-squared with the "r.squaredGLMM" function of the R package "MuMin" (version 1.15.6;Barton 2016).Plot regression models were obtained through the "sjPlot" package (Lüdecke et al. 2021).All differences and models were statistically significant at the significance level p < 0.05.

Riparian forests description
The surveyed stands captured forest stand composition variability within and among different rivers (SI 01).On the basis of the proportion of the tree layer basal area, we found two kind of riparian woodlands, i.e., those from Merse and Ombrone (referring to the Populion albae phytosociological alliance) rivers and those from Tevere (referring to the Alnion glutinosae phytosociological alliance) river.Sporadic species, sometimes also relatively Ombrone river is characterised by significantly lower tree density and higher mean diameter with respect to the other two river stands, probably due to the different water flow and sediment dynamics.The high value of diameter quadratic means highlights the remarkable variability of tree size within a stand.Overall, diameter distribution frequency (Fig. 2) shows a clear dominance of low diameter trees (< 20 cm, at 15 cm above the ground), with most trees under the threshold of 12 cm.A slightly positive relationship between tree DBH and distance from the riverbank (SI 02a), emphasises the tree spatial arrangement of surveyed riparian woodlands, where largest trees are mainly located in the farthest part of the riverbank.This stand structure can be explained, especially in the case of Salicaceae, with the morphologic variation of the riverbank (i.e., increasing elevation from riverbank to the innermost portion of riparian forest) that determines the physical and ecological conditions for tree species recruitment (Mahoney and Rood 1998).

Impacts on woody vegetation
Riparian forests showed beaver gnawing signs up to 12.5 m from the riverbank (Table 2, Fig. 3); beaver activity was recorded almost entirely within the first ten metres from the riverbank.The highest diameter (h = 15 cm above the ground) of a gnawed tree was 32 cm; most of the gnawed trees fall under 12 cm diameter class (Fig. 3, SI 02b).Mean diameter values varies among stands (SI 03) and significant difference between diameter means of selected and unselected trees is highlighted (Fig. 4).Up to 5 tree species (i.e., Tevere river) were impacted by beaver activity (Table 2).Even if most of the selected trees fall into the taxonomic family of Salicaceae (Salix spp., Populus spp.), other less frequent species can Fig. 2 Tree diameter (15 cm above the ground) frequency distribution in the surveyed rivers be sporadically impacted (SI 04a).Among selected trees, Salix alba and Populus nigra had the highest diameter mean values (SI 04b).No significant diameter difference among groups of species (SI 04c) was found for selected trees, even if significant differences were found only for unselected ones.  1 3 Trees were entirely harvested for more than 90% of the cases; it is noteworthy a significant presence of resprouting stumps of Salicaceae (see SI 05).In a very few cases (5% of the total selected trees), trees were gnawed just in some portions of the stem, with partial damage to the wood and phloem or to the bark.

Influence of factors on the probability of tree selection by beavers
GLMM showed a strong influence of specific taxonomic groups (Salix spp.and Populus spp.) on beaver tree selection (χ 2 (2) = 16.15,p < 0.001).A significant influence is driven even by the interaction of distance from the riverbank and diameter size class (χ 2 (1) = 5.49, p < 0.05), and diameter size class per sé (χ 2 (1) = 4.47, p < 0.05).Interestingly, there is no influence of the distance factor per sé.The overall model explains 44% of the beaver choice (R 2 marginal = 0.44).
Interestingly, the probability curves showed two different trends on the basis of diameter size classes (i.e., middle-low: ≤ 6 cm; middle-high: > 6 cm: Fig. 5); for the middle low class of Salicaceae a slight probability of impacted trees with increasing distance was detected.This effect is stronger for Alnus and other species, with respect to willows and poplars.Opposite trends can be assessed for the middle-high diameter size class, for which the probability of impacted trees with distance strongly decreases for Salicaceae.

Discussion
Mediterranean riparian forests are characterised by high environmental suitability for beavers, but they are critical as deeply anthropized (Hupp and Rinaldi 2007).In over 500 years, rivers of the Mediterranean basin have changed their flowing and the ecosystems they host (e.g., Bellotti et al. 2004;Tarragoni et al. 2011), mainly due to human interventions.The Eurasian beaver is widely considered as the quintessential ecosystem engineer, as showing impacts upon geomorphology (including canal and den excavation, increase of woody debris, dam and lodge building), hydrology, aquatic ecology and nutrient cycling (Rosell et al. 2005;Brazier et al. 2021).Therefore, the occurrence of beavers in Central Italy after centuries of absence (Höfle et al. 2014;Juhász et al. 2020;Salari et al. 2020) may lead to additional ecosystem modifications which can affect the local biodiversity (e.g., Brazier et al. 2021;Viviano et al. 2022).With this study, we tested the effectiveness of a novel and easily replicable methodological approach, specifically designed for Mediterranean riparian woodlands recently re-colonised by Eurasian beavers.To compare harmonised data among other Mediterranean sites, we provide a step-by-step protocol which can be implemented in other riparian woodlands (SI 06).With respect to other protocols (e.g., Zwolicki et al. 2019;Piton et al. 2020;Pejstrup et al. 2023), it gives priority to collect data in a fewer number of small plots, but clustered for homogeneous and representative stands scattered in different river sites.
Our results show that the main changes on overall stand composition and structure of the riparian forest can be assumed for the first ten metres from the riverbank; standrelated variations are based on diameter frequency classes and tree spatial distribution.The lower mean diameter values of selected trees was consistent with some studies from central-northern Europe and North America (Janiszewski et al. 2017;Jackowiak et al. 2020;Mahoney and Stella 2020).Compared to those results, our work emphasises that, regardless of the available variability of stem diameter classes, the diameter class size of gnawed stems did not vary significantly among stands and sites.Most of the gnawed trees fall into a little lower diameter class with respect to Jackowiak et al. (2020) and Juhász et al. (2022).This evidence should be tested also in rivers where flooding disturbances are less frequent, with different tree regeneration dynamics (especially for Salicaceae, see Karrenberg et al. 2002); indeed, riparian forest structure can differ (e.g., fewer trees in lowest diameter classes) between regulated and semi-natural rivers (like those of our study, characterised by a predominance of free-flowing dynamics).Differently from other studies (Haarberg and Rosell 2006;Dvořák 2013), we found sporadic trees over 25 cm gnawed in all of the surveyed sites within the first ten metres from the riverbank.This evidence is of particular importance for possible structural changes of riparian forests and the management of river basins.
It's confirmed beaver tree preferences for Salicaceae species (Salix spp., Populus spp.) as reported for several European studies (Juhász et al. 2022;Jackowiak et al. 2020;Mikulka et al. 2022;Piton et al. 2020).This may be due to the abundance of this taxonomic family in the surveyed riparian woodlands, but different forest types can change beaver attitudes (Haarberg and Rosell 2006;Janiszewski et al. 2017;Piętka and Misiukiewicz 2022).Due the predominance of harvested stems compared to other type of damages and the expected growing beaver population size (Misiukiewicz et al. 2016), we can assume significant changes on river ecosystem structure and functionality (King et al. 1998;Mahoney and Stella 2020) in the mid-term.Interestingly, gnawing activities will tend to replace Salicaceae high trees with their shrubby growth forms.Further, Hall (1960) showed how consumption of overgrazed poplars is continuously replaced by proportionally higher consumption of willows, due to the slower resprouting of poplars in comparison to willows.The regeneration from resprouting stumps could also result in a change in the genotypic composition of a dominant species and hybridising complex, acting beavers as agents of natural selection (Bailey et al. 2004).
Our results showed that sporadic species can be selected even if not so common within the stand.This can be due to beaver diet that tend to avoid deficiencies problems (Nolet et al. 1994) as well as beaver preferences in terms of tree size and distance from the riverbank.Indeed, the results of regression analysis support the evidence that, with the increasing distance from the water, the probability of beaver felling on large woody plants (i.e., with diameters > 12 cm) decreased, with some different trends among taxonomic groups (Haarberg and Rosell 2006;Mikulka et al. 2022).Interestingly, the probability curves show two different trends that are diameter class-related, with similar tendencies among species groups.This supports the central-place foraging theory (Schoener 1979) which were tested for beavers; the theory highlights that the species should become increasingly more selective the farther the trees are found from the central place, i.e. the riverbank (Fryxell and Doucet 1991;Haarberg and Rosell 2006;Mahoney and Stella 2020).Jenkins (1975) described beaver as from 'choosy generalist' to simply 'generalist', on the basis on the different riparian stand characteristics; our study evidences the overall "opportunistic" attitude of the species, which tend to prioritise the characteristics of the tree (i.e., diameter and distance from the river) rather than the species per sé.Selective foraging on softwood genera may increase local forest regeneration (Fustec et al. 2001;Jones et al. 2009) and, on the other side, may also accelerate the spread of alien invasive species (IAS), due to new forest gaps that are created by gnawing activities (see Kurokochi et al. 2010).Otherwise, eventual gnawing activities on invasive trees may favour new sprouting roots of Robinia pseudoacacia, even if this species seems to be less used by beavers (Juhász et al. 2020(Juhász et al. , 2022)).Therefore, monitoring beaver distribution should be required also to control the dispersal-related processes of IAS (Catford and Jansson 2014), as well as to implement strategies for enhance the regeneration of native riparian trees (Juhász et al. 2020).Further research can implement our protocol in riparian forests where the presence of IAS (see Juhász et al. 2022) is more consistent due to different riparian conditions (Aguiar and Ferreira 2013;Catford and Jansson 2014).
A recent species distribution model suggests that beavers may further expand their distribution in the Italian and Iberian Peninsula (Serva et al. 2023).Currently, the main predator of the Eurasian beaver, i.e. the grey wolf Canis lupus, even if occurring in our study area, is not preying this large rodent yet (Mori et al. 2022b).Therefore we expect a further range expansion of beavers in Italy.Despite their historical presence, the return after 500 years of absence seems to be almost unlikely (Pucci et al. 2021;Capobianco et al. 2023).If on one side the Eurasian beaver is protected by the Habitats Directive, the species should anyway be treated as alien if released by humans without a formal authorization.However, removal strategies for unofficially released Eurasian beavers have been both ineffective and expensive (e.g., in Spain: Calderón et al. 2022).Moreover, Viviano et al. (2023) showed that, in Central Italy, a widespread positive attitude towards beaver presence is occurring; imposing strong management actions should thus include a detailed communication campaign, as it could trigger a chronic and expensive problem.Therefore, in our study areas, intensive monitoring coupled with structural measures may prevent conflicts and avoid drastic measures aimed at beaver removal (Swinnen et al 2017;Campbell-Palmer et al. 2021).Our work details new insights on the impact of beavers on woody species in the Mediterranean river basins, and the results deserve as base for further studies to predict structural and compositional changes of Mediterranean riparian forests in the mid-term (Brazier et al. 2021;Mikulka et al. 2022).The collection of harmonised data will facilitate clear communication by the scientific community to better understand potential conflicts with human needs (Auster et al. 2020;Campbell-Palmer et al. 2021;Viviano et al. 2023) as well as new ecosystem services provision by the species (Auster et al. 2020;Blewett et al. 2021).
Italian Ministry of University and Research funded by the European Union -NextGenerationEU; Project code CN_00000033, Concession Decree No. 1034 of 17 June 2022 adopted by the Italian Ministry of University and Research, CUP B83C22002930006, Project title "National Biodiversity Future Center-NBFC".

Fig. 1 a
Fig. 1 a Climate diagram according to Walter and Leith and location of study sites within the boundaries of Toscana and Umbria administrative region.In the map, stars with different colours represent the three river sites where the stands are located.The black square and triangle represent the position of the weather stations.In the climate diagrams, the red line represents the variations of the mean monthly temperatures, the blue line represents the variations of the monthly precipitation, the blue hatched area the humid period, the red dot area the dry period and the blue filled area represents the wet period.The coloured boxes along the x axis show months where frost is likely (cyan) or definite (blue).Validated and continuous temperature and precipitation data were provided by the Regional Hydrological Sector of Tuscany (SIR) and cover the period 1993-2022 (https:// www.sir.tosca na.it/ consi stenza-rete).For our purposes we selected two SIR meteorological stations located near the study plots and representative of the homogeneous areas.The weather station of Sansepolcro (code T0S11000039, Lat 43.559,Long 12.097) was selected as representative of Tevere river whilst Buonconvento (code T0S11000067, Lat 43.092,Long 11.439) was selected for Merse and Ombrone river.The climate diagrams were created in R using the "climatol" package(Guijarro 2019).b Sampling design exemplification.The picture shows an exemplificative site where two homogeneous stand are surveyed collecting data in two representative plots for each stand.Minimum required distances between plots and stands are also reported

Fig. 4
Fig. 4 Differences between diameter (h = 15 cm above the ground) means (unselected and selected trees).Data are not normally distributed, thus a non-parametric test is used.Medians are shown

Fig. 5
Fig. 5 Probability of beaver selection on woody plants with increasing distance from the riverbank for the three different taxonomic groups.Probability curves are based on the GLMM, considering the first ten meters of distance (i.e., where beaver activity was almost entirely observed) from the riverbank.Trunk diameters were transformed into two categorical classes based on the median value of gnawed trees (i.e., middle-low: ≤ 6 cm; middle-high: > 6 cm)

Table 1
Type of data collected and its explanation With 'Selected standing element' we meant standing elements with whatever type of tree damage caused by beaver gnawing activity Robinia pseudoacacia L.), are present both in Tevere and Merse river.

Table 2
Extent and characteristics of trees selected by beaver within each river stand