Temporal changes in boreal vegetation under 70 years of conservation

Biodiversity conservation through protected areas (PAs) is often based on the idea that biodiversity is relatively static. This assumption is increasingly being challenged as species and communities shift their distributions in response to changing environmental conditions. Empirical evidence on the performance of PAs over decades is still sparse or lacking from several environments, although it is needed to understand species dynamics, support modelling of PA performance, assist PA management and ultimately, to achieve global biodiversity conservation goals. In 2021, we resurveyed vegetation of five boreal habitat types (heath forests, paludified forests, sun-exposed sites, mires and eulittoral sites) in Rokua National Park in Finland, where one of the conservation targets is to preserve the flora characteristic of the area. The study sites were originally surveyed in 1945-49, just before the National Park was established. Study sites have also remained free from the disturbances (forest fires and reindeer grazing) typical of boreal regions. We show that the compositional similarity of plant communities between habitat types has increased over time and is associated with the increase of forest species in several habitat types and the loss of many habitat-specific species. Drivers of change were most often linked to ongoing succession (understory closure) and changes in moisture conditions. Our results suggest that without natural disturbance or appropriate management efforts, the original conservation targets may be compromised over the decades. Our study demonstrates that resurvey of historical vegetation data can be effectively used to estimate long-term PA performance, helping to fill in missing temporal evidence.


Introduction
Conserving species, habitats, and landscapes through the establishment of protected areas (PAs) dates back over a century (Watson et al. 2014).The aim of PAs is to reduce or exclude direct human-induced disturbance and stress that may be harmful for conservation objectives.PAs continue to be a key conservation strategy, with new areas frequently being established to conserve nature and fight against biodiversity loss (Watson et al. 2014;Gray et al. 2016).In December 2022, it was agreed at the United Nations Biodiversity Conference (COP15) in Montreal that at least 30% of the Earth's land, sea, and freshwater ecosystems will be effectively conserved and managed by 2030 (Dinerstein et al. 2019).However, the assumption that biodiversity remains relatively static under conservation has been challenged (e.g., Pressey et al. 2007;Tingstad et al. 2020), because species and habitats change over time for many reasons.For example, the accelerating rate of current climate change has raised concerns whether PAs are still able to meet their goals as species and entire communities shift their distributions in response to changing environmental conditions (Thomas and Gillingham 2015;Hoffmann 2022).It is therefore of great importance to gain empirical evidence on how PAs manage to conserve biodiversity in the long term.
While the human-induced disturbances are often excluded, reduced, or otherwise controlled within PAs, they do not have control over all the drivers that shape species distributions and community composition over time.Because species respond differently to environmental drivers such as climate change (Parmesan 2006), or are protected during a successional stage that may change over time (Ranius et al. 2023), the establishment of PAs alone may not protect all species or habitat types effectively, equally, or as originally intended in the long-term.Implementing management measures (e.g., monitoring and restoration measures) is often critical to maintaining the desired state of conservation targets, but the lack of resources often prevents these efforts within PAs (Watson et al. 2014;Hoffmann 2022).Total exclusion of disturbances without active management may be problematic within PAs, as environmental changes may lead to faster or unexpected changes in species and communities.
Compared to our growing knowledge of global biodiversity loss, there is relatively little empirical evidence of temporal biodiversity changes in PAs over several decades (Thomas and Gillingham 2015).Empirical evidence would be highly useful in updating and adjusting policies, planning and management of PAs (Heller and Zavaleta 2009) and modelling of PA performance (Thomas and Gillingham 2015).Different frameworks and methods have been proposed to identify areas that will be important for species and biodiversity also in the future.They can be divided into ones that focus on protecting the "actors" or species by forecasting their potential distributions (e.g., ecological niche modelling), or the ones that protect the "stage", or the abiotic landscapes and features that are expected to offer a variety of refugia for different species to co-exist (e.g., geodiversity analysis) (Zhu et al. 2022).Empirical evidence on the long-term changes in species and biodiversity inside PAs offer valuable insights into the application of these methods.Long-term empirical evidence is also needed to improve understanding of the baseline conditions.This is essential to avoid the Shifting Baseline Syndrome, where insufficient knowledge of baseline conditions can lead to current conditions being regarded as "normal" (Soga and Gaston 2018).This could be particularly problematic in PAs, which are often considered as ecological reference areas (Josefsson et al. 2009).
Long-term empirical evidence of biodiversity changes is especially sparse or lacking from northern PAs (but see for example, Lehikoinen et al. 2019), and especially for their plant communities and habitat types.Northern PAs are experiencing particularly intense climate warming (Rantanen et al. 2022) but are also considered as future refuges for biodiversity following suitable climatic regimes (Berteaux et al. 2018).Boreal biome, the largest biome on Earth, is dominated to a large extent by coniferous forests, which are the backbone of the biome's well-being (Gauthier et al. 2015).However, the biome encompasses also other important habitat types, such as peatlands and freshwater habitats, which are critical for multiple ecosystem functions and services, yet are remarkably impacted by various global change effects (Beaulne et al. 2021;Heino et al. 2021).Boreal vegetation has been historically dependent on and is affected by several natural disturbances, such as wildfires, that have a vital role in the functioning of the boreal ecosystem, for example, by driving the understory composition in forests (Niklasson and Granström 2000;Nilsson and Wardle 2005;Gauthier et al. 2015).Grazing by large herbivores is another source of disturbance that strongly modifies ecosystem functions, community composition and diversity in the boreal biome (Väre et al. 1995;Happonen et al. 2021).However, little is known about the long-term trajectories of protected boreal habitat types in the absence of disturbances typical of the biome (fires, grazing) and concurrent climate warming.Succession can be expected to play a key role (Nilsson and Wardle 2005), but how this is reflected in biodiversity changes in the long-term is not fully understood.
Here, we provide empirical evidence of how vegetation has changed in five boreal habitat types (heath forests, paludified forests, sun-exposed habitats, mires, and eulittoral habitats) over 70 years of conservation.We use vegetation resurvey data from Rokua National Park in Finland, where one of the key conservation targets is to protect vegetation and habitat types on distinct geological formations.Our resurvey data is exceptional in including not only vascular plants but also non-vascular cryptogams (bryophytes and lichens).These are often underrepresented, even though they have important functions for several ecosystem processes (Turetsky et al. 2012) and typically constitute a prominent share of understory plant diversity in boreal habitats (Maliniemi et al. 2019;Maliniemi and Virtanen 2021).In addition to being free of land use disturbances throughout the study period, the study sites have remained unaffected by forest fires and reindeer grazing, which is restricted in the area.To investigate how well conservation has succeeded in maintaining the flora and habitat types of the area, and how the vegetation has responded to the exclusion of major disturbances over the past 70 years, we ask 1) how the composition and diversity of plant communities has changed (a) across the protected area and (b) within the habitat types that include both common (e.g.mires) and threatened (e.g.sun-exposed environments) boreal habitat types, and 2) what environmental drivers are linked to the vegetation changes and do they depend on habitat type?

Study area and the floristic surveys
Study sites are located in Rokua National Park in Finland (64.56 o N, 26.40 o E) in an area of approximately five km 2 (Fig. 1a).The national park was established in 1956 to conserve the distinct esker and dune formations of the area, which together with their vegetation create a distinct entity compared to the surrounding area.The national park is also part of Rokua UNESCO Global Geopark that was founded in 2010 owing to its geological history.Climatic conditions of the area are continental.Mean annual and summer (June, July, August) temperatures from the period of 1991-2021 are 2.5 o C and 14 o C, respectively, and annual and summer precipitation is 600 mm and 230 mm, respectively (Jokinen et al. 2021).Temporal climate trends from 1950 to 2020 suggest an increase of approximately 2 o C in annual mean temperatures and 0.7 o C in summer mean temperatures, whereas annual precipitation increased approximately 100 mm and summer precipitation 60 mm (Appendix S1).Moose (Alces alces) is the only large grazer in the study area, with a 10-year average of 3.4 individuals per 1000 hectares (Natural Resource Institute Finland, 2022).Reindeer (Rangifer tarandus) husbandry in the area ceased in 1800s and since then there have been no domestic reindeer (Hoppa et al. 2008).No major forest fires have occurred in the study area after 1860, when almost the entire area burned down (Hoppa et al. 2008).No signs of more recent, minor fires were found from any study site.Some fallen trees due to strong winds were observed in some of the study sites, but there were no signs of recent major wind throws in any site.Groundwater level fluctuations are typical to the area, but a general, slow decline has occurred during the recent decades (Karjalainen et al. 2013).
The study area belongs to the boreal vegetation zone.More precisely, it is located in the middle boreal vegetation zone where the southern boreal flora gradually changes into the northern boreal flora (Ahti et al. 1968).Species richness in the area is relatively low due to the dominance of sandy and nutrient-poor soils (Jalas 1953).There are less than 20 threatened plant and insect species in the area, many of which are found in old forests and sun-exposed habitats (Hoppa et al. 2008).Main habitat types in the study area include old pine-dominated mesic, sub-xeric, xeric, and barren heath forests, transitional and wooded mires and mire complexes, esker ponds and small humic ponds and lakes (Hoppa et al. 2008).Most vegetation types in the study area have well defined boundaries and are thus clearly distinguishable.The study area is one of the most important protected areas for barren heath forests in Finland (Hoppa et al. 2008), which in Rokua are known for their pristine and extensive lichen cover (Fig. 1b).Other specific habitat types to the area include sun-exposed esker and dune habitats (Fig. 1c).Both barren heath forests and sun-exposed habitats are classified as threatened habitat types according to the IUCN classification (endangered and vulnerable classes, respectively; Kontula and Raunio 2019).The reduction of large-scale disturbances, including natural forest fires and controlled burning, and the early stages of natural forest succession are the main reasons why these habitat types are threatened and are also considered future threats (Kontula and Raunio 2019).Some of the sun-exposed habitats in the area are managed by keeping them open.However, the sun-exposed sites of this study have remained purposefully unmanaged to follow their overgrowth and succession (Hoppa et al. 2008;pers. comm. Päivi Virnes).
The original vegetation surveys were conducted in 1945-49 by Jaakko Jalas (Jalas 1950(Jalas , 1953)).His purpose was to describe the flora and the vegetation types of the study area before the establishment of the national park (Jalas 1953).He inventoried a total of 60 sites including different types of forests, mires, eulittoral habitats and meadows.Of these sites, 43 sites remained within the boundaries of the national park.Due to the species poor vegetation in the area, relatively low number of study sites was considered to be sufficient to describe the flora of different habitat types.In each habitat, sites were established on a homogeneous stand of vegetation that was as comprehensive representation of the community as possible.The size of a site varied from 25 to 250 m 2 , with the most common size being 100 m 2 .The vegetation of a site was recorded generally from ten 1 × 1 m squares by identifying all vascular plant, bryophyte and lichen species and estimating their percentage covers.If the site was smaller than 100 m 2 , 4 to 7 squares were made.The original data (Jalas 1953) includes the frequencies and mean cover of species (calculated from vegetation squares) per site.Also, in his other work from 1950, Jalas mapped dune and esker vegetation from a large geographical area in Finland, including nine sun-exposed sites from Rokua.In this dataset, the size of a site was 100 m 2 , from which presences of all vascular plant, bryophyte and lichen species were identified.This data was included into our study to complement the study from 1953 that included limited information on sun-exposed habitats.Because no percentage cover information exists in this data from 1950, species frequencies are primarily used in the analyses of this study and complemented with cover data where possible.
In summer 2021, we relocated 42 sites consisting of a total of 293 vegetation squares within the national park.Of these, 11 were heath forest sites (94 squares of 1 × 1 m), nine paludified forest sites (84 squares), nine mire sites (59 squares), six eulittoral sites (49 squares) and seven sun-exposed sites (Fig. 1).Vegetation within each site were resampled using the same methods and the same number of vegetation squares per site that Jalas (1953) used over 70 years ago.The relocation of sites was based on detailed descriptions of their location.Thus, sites are considered as quasi-permanent (Kapfer et al. 2017).Most sites were easily relocated (but see the exception of a paludified forest subtype in Appendix S2) and a few sites with inexact or uncertain description were left out from this study.Inevitably, some relocation error may exist among the pairs of original and resurveyed sites.However, this error is minimized by the fact that all the habitat types are relatively species-poor, and the vegetation is generally homogenous within habitats.In statistical analyses, we also avoided pairwise comparisons of individual sites and concentrated on habitat-level comparisons.Because of the differences in the original sampling between the datasets from 1950 to 1953 and uneven number of plots and sites between habitat types, measures such as species richness are not directly comparable between habitats.However, we were interested in general trends in vegetation within the area and habitats, which can be examined with quasi-permanent data and applied to habitat conservation issues (Alfonsi et al. 2017).Species' taxonomy was harmonized between the datasets for temporal comparisons.A few closely related Sphagnum species were combined.In addition, the smallest liverworts (Hepaticae) and cupand horn-forming Cladonia species (i.e., other than the "reindeer-lichens") were treated at the generic level due to some changes in their taxonomy and because the same was done in the original survey due to difficulties in separating species (species list in Appendix S3).

Indicators of environmental change
To identify possible drivers of floristic changes, we used species-specific Ellenberg indicator values (EIVs) for vascular plants, which were adjusted to Northern European conditions by Tyler et al. (2021).These indicators were used to represent changes in environmental conditions, when no historical environmental data on, for example, soil pH or nutrients was available.Altogether nine relevant continuous indicators were chosen as predictor variables to explore drivers underlying changes in species occurrences (Table 1).Values of all indicators were found for all vascular plants in the data.We also considered modelling temporal bryophyte occurrences by their species-specific indicator values (Hill et al. 2007), but this was not meaningful due to large multicollinearity (> 0.8) between several indicators across all sites and within most habitat types.

Statistical methods
Statistical analyses were carried out using R version 4.1.0(R Core Team 2022).Temporal changes in community composition and diversity were analysed with R package vegan (Oksanen et al. 2022).To test if compositional shifts had taken place in relation to original community compositions across the study sites and within the habitats, we used permutational manova (PERMANOVA; Anderson 2001) with 999 permutations.Due to the repeated measures, permutations were not allowed within paired sites.We further tested whether community compositions had become increasingly similar across all sites and within habitats over time, i.e. if communities had undergone biotic homogenization.This was done based on the homogeneity of multivariate dispersions (PERMDISP; Anderson et al. 2006) and by calculating distances from each site to the spatial median of sites in a multivariate space during both surveys.Compositional shifts were illustrated with nonmetric multidimensional scaling (NMDS) ordination using species 0/1-data and Jaccard distance.To complement the interpretation of compositional changes and to see if changes in dominance structure (rather than in species identities) were driving compositional shifts, we made additional NMDS ordinations and corresponding analyses based on Bray-Curtis distance for those habitats with available cover data (heath forests, paludified forests, mires and eulittoral sites).Complete species data, i.e., vascular plants and non-vascular cryptogams, were used throughout.
We used richness-based species exchange ratio (SERr; Hillebrand et al. 2018) to explore how species identities had changed between the surveys.SERr applies presence-absence data and is calculated by dividing the sum of immigrations and extinctions by total number of species found across the surveys and is therefore a complement of Jaccard's similarity index.We calculated SERr at the study area level and habitat level and thus, avoided sitespecific comparisons.In addition, we calculated abundance-based species exchange ratio (SERa; Hillebrand et al. 2018) for those habitats with available cover data.SERa is closely related to Simpson's diversity index and is not that sensitive to species richness or rare species as SERr.SERa is calculated using species' proportional abundances and it characterizes changes in both identities and abundances of most common species.Both SERr and SERa get values between 0 and 1, indicating unchanged identity structure (SERr) or identity and dominance structure (SERa) at 0 and complete change at 1. Mean species richness was calculated for all species (total), vascular plants and non-vascular cryptogams (bryophytes and lichens) across all sites and within the habitats during both surveys.Because of the different number of sites and vegetation squares used per habitat, species richness is not directly comparable between the habitat types.However, data allow temporal comparisons in species richness, which were analyzed using paired t-tests.
To analyse if temporal changes in vascular plant frequencies were related to their associated EIVs (Table 1), we run binomial GLMMs using the package glmmTMB (Brooks et al. 2022).Models were made across all sites to explore general trends and run separately for each habitat type to explore whether temporal changes were related to the same or different indicators.Presences and absences of each species within each site were used as a response variable.Specifically, the response variable in each model consisted of the information on how many times each species was present and absent (success-failure) during each survey.Survey (original vs. resurvey), EIV and the interactions between survey and each EIV were used as fixed effects.Species identity was added as a random effect.Of the fixed effects, our focus was on the interactions, which revealed whether temporal changes in species frequencies were related to certain environmental indicators.Main effects of survey and each ecological indicator themselves were not interesting in our case.The main effect of survey only informed if temporal changes occurred in the presences and absences of all species and the main effect of each ecological indicator how its values are distributed across both surveys (see similar approach in Tyler et al. 2018).A model described above was first made for all the study sites and then separately for each habitat type.All indicators were included in the initial models, as there were not pairwise correlations larger than 0.7 (Appendix S4).Using a backward selection, non-significant interactions (p > 0.05), were removed.Main effects were also removed upon the removal of an interaction due to reasons described above.This procedure was continued until only significant interactions were left (removed variables in Appendix S5).Residual diagnostics were checked for each model.

Temporal changes in species composition and diversity
Species composition had changed across all sites and within each habitat type in 70 years (Fig. 2, PERMANOVA in Table 2).A general temporal trend, as described by the first axis of the NMDS ordination (Fig. 2), was that originally driest habitats (sun-exposed sites and heath forests) shifted towards moister types (paludified forests, mires, eulittoral sites), which in turn, shifted towards drier types.Thus, compositional extremes diminished over time.However, the increasing compositional similarity in terms of biotic homogenization was not confirmed across all sites (PERMDISP in Table 2).Of the different habitat types, biotic homogenization took place only within mires (Table 2).The results of the corresponding PERMANOVA and PERMDISP analysis made with cover data (Appendix S6) were in line with the results described above.The only exception when using cover data was that biotic homogenization had taken place within paludified forests.A closer look at species-specific changes revealed a general increase of late-successional vascular plants (e.g., ericaceous species) and bryophytes (e.g., Pleurozium schreberi, Dicranum polysetum) across the study sites (Appendix S3, S6).A total of 124 different species were recorded from the study sites in the original survey, while 96 species were recorded during the resurvey, a decrease of 23%.The temporal turnover, quantified as SERr, was 0.37 across the study sites, largely due to loss of species (Fig. 3a).Notably, none of the gains (Appendix S3) were entirely new species to the study area but were observed outside the study sites during the original sampling (Jalas 1953).Losses dominated the forested and sun-exposed sites while gains were more common in mires and eulittoral sites (Fig. 3a).Relating SERr to SERa indicated that in heath forests, rare species were replaced without pronounced changes in the dominance structure (i.e., no large shifts in the abundances of dominant species took place) (Appendix S7a).In eulittoral sites, simultaneuous shift in both species identities and dominance structure took place.In mires and paludified forests, moderate turnover took place in species identities and dominance structure.Mean total richness and mean cryptogam richness calculated across all sites increased slightly over time (Fig. 3b).No changes in mean species richness took place in heath forests and paludified forests while mean richness of all groups increased in mires and eulittoral sites.In sun-exposed sites mean cryptogam richness increased.In general, the habitat-specific gains were often species typical to the surrounding heath forests (Appendix S3).

Indicators of floristic changes
According to the final binomial GLMM across all sites, species indicating soil pH conditions (χ 2 = 25.44,p > 0.001) and phosphorus availability (χ 2 = 11.10,p > 0.001) had changed in frequency over time.Compared to the original survey, species preferring strongly acidic soils were found more often in the resurvey, while those favoring more neutral soils were found less often (Fig. 4a).Simultaneously, species preferring higher phosphorus availability became less common over time (Fig. 4b).

Discussion
We studied how boreal vegetation has developed after 70 years of protection and in the absence of major disturbances or direct management actions.We analyzed temporal changes in plant community composition and diversity across and within five habitat types (heath forests, paludified forests, sun-exposed habitats, mires, and eulittoral habitats).By analyzing differences in community composition between past and present states, and the drivers of these changes, we provided long-term empirical evidence of PA performance from a northern environment undergoing rapid climate warming.Even though the plant communities of the studied habitat types had preserved many characteristics of their past states, the compositional variation was clearly decreased among the present plant communities.Empirical evidence from our study supports the view that without appropriate biodiversity management, species composition and diversity may not correspond to the state that was originally conserved.This will become increasingly important under climate change (Ranius et al. 2023), and highlights the importance of local conservation efforts in achieving global conservation goals (Dinerstein et al. 2019).
Our results revealed clear shifts in community composition across the study sites and within habitat types.These changes were largely related to the increase of generalist forest species in several habitat types and loss of many habitat-specific species.However, evidence for across-habitat biotic homogenization was statistically weak (Table 2).Drivers of vegetation change, assessed as temporal changes in EIVs, were limited to a rather small set of indicators.They suggested that the observed changes in vegetation were related to ongoing succession (understory closure), which was expected given the general absence of disturbances and management actions.Changes in moisture conditions were also often associated with changes in vegetation.These results add to the existing evidence that compositional trajectories of different undisturbed boreal habitat types are not parallel in the long-term (Maliniemi et al. 2019;Kolari et al. 2021).We did not record any new species to the study area (e.g., no invasive species or range shifts of southern species), which is often one of the key goals in conservation.This is contrary to what could be expected under long-term climate warming yet may change in the future (Berteaux et al. 2018), highlighting the importance of more frequent monitoring in northern PAs.Instead, we recorded several local losses within habitat types.These losses were due to the local disappearance of species requiring open space (such as Antennaria dioica, Eleocharis sp.).These losses are clearly indicating undesired trends in local plant diversity and are concerning in a relatively species-poor environment.
Apart from mires, within habitat biotic homogenization had not taken place (Table 2).This is likely, at least partly, related to the appearance of species into habitat types from the surrounding forested sites, which contributes to maintain the compositional variation within plant communities.The appearance of these species had also increased species richness in non-forested habitat types (Table 3b, Appendix S3).Nevertheless, more frequent monitoring is needed to determine whether such trend is transitional (i.e., gains will increase towards dominance in the future that could possibly advance biotic homogenization) or more stable (i.e., original species will dominate, and gains will remain as subordinate species).The general characteristics of both lost (requirement of open space) and gained species (from forested sites) indicate ongoing succession.
A general temporal trend in EIVs of vascular plants across the sites was a decrease in species that demand high soil pH and P availability (Fig. 4).These changes are most likely linked to secondary succession that has been ongoing without any frequent (grazing) or intensive (forest fires) disturbances during the study period.In the absence of fires, soil pH and nutrients generally decrease over time (Hart and Chen 2006).The increase of poor-quality and more acidic litter due to increase of ericaceous species have likely led to more nutrient poor soil conditions and lower soil pH (Mallik 2003).A well-developed bryophyte mat absorbs phosphate effectively and thus reduces its availability for vascular plants (Chapin et al. 1987).Thus, the increase in species favoring soils with less P likely reflects the closure of understory vegetation.Interestingly, there were no distinct climate warming effects in terms of change in climatic EIV's (heat or cold requirement of species), despite the increase in temperatures and the location of the study area in the transition zone between southern and northern boreal floras.However, it is possible that climate warming increases productivity of the already dominant ericaceous species and thereby boosts overgrowth (as also discussed by Kolari et al. 2021 in relation to boreal mire vegetation).At the habitat-level, changes in EIVs related to succession were common, but changes in moisture conditions were also often observed.This can be related to earlier snowmelt or drier summers, yet these changes are likely linked to groundwater level fluctuations in the study area (Karjalainen et al. 2013).In any case, our results suggest that succession and changes in moisture conditions are key drivers of long-term vegetation change in the study area.
In heath forests, the most distinct change was the closure of understory vegetation.The grounds of heath forest sites were originally more open and characterized with sandy spots (Jalas 1953) that were absent during the resurvey.We observed an increased dominance of ericaceous dwarf shrubs (Appendix S6a), which according to Nilsson and Wardle (2005) occurs after long-term fire exclusion and may further reduce forest regeneration.Indeed, no closing of the canopy occurred during the study period (Appendix S8) nor changes in relation to light EIVs took place.Temporal constancy in EIVs of vascular plants (no indicators were left in the final model for heath forests) point towards relatively stable communities, however, EIV models were made for presence-absence data that do not account for cover changes of dominant ericaceous species that were already present in the original survey.
Understory community composition in paludified forests had shifted towards heath forest communities, and the shift was linked to more acidic soil pH and drier conditions (Fig. 5a, b).Notably, a subtype of paludified forests (seasonally hydrophilous forests), had completely disappeared (Appendix S2).Increased forest species abundance, especially ericaceous dwarf shrubs with acidic litter, likely lowered soil pH over time (Mallik 2003).These sites were originally comparably small in area, allowing the surrounding forest species to take over these habitats relatively efficiently.Even though we cannot address the exact reasons for drier conditions, they are likely linked to a general long-term decline in groundwater levels in the area (Karjalainen et al. 2013), rather than annual precipitation that showed long-term increase in the study area.From the conservation perspective, our results indicate that boreal habitat types that are small in area may be especially prone to disappearance as a result of long-term succession.
Our results add to the existing evidence showing how succession, as indicated by the decreased soil pH and nutrient availability, leads to closure of understory vegetation in sun-exposed habitats (Tukia et al. 2015).Species favoring openness and disturbances had greatly declined at the expense of increased forest species abundance.The decrease of species indicating high N availability in soil is likely related to the increase of N-fixing nonvascular cryptogams which form recalcitrant litter and further immobilize plant available N in soils (DeLuca, 2002;Nilsson and Wardle 2005).This effect is likely emphasized in this habitat type that was characterized by large areas of sandy and vegetation-free ground during the original survey (Jalas 1950).Losing plant species and habitat characteristics that are special to this habitat type will further negatively affect other aspects of biodiversity, such as several insects that may decline or disappear due to loss of their host plants or habitat preferences (Toivanen and Kotiaho 2007).
Long-term temporal evidence from mires and eulittoral habitats is very few to date.Of the studied habitat types, mires were the only type that had undergone biotic homogenization over time.Mires seem to be drying out to some degree, as found also by Kapfer et al. (2011) in southern Sweden and Kolari et al. (2021) in Eastern Finland, but some mechanisms also maintain locally wetter conditions.The increase in the abundance of species typical to drier wooded mires (Pinus sylvestris, Sphagnum fuscum) but also to wetter mire types (e.g., Rhynchospora alba, Sphagnum cuspidatum) indicate complex temporal dynamics and underline the importance of local moisture conditions to vegetation changes in mires.Most likely the groundwater level (water table depth) can be locally high, despite the general temporal decreasing trend in groundwater level in the study area (Karjalainen et al. 2013).In addition to mires, the local variation in groundwater levels had most likely also influenced the compositional trajectories of eulittoral communities.However, successional changes were also clear in this habitat type.The eulittoral zone was still characterized in places with tall sedges, yet open sand and related species (e.g., Eleocharis sp.) were no longer found.In many places the zone was paludified and vegetal invasion had taken place.More frequent monitoring in this habitat type is of great importance as this type is data deficient in national habitat type evaluation (Kontula and Raunio 2019).
Our study demonstrates that resurveying historical vegetation data, if available, can be effectively used to estimate long-term PA performance.Such efforts are extremely few to date but would have great potential to fill in the lack of temporal evidence on the perfor-mance of PAs, which is in demand globally and of increasing importance (Thomas and Gillingham 2015).Resurvey data can also continue to be used for biodiversity monitoring, as once old sites have been relocated, sites will be available for more frequent monitoring in the future.As our study sites are also located within a UNESCO Global Geopark area designated for geoconservation, data from this study can be further used to explore the relationships between abiotic characteristics and vegetation dynamics over time.This has rarely been done, but the different aspects of the environment (both actors and the stage; Zhu et al. 2022) should be considered, when assessing the resilience of PAs and their ability to conserve biodiversity under climate change (Tukiainen and Bailey 2022).
Our study provides insights into the long-term outcomes of excluding disturbances, such as wildfire, which are important for maintaining biodiversity in boreal environments (Thom and Seidl 2016).Our results suggest that effective management actions and further monitoring are highly important in PAs, especially when natural disturbances are completely excluded.In Rokua National Park, active management and restoration efforts to inhibit understory overgrowth are carried out outside the sites of this study (Hoppa et al. 2008).These include tree removal and controlled fires in heath forests, and disturbance of understory vegetation (such as moderate trampling and small-scale restoration fires) in sun-exposed sites where species require frequent disturbance to thrive (Tukia et al. 2015).In addition, smaller-scale drivers such as windthrows and insect outbreaks can be important local drivers that trigger re-vegetation (Angelstam and Kuuluvainen 2004).Long-term decline in groundwater level likely contributes to the successional development in paludified forests, mires and eulittoral zone.The management of groundwater dependent systems can be more complex (Eamus et al. 2006;Kløve et al. 2014), but even small-scale rehabilitation of lake shores, for example, may help to get certain species to come back (Geist and Hawkings, 2016).
Appropriate and effective management of PAs is essential to achieve the United Nations COP 15 target of protecting biodiversity (Dinerstein et al. 2019).In addition to climate change impacts, the long-term influence of local environmental drivers and disturbances (e.g., groundwater fluctuations and wildfires in our study area) should be extensively considered in the management of PAs, as they may be even more or equally important to climate change impacts, which are often prioritized.Local drivers and disturbances take part in maintaining long-term biodiversity and original conservation objectives, but are highly area-specific and dependent on, for example, habitat type and the underlying ecological processes.Thus, although many nature conservation targets, such as biodiversity protection and increasing the area of PAs, are global, the knowledge and actions that underpin them are ultimately local (Wyborn and Evans 2021).
Funding TM, HT and JA were funded by Academy of Finland project number 322652, TM and LM by Finnish Cultural Foundation and KH by The Ella and Georg Ehrnrooth Foundation and Jenny and Antti Wihuri Foundation.Open Access funding provided by University of Oulu (including Oulu University Hospital).

Fig. 1
Fig. 1 Map of the study area and focal habitat types.(a) Location of Rokua National Park and the study sites.Gray line indicates the border of the national park.(b) Barren heath forests (IUCN class endangered), characterized by extensive and undisturbed lichen cover.(c) Sun-exposed habitats (IUCN class vulnerable) harbor unique plant communities and are important for many specialized insects.(d) Lowsedge and flark fens (IUCN class least concern) form mire vegetation that stand out from the surrounding environment.(e) Eulittoral habitat dominated by tall sedges (IUCN class data deficient) forms a strip between the open water and the forest.IUCN-classes are the statuses in Finland (Kontula and Raunio 2019).Photographs: Tuija Maliniemi, Helena Tukiainen

Fig. 2
Fig. 2 Temporal changes in species composition.Non-metric multidimensional ordination (NMDS) based on species presences and Jaccard dissimilarity for each habitat type with 95% -confidence ellipses.Indented image is the same ordination with ellipses drawn across all sites during both surveys.Habitatspecific ordinations based on species' covers and Bray-Curtis dissimilarity are shown in Appendix S6

Fig. 4 Fig. 3
Fig. 4 Probability of occurrence of vascular plants related to environmental indicators across all study sites during the original survey (t1) and the resurvey (t2).Probabilities are model effects with 95% confidence intervals derived from the generalized linear mixed effects model and are shown for the indicators that were left in the final model.(a) Soil pH; from strongly acidic (1) to circumneutral (6) and (b) Phosphorus; from avoidance of soils with high P (1) to confinement of soils with high P (5).Detailed classification in Tyler et al. (2021)

Fig. 5
Fig. 5 Probability of occurrence of vascular plants related to environmental indicators in paludified forests (a-b), sun-exposed habitats (c-f), mires (g) and eulittoral habitats (h-i) during the original survey (t1) and the resurvey (t2).Probabilities are model effects with 95% confidence intervals derived from the generalized linear mixed effects models and are shown for the indicators that were left in the final models.Moisture: from mesic (4) to permanently under water (12); Soil pH: from strongly acidic (1) to circumneutral (6); Soil disturbance: from species colonizing and successfully competing in already established vegetation (1) to species reproducing only on naked or disturbed soils yet may persist for decades upon vegetation closure (6); Phosphorus: from avoidance of soils with high P (1) to confinement of soils with high P; Nitrogen: from very N poor (1) to N rich (6) soils.Detailed classification in Tyler et al. (2021).

Table 2
Summary of analyses on compositional (PERMANOVA) and beta diversity change (PERMDISP) across all sites and within habitats over time.Analyses are based on species presences and Jaccard dissimilarities.Delta (Δ) indicates the direction of beta diversity change, i.e., how the average distance of sites to the spatial median of species composition in a multivariate space has changed between the surveys, decrease indicating biotic homogenization and increase biotic heterogenization