*3.1. The Problem and the Setup*

We exemplify the potential contribution of focal-species modeling based on a recent decision in Estonia to include into national reserve network some forest types for which natural areas have been largely lost. The situation that the most productive forests are underrepresented in reserves and ecologically impoverished outside due to intensive use is common in developed countries [154–156]. Research has shown that some forest structures can spontaneously recover within a few decades (e.g., [157,158]), but demanding species re-colonize with a delay [159–161], and it is unclear whether these processes would benefit from active restoration [72,162,163].

The Estonian case followed from an analysis made in 2002 for the national forestry development plan, which identified forest protection gaps for old-growth biodiversity by site type [164]. A 2016 ministerial review set the remaining gaps as quantitative targets, prioritizing to set aside additional eutrophic (149 km2) and meso-eutrophic (147 km2) forests. It was clear that reasonably large patches of such areas only existed in impoverished states, but emerging research suggested that their protection might still pay off in the long run [165,166]. In 2017–2018, the Ministry of the Environment, the State Forest Management Centre, researchers, and environmental NGOs collectively identified a cost-effective selection of state lands that would cover most of the gap. A total of 286 km2 (1.2% of Estonian forest land) was set aside as a result, mostly by a single governmental decision in February 2019 (58 new strict reserves; 267 km2; including 25% meso-eutrophic, 36% eutrophic, and 16% eutrophic-paludified types). Here, we use spatial modeling to analyze the reversibility of the most vulnerable and degraded ecological conditions significant for biodiversity. The models will be used as a basis to assess restoration potential in these new reserves.

The reserves comprise 106 distinct patches all over the country, on average 2.7 km2 (range 0.03–28.5 km2) in size, and with a heavy management footprint. According to historical maps, the area has had >100 km of natural watercourses; ca. 30% remains in its natural streambed, but 56% has been straightened, and 14% is lost due to forestry drainage. Until the mid-20th century, 8% of the forest area was under some agricultural use (arable land; pastures; wooded grasslands). Most other forests have been converted by production forestry into mosaics of forest stands in various successional stages interspersed by networks of drainage ditches and forest roads. Stands >100 years old, which are the ecosystem targets of strict protection [164], cover only 7% of the area. The rest are recent clearcuts (11%), stands <20 years (28%), 21–60 years (28%), or 61–100 years old (27%). Artificial regeneration has been used on 45% of all forest land, mostly with Norway spruce (*Picea abies*), but the current share of the planted component varies widely among stands. Besides clearcutting, pre-commercial thinning (9% of forest land), and thinnings and sanitation cuttings (24%) have been used within the last 20 years. Of stands >20 years old, 47% are mixed, 30% conifer, and 23% deciduous forests. The main tree species are *P. abies* (42%), *Betula spp*. (25%), *Pinus sylvestris* (17%), *Populus tremula* (10%), *Alnus incana* (2%), and *A. glutinosa* (2%). Nemoral hardwoods, characteristic of such natural forests (*Quercus robur*, *Tilia cordata*, *Fraxinus excelsior*, *Acer platanoides*, *Ulmus spp.*), now only occur at small frequencies.

#### *3.2. The Modeling Approach and Inference*

We defined major threats as distinct empirically supported habitat dimensions, along which production forestry can reduce natural species pools of eutrophic and meso-eutrophic forests. The expert-based process (including two meetings) involved lead forest biodiversity experts in the country, with knowledge of multiple taxon groups. The forest types under question are well defined by topographic and soil conditions and have diverse species pools [167,168]. In a natural state, these forests have complex uneven-aged or all-aged structure created by gap-dynamics and, depending on moisture, rare stand replacements (mostly due to storm or pathogens) [164,169,170]. This structure is greatly simplified by clearcutting based forestry that uses 30–70-year rotations and a few selected tree-species, notably pioneer deciduous trees and planting of *Picea abies* over large areas [170–173]. However, rapid tree growth also accelerates structural recovery after abandoned management or long rotations based on natural regeneration [170,173,174]. Many well-dispersing old-forest species can colonize such forests at longer rotations [92,175], while others remain excluded due to the absence of old-forest structures [93] or (poor dispersers) lack of local refugia [176].

Considering these patterns, we defined seven main 'ecological profiles' of focal specialist species and their management-affected limiting factors in priority order (Table 2). The factors were then formalized as decision trees and parameterized based on requirements of the focal species and using practical habitat proxies (available in GIS). The output was designed as eleven threat-related habitat quality scores on the ordinal scale, which can be grouped qualitatively from non-habitat to quality habitat (Figure 2). The 'ecological profiles' were (Table 2): D1, a poorly dispersing perennial plant of gap-dynamic eutrophic forest, vulnerable to continuity disruption and dense shade in monoculture stands; D2, a poorly dispersing saproxylic species of natural *Picea* forests, which is vulnerable to the disrupted continuity of large downed trunks in moderate shade; D3, a rare species inhabiting senescent or dead *Populus tremula* in mid-succession, which is lost both in intensively managed forests and old-growth without *P. tremula* recruitment; D4, an epiphyte on old nemoral hardwood trees that are characteristic in natural forests (see above) but suppressed by production forestry; D5, an area-sensitive vertebrate in vertically well-structured stands, vulnerable to structural simplification and patch fragmentation; D6, a terrestrial invertebrate on stable moist ground that suffers from stand-continuity loss and unpalatable litter of forestry-favored conifer and *Betula* trees; D7, a (semi)-aquatic species of small forest streams, which is threatened by loss of microhabitats due to dredging of stream channel and upstream pollution from agriculture and drainage systems.

