Relics of beavers past: time and population density drive scale-dependent patterns of ecosystem engineering Research

effects can be scale-dependent, indicating researchers should evaluate the ecological impact of engineers across diverse spatiotemporal scales to fully understand their functional roles in ecosystems.


Introduction
A central tenet of ecology concerns understanding how ecological patterns and processes scale and vary across space and time. Because there is no 'correct' scale that appropriately describes the variation or predictability of all ecological systems (Levin 1992), ecologists often incorporate or account for scaling factors in research (Wiens 1989, Jackson andFahrig 2014). For instance, changing the spatial extent or resolution at which landscape and community dynamics are evaluated reveals scale-dependent patterns in indices such as species diversity and richness , Schneider 1994, Chase and Leibold 2002, Rahbek 2004, due in large part to spatial heterogeneity within landscapes (Turner , 2005. Further, ecological processes often generate patterns at spatial scales that are different from scales at which the processes act (Levin 1992, Chave 2013. Numerous ecological concepts, such as metapopulation and metacommunity dynamics (Hanski 1998, Leibold et al. 2004, focus on evaluating how finer-scale processes are linked to broader-scale patterns (Wiens et al. 1993).
Natural disturbance dynamics provide a good example of how finer-scale processes can generate patterns at broader spatial scales (White 1979). Disturbance is recognized as a key ecosystem component that can structure ecological systems by driving community assembly (Levin andPaine 1974, Paine andLevin 1981), increasing landscape heterogeneity (Sousa 1984, Turner 2010 or even causing dramatic shifts in ecosystem stable states (Scheffer et al. 2001, Scheffer and Carpenter 2003, Ratajczak et al. 2018). Both small-and large-scale disturbances can contribute towards establishing a shifting mosaic of landscape patches with plant communities at varying successional stages (Bormann andLikens 1979, Pickett andWhite 1985). The degree of spatial heterogeneity within landscapes is inherently a function of the spatial and temporal scale of the disturbance events (Turner 2010) and the scale at which the processes are examined (Levin 1992).
Ecosystem engineers often create patch disturbances on the landscape that alter ecological processes across diverse spatial and temporal scales (Jones et al. 1994, Hastings et al. 2007), but little is known about how engineer-mediated effects scale or vary across space and time. We explore how ecosystem engineer-created patch disturbances affect ecological processes across space and time using beavers Castor canadensis as a study organism. Beavers are prolific engineers that impound and store surface water by building dams (Naiman et al. 1988, Gurnell 1998), a process recognized as a patch disturbance because the characteristics of the dammed site are radically altered (Johnston and Naiman 1987, Wright et al. 2004, Nummi and Kuuluvainen 2013, Hood and Larson 2015. Beavers are often the primary drivers of wetland creation in boreal forest ecosystems (Johnston and Naiman 1990a, Snodgrass 1997, Hood and Bayley 2008, Gable et al. 2020, and the accumulation of ponds (patch disturbances) through time influences the composition and dynamics of landscapes Naiman 1990a, b, Kivinen et al. 2020). The ecological effects that beaver ponds have on plants, animals and ecosystem processes at both the pond and landscape scales are well-documented (Naiman et al. 1994, Wright et al. 2002, Rosell et al. 2005, Johnston 2017, Willby et al. 2018.
Beavers live in multi-generational family groups ('colonies') that periodically occupy and abandon ponds through time, anywhere from one year to several generations (Hood 2020, Kivinen et al. 2020. The periodic (re)occupation of ponds, combined with the creation of new ponds, results in a mosaic of patches at varying stages of ecological succession (Naiman et al. 1988, Johnston andNaiman 1990b) that increases landscape heterogeneity (Remillard et al. 1987, Gable et al. 2020, Kivinen et al. 2020. Previous research has evaluated how this beaver pond mosaic affects organisms (Schlosser et al. 1998, Snodgrass and Meffe 1998, Schlosser and Kallemeyn 2000, Hood and Larson 2014, Bush et al. 2019, and there is evidence that both local and landscape drivers influence how beaver engineering affects species richness and abundance (Wright et al. 2002, 2003, Nummi et al. 2019b. Less is known about how the beaver pond mosaic influences ecological processes across multiple temporal and spatial scales (Brazier et al. 2021, Wohl 2021) and, in particular, how abandoned ponds continue to influence ecosystems (Gurnell 1998, Levine and Meyer 2019, Larsen et al. 2021. Beaver ponds can be remarkably long-lasting features (Johnston 2015) and may continue to influence riparian ecosystems when abandoned -for instance, by storing sediment and water upstream of the abandoned dam (Butler and Malanson 2005, Levine and Meyer 2014. Additionally, the long-term occupancy dynamics of beaver ponds at large spatial scales remains under-explored (but see Johnston 2015, Kivinen et al. 2020. Assessing long-term patterns in how occupied and abandoned ponds influence ecological processes will provide a more comprehensive understanding of the beaver's keystone role within ecosystems. In the present study, we investigated how beavers affect patterns in surface water dynamics across multiple spatial scales (pond, watershed and regional scales) over a 70-year period (1948-2017) within five watersheds in northeastern Minnesota, USA. We used aerial imagery to evaluate 1) temporal patterns of surface water area associated with beaver ponds, 2) the influence of beaver population density on surface water area, 3) whether the influence of population density on surface water area varies through time, 4) how surface water area is influenced by pond occupancy status (occupied/abandoned) and 5) whether patterns in surface water area vary across spatial scales. We focused on the influence of time and population density on surface water area within beaver-created impoundments because a recent large-scale analysis from Minnesota demonstrated changes in population density were predominantly driven by densitydependent mechanisms whereas environmental factors (precipitation, temperature, habitat quality) had minimal effect (Johnson-Bice et al. 2021a). Here, time is a representative measure for several underlying mechanisms of beaver pond dynamics (pond availability, suitability and rates of decay/ recovery) primarily governed by landscape and habitat attributes that may be independent of population size. Our work, given the spatiotemporal scale (five watersheds over 70 years), allowed us to evaluate general patterns of beaver-created surface water dynamics in a low-gradient system. Furthermore, this work is unique because it encompasses the recolonization, recovery and stabilization of beaver populations within these watersheds following decades of near-eradication, providing an excellent opportunity to evaluate how the recovery of an ecosystem engineer influences ecological patterns and processes across space and time.

Study area
We conducted our study along the northern shore of Lake Superior in northern Minnesota, USA (hereafter referred to as the 'North Shore'; Fig. 1). The Minnesota North Shore comprises a 5651 km 2 area and has topography that varies in elevation from 130 to 700 m. The lower segments of North Shore streams (as they enter Lake Superior) are bedrock-confined and characterized by steep slopes and ravines, so most beaver ponds are in the low-gradient headwaters where there is less topographic relief.
The North Shore includes 19 watersheds that drain separately into Lake Superior. We selected five of these watersheds comprising ~18% of the total North Shore area for our analysis: Cascade River (289 km 2 ), Kadunce River (122 km 2 ), Knife River (224 km 2 ), Manitou River (254 km 2 ) and Split Rock River (113 km 2 ) watersheds (hereafter referred to as the 'study watersheds'; Fig. 1). The study watersheds were selected specifically because they capture the variation in topography and habitat characteristics found throughout the North Shore. We assumed these watersheds were independent because each watershed drains separately into Lake Superior and the geographic separation between watersheds means there was likely little-to-no beaver dispersal among them.
Annual precipitation in the North Shore averages 70-75 cm, with nearly half occurring during the growing season (May-September; generally, 121-135 days). Terrestrial vegetation is largely composed of aspen Populus spp. and paper birch Betula papyrifera, though some sections of conifer bogs and mixed hardwood-pine forest remain despite extensive logging activity.
Like most beaver populations in North America, North Shore beaver populations were decimated during the fur trade era and the decades following (approximately between the years 1600 and 1850). Trapping was prohibited in Minnesota starting in 1909 (Longley and Moyle 1963)

Aerial imagery products
Historical aerial imagery has frequently been used to evaluate landscape-scale effects of beaver ecosystem engineering. We   Figure 1. Study area depicting the North Shore region and the five smaller watersheds that we selected for our analysis along the shoreline of Lake Superior in northern Minnesota, USA. compiled aerial photos from 1948, 1961 and 1980-1982 for each of our study watersheds. We also compiled photos for the Cascade and Kadunce River watersheds from 1934 and the Knife River watershed from 1939. Imagery pre-1948 was used only for context on the historical presence of beavers, while imagery from 1948 onward was used in our analysis of surface water dynamics (1948 is the earliest year imagery was available for all study watersheds; Supporting information). In total, 780 aerial photographs spanning 1934-1982 were scanned and georeferenced. We also used digital orthoimagery from five additional years (1991/1992, 2003, 2008, 2013, 2017) at a 1-m resolution scale for our analysis. We obtained the 1991/1992 imagery from the U.S. Geological Survey (black/white), while the imagery from 2003 to 2017 was acquired from of the U.S. Department of Agriculture National Agriculture Imagery Program (3-and 4-band color). All imagery products were acquired during late spring or summer of their collection year.

Digitizing surface water associated with beaver activity
We digitized all surface water features associated with beaver activity, whether the ponds were full of water or not. We considered a surface water feature to include both 'open water' and 'deep marsh' wetland classes (Eggers and Reed 2015) (i.e. where water could be visually detected from imagery; Supporting information). Additional details on digitizing methods and field validation of beaver features are provided in the Supporting information. Although this approach captures temporal changes in surface water storage, it underestimates the total amount of surface water stored within beaver ponds, as water is often within emergent wetland classes. Thus, our surface water values should be considered minimum estimates, while the trends should be accurate because all features were digitized using a consistent method by a single observer (the lead author). Due to the strong correlation between pond surface area and pond volume (Karran et al. 2017 reported an R 2 of 0.98), we assumed that trends in surface water area reflect trends in surface water storage. Beaver dams frequently alter the geomorphology of streams (Gurnell 1998), leaving modifications that can persist beyond the destruction of dams (e.g. widening of a stream; Supporting information). We likewise included surface water features associated with beaver-created changes to stream geomorphology. Surface water features along continuous stream reaches were bounded at the location where the stream narrowed to its normal width, or the first major bend if the stream did not noticeably narrow. Features within each pond were digitized manually, and the surface water area was calculated in ArcMap 10.6.1.

Beaver colony density index
Due to the absence of population survey data within all study watersheds, we derived an index of beaver population density from aerial imagery similar to the method used by Snodgrass (1997). We first identified all ponds that were 'occupied' by beavers as evident from aerial imagery. Ponds were considered occupied if they were completely full of water (i.e. water across the full length of the dam), as fully intact dams indicate that beavers are actively maintaining the dam (Snodgrass 1997, Johnston andWindels 2015). Occupied ponds that were separated by 50 m of unimpounded (free-flowing) stream reaches/segments were considered separate colonies. When > 3 adjacent occupied ponds were not separated by > 50 m of stream, we assumed that the number of colonies was half of the number of occupied ponds and we rounded down to the nearest whole number (e.g. five adjacent ponds = two colonies which collectively maintained five ponds; six adjacent ponds = three colonies which collectively maintained six ponds). We considered the presence of any unoccupied pond between occupied ponds to be a boundary between separate colonies. We used the number of identified colonies per km 2 within each watershed-year as our beaver density index for each imagery year. This approach likely does not reflect the 'true' density of beaver colonies, as beavers often occupy partially filled ponds or reside in free-flowing stream sections. Rather, this method provides an index of beaver colony density that was consistent across imagery years and between watersheds, and leverages the intricate relationship between beavers and the aquatic habitats they create. We use the terms colony density and population density interchangeably here, as the former has been repeatedly shown to be a good proxy for the latter (Johnson-Bice et al. 2021a).

Beaver pond site longevity
To better understand the long-term effects of beaver-created surface water dynamics, we investigated the number and proportion of surface water features identified from 1948 and 1961 aerial imagery that were occupied by beavers > 48 years later. Using aerial imagery from 2009, 2010, 2013, 2015 and 2017, we evaluated whether each feature from 1948 and 1961 imagery was occupied by beavers at the same general impoundment site at least once from 2009 to 2017 using the same criteria as our density index (i.e. ponds completely full of water).

Precipitation data
Previous research has demonstrated precipitation has little influence on total surface water area within beaver ponds (Hood and Bayley 2008) or rates of beaver population change (Johnson-Bice et al. 2021a). Nonetheless, to ensure precipitation did not influence total surface water area within beaver ponds at the watershed-scale we incorporated annual precipitation values into these models (e.g. during years with greater precipitation, more water might be stored within ponds across the landscape). We acquired monthly precipitation raster data (mm) at a 4-km scale from the PRISM Climate Working Group (PRISM Climate Group 2014). For each watershed and month, we calculated the average monthly precipitation value using the Zonal Statistics as Table tool in ArcGIS. Since most aerial photos were taken in late spring, we calculated annual precipitation as the sum of precipitation from June of the previous year to May of the year imagery was acquired. We evaluated both the 1-year and 2-year average precipitation, the latter of which was the most influential precipitation value in a similar analysis (Hood and Bayley 2008).

Statistical analysis
We fit our surface water data to four linear mixed-effects models implemented in the 'lme4' package (Bates et al. 2015) in R (ver. 3.6.3; <www.r-project.org>). We examined pondscale patterns by evaluating the influence of time and colony density on 1) the average area of surface water stored within all ponds (both occupied and abandoned), and 2) the average area of surface water stored within only occupied ponds.
To meet assumptions of normality, we log-transformed these data for the pond-scale analyses. We then examined watershed-scale patterns by evaluating the influence of time and colony density on 3) the density of ponds storing surface water (features km −2 ) within each watershed and 4) the total area of beaver-created surface water within each watershed (adjusted by watershed area; ha km −2 ). To evaluate whether beaver management in the Knife River watershed influenced the total amount of surface water, we performed a post hoc analysis by fitting the total surface water data to models excluding 4a) all data from the Knife River and 4b) data from just the last four time periods (2003-2017) after management began. To account for the potential influence of precipitation on total surface water area in beaver ponds, we also included precipitation as a covariate in models 4, 4a and 4b, and fit separate models for 1-year and 2-year average precipitation terms (i.e. 6 total fitted models covering a range of possible ways that beaver management and precipitation could have influenced total surface water area). Within each model, we included an interaction term between 'year' and 'density' to evaluate whether the influence of density varied through time -for instance, whether colony density influenced surface water dynamics differently as ponds accumulated on the landscape. We centered the data around a mean of 0 for each model and used 'watershed' as a random intercept term to account for variation among watersheds (exploratory analyses revealed including random slopes did not significantly affect parameter estimates). Although numerous ponds were measured in multiple imagery years, we did not include a random effect term for 'pond ID' due to inherent difficulties in assigning consistent and statistically/ biologically meaningful ID values for each pond. Ponds in our study area constantly shifted between periods of occupation and abandonment and were often divided up via new beaver dams or combined with other ponds through time (Supporting information). Any rules we developed to consistently assign pond IDs would add subjectivity and bias to our models and thus we avoided this approach.
We assessed collinearity by calculating variance inflation factor (VIF) values for each variable; VIF values were < 2 for all variables, indicating no collinearity (Dormann et al. 2013). We considered modeled trends to represent regional-scale patterns in surface water patch dynamics, since our study watersheds were selected specifically because they characterize the environmental variation found throughout the North Shore.

Results
We digitized 11 609 surface water features associated with beaver engineering from aerial imagery spanning 1948-2017 across our five study watersheds and 92 features from imagery pre-1948 that were used as context for historical beaver presence.
Based on our index of beaver colony density, few beavers were present in the North Shore in the 1930s (Supporting information). Although patterns were variable among watersheds, beavers were generally recolonizing watersheds by 1948, and populations had started to stabilize by 1961 and 1980-1982. Beaver populations throughout the North Shore region appear to have remained stable since then, albeit with asynchronous density fluctuations among watersheds (Supporting information).

Beaver pond site longevity
The general pond sites from approximately half of all surface water features digitized from aerial imagery in 1948 (45.9%, n = 241) and 1961 (50.7%, n = 561) were occupied and maintained in at least one imagery year between 2009 and 2017 (Supporting information).
Despite the increase in surface water features within watersheds, the cumulative area of beaver-created surface water (adjusted by watershed area; ha km −2 ) remained remarkably stable through time (β = 0.000, CI = [−0.001, 0.001]) ( Table  1, Fig. 3d). That is, surface water features decreased in average area but increased in density, resulting in little change in total surface water storage through time. Beaver colony density was strongly related with total surface water area (β = 0.854, CI = [0.750, 0.958]) (Table 1, Fig. 3e), suggesting the region-wide stability of surface water was influenced by variation in beaver colony densities within and among watersheds through time. The interaction between year and density was not significant (β = −0.001, CI = [−0.006, 0.004]). Results from the post hoc models evaluating the influence that beaver management in the Knife River watershed had on watershed-scale surface water area indicated that management practices did not affect the overall (i.e. regional-scale) trend (Supporting information). Precipitation (both 1-year and 2-year average) did not influence total surface water area associated with beaver ponds in neither the full nor any post hoc models (Table 1, Supporting information).
The proportion of cumulative surface water that was stored within abandoned ponds (i.e. ponds that were not full of water) increased through time (Fig. 4), such that 54.0% of beaver-created surface water area in 2017 was stored in abandoned ponds compared to 24.6% in 1948 (pooled across all watersheds).

Discussion
By examining beaver-created surface water dynamics over 70 years across multiple spatial scales (pond, watershed and regional scales), we demonstrate how the interrelated effects of time and population density result in scale-dependent patterns of beaver ecosystem engineering. These scale-dependent patterns appear to be driven by a combination of densitydependent population dynamics and underlying temporal patterns of pond availability, suitability and decay/recovery. Our results suggest that the legacy effects (i.e. ecological effects that last beyond the lifetime of the organism) left by beaver ecosystem engineering play a pivotal role in beaver-created surface water dynamics. Given the spatiotemporal extent of the data, our study provides novel insight into how ecosystem engineers can influence large-scale ecological processes across space and time during and after population recovery.

Beaver-created surface water dynamics across space and time
Several temporal patterns in beaver-created surface water dynamics at both the pond-and watershed-scale were similar across all study watersheds. The influence of time on beaver-created surface water dynamics was likely driven by the availability, suitability and rates of decay/recovery of beaver ponds, which were likely also driven by beaver population dynamics (Wright et al. 2004). As expected, the number of ponds storing surface water increased as North Shore beaver populations recovered. The sustained increase in the density of ponds storing surface water reflects two important aspects of beaver pond dynamics: 1) the ability of abandoned ponds to store surface water, and 2) the sustained creation of new ponds. As beaver ponds are created and abandoned, there is a steady accumulation of abandoned ponds on the landscape that store surface water (Fig. 3a), which effectively lowers the average area of surface water across all ponds (Fig. 2a).
The fact that beaver ponds continue to be constructed several decades after reaching population stabilization indicates that either the rate at which beavers deplete resources exceeds the recovery rate of abandoned ponds, that beavers prefer unmodified sites, or both (Wright et al. 2004). This also likely explains why surface water features accumulated faster at high beaver densities (Fig. 3c) (i.e. with more beavers present, more ponds are created). The decline through time in the average area of surface water within occupied ponds reflects another reported aspect of beaver pond dynamics -their preference during colonization for optimal sites with a large surface water area to dam length ratio Naiman 1990a, John et al. 2010). Once most optimal sites have been colonized, beavers that choose to construct new dams must do so at progressively less optimal sites that, on average, store less surface water (Fig. 2c). In other words, the longer engineers (beavers) occupy the landscape, the more the structures created by their predecessors constrain the ecological impact that existing engineers have. The influence of beaver colony density on surface water dynamics depends on the scale at which it is assessed. Beaver density clearly positively influenced total surface water area at the watershed scale (Fig. 3d), confirming beavers are biotic agents that increase surface water across landscapes (Hood and Bayley 2008, Puttock et al. 2017, Jones et al. 2020. However, at the pond scale, beaver density was inversely related to average surface water area (Fig. 3e)  8 interaction between time and density on the area of surface water within occupied ponds was not significant, the rate of decrease in pond-scale surface water area tended to be greater at higher densities (Fig. 2f ). By synthesizing results from our study with findings from previous research, we can provide a comprehensive narrative of beaver-created surface water dynamics through time within lower gradient systems. Beavers preferentially select larger pond sites during colonization Naiman 1990a, John et al. 2010) that store more surface water on average. But as these early ponds degrade and beavers continue to recolonize, beavers must build dams at sites that store less water (i.e. ponds get progressively smaller; Fig. 2d) -even building dams within large ponds, which creates multiple smaller ponds (Fig. 5, Supporting information). Competition for pond sites exhibits a density-dependent response such that beavers in higher-density populations occupy ponds that store less water on average (Fig. 2e), a response that may increase through time as beaver populations stabilize (Fig. 2f ). However, despite our finding that occupied ponds stored less water through time, the influence that beaver density had on total surface water area at the watershed scale (Fig. 3e) did not change across the study period (Fig. 3f ). If beavers occupy progressively smaller ponds through time, how can total beaver-created surface water remain stable if beaver density is held constant? The answer lies in the accumulation of abandoned ponds, which through time store a greater proportion of the total beavercreated surface water (Fig. 4). These relics of past beaver activity (which accumulate at different rates depending on beaver density; Fig. 3c) appear to offset the smaller size of occupied ponds, which results in the temporally constant relationship between beaver density and watershed-scale water storage. This study provides an excellent example of how organism-created patches can have legacy effects on ecosystems by continuing to affect ecological patterns and processes beyond the lifetime of the engineer (Hastings et al. 2007, Cuddington 2011. Continued long-term monitoring will provide insight into whether the trajectory of these trends in beaver pond dynamics will continue or reach an asymptote at some point. Several methods have been used to estimate beaver density from metrics derived from aerial imagery when other population information is absent. For instance, beaver dam counts have been used as a proxy for population density from studies in Michigan and Alaska, USA  Figure 3. The influence of time and beaver colony density on watershed-scale patterns of beaver-created surface water dynamics. An asterisk next to the panel letter indicates the 95% confidence intervals of the coefficient estimate did not overlap 0. (a) and (b) depict the influence of time (year) and colony density on the density of ponds storing surface water, while (d) and (e) depict the influence of time and density on the cumulative area (adjusted by watershed area; ha km −2 ) of beaver-created surface water. Colored lines represent the average linear slope for each watershed, while the black dashed line and grey ribbon is the linear mixed effects model slope and 95% confidence interval. (c) and (f ) depict the estimated interaction between time and beaver density on the density of ponds storing surface water (c) and the cumulative area of beaver-created surface water. The interaction is shown at three density levels -low, medium and high -where medium represents the mean average density (0.33 colonies km −2 ) and low/high represent ± one standard deviation (SD = 0.16). (Martin et al. 2015, Jones et al. 2020). The method we used assumed ponds that were full of water were occupied by a beaver colony that was actively maintaining the pond's dam and therefore impounding water. Our density index results (Supporting information) are consistent with other research from Minnesota that demonstrated rapid population recovery and recolonization from the 1940s into the 1960s/1970s followed by a prolonged period of population stabilization (Johnston and Naiman 1990a, Johnston and Windels 2015, Johnson-Bice et al. 2021a, indicating this method was appropriate for our study area. However, one potential issue with the index we used is that it would not account for the possibility that individual beaver colonies may have occupied more ponds, on average, through time as occupied pond size declined. To assess this issue, we evaluated whether the proportion of occupied ponds comprised of single-pond colonies changed through time but found no change throughout the study period (Supporting information). These results, in combination with the agreement in density trends between our study and elsewhere in Minnesota, indicate our index of beaver colony density was reliable. Although this method relies on information partly related to surface water metrics, this is a byproduct of the inherent relationship, and resultant feedbacks, between ecosystem engineers that create their own habitat and the ecosystem (Jones et al. 2010). We note that our study was conducted within a low-gradient area that generally experiences little-to-moderate variation in precipitation, and that our results may not be directly transferrable to other beaver-modified environments. Beaver dams generally fail and degrade quicker in higher gradient systems (Gurnell 1998), suggesting that the legacy effects from abandoned dams may not last as long in these systems. Additionally, areas with large variation in precipitation may exhibit greater variation in landscape-scale water storage than we found here. Beaver recovery is ongoing in many mountainous and dry regions, and the long-term monitoring of beaver engineering from these other systems will provide greater insight into how generalizable our results are throughout the beaver's range.   . Visual representation of general beaver-created surface water patterns through time in low-gradient systems. During recolonization (a), beavers selectively build dams in areas that result in larger ponds that store more surface water. As beavers expand on the landscape (b), they progressively create smaller ponds as availability of larger sites declines. At this stage, beavers also often build dams within large ponds (top pond, b), creating multiple, smaller ponds in the process. Although beavers may abandon ponds, dams at abandoned sites nonetheless often store surface water (middle pond, b). Once beaver populations have stabilized (c), density-dependent mechanisms and patterns in pond site availability and decay/recovery drive beavers to build still smaller ponds through time. The steady accumulation of ponds generates a mosaic of patches at various stages of ecological succession, ultimately resulting in a steady increase in the total number of ponds storing surface water and an increase in the proportion of beaver-created surface water that is stored within abandoned ponds.

Implications for rewilding, conservation and management efforts
Restoring beavers to areas they have been extirpated from has received attention as a natural framework for achieving ecological restoration objectives. Beaver-assisted restoration has been touted as a tool for repairing incised and degraded stream systems throughout the western USA and parts of Eurasia (Burchsted et al. 2010, Pollock et al. 2014), while others have advocated for using beavers to rewild ecosystems (Law et al. 2017, Willby et al. 2018) and facilitate wetland biodiversity (Nummi and Holopainen 2014, Law et al. 2019. We suggest beaver engineering could be used to increase surface water storage and promote freshwater conservation efforts in forested ecosystems. If given the opportunity to recolonize landscapes, beaver-created surface water can reach capacity levels within a relatively short period of time (< 20-40 years), and moreover, remain important aquatic features on low-gradient landscapes that often persist for decades. Our study supports previous work that demonstrated precipitation had little influence on surface water stored in beaver ponds at the landscape-scale (Hood and Bayley 2008). Even if beaver densities remain stagnant or decline through time, the accumulation of abandoned ponds that continue to store surface water can preserve important ecological functions like mitigating drought effects (Hood and Bayley 2008), offering refugia during wildfires (Fairfax and Whittle 2020) and providing essential aquatic habitats to the numerous other species that thrive in beaver-altered environments (Cunningham et al. 2007, Karraker and Gibbs 2009, Nummi and Holopainen 2014, Nummi et al. 2019a. Slowing or retaining water on the landscape is also a common goal of ecological restoration efforts designed to reduce flooding and export of sediment and nutrient into receiving water bodies (LCCMR 2008). Comparatively less is known about the degree to which abandoned dams in other systems (e.g. mountainous or arid environments) are able to similarly preserve ecological functions.
The influence of time and population density on beavercreated surface water dynamics has several conservation and management implications. Because beaver density is positively related with surface water area, reducing beaver density via removal or control efforts will reduce surface water storage across landscapes (Hood and Bayley 2008). However, the extent of surface water reduction appears related to the duration that beavers have occupied the landscape: the longer beavers occupy the landscape, the less of a disturbance removing one colony may have on surface water dynamics since beavers occupy progressively smaller ponds through time. The reduced impact that beaver control has on regional-scale surface water area through time may explain why control efforts in the Knife River watershed had little influence on surface water dynamics. But focusing solely on the impact of management efforts on surface water area overlooks the longevity of beaver-driven surface water patches in low-gradient systems, many of which undergo occupation and maintenance even > 48 years after their creation (Supporting information). A smaller-scale study from the Upper Peninsula of Michigan revealed beavers may occupy ponds even 150 years after creation (Johnston 2015) but our multi-watershed analysis demonstrates the occupancy of old pond sites is pervasive among beaver populations. For areas where invasive beavers are drastically altering local environments, such as in Patagonia (Anderson et al. 2006, Anderson and Rosemond 2007, Huertas Herrera et al. 2020, Jusim et al. 2020, our study indicates that extensive dam removal efforts may be required to restore ecosystem function after beavers are eradicated, since the legacy dams left behind will likely continue to influence ecological processes.

Scale-dependent patterns of ecosystem engineering
The classic perspective on patch disturbance dynamics suggests that natural disturbances at one spatial scale may produce patterns at another scale (Bormann andLikens 1979, Pickett andWhite 1985). Individual patch disturbances may vary in their successional stage from one another, but at larger spatial scales the assemblage of patches results in a disturbance mosaic that is a key driver of landscape heterogeneity (Turner 2010). Although most commonly evaluated with natural disturbances such as fire or windthrow events, beaver-created disturbances are also recognized as important drivers of spatiotemporal landscape heterogeneity (Remillard et al. 1987, Kivinen et al. 2020. Here, we found that the increased heterogeneity provided by the beaver pond mosaic can also generate patterns in ecological processes at larger spatial scales. Despite large variations in watershedscale surface water, the total amount of surface water area within beaver ponds across all watersheds remained remarkably stable through time (Fig. 3d). By evaluating pond-and watershed-scale patterns in relation to time and population density, we conclude that this regional-scale stability was driven by the interaction of density-dependent mechanisms at the pond and watershed scales and legacy effects from the steady accumulation of abandoned ponds. Thus, organismcreated small patch disturbances not only conform to classic patch disturbance dynamics but it appears that the organismcreated disturbance mosaic can also generate stability in ecological processes at larger spatial scales.
A comprehensive understanding of how organisms modify their environment requires evaluating whether their ecological effects scale or vary across space and time. The effects of ecosystem engineering may vary between sites Jones 2006, Baker et al. 2012) or may change in magnitude in response to environmental variation (Phillips et al. 2019). Ecosystem engineers may also have temporal or spatial scaledependent effects on other organisms. For instance, the impact of ecosystem engineering by Arctic foxes Vulpes lagopus on other wildlife depends on whether the engineering pathways operate on short or long time scales (Zhao et al. 2021), while elephant Loxodonta africana foraging affects patch-scale site selection of small vertebrates but landscape-scale habitat selection of large vertebrates (Pringle 2008, Valeix et al. 2011. The fact that ecosystem engineers can have context- (Padilla 2010) and scale-dependent effects (Ferry et al. 2020) is an important yet often over-looked aspect of the ecosystem engineering concept. Our study clearly demonstrates that, for organisms with keystone roles within ecosystems, our understanding of the patterns and magnitude of their ecological impact is enriched when evaluated across time and multiple spatial scales.

Supporting information
The supporting information associated with this article is available from the online version.