Local disturbance by muskrat, an ecosystem engineer, enhances plant diversity in regionally-altered wetlands

. Biotic ecosystem engineers are increasingly recognized as important drivers of biodiversity, structure, and function in many ecosystems. By regulating physical processes and creating local disturbances, ecosystem engineers can serve as important elements of passive habitat restoration, as they continuously alter and shape their environments. Native to North American temperate wetlands, the common muskrat ( Ondatra zibethicus ) is an ecosystem engineer that alters its environment through herbivory, house and structure building, and creation of activity networks. While the physical consequences of muskrat disturbances are well-known, few studies have quanti ﬁ ed their impact on the wetland biotic community. We conducted a ﬁ eld study to investigate the effects of muskrat herbivory and structure building on plant biodiversity in wetlands along the upper St. Lawrence River (New York, USA) that have been extensively modi ﬁ ed by long-term water regulation and non-native cattail invasion. Plant species richness and diversity were greater in muskrat-disturbed areas compared to at-large reference plots within the same wetlands. Soil saturation mediated the biodiversity impacts of muskrats, with intermediate moisture levels resulting in the highest species richness. These interacting drivers, muskrat activity and hydrology, had a compensatory effect on plant biodiversity loss associated with non-native cattail invasion. Multivariate analysis indicated a distinct plant community associated with muskrat disturbances. Analysis of spatiotemporal patterns of house locations suggests that some muskrats reuse house locations in sequential years, likely amplifying the intensity and duration of their impacts. This study demonstrates that muskrat disturbance in ﬂ uences wetland plant diversity at scales relevant to regional drivers of plant diversity.


INTRODUCTION
Ecosystem engineers are organisms that alter their surrounding environment to affect the quantity and quality of available resources (Jones et al. 1994, Gutierrez et al. 2003. By changing the amount and distribution of available resources in an area, ecosystem engineers often create unique habitats otherwise absent in the environment, thus opening up additional ecological niches that can support a more diverse species assemblage (Gutierrez et al. 2003, Law et al. 2017. In wetlands, the presence of ecosystem engineers can speed up, slow down, or reset ecological transitions, by creating disturbances of varying intensity and duration. In systems where physical drivers have been strongly modified, ecosystem engineers may enhance restoration efforts and ecological recovery (Curran andCannatelli 2014, Law et al. 2017). This can be especially important where altered physical conditions have decreased the competitive advantages of native species and facilitated biological invasions (Hecky et al. 2010, Kovach et al. 2015. In the wetlands of North America, the North American beaver (Castor canadensis) and the common muskrat (Ondatra zibethicus) are two prominent mammalian ecosystem engineers (Higgins and Mitsch 2002, Mott et al. 2013, Curran and Cannatelli 2014, Rozhkova-Timina et al. 2018. A highly adaptable species, muskrats are relatively common throughout their native range in temperate North American wetlands, as well as in Europe and Asia where they are invasive and modify native ecosystems and human water control infrastructure (Cotner andSchooley 2011, Rozhkova-Timina et al. 2018).
Muskrats are medium-sized, semi-aquatic rodents, which primarily feed on robust emergent hydrophytes and occasionally on small animals such as freshwater mussels (Willner et al. 1980, Allen and Hoffman 1984, Ecke et al. 2014, Edelman et al. 2015. Like most rodents, muskrats can reproduce quickly and can have two to three litters per season with 4-8 kits each (Alexander 1956, Willner et al. 1980. They build a variety of structures using locally available plant material and substrate, including free-standing houses and bank burrows, which serve as refuge from predators and exposure to the elements (Dozier et al. 1948, Danell 1978.
Under optimal conditions, muskrat herbivory can remove large quantities of vegetation in a short period, creating extensive impact areas called "eat-outs" (Allen andHoffman 1984, Higgins andMitsch 2002). Muskrat herbivory can drastically reduce the biomass and cover of dominant emergent hydrophytes such as cattails (Typha spp.), because the animals remain active throughout the winter, continuously foraging for succulent roots and shoots of perennial plants (Alexander 1956, Willner et al. 1980, Toner et al. 2010. Ecological studies have documented the effects of muskrat activity on enhancing habitats for waterbirds (e.g., Chlidonias niger; Kiviat 1978), altering water chemistry (Connors et al. 2000), changing macroinvertebrate assemblages (De Szalay andCassidy 2001, Nummi et al. 2006), and increasing biomass of non-dominant plants atop their houses (Kangas and Hannan 1985). However, few quantitative studies have investigated how muskrat disturbance interacts with other wetland ecosystem drivers to influence the biotic community as a whole (Danell 1978, Kadlec et al. 2007, Skyrien _ e and Paulauskas 2013.
Muskrats' capacity to remove large areas of robust emergent vegetation may create new ecological niches by changing the abiotic and biotic conditions in their habitats. Removal of large, canopy-forming, robust emergent hydrophytes would effectively reduce the interspecific competition among wetland plants in conjunction with increased availability of light and space. The construction of houses by muskrats creates topographic heterogeneity in the landscape, specifically locally elevated habitats that may allow less flood tolerant and more facultative wetland plants to occur. Aside from the direct impact on the wetland plant assemblage, these changes in environmental qualities could have emergent effects on the wetland community as a whole through the long-term and cumulative local disturbances present in wetlands with high muskrat densities.
Despite muskrats' ubiquity and ecological importance, few studies have investigated and quantified their capacity as ecosystem engineers to influence the broader ecological community, including wetland plant assemblages. As such, muskrat-mediated changes in plant species composition and species diversity are poorly understood, particularly in respect to how muskrat disturbance interacts with other ecological drivers such as hydrology and biotic invasions. This question is particularly pressing with the introduction and spread of invasive cattails throughout large parts of the muskrats' natural habitat in the North American Great Lakes region (Jansson et al. 2000, Wilcox et al. 2007, Tuchman et al. 2009, Freeland et al. 2013. One area of concern in North America is the Thousand Islands (TI) archipelago of the St. Lawrence River (SLR), where decades of water level stabilization in conjunction with the expansion of the invasive cattails has led to degradation of wetland habitat throughout the region (Jansson et al. 2000, Wilcox et al. 2007, Keddy 2010, Roseman et al. 2014. The resulting loss of native sedgedominated wetlands in this biome is especially concerning as these habitats support high biodiversity, primarily through the creation of heterogeneous microhabitats (Peach and Zedler 2006).
In this context, muskrats' ability to remove large amounts of vegetation, primarily cattails, allows them to act not only as a biological control agent, but also as an ecosystem engineer. Ecological theories of disturbance and community assembly suggest that the hydrological and ecological modifications to wetlands by muskrats would increase biodiversity and form a more resilient plant community (Brown 1995, D eri et al. 2011, Larkin et al. 2016. While incorporation of ecosystem engineers in habitat restoration plans is relatively new, recent studies have shown their stabilizing effect on restored habitats (Law et al. 2017, Bakker andSvenning 2018).
In this study, we investigated how local disturbance associated with muskrat activity influences wetland plant communities. We took a comparative approach, contrasting areas around muskrat houses with local reference sites in five wetlands of the Upper SLR system. Our principle objectives were to: (1) quantify the differences in plant community diversity and species composition between wetland sites used by muskrat and adjacent undisturbed areas; (2) model the complex relationships between disturbance, resource supply, and biotic competition via interaction effects of muskrat activity, wetland hydrology, and cattail litter on biodiversity; and (3) analyze the spatial and temporal distribution of house locations to assess the geographic scope and potential long-term modifications to the ecosystem.
We hypothesized that muskrat disturbances would increase plant biodiversity by reducing the extent of the dominant canopy-forming species (in particular Typha spp.), thereby increasing habitat heterogeneity and available growing space for other species. A further question was whether muskrat disturbance would shift the plant community toward more wetland obligate or alternatively terrestrial species, and how local conditions of hydrology may mediate that shift (Reed Jr. 1988). Finally, we hypothesized that high densities of muskrat, whether due to natural site factors or due to wetland water level manipulation, should result in landscape-scale changes. We tested this idea by evaluating the spatiotemporal persistence of muskrat houses in a managed wetland with water levels manipulated to prevent artificial winter drawdown of the SLR hydrosystem.

Study sites
The TI region of the Upper SLR is a boundary water that separates northern New York State (NY) from Canada's Ontario and Quebec provinces and spans approximately 78 km of the river's length from the Lake Ontario outlet to Morristown, NY (Fig. 1). This reach of the SLR forms part of the St. Lawrence Seaway, which spans the waters between the Gulf of St. Lawrence and Lake Superior. The water levels of the TI region and the seaway are regulated by policies of the International Joint Commission (IJC) via the Moses-Saunders Power Dam at Massena (Roseman et al. 2014).
Currently, multiple anthropogenic stressors impact the TI region including agriculture runoff, human alterations to the SLR hydrologic regime, non-native species invasions, and pollution from industry and maritime trade (Roseman et al. 2014). River regulation following the 1958 completion of the Moses-Saunders Power Dam has dampened annual and interdecadal hydrologic fluctuations and lowered the river's average winter water level , Roseman et al. 2014. While a new IJC water regulation plan (International Joint Commission 2014) is intended to restore some of the region's natural hydrologic variability, over five decades of water level stabilization in the SLR has facilitated the expansion of invasive cattails throughout the region's coastal wetlands (Wilcox et al. 2008). The spread of hybrid Typha 9 glauca Godr., which is a cross of native broad-leafed cattail (Typha latifolia) with the non-native, invasive narrow-leafed cattail (Typha angustifolia), has led to a substantial decrease in ecologically important sedge meadow habitats (Galatowitsch et al. 1999, Jansson et al. 2000, Wilcox et al. 2008, Olson et al. 2009, Freeland et al. 2013, Szabo et al. 2018). T. 9 glauca exhibits hybrid vigor, which allows it to outcompete and displace T. latifolia throughout much of its native range (Pieper et al. 2018, Szabo et al. 2018. v www.esajournals.org Five sites were surveyed for this study, spanning two broad wetland systems, French Creek and Cranberry Creek (Fig. 1). Both wetlands have drowned river mouth geomorphology connected to the SLR and share a similar ecological community structure and invasive Typha dominance . The ubiquitous effect of the IJC-managed regulation influences nearly 7000 ha of coastal wetlands and includes Cranberry Creek and other sites in French Creek . Historically, Cranberry Creek and part of French Creek (Carpenters Branch) had local water control weirs as part of a field experiment by Farrell et al. (2010). While the water control structure at Cranberry Creek had been removed in 2003, the local water level at Carpenters Branch in French Creek is still experimentally managed via a passive weir to prevent artificial fall drawdowns and maintain higher winter water levels relative to the IJCmanaged Moses-Saunders Dam-regulated water levels of the SLR and Lake Ontario. The locally water-controlled wetland has been shown to maintain a higher muskrat house density than the other sites , Toner et al. 2010, Kua 2019) and allows assessment of disturbance effects relative to reference areas across a greater range of muskrat activity levels and hydrologic conditions.

Field sampling
To assess the local influence of muskrat activity on the wetland plant community, we conducted vegetation surveys in the TI region in v www.esajournals.org July and August 2018 during the peak growing season, comparing diversity and community composition in the vicinity of muskrat houses, which are disturbed areas, to relatively undisturbed reference areas (Table 1). To establish a baseline condition, a site-wide assessment of wetland plant diversity ("reference" plots) was conducted using variable length transects that spanned an elevation gradient from open water to the upslope wetland boundary. The sampling approach was designed to encompass the range of environmental heterogeneity within the wetlands. Plots were sampled with individual quadrats (1 m 2 ) set along each transect at fixed elevation increments (0.5 m). Elevations as specified by the International Great Lakes Datum (1985) were determined from known-elevation benchmarks established using Trimble GeoXH 5800 at each plot location . The percent cover of all plant species observed in the quadrats was estimated to the nearest percent by trained field technicians. Wetland plants were identified to the genus level at minimum, and species richness and Shannon's diversity index were calculated for each quadrat. In addition, we characterized substrate moisture status as either "Dry," "Saturated," or "Flooded." Plots dry to touch were considered "Dry," plots wet to touch and/or with standing water <10 cm were categorized as "Saturated," and plots inundated to >10 cm were classified as "Flooded." Live Typha cover and dead litter cover were estimated (Kua 2019). Across both French Creek and Cranberry Creek systems, a total of 133 reference plots were surveyed along 26 transects (Table 1).
To assess muskrat impacts, the local area around selected houses ("muskrat" plots) was surveyed with quadrat sampling along 5-m length transects. Five meters corresponds with a radio-telemetry study on winter muskrat movement (MacArthur 1978). Twenty-eight muskrat houses were selected at random from a list of active houses detected during the prior winter house count survey (January-February 2018). The average area of the muskrat houses (as ellipses) was 1.9 m 2 AE 0.9 (AESD) with an average height of 61.8 cm AE 17.4. Using a square PVC frame, three 1 m 2 quadrats were sampled at 1 m spaced intervals along each transect with an additional quadrat centered on the muskrat house (Appendix S1: Fig. S1). Muskrat activity at each location was verified from field observations of well-maintained house structure, herbivory (i.e., signs of chewed plant material or cut plant stems) or cleared rivulets which signified muskrat runs. Data collection protocols were identical to those for the reference quadrats. Across both wetland systems, a total of 184 muskrat plots were surveyed at 28 separate house locations (Table 1).

Statistical analyses
Muskrat disturbance and habitat structure, plant diversity, and community composition.-Typha spp. live cover between muskrat and reference plots were compared to test the difference in Typha live cover (abundance) between these two habitats. Species richness and Shannon's diversity indices from vegetation plot data were compared to test for differences between muskrat and non-muskrat areas (i.e., a disturbance effect), and a multivariate community analysis was performed to explore the compositional dissimilarity between muskrat-disturbed and reference transects.
We used a mixed modeling approach to estimate the effects of active muskrat disturbance Note: Reference sampling units were variable length transects with plots located at fixed elevation increments of 0.5 m from the mean water level datum (IGLD 85). Muskrat transects indicate number of house locations surveyed. For both muskrat and reference transects, plots consisted of 1 9 1 m square quadrats. v www.esajournals.org and two environmental covariates, quadrat moisture class (moisture), and Typha spp. litter cover (litter), on vascular plant species diversity and richness. Wetland site was nested within wetland system (French Creek or Cranberry Creek) as a random effect to account for potential spatial covariance across all of the surveyed sites. We fitted a generalized linear mixed model to the species richness data. The richness counts were right-skewed as demonstrated by a Kolmogorov-Smirnov test of normality (P < 0.001; KS test); therefore, a negative binomial model with log link function was applied to the data after comparison with other distributions (e.g., Poisson). For the species diversity model, which also accounts for species evenness, we used a linear mixed model, as the KS test of the data distribution found no significant deviation from normality (D = 0.046; P = 0.11).
Typha spp. live cover differences were compared through similar mixed modeling approaches described above. Wetland site nested within wetland system was used as the random effect with muskrat disturbance and quadrat moisture class (moisture) as the fixed effects. KS test suggested that Typha spp. live cover was right-skewed and not normal (P < 0.001). Therefore, a Gamma model with inverse link function was fitted after comparison with other distributions (e.g., Beta).
For all three response variables (Typha cover, species richness, and evenness), a full model containing main effects and their two-way interactions was reduced and optimized through a candidate model selection process using the Akaike information criteria values corrected for small sample sizes (AIC c ). Following Burnham and Anderson (2004), we considered models with DAIC c > 2 to indicate substantial differences in support. Tukey-adjusted estimated marginal means (EM means; Lenth 2019) were used to test the pairwise differences between groups while accounting for the unequal plot sample sizes for species richness and evenness.
We used non-metric multidimensional scaling (NMDS) to explore the community differences between muskrat-disturbed and reference areas. Species percent cover estimates in each quadrat were averaged at the transect level for reference plots and by house location for muskrat plots and subsequently standardized using the Wisconsin double standardization (Bray and Curtis 1957). Bray-Curtis dissimilarity indices were calculated, and the effect of muskrats was quantified through a permutational analysis of variance (PERMANOVA). The variance (i.e., dispersion) between disturbance groups was compared through a permutation test of multivariate homogeneity following procedures from Anderson (2006). Additionally, a Kruskal-Wallis test (KW test) was used to test for differences in standardized plant species cover values between muskrat-disturbed and the reference sites.
The NMDS and subsequent community analyses were repeated after excluding the muskrat house quadrats to help us better understand the influence of the house structure on the local plant community assemblage.
Muskrat house density and distribution.-Muskrat house counts were conducted in Carpenters Branch (the water-controlled site) during the winters of 2017, 2018, and 2019 to better understand the impact and effect size of muskrat activity on regulated wetlands. We were interested in both the density of houses and their relative distribution from year-to-year, as indications of the spatial scope and longer-term (i.e., interannual) impacts of muskrat disturbance on their local environment. These data were pooled with historical house count surveys conducted between 2001 and 2012 (Toner et al. 2010;J. M. Farrell, unpublished data). While the historical data had irregular lapses in annual surveys, they were useful to compare house density trends over time (Fig. 2).
We used these pooled data to analyze the spatial-temporal correlation of muskrat house distribution in Carpenters Branch, for eight survey years between 2001 and 2019. A Knox test was conducted on muskrat house locations with a critical spatial distance of 0.0001 decimal degrees (approximately 7.87 m at 45°N) and recurrence lag of one year to quantify and compare the observed spatiotemporally "close" house pairs with their expected value, under the null hypothesis of no space-time interactions (Knox and Bartlett 1964). The expected value was calculated through a Monte Carlo (MC) permutation approach that compared distances for all possible pairings of houses in the data set (Mantel 1967). With this approach, a significant difference between observed and expected house pairs v www.esajournals.org would indicate that a spatial-temporal correlation for muskrat houses exists, as would be expected if houses were reoccupied in sequential years.
An alpha value of 0.05 was used to determine significance in all tests apart from the Kruskal-Wallis test (a = 0.01). All analyses were conducted using R 3.6.0 (R Core Team 2019) with packages car (

Muskrat disturbance and habitat structure
The best model selected to predict Typha spp. live cover included both muskrat disturbance and moisture as fixed effects without the interaction effect (Table 2). Results from the EM means comparison showed that Typha spp. live cover was primarily influenced by substrate moisture. In particular, Typha spp. cover in quadrats with "Saturated" moisture had significantly higher percent live cover (Fig. 3). While the absolute percent of live Typha spp. cover was lower in muskrat-disturbed quadrats compared to reference areas, the EM means comparison indicated that the differences were not statistically significant (t = 1.899; P = 0.058; Fig. 3).

Muskrat disturbance and plant diversity
Plant species richness and diversity were influenced by both muskrat disturbance and the environmental factors of substrate moisture and Typha litter cover. The best model predicting species richness included the three main effects for disturbance, litter, and moisture and the v www.esajournals.org interactions of Typha spp. litter 9 disturbance and moisture 9 litter (Table 3). Comparison of the EM means showed that species richness was significantly higher in the habitat around muskrat houses than in reference sites for both "Saturated" and "Dry" quadrats but was similar in "Flooded" quadrats ( Fig. 4). Typha spp. litter cover had a negative effect on species richness (Fig. 5). The disturbance 9 moisture interaction showed a significant difference in magnitude of the effect of disturbance under different moisture classes, whereas the disturbance 9 litter interaction showed a difference in the loss rate of species richness with increased litter under each disturbance (Fig. 5).
The best model predicting species diversity included disturbance, moisture, and their interaction with litter included as a covariate (Table 3). Estimated marginal means showed that within disturbances, there were no significant differences in Shannon's diversity. However, Shannon's diversity was higher around muskrat houses in both "Dry" and "Saturated" quadrats but was not significantly different in "Flooded" quadrats ( Fig. 4). Similar to species richness, Shannon's diversity decreased as Typha spp. litter cover increased (Fig. 5). The interaction between disturbance and moisture suggests a difference in the contrasts between disturbances in each moisture class (Fig. 5).
For all three models, there was no obvious deviation from normality and homoscedasticity based on visual inspection of the model residual plot.

Muskrat disturbance and wetland plant community composition
The PERMANOVA showed that community composition of wetland plants was significantly Note: Candidate models with different predictor combinations of disturbance (muskrat or reference) soil moisture class, and litter cover, plus their interactions, were compared using the Akaike information criterion corrected for small sample sizes (AIC c ). Coefficient values are shown for continuous variables; "+" indicates inclusion of categorical variables; and "NA" indicates omission. We considered values of DAIC c > 2.0 to indicate substantially better models, as per Burnham and Anderson (2004). The selected model is highlighted in bold. "Dist" indicates disturbance. Global model = Typha CoverD isturbance 9 Moisture. Fig. 3. Estimated marginal (EM) means of cattail live cover among quadrat moisture classes (Dry, Saturated, Flooded) by disturbance (Muskrat and Reference). Differences in cattail live cover between muskrat-disturbed and reference areas were not significant, but among moisture classes, quadrats with "Saturated" moisture had significantly higher live cover. Points indicate the EM mean, and error bars indicate the 95% confidence interval of the mean. Means sharing a letter are not significantly different under Tukeyadjusted EM means comparisons (P < 0.05).
We observed differences in distribution for nine plant taxa between muskrat and reference quadrats, as indicated from the KW test. Of these, six taxa had significantly greater abundance in muskrat-disturbed quadrats, whereas three taxa were significantly more abundant in reference transects (Fig. 7). Among the species that were more abundant around muskrat-disturbed areas, all are obligate wetland plants except for the facultative wetland species in the genus Galium. In reference transects, only Carex lacustris is an obligate wetland plant, whereas Phalaris arundinacea and Spirea alba are facultative wetland plants.
When quadrats over muskrat houses were removed, the PERMANOVA suggested that there were no significant differences between the wetland plant community assemblages (Pseudo-F = 1.3; R 2 = 0.02; P > 0.1). There were no significant differences in the dispersion between muskrat and reference groups as well (F = 0.32; P > 0.5). Consequently, there was no clear separation between the two clusters (N = 54; stress = 0.18; Fig. 6; Appendix S1: Fig. S3). Additionally, there were no significant differences in the distribution of plant taxa.

Muskrat house distribution
Muskrat house density at the Carpenters Branch site increased sixfold between 2001 and 2019, from an average of 8.33 houses (0.61 houses per hectare [ha]; n = 3 yr) prior to the installation of the water control weir, to an average of 50.8 houses (3.72 houses per ha; n = 5 yr) in the years afterward. The Knox test showed there was a significant spatiotemporal correlation (P < 0.0001) for muskrat houses in Carpenters Branch during the survey periods; the number of house pairs occurring within 0.0001 decimal degrees and 1 yr was significantly higher than expected (Appendix S1: Table S2). Based on the Monte Carlo permutations, the expected number of close pairs was 78, which was smaller than the observed 104 close pairs (Appendix S1: Fig. S4). Note: Candidate models with different predictor combinations of disturbance (D; muskrat or reference) soil moisture class (M), and litter cover (L), plus their interactions, were compared using the Akaike information criterion corrected for small sample sizes (AIC c ). Coefficient values are shown for continuous variables; "+" indicates inclusion of categorical variables; and "NA" indicates omission. We considered values of DAIC c > 2.0 to indicate substantially better models, as per Burnham and Anderson (2004). The selected model is highlighted in bold. Global model = Response~Disturbance 9 Moisture 9 Typha litter.

DISCUSSION
Our results suggest that muskrats act as ecosystem engineers in Typha-dominated wetlands of the Thousand Islands region and that they generate a disturbance regime with cascading effects on the biotic community, specifically on plants. Muskrat activities centered in these North American temperate wetlands create a more heterogeneous physical habitat and facilitate a more diverse and distinct local plant assemblage relative to areas with lower disturbance. We found lower (absolute) percent live cover of cattails in muskrat-disturbed areas; however, there were no statistical differences in overall Typha live cover in our surveyed locations. The wetland plant assemblage around muskrat houses had a greater biodiversity and a distinct community composition relative to the more homogeneous, adjacent wetland areas dominated by invasive Typha. Effects on biodiversity were influenced by abiotic factors (moisture) and biotic interactions (Typha litter), and these driver interactions were non-linear. Most of the impact of muskrat disturbance on the plant community centers on the house structure. Multiyear data on house distributions revealed that at high densities, muskrat houses are annually maintained or rebuilt in similar locations, which may amplify the intensity and scale of their impacts on the regional plant community and larger wetland ecosystem.

Muskrat disturbance increases wetland plant biodiversity
The higher levels of plant diversity and species richness observed around muskrat houses indicate that this local disturbance regime has a strong influence on the local wetland plant community. The distinct plant community composition around muskrat houses revealed in the multivariate analysis is consistent with the observed increased species richness and diversity relative to reference plots. We suggest that muskrat disturbance induces a type of patch dynamics, in which heterogeneity in microtopography and soil quality caused by house construction, in conjunction with herbivory on cattails, induced a competitive release for non-dominant plants in the wetland. Muskrat houses typically stand around 50 cm tall and create a distinct in situ elevation gradient within the surrounding habitat. The muskrat houses sampled in this study averaged approximately 61 cm tall, with most primarily constructed from cattails (in combination with some substrate and other plant materials). The difference in height between muskrat houses and the surrounding wetland creates a microtopography difference which may allow different species to establish and persist in the habitat. Moreover, Kangas and Hannan (1985) have shown that nondominant plant species were more abundant on v www.esajournals.org muskrat houses. Consequently, the observed community differences suggest that the area around active muskrat houses are distinct habitat patches in the wetlands.
The combination of muskrat structure building and foraging activities initiates a disturbance regime and diversity of habitat (i.e., removal of Typha spp., deeper water in rivulets, and drier habitats atop muskrat houses) that distinguishes the local environmental conditions from the surrounding wetland matrix (Fig. 8). As expected, the environmental factors of substrate moisture and cattail litter cover had a significant effect on species richness and diversity. However, these contrasts were more pronounced under muskrat disturbance and the effect was most apparent in the quadrats of intermediate moisture levels (i.e., the "Saturated" group; Fig. 4). Our observations suggest that the maintenance of rivulets and herbivory around muskrat houses create an active transitional area between the most intensely disturbed areas and the surrounding wetland matrix. These activities create a higher level of spatial-environmental variation in the wetlands which would likely increase both microtopographic complexity and overall surface area available for plant colonization, similar to other ecological analogs such as tussock mounds in the Great Lakes wetlands (Peach and Zedler 2006). While this pattern seems consistent with predictions under the intermediate disturbance hypothesis (Connell 1978), where some studies validate the idea that levels of intermediate disturbance support the greatest species diversity in ecosystems with moderate levels of productivity (Petraitis et al. 1989, Kondoh 2001, many counter examples also exist and there is not general agreement on the mechanisms that underlie these patterns (Fox 2013).
In addition, disturbances caused by muskrats through the removal of the dominant emergent plant, Typha spp., likely contributed to the higher species richness and diversity in these wetlands. Removing cattails through herbivory and constructing structures benefit other plant species, not only by reducing interspecific competition of a dominant invasive, but also by recycling some of the nutrients and energy through decomposition of feces and processed plant material (Kangas and Hannan 1985, Berg 1987, Wainscott et al. 1990. Other studies have shown that muskrat activities can lead to changes in abiotic conditions beneficial to plants, including increases in net nitrogen mineralization and nitrification rates Fig. 6. Non-parametric multidimensional scaling (NMDS) plot mapping three dimensional (k = 3) axes (black lines) of wetland plant species community by transects including all sampled quadrats (N = 54; stress = 0.17; A) and excluding the muskrat house center quadrats (N = 54; stress = 0.18; B). Note the 95% confidence ellipses drawn to highlight the different clusters for muskrat transects (red points) vs. reference transects (blue points). The clusters were significantly different when all sampled quadrats were included (left; P < 0.05). (Connors et al. 2000) as well as increases in temperature, pH, and dissolved oxygen (De Szalay and Cassidy 2001).
Although this study did not show significantly lower abundance of cattails in muskrat-disturbed areas, previous studies in the region demonstrate the capacity of muskrats to reduce Typha spp. biomass (Connors et al. 2000, Toner 2006. Reducing the abundance of the most dominant (invasive) species, specifically T. 9 glauca, would presumably allow other displaced native species to reestablish a foothold in the ecosystem. This idea is also supported within disturbance ecology where higher levels of disruption in the ecosystem handicap stronger competitors and give the advantage to better-adapted colonizers (Grime 1973, Kondoh 2001. While this study clearly shows the positive impact of muskrat activity on the wetland plant community, this effect is very likely dictated by greater seasonal hydrological patterns. This notion is supported by a study in the same region which showed that water levels were positively correlated with cattail stem density at half of the surveyed sites . It is important to note that muskrats depend heavily on seasonal hydrological patterns and are likely to occupy optimal habitats in the wetland that allow them to overwinter. In this context, our results help shed some light on the complex interactions between muskrats, their surrounding wetland community, and the local and Great Lakes regional hydrology (Fig. 8). v www.esajournals.org Supporting a distinct plant community with more wetland obligate species The plant community in muskrat-disturbed areas was significantly different compared to the non-disturbed areas, with greater presence and coverage of some obligate wetland plants compared to the rest of the wetland. Muskrat house construction uses locally available material, typically a combination of soil, submerged aquatic hydrophytes, and robust emergent hydrophytes, which creates some degree of heterogeneity within the Typha-dominated wetland (Alexander 1956). This structure-building behavior creates microtopographic differences in the wetland and allows facultative wetland plants to occur and was shown to be the key driver in generating the community differences found in muskratdisturbed sites.
Close to 70% of the taxa that were found in higher abundance around muskrat houses are classified as obligate wetland plants (Lichvar et al. 2016). Of the non-obligates, Galium spp. were typically found at the base of cattail stems or growing out from the house, while the wetland obligate, Bidens spp., was also largely found growing on cattail mats (Z. X. Kua, personal communication). Among the species identified to occur more in muskrat-disturbed transects, three taxa (i.e., Lysimachia spp., Scutellaria spp., and Cicuta bulbifera) were plants associated with sedge tussock mounds (Johnston and Zedler 2012). In this case, muskrat house structures may create similar physical structures which allow these tussock-associated plants to occur more frequently around muskrat houses. Moreover, the higher abundance of Lemna minor around muskrat-disturbed areas may be associated with the habitat type associated with muskrat activity. The lesser presence of C. lacustris around muskrat houses, alongside P. arundinacea, was also identified as wetland plants associated with tussock mounds (Johnston and Zedler 2012). However, inferences made based on these observed differences remain speculative and will require further studies and observations to help properly understand these phenomena.

Regional scale of muskrat disturbance
The results from the spatiotemporal analysis strongly suggest that muskrat house locations were not randomly distributed in space and time, but were more likely to occur in the same location from year-to-year. Although muskrat houses are typically transient in temperate zones, with spring floods flushing away the old houses, the results from this study suggest that under relatively stable hydrologic conditions, houses may be used for more than one year. Yearly winter surveys of muskrat house locations from Carpenters Branch suggested that houses, or at the very least house bases, were being reused interannually. Similarly, Alexander (1956) documented the reuse of house bases for several years in the marshes of central New York. Reusing whole or portions of houses whenever possible makes sense as muskrat density increases in a wetland Fig. 8. Conceptual diagram hypothesizing the relationships between wetland water level management, muskrat activity and wetland plant community dynamics. Solid arrows indicate direct relationships, and dashed arrows indicate indirect relationships. and as prime space and limiting resources decrease. Due to the geography and water management in Carpenters Branch, spring melts are mitigated and dampened locally, which extends the lifespan of muskrat structures and allows reuse of houses interannually. Indeed, the structural permanence of houses in certain marshes has required researchers to develop alternative methods for estimating muskrat populations than simple house counts (Engeman and Whisson 2003). Thus, verifying the active status of muskrat structures is important for ensuring accurate population estimates.
A comparison of the muskrat house counts across all surveyed sites in the TI archipelago over the three periods spanning 17 yr (2001-2005, 2012-2013, and 2016-2018) yielded the highest densities in sites with local water control structures (Kua 2019). The highest absolute number of houses recorded was from Cranberry Creek in 2001 when it was under water level management, while the highest muskrat house density recorded was in Carpenters Branch in 2017 at 6.88 houses per ha (Kua 2019). Though the numbers fluctuate in relation to annual water levels throughout the region, there is a clear difference in house densities between sites with water control weirs (3.39 house per ha) and unmanaged wetlands (0.53 house per ha; Kua 2019). These results agree with Toner et al. (2010) where muskrat house density was highest in Carpenters Branch and had a strong relationship with winter water levels.
These spatiotemporal patterns in conjunction with high house densities suggest that the cumulative impacts of local muskrat disturbances at a landscape scale can be substantial. While published investigations of large-scale ecological changes caused by muskrat activity are rare, studies by Alexander and Radway (1951) in Montezuma Wildlife Refuge, New York, and Danell (1978) in northern Sweden documented extensive changes in wetland vegetation associated with high muskrat density. Therefore, it is likely that our findings of strong local effects of muskrat disturbance on wetland plant diversity and community composition were amplified at a landscape level in the wetlands with high muskrat densities and long-term longevity. Many wetland systems with regional hydrologic alteration such as the flow-regulated SLR have depressed muskrat populations, and thus, observations of potential landscape-scale muskrat impacts are largely confined to wetlands with locally managed water levels and corresponding high muskrat densities (e.g., Carpenters Branch and Montezuma Wildlife Refuge).
Visual comparisons of house distributions across our surveyed sites suggested that in wetlands without local water control structures (i.e., subject to the SLR system-wide flow regulation), muskrats tend to occupy areas adjacent to main channels. In wetlands where passive weirs were constructed to raise water levels, they colonize areas further into the wetlands, but only in years in which the weirs were operational. This observed expansion into the water-controlled wetland could be due to the deeper average winter water level, which effectively expanded the available habitat for overwintering animals. The increased area of high muskrat activity would likely amplify the effects of herbivory and house building on the ecological community throughout the wetland.

CONCLUSIONS
The high intensity of human development has led to a dramatic decline of wetlands in recent decades, with an estimated reduction in total area of 54-57% globally (Davidson 2014, Hu et al. 2017. Combined with climate change and other stressors such as systematic water regulation and species invasions, the decline of natural wetlands highlights the need to understand potentially mitigating effects of ecosystem engineers, including muskrats.
This study showed that muskrat disturbance has a positive effect on enhancing species richness and diversity in wetland ecosystems and supports a distinct plant assemblage. These effects are mediated through non-linear biophysical feedbacks with hydrology and with herbivory. Species richness peaked at intermediate moisture levels in muskrat-disturbed areas and biodiversity gains were greatest for wetland obligate plants.
These results improve our understanding of the ecosystem-level effects that muskrats create as an ecosystem engineer in temperate North American wetlands. Because they are sensitive to changes in water level, muskrats hold great importance as an integrative indicator of wetland health and, as such, have been used to inform international water regulatory policy (Toner et al. 2010, International Joint Commission 2014. Results from the spatial-temporal case study in Carpenters Branch suggest that water level management in wetlands plays a role in amplifying the impacts of muskrat disturbance across a greater geographic scope within Typha-dominated wetlands. Restoring wetland hydrology by managing water levels to enhance muskrat activity may yield system-level benefits, especially in areas dominated by invasive plant.