We illustrated the mapping approach by predicted changes from 10-years past (2009; based on real data) to current (2019) and 10-years future (2029; predicted by individual variables of the decision trees; Figure 2). We selected these relatively short symmetric time-frames here to assess both the delayed establishing of the reserves (compared to identifying their necessity) and planning and implementing a restoration program (compared to natural succession). Such short time frames also allow us to use a simplified approach to model uncertainty estimation. We ran the decision trees on openly available GIS sources and some critical elements digitalized for this project (notably forest continuity and stream channel changes from historical topographic maps). The basic spatial unit was the forest stand, as defined in the national forest registry. However, since other spatial data divided stands and subsequent forest surveys changed their borders, we performed areal calculations by rasterizing the maps (20 m grid). Certain subjective decisions were made (e.g., we applied the effect of thinning in a 20-year time frame), but the ordinal scale used appeared relatively robust to that. We used QGIS 3.10.2 [177] and R packages, dplyr [178] and sf [179], for the spatial analyses.

**Figure 2.** Decision trees for spatial modeling of the focal taxon groups D1–D3 and D5 in meso-eutrophic and eutrophic forests in Estonia (cf. Table 2): **(a)** poorly dispersing perennial plant of eutrophic gap-dynamic forests; (**b**) saprophytic fungus inhabiting continuous supply of large downed trunks of *Picea abies*; (**c**) rare fungus inhabiting senescent or dead *Populus tremula* in mid-succession; (**d**) area-sensitive vertebrate of vertically structured stands. The decision order is from left to right, and from top to bottom; the bottom row comprises habitat-quality scores (0...10; colors referring to broad classes). The parameters marked with asterisk (\*) were modeled as dynamic in the 10-year future scenario. If not specified, the tree variables refer to the 1st layer. The site type codes in (a): ND, *Aegopodium*; SL, *Hepatica*; kSJ, drained *Dryopteris*; JK, *Oxalis*; AN, *Filipendula*.




scores 0–2; High, scores 7–10. Changes in High (HighC) are percentage points relative to the 2019 level. References are for regional justification of the species and the model variables; the process also used unpublished expert knowledge.

We found that (i) although the areas had degraded age and tree-species structure (see above), the reserve selection had been generally successful. Net habitat loss in the last 10 years was only apparent for the perennial plant (D1); it will recover in a decade (Table 2). (ii) Regarding restoration potential, some popular stand-structural targets of active forest restoration (diversification of stand structure; dead wood creation) [162] were likely to be met at reasonable rates also by protection (Table 2; Figure 3b–d). (iii) In contrast, the pronounced lack of old nemoral hardwood trees for D4 will not be healed (Table 2). Since it cannot be rapidly addressed by restoration, too, another perspective is needed—perhaps protecting residual trees in the surrounding landscapes for the long term [200–202]. The same factor has degraded the habitat of litter-dwelling invertebrates (D6), but we expect their habitat quality to recover sooner along with undergrowth development. (iv) Another issue revealed was that even though 10% of the area was currently non-habitat for every focal species defined, we expect considerable passive recovery from that status (Figure 3e). Even fewer of such universally degraded stands appear concentrated enough to allow cost-efficient restoration. (v) The forest area containing quality habitats for lotic invertebrates is very small and, thus, potentially vulnerable to occasional disturbance. To sustain this part of biodiversity, we need a better basic understanding of its functioning in degraded forests and in relation to protection regimes.

Uncertainty of our models contains three major components. First, the priority order of the variables (sequence of decision nodes), which can be assessed by field-checking alternative decision trees. Second, parameter values at decision nodes to be analyzed for sensitivity. Third, parameter accuracy in the GIS sources that can be addressed by combining different sources. To exemplify, we report sensitivity of two models to high-priority nodes of tree age: 'oldest trees >60 yr' in D2 (Figure 2b) and four classes of 'mean tree age' in D5 (Figure 2d). Model projections for these threshold values changed by ±5 years did not yield abrupt changes in habitat quality distributions (quality habitat areas were not affected at all in D2). Trends were least sensitive: only one projection was affected by >1% percent point. Thus, a +5-yr threshold in D2 predicted a 10% decrease in non-habitat by 2029 instead of a 7% decrease at the original threshold and 6% for the 5-yr threshold.

**Figure 3.** Predictive habitat modeling of new reserves for eutrophic and meso-eutrophic forests in Estonia. (**a**) Locations of the reserves. (**b–d**) Predicted distribution of quality habitats (scores 7–10) in 2029 for three focal taxa, zoomed in for a selected reserve (**a:** red box). The colors refer to 2019 habitat quality (cf. Figure 2) and reveal: (**b**) moderately favorable, but only slowly improving situation for the perennial plant (D1); (**c**) poor, but rapidly improving, situation for the *Picea*-inhabiting old-forest fungus (D2); (**d**) favorable and further improving situation for the old-aspen inhabiting fungus (D3). (**e**) Shrinking of non-habitat (score 0–2) for any terrestrial focal species (D1–D6) by 2029 (black) from its current distribution (red).
