Effects of Long-Term Habitat Protection on Montane Small Mammals: Are Sorex araneus and S. minutus More Sensitive Than Previously Considered?

Protection of natural areas by restricting human activities aims to preserve plant and animal populations and whole communities, ensuring the conservation of biological diversity and enhancement of ecosystem services. Therefore, it is expected that the longer the protection, the stronger the desired effects. We evaluated the responses of small mammals at the population and community levels under protection in the southern Carpathian Mountains. We surveyed small mammals for five years in sites with long- and short-term protection and non-protected. Besides protection status, we included elevation, habitat heterogeneity, and the month of survey as predictors in our models. As response variables, we considered abundance, presence, species composition and species richness. Community abundance responded to all four predictors and species composition was influenced by protection status and month of study. The shrews Sorex araneus and S. minutus had positive responses to protection, both in terms of abundance and relative abundance (their ratio within the community). Our results suggest that overall, montane small mammal communities respond positively to long-term protection, especially S. araneus and S. minutus. These shrew species are considered habitat generalists, but they appear to be in fact sensitive to the habitat quality enhanced through protection.


Introduction
The primary causes of decline in global biodiversity are the destruction and degradation of natural ecosystems [1,2]. Forests have one of the longest histories of aggressive human intervention, especially in densely populated areas [3]. For a long time, forests were completely removed because of the need for timber [4] and for expanding pastures [5], agricultural fields and settlements [4], while later management practices included a variety of actions for enhancing forest regeneration [6]. Consequently, most contemporary forests of Europe are monoculture plantations of intensive commercial forestry [7]. Although various nature-friendly management methods have been developed and applied within organized forestry in the last 300 years [7], only a small proportion of the European woodland area remains spared of human impact, mostly in montane areas. Because of their inaccessibility, low human population densities and economic underdevelopment, the Romanian Southern Carpathian Mountains still host an important part of the European virgin forests [7], but these lack strict protection [8].
Protected areas are major evidence of humanity's efforts toward maintaining biodiversity [9]. Most of the time, a well-managed protected area reduces habitat loss and maintains Diversity 2022, 14, 38 2 of 15 population levels, even for threatened species [10]. For example, Gray et al. [11] found that species richness is 10.6% and abundance is 14.5% higher in protected areas than in non-protected ones, this being the result of differential land use between the two types of areas. Recently, protected areas have been created not only to conserve unique landscapes and ensure habitats for endangered species but also to contribute to sustaining human communities and, among other functions, to adapt to climate changes [10].
Studies focused on the response of small mammals to a particular human impact, such as the Peru LNG pipeline [29] or the clearings in managed Central European forests [12,[30][31][32][33], have shown that small mammals may also be useful indicators of environmental change, responding rapidly to changes in woodland management [29,30]. Therefore, it is expected that habitat protection would positively affect small mammal communities, with differential responses of various species, depending on their ecological requirements.
Numerous studies on small mammals have been conducted in various protected areas, such as wet meadows or riparian forests [34], wooded urban areas [35] and forests [36]. Some studies aimed to compare small mammal abundances and diversities within and outside protected areas [37,38]. However, to the best of our knowledge, no studies have tested the effect of long-term protection on small mammals in mountain landscapes.
We sought to assess the response of small mammals to the duration of habitat protection in mountain landscapes. We assessed small mammal responses to short-and long-term protection in montane habitats at community and population levels. Because natural habitats tend to be more complex than human-altered ones, the lack of protection may affect animal communities by diminishing habitat heterogeneity; therefore, we included this variable in our models. We also considered elevation and the month of study because we expected a significant elevational pattern and seasonality in the small mammal community parameters. Thus, we examined relationships between species abundance or presence, total abundance, species composition, and richness of small mammals as response variables and protection duration, habitat heterogeneity, elevation and month as explanatory variables. We predicted an overall positive effect of duration of forest protection on small mammal abundance and richness, independent of habitat heterogeneity.

Description of the Study Areas
We conducted our study in four areas of the Carpathian Mountains (Romania): Apuseni Mountains, Lotru Mountains, Râu S , es River Basin and Retezat Mountains ( Figure 1).
The Apuseni Mountains, the northernmost of the study areas, are characterized by low elevations (the highest peak is 1849 m a.s.l.). Covered in forests, with no subalpine shrubs or meadows and with scattered settlements, these mountains have a long history of logging, with many plantations and artificially thinned forests, clearings, hayfields and pastures. Although the study sites were part of the Apuseni Nature Park (established in 1990) at the time of survey, the recent and ongoing human disturbance (mainly forest exploitation) was very intense throughout the park, including in the near vicinity of some of our trapping sites.
The Lotru Mountains have been intensely logged during the last century. As a result, most forest plots are secondary forests, either planted or naturally regenerated. In some places, spruce forests were planted below their natural elevation. During the study period, no logging occurred in the trapping sites or nearby. peak is 2509 m a.s.l.), with extensive high ridges and peaks and include Retezat National Park, the oldest protected area in Romania, founded in 1935. In 1979, it was declared a Biosphere Reserve and included in the UNESCO World Heritage. After the change of the political regime in 1989, the number of protected areas in Romania increased greatly, and in the 1990s, the size of Retezat National Park was also increased by including areas adjacent to the limits of the early park. Figure 1. Location of the trapping sites. Red dots represent sites in non-protected areas; blue dots represent sites in newly protected areas; green dots represent sites within the limits of the early Retezat National Park. The map was made in QGIS [39] with a base map from Natural Earth [40].
In these areas, we surveyed forested habitats. At lower elevations, there are mainly beech (Fagus sylvatica) forests, usually with a dense tree canopy but poor herbaceous layer and understory. Mixed forests, composed of differing proportions of beech and Norway spruce (Picea abies) with scattered silver fir (Abies alba) and sycamore (Acer pseudoplatanus), are characteristic of mid-elevations. A Norway spruce forest belt reaches up to the timberline, which usually is present at elevations between 1600 and 1800 m, depending Figure 1. Location of the trapping sites. Red dots represent sites in non-protected areas; blue dots represent sites in newly protected areas; green dots represent sites within the limits of the early Retezat National Park. The map was made in QGIS [39] with a base map from Natural Earth [40].
The Râu S , es River Basin (including parts of T , arcu and Godeanu Mountains) comprises a mix of virgin, natural and planted forests, and subalpine shrubs and meadows. Old exploitations affected most of the spruce forests, which were clear-cut and then replanted, but forests at the timberline were usually kept intact. Beech and mixed forests were exploited by clear-cut or felling of mature trees. Tree trunks were removed while branches were usually left in place, representing most of the coarse woody debris in logged forest plots.
The Retezat Mountains are one of the tallest mountain ranges in Romania (the highest peak is 2509 m a.s.l.), with extensive high ridges and peaks and include Retezat National Park, the oldest protected area in Romania, founded in 1935. In 1979, it was declared a Biosphere Reserve and included in the UNESCO World Heritage. After the change of the political regime in 1989, the number of protected areas in Romania increased greatly, and in the 1990s, the size of Retezat National Park was also increased by including areas adjacent to the limits of the early park.
In these areas, we surveyed forested habitats. At lower elevations, there are mainly beech (Fagus sylvatica) forests, usually with a dense tree canopy but poor herbaceous layer and understory. Mixed forests, composed of differing proportions of beech and Norway spruce (Picea abies) with scattered silver fir (Abies alba) and sycamore (Acer pseudoplatanus), are characteristic of mid-elevations. A Norway spruce forest belt reaches up to the timberline, which usually is present at elevations between 1600 and 1800 m, depending on slope exposure and other geomorphlogical characteristics of the site. The cover of shrub and herbaceous layers varies greatly in spruce forests, the understory being composed mainly of spruce saplings, often with blueberry (Vaccinium myrtillus) bushes and a thick moss layer. The tree canopy cover of spruce forests decreases near the timberline, but the herbaceous layer is also reduced due to rocky outcrops and surfacing stones. At the timberline, mountain-ash (Sorbus aucuparia), stone pine (Pinus cembra) and juniper (Juniperus communis) shrubs are interspersed among dwarf spruce trees. This natural elevational succession of forest habitats is sometimes altered by temperature inversions or past logging and reforestation, which artificially lowered the lower limit of spruce forests.

Protection Status and Habitat Variables
We coded the protection status of the trapping areas as an ordinal variable, based on the duration of protection as: 2-sites benefitting from long-term protection (within the borders of the early Retezat National Park), 1-sites with short-term protection (in areas included in the park after 1989), and 0-sites in unprotected areas (here we also included the Apuseni Mountains) (Figure 1).
Because trapping was conducted in forested habitats, we expressed habitat heterogeneity based on five forest features evaluated for each transect using ordinal variables on a scale from 1 to 4, as follows: for tree canopy diversity: 1-only one species, 2-one dominant species (over 75% cover), 3-two dominant species (over 40% cover each) and 4-three or more codominant species. For cover of the shrub canopy: 1-no shrubs, 2-<20% shrub cover, 3-21-40% shrub cover and 4-41-60% shrub cover. For cover of the herb layer: 1-0-25% herb cover, 2-26-50% herb cover, 3-51-75% herb cover and 4-76-100% herb cover. For abundance of coarse woody debris: 1-absent, 2-scarce, 3-moderate and 4abundant and rocky outcrops. For surfacing stones: 1-absent, 2-scarce, 3-moderate and 4-abundant. We calculated habitat heterogeneity as the geometric mean of these five variables. We chose the geometric mean over the arithmetic mean to give more weight to more equal contributions of the five elements determining heterogeneity. We also considered the elevation (between 820 and 1910 m a.s.l.) of each surveyed habitat both as a quantitative variable (expressed in km for similar scaling with the other variables) and a factor with three levels (low-<1150 m, medium-1150-1450 m, high->1450 m) because in some response variables, we expected a unimodal response with a maximum at mid-elevation.
Small mammal trapping was conducted during the warm season (from June to September), and the month of survey was considered an ordinal variable: 1-June, 2-July, 3-August, 4-September.
To avoid multicollinearity, we tested the correlation between our four predictors prior to model construction, using the non-parametric (because of the ordinal variables) Spearman correlation. The only two variables that showed a significant correlation (p = 0.008) were protection status and elevation, but the correlation coefficient was low (rho = 0.3); therefore, we kept all the predictors in the analyses. In addition, we tested for multicollinearity of the multiple models using the variance inflation factor, but none was higher than 2, confirming no multicollinearity of predictors.

Small Mammal Trapping
Between 2003 and 2007, we conducted several trapping sessions using artisanal wooden and plastic box-traps (18 × 10 × 8 cm) set in transects. Transects included 30 to 40 live traps set 15 m apart, along the contour lines and parallel to the closest watercourse, forest edge, road or trail, within homogenous habitats. Trapping sites were selected along several valleys to cover the elevational gradient of forested habitats in the Retezat National Park and at similar elevations in the non-protected areas. In total, we surveyed 39 habitats, 25 in non-protected areas (1 broadleaved forest, 11 mixed, 13 spruce), 6 in newly protected areas (2 broadleaved, 2 mixed, 2 spruce) and 8 in the early park (1 broadleaved, 2 mixed, 5 spruce). Habitats were randomly surveyed between one and five times (once each year), resulting in 75 transects (visits), 34 in non-protected areas (1 broadleaved, 18 mixed, 15 spruce), 24 in newly protected areas (7 broadleaved, 7 mixed, 10 spruce) and 17 in the early park (3 broadleaved, 2 mixed, 12 spruce).
Traps were baited with sunflower seeds, apple slices and meat pasta. Hay was provided as bedding material. Traps were checked at dawn and dusk for two or three consecutive days, resulting in a trapping effort that varied between 60 and 120 trap nights (TN) per transect. Many traps were stolen, disturbed by animals (e.g., shepherd dogs, foxes, wild cats, squirrels) or closed by wind, heavy rain or invertebrates, mainly large coleopterans and slugs attracted by the bait. Therefore, to calculate the effective TN, we subtracted from the total number of TN those that were non-functional. The effective trapping effort varied between 20 (we eliminated the few transects where the effective effort was lower) and 120 TN per transect, with a total of 3657 TN. We identified captured animals to species based on morphological traits, and after marking, we released each animal at its trapping site. Because we were not interested in recaptures from one year to the next and only needed to identify recaptures during the same trapping session, we used temporary marking by fur clipping. We did not include recaptured individuals in any analysis.

Data Analysis
We used the following response variables: abundance, presence, species richness and species composition of small mammals. As a proxy for abundance, we used the number of captured individuals, but because the trapping effort strongly influences capture success, we also considered the number of effective TN. This was included in the models or used to calculate the capture index, as the number of captured individuals per 100 effective TN. As measures of species richness, we used the number of species captured per transect and the Chao estimator [41], calculated in the vegan package [42] in R version 3.6.1 [43]. To compare the overall diversity, we constructed the rarefaction curves for each protection status using the function iNEXT in the iNEXT package [44]. Because of the variation in sample size (different trapping efforts among transects), we used the individual-based approach for constructing rarefaction curves.
Some transects were surveyed repeatedly from year to year, but the strong temporal fluctuations in species abundance and spatial synchrony (at least at a local scale) [45] caused the effect of the site to be non-significant. Therefore, we accounted only for temporal dependency among transects, including year as a random factor when its effect was significant, using mixed-effects models in the lme4 package [46] in R. When the effect of year was not significant, we constructed (generalized) linear models ((G)LMs).
We modeled the total abundance and the abundance of the three dominant species of small mammals. We included the number of captured individuals as the response variable and the trapping effort as offset. Because overdispersion was significant, we used negative binomial generalized linear mixed models (GLMMs). For less abundant species (but only those captured in more than five transects), we constructed models of their presence using binomial GLMs (probability of presence did not vary significantly among years). Species richness (number of captured species and Chao estimator) followed the Gaussian distribution and did not show significant interannual variations; thus, we used LMs. To compare the models and find the most parsimonious one, we used the stepwise forward selection procedure (starting from the null model) with the likelihood-ratio test, which assesses the goodness of fit of two competing nested models based on the ratio of their likelihoods. In the GLMMs, the explained variation for the best model was expressed by the conditional and marginal pseudo-R 2 statistics based on Nakagawa et al. [47], computed in the MuMIn package [48]. The marginal R 2 represents the variance explained only by the fixed part of the model, while the conditional R 2 is interpreted as the variance explained by the entire model, including both fixed and random factors. For the GLMs, we calculated the explained deviance, and for the LMs, the adjusted coefficient of determination (R 2 ). The effects of interactions among predictors were also tested, but none were significant.
The effect of protection status and other variables at the community level was analyzed using Canoco 5.12 software [49]. We first ran a detrended correspondence analysis (DCA), and the length of the gradient was 3.2, suggesting the use of the unimodal method, canonical correspondence analysis (CCA). In the CCA, the response variables are the species' ratios within the community (the relative abundances of species), which we refer to as species composition hereafter. In the partial CCA, the effect of the considered predictors on species composition is evaluated after removing the effect of a covariate, in our case, the year. Rare species (captured in less than five transects) were not included in the analyses. Response variables were log-transformed by the expression y' = log(y + 1). The interactive forward selection was applied to choose the most parsimonious set of predictors for the partial CCA. Probabilities were adjusted (p adj ) to correct type-I error inflation caused by multiple testing, using the false discovery rate values [50]. The significance of ordination axes was tested by the Monte-Carlo permutation test with 999 permutations per test. Permutations were restricted to blocks defined by the year of survey.
The significance of the response of individual species to the selected predictors was illustrated through t-value biplots, which approximate the t-values of the regression coefficients of a multiple regression with the particular species as the response variable and all the habitat types as predictors, revealing statistically significant pairwise relationships between each response variable and each predictor [50].

Effect of Long-Term Protection on Montane Small Mammals
The overall capture index increased slightly with the length of protection (Figure 2a) due to the higher abundances of M. glareolus (Figure 2b) and S. araneus in the protected areas. The median capture index of S. araneus in the long-term protected areas was zero because in most of the transects it was not captured, but when present, abundances were higher than in other areas (Figure 2c). Apodemus flavicollis had similar values in non-protected and short-term protected areas and lower values in long-term protected areas (Figure 2d).
Muscardinus avellanarius and Sorex minutus were captured in 8 (10.6%, SE = 3.5%) and 5 (6.7%, SE = 2.8%) transects, respectively, while the other 8 species were found in less than 5 transects. Trapping success varied between 0 and 92.6 individuals/100 TN, with a mean of 11.32 individuals/100 TN (SE = 1.7 individuals/100 TN). Species richness varied between zero and five species per transect, while the estimated number of species (Chao estimator) varied between zero and six.

Effect of Long-Term Protection on Montane Small Mammals
The overall capture index increased slightly with the length of protection (Figure 2a) due to the higher abundances of M. glareolus (Figure 2b) and S. araneus in the protected areas. The median capture index of S. araneus in the long-term protected areas was zero because in most of the transects it was not captured, but when present, abundances were higher than in other areas (Figure 2c). Apodemus flavicollis had similar values in nonprotected and short-term protected areas and lower values in long-term protected areas (Figure 2d). When accounting for interannual variations and the effect of the other predictors (habitat heterogeneity, elevation and month of survey), in the mixed models with year as random factor, the protection status had a significant positive effect on the total abundance and of two dominant species, M. glareolus and S. araneus, but it did not influence A. flavicollis (Table 2).
Interannual variations in total abundance were very strong, with the pseudo-R 2 for fixed effect component 0.204 of the 0.664 for the whole model. Total abundance also proved to be significantly affected by the other three predictors. There was an increase with habitat heterogeneity, from June to September and a decrease with elevation ( Table 2). Myodes glareolus also showed significant interannual variation and a significant unimodal elevational pattern, with maximum abundances at medium elevations between 1150 and 1450 m a.s.l. and minimum abundances above 1450 m a.s.l ( Table 2). Sorex araneus showed a significant seasonal pattern, with increasing abundances from June to September but no significant interannual variation. Variation among years was the strongest in A. flavicollis, for which year explained more variation in abundance than elevation and habitat heterogeneity, as the pseudo-R 2 for the fixed effect component of the model was only 0.122 of the 0.772 for the whole model. Unlike M. glareolus, A. flavicollis had a linear, negative response to elevation and a positive response to habitat heterogeneity (Table 2). Table 2. Best univariate models for the parameters of the small mammal community, considering the protection status (Prot), month of survey (Month), habitat heterogeneity (Hab.heterog) and elevation (Elev) as predictors. In the model for M. glareolus, elevation was included as a factor with three levels (high, low, medium). S-captured number of species and Chao-the Chao estimator of species richness. In the mixed models (indicated by *), year of survey was included as a random factor. Probability distributions of the models' random components are given in parentheses. The presence of the less abundant species, M. avellanarius and S. minutus, did not show significant interannual variations and were significantly affected by protection status ( Table 2). The presence probability of S. minutus was positively affected by the duration of protection, and the explained deviance was 0.167. For M. avellanarius, the capture probability decreased with the duration of protection, but its effect was only marginally significant; therefore, the model fit was low ( Table 2).

Model (Predictors)
In contrast, species richness (observed and estimated) was not affected by the protection status but showed a significant seasonality, increasing from June to September. In addition, the captured number of species also showed significant interannual variation.
However, the variation in species richness explained by these models was low (Table 2), which indicates that species richness is more dependent on other variables that were not included in the analyses, being beyond the scope of this paper.
Rarefaction curves (Figure 3) indicate that recently protected areas had the highest overall species richness, while in the early park, the species richness was the lowest. However, the large confidence intervals, caused by the high number of singletons and doubletons, show that a much larger sample size is needed to compare the species richness of areas with different protection statuses. status ( Table 2). The presence probability of S. minutus was positively affected by the duration of protection, and the explained deviance was 0.167. For M. avellanarius, the capture probability decreased with the duration of protection, but its effect was only marginally significant; therefore, the model fit was low ( Table 2).
In contrast, species richness (observed and estimated) was not affected by the protection status but showed a significant seasonality, increasing from June to September. In addition, the captured number of species also showed significant interannual variation. However, the variation in species richness explained by these models was low (Table 2), which indicates that species richness is more dependent on other variables that were not included in the analyses, being beyond the scope of this paper.
Rarefaction curves (Figure 3) indicate that recently protected areas had the highest overall species richness, while in the early park, the species richness was the lowest. However, the large confidence intervals, caused by the high number of singletons and doubletons, show that a much larger sample size is needed to compare the species richness of areas with different protection statuses.

Responses of Small Mammal Species Composition
The variables that we considered were significant predictors of small mammal species composition (pseudo-F = 3.1, p = 0.002), explaining 19.5% (13.2% adjusted) of the partial variation in the community structure, i.e., after we removed the effect of year. All variables except habitat heterogeneity were significant predictors of community composition when considered separately, but their effects were partially overlapping. The best model, including protection status and month of study (pseudo-F = 4.7, p = 0.002), accounted for 15% (11.9% adjusted) of the residual variation in species composition after

Responses of Small Mammal Species Composition
The variables that we considered were significant predictors of small mammal species composition (pseudo-F = 3.1, p = 0.002), explaining 19.5% (13.2% adjusted) of the partial variation in the community structure, i.e., after we removed the effect of year. All variables except habitat heterogeneity were significant predictors of community composition when considered separately, but their effects were partially overlapping. The best model, including protection status and month of study (pseudo-F = 4.7, p = 0.002), accounted for 15% (11.9% adjusted) of the residual variation in species composition after removing the effect of interannual fluctuations. The two predictors were mostly independent in the constrained ordination space. Duration of protection was positively associated with the relative abundance of S. minutus and S. araneus and negatively associated with that of M. avellanarius and A. flavicollis (Figure 4). The responses of S. araneus and A. flavicollis were significant, as shown by the t-value biplots ( Figure A1a). The relative abundance of the two shrews increased from summer to autumn, while the relative abundance of the rodents decreased ( Figure A1b), but only the response of S. araneus was significant ( Figure A1b).
The interaction between protection status and month of survey had no significant effect on species composition (pseudo-F = 1.6, p = 0.136). associated with that of M. avellanarius and A. flavicollis (Figure 4). The responses of S. araneus and A. flavicollis were significant, as shown by the t-value biplots ( Figure A1a). The relative abundance of the two shrews increased from summer to autumn, while the relative abundance of the rodents decreased ( Figure A1b), but only the response of S. araneus was significant ( Figure A1b).
The interaction between protection status and month of survey had no significant effect on species composition (pseudo-F = 1.6, p = 0.136).

Discussion
We surveyed small mammal communities in protected and unprotected areas, in forest habitats differing in elevation and heterogeneity. Our study is the first to evaluate the response of small mammal communities to the duration of habitat protection in a montane landscape along an elevational gradient, at least in central Europe. During our survey, we identified 13 species. Myodes glareolus, Apodemus flavicollis and Sorex araneus were dominant. All four predictors that we tested (protection duration, habitat heterogeneity, elevation and month of survey) affected one or several of our response variables (species and total abundance, species presence, species richness and composition), but none had a significant effect on all of them.

Effects of Habitat Protection Status
The abundances of M. glareolus, S. araneus and S. minutus increased significantly with protection duration, as did the total abundance ( Table 2). The usual ratio of shrews within small mammal communities is between 8-15% [34]. This value is low compared with the overall relative abundance of S. araneus in our study (23.6%) and especially with the value for the early park (35.5%). These values are uncommonly high for S. araneus [51], especially considering that we used box-traps for sampling, which usually underestimate the abundance of shrews [52,53]. The increase of S. araneus population abundances due to

Discussion
We surveyed small mammal communities in protected and unprotected areas, in forest habitats differing in elevation and heterogeneity. Our study is the first to evaluate the response of small mammal communities to the duration of habitat protection in a montane landscape along an elevational gradient, at least in central Europe. During our survey, we identified 13 species. Myodes glareolus, Apodemus flavicollis and Sorex araneus were dominant.
All four predictors that we tested (protection duration, habitat heterogeneity, elevation and month of survey) affected one or several of our response variables (species and total abundance, species presence, species richness and composition), but none had a significant effect on all of them.

Effects of Habitat Protection Status
The abundances of M. glareolus, S. araneus and S. minutus increased significantly with protection duration, as did the total abundance ( Table 2). The usual ratio of shrews within small mammal communities is between 8-15% [34]. This value is low compared with the overall relative abundance of S. araneus in our study (23.6%) and especially with the value for the early park (35.5%). These values are uncommonly high for S. araneus [51], especially considering that we used box-traps for sampling, which usually underestimate the abundance of shrews [52,53]. The increase of S. araneus population abundances due to habitat protection [35] and the decrease due to forest management [54] were previously reported.
Sorex araneus and S. minutus are euryoecious species populating both lowlands and high areas [55]. The abundance of invertebrate communities is one of the most important microhabitat characteristics for shrews [23]. In harsh conditions, such as montane habitats, the spatial and trophic niches of S. araneus and S. minutus overlap [56], which is probably why the two species showed similar responses to protection.
The positive association between duration of habitat protection and the abundance of S. araneus and S. minutus in our study suggests that, despite their ability to inhabit a wide range of habitats, these species are sensitive to habitat quality, which may be enhanced by long-term habitat protection.
In contrast, Muscardinus avellanarius showed a marginally significant negative response to habitat protection status. M. avellanarius is a European protected species, listed under Annex IV of the European Commission Habitats Directive, and was the only species of Community interest captured during our study. Therefore, its negative response to habitat protection is, at least apparently, contrary to expectations. However, being a strictly arboreal rodent [57], the low capture rate with box-traps set on the ground does not reflect its real abundance in the protected areas because this survey technique, although it is the standard monitoring method in Romania, is inappropriate (and therefore not effective) for assessing M. avellanarius. Individuals captured in box-traps are those forced to forage or move on the ground due to reduced food resources and movement restrictions at the canopy level in suboptimal habitats.
Ancient woodland habitats are usually considered the primary habitat for M. avellanarius and other dormice [58]. However, regular forest thinning for regeneration and clear-felling seems beneficial for M. avellanarius, as long as these activities are well managed [59]. Goodwin et al. [59] highlight the importance for this species of traditional forest management and especially coppicing, a forestry practice that is seldom used in Romania.

Other Patterns
Heterogeneous habitats have a wide range of microhabitats, which provide optimal resources for several species, promoting local-scale species richness [60,61]. Therefore, we expected a positive effect of habitat heterogeneity on species richness, but we did not find it to be significant. Some authors (e.g., [62]) contend that the failure to confirm the positive effect of habitat heterogeneity on small mammal species richness is due to poor choices of structural habitat variables.
A wider range of microhabitats may benefit not only more species but also individuals, providing more resources, i.e., food and especially shelter, leading to increased population density, hence the overall positive effect of habitat heterogeneity on total abundance ( Table 2). Among the species, only A. flavicollis had increased abundances in more heterogeneous forests, possibly because of its mainly granivorous diet. In forests with a homogenous tree canopy, food is abundant during masting years but otherwise scarce, while diverse vegetation (especially woody) ensures more stable feeding conditions. The diet of A. flavicollis also includes animal food items (i.e., small invertebrates), which are more diverse [62] and abundant [63] in heterogeneous habitats, which is another potential mechanism behind the positive effect of habitat heterogeneity on A. flavicollis abundance. This effect should be even stronger in insectivorous species, but it was not confirmed by our results, probably because of a stronger response of these species to other habitat features enhanced by long-term protection of forests.
Although species richness did not respond to elevation, we did find a significant elevational pattern for total abundance and for the two dominant rodents (Table 2). Due to its more thermophilous character [64], A. flavicollis showed a strong negative response to elevation. Myodes glareolus showed a unimodal pattern of elevational distribution, with maximum abundances between 1150 and 1450 m a.s.l. At low elevations, M. glareolus may be outcompeted by the larger, more agile and aggressive A. flavicollis. In Britain, where the ranges of the two species overlap, M. glareolus is less abundant where A. flavicollis is also present [64]. Because of its strong positive response to tree cover [65], M. glareolus population densities decreased at high elevations, being negatively affected by forest rarefaction toward the timberline.
Total abundance and species richness increased significantly from June to September (Table 2). However, the higher observed diversity in autumn is probably an artifact, given that the number of captured species increases with the number of individuals captured. This is also supported by the fact that estimated species richness showed only a marginally significant seasonality, probably an indication that we undersampled rare species.
The total densities of small mammal communities are usually much higher during the autumn than during the breeding season [66] because of the recruitment of new individuals into the population. Furthermore, the breeding season becomes shorter with colder temperatures, corresponding to the decreasing length of the favorable season [67]. Therefore, the breeding season ends earlier at higher elevations, and we infer it was already over in most of our study sites in September, yielding maximum densities then.

Future Research Directions
Further studies aiming to evaluate the effect of long-term protection on non-volant small mammal communities should be conducted by using multiple survey methods. The use of nest boxes [68] and the inventory of nests along transects [69] are the most suitable methods for evaluating the abundance of dormice, while combined quadrate trapping with snap-and pitfall traps is generally considered the most suitable method for capturing shrews [34]. However, for ethical reasons, especially in protected areas, live trapping with sensitive box-trapping should replace snap-trapping.
A large-scale, long-term study in Great Britain revealed that although artificially established forests become gradually better over time for woodland species of small mammals, new plantations cannot act as a replacement for ancient woodlands, because the value of secondary forests for biodiversity does not match that of older, larger, undisturbed woodlands [12]. To test whether or to what degree the long-term protection of natural habitats can lead to small mammal communities similar to those found in virgin forests, comparative studies are needed.
Further research should also consider, besides species richness and abundance of small mammals, their functions in the ecosystem (e.g., interspecific relationships, seed and fungus dispersal, soil aeration) to evaluate whether long-term protection management could enhance the ability of small mammals to provide various ecosystem services.  Data Availability Statement: The dataset used for the presented analyses are available upon request to the corresponding author. Figure A1. The t−value biplots displaying the first two partial CCA axes with van Dobben circles drawn for: (a) protection status and (b) month of study. Red triangles represent the two predictors. Species indicated by arrows ending inside the pink circle show a significant positive response to that variable at p = 0.05, having increased relative abundance, while the species indicated by arrows ending inside the blue circle show a negative response, significant at p = 0.05, the response being stronger when arrows are shorter. For the other species, relative abundance is not significantly related to protection status or month of study after removing the effect of year.