*2.1. Study Area*

The Moscow Region is located in the central part of the East European (Russian) Plain: 35◦10 –40◦15 East, 54◦12 –56◦55 North, and covers an area of 4.69 million hectares (including Moscow: 0.26 million hectares) (Figure 1). The population is 20.4 million people (Moscow: 12.7 million people). The average population density is 2.29 people per square kilometer. The region has several important natural and phytogeographical boundaries. The significant border is the south edge of the natural range of the spruce forests in the broad-leaved spruce forests zone.

**Figure 1.** Study area and field sample plots. The locations of 1684 sample plots are shown on the figure along with the formation identity of each sample plot (color). A: Spruce, B: Spruce—aspen/birch, C: Pine—spruce, D: Pine, E: Oak—Spruce, F: Broadleaf—spruce, G: Linden, H: Birch, I: Aspen, J: Grey alder, K: Black alder.

Due to the high soil fertility, the forest cover of the region has experienced an extensive anthropogenic impact (felling, plowing) for several centuries. At the beginning of the 20th century, and especially after the Second World War, there was a significant change in the direction of impact: the use of forestry practices on the site of the former arable land. In 1947, all forests of the Moscow Region were recognized as a green zone, including a ban on industrial felling. However, up to the end of the 20th century, there was an increase in the pace of industrial development: the construction and operation of machine-building plants and related infrastructure, including enterprises of the energy complexes and oil refining complexes [35]. This has led to an increase in emissions of pollutants into the atmosphere and hydrosphere. At the beginning of the 21st century, the direction of influence changed to post-industrial. The growth of the population along with the construction of housing stock and transport infrastructure started in connection with the development of the financial sector of the economy [36,37]. In addition, deterioration in the volume of forestry activities aimed at maintaining the sustainability of forest plantations was recorded. The nature protection regime violations are as follows: unauthorized felling and household/industrial waste dumps. These violations result in large pest outbreaks, forest fires and degradation of species composition [17].

#### *2.2. Design of the Study*

1684 sample plots (including 1494 forest sample plots) were collected during 2010–2019 (Figure 1). The sample plots' locations selection is based on the principle of representativeness for the main

species prevailing in the forest stand, forest types, taking into account age groups, as well as the ecology of the habitat. The sampled vegetation communities are homogeneous in terms of the general floristic composition, the composition of the dominants of each layer, physiognomy (aspect, community structure) and habitat conditions. According to the methodology, the sample plots were limited to 20 × 20 m and located at a distance of at least 200 m from each other. When localizing the sample plots, the representation from the main types of communities was taken into account. Lack of road infrastructure was a serious limitation in the uniform laying of test plots. The following properties of vegetation communities are recorded at each sample plot:


The previous study of Chernenkova [38] has shown that the ecological-phytocenotic classification is best for analysis of communities in this region. The techniques allows for great differentiation of communities across the region [38]. The use of units of ecological-phytocenotic classification in the rank of formation and association groups is explained by a number of reasons: (1) good correspondence of typological and mapped units, (2) compliance with Russian units of forest typology and (3) consideration of rare types of forest areas, as well as secondary derivative communities, which is important from an environmental point of view.

Formations are identified according to the dominant forest species. Each formation is represented by communities with a different combination of common dominants in the lower stories. The classifications of syntaxons at the level of association groups is carried out according to the dominant ecological and morphological groups of plants of subordinate stories [38]. The association groups are connected with the features of the herb-dwarf shrub layer, trophic conditions and moisture conditions. Comparison of the formations' composition and association groups allows one to assess the direction of successional development, the degree of anthropogenic impact and the stability of forest communities.

For mapping the spatial structure, the maximum entropy is chosen [39]. The choice of method is based on the nature of the field data collection. There are two main ways to collect data on the distribution of natural objects [40], which require the use of various modeling algorithms: One is collection of only occurrence points. With this method of analysis (it is called "presence only" or "presence/background" presence analysis), one can use GPS data, materials of collections and herbaria, publications describing the places of species registration. Another is separate collection of occurrence and absence points. This is called a "presence absence" analysis. The possibility of collecting correct information about the absence of a natural phenomenon in any territory is debatable and may turn out to be false. The maximum entropy method is designed specifically for processing "presence/background" data [9].

SDMtoolbox software includes additional tools, namely calibration of multiple models, their testing and final selection of the best model using the quality indicators [41]. An important part of the tool is the spatial jack-knifing (or geographically structured k-fold cross-validation) [42]. This tool tests evaluation performance of spatially segregated, spatially independent localities. The tool also allows for testing different combinations of model feature class types (linear, quadratic, hinge, product and threshold) and regularization multipliers to optimize MaxEnt model performance.

Landsat 8 and Landsat 5 spectral reflectances and spectral indices as well as Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) and morphometric variables are used as the environmental variables [43,44]. Palsar-2 25 m resolution images are also included: ortho and slope-corrected backscattering coefficient (horizontal–horizontal (HH) and horizontal–vertical (HV) polarization) for 2019 [45]. To manage the autocorrelation, the variables that have 95% correlation

with other variables are removed, and 54 of 83 variables are left, the list of variables is provided in Appendix A Table A1. The Global Forest Watch dataset is utilized to prepare forest/non-forest masks as well as the loss year mask and water mask [46].

The modeling of forest cover is performed within two hierarchical levels: (1) forest formations and (2) association groups. Forest formations are aggregated syntaxons that are warranted by statistically sufficient and homogenous training samples. This makes forest formations more robust for spatial modeling. Association groups are more divisional, in that they have heterogenous training samples and thus they are more sensitive natural objects for modeling. The overall quality of modeling is evaluated by two confusion matrices.

On the first upper level, forest formations are modeled. Using the Global Forest Watch dataset, the territory is divided into two strata based on 30% forest cover threshold and forest loss/gain data: forest stratum and non-forest stratum. Raster masks are used to apply strata during mapping. Multiple approaches are tested for selection of bias layers [40–42]: no bias, and 10 and 25 km biases around occurrence points. The best results (highest area under the curve and 1-omission error) are obtained using forest layer as bias for forest formations and non-forest layer as bias for open lands and agriculture.

The systematic sampling approach that showed the best results along with bias layer, clustering, splitting and background restriction approaches is utilized [21]. To provide balance between sample sizes and sample equality, the sample size for each formation is systematically decimated to around 100 sample plots, because MaxEnt allows to include all non-linear feature class types only when the size is over 80 samples. To achieve 100 sample plots as well as to equalize sample sizes of formations and to reduce spatial autocorrelation, the following compromising approach is used. The spatial rarefication of occurrence data is applied.

Spatial jack-knifing is performed by three regions, five model feature class types are used (linear, quadratic, hinge, product and threshold) and three regularization multipliers (0.5, 1, 2). Final models of formations are integrated into one map by the method of highest position in ArcGIS. The confusion matrix is calculated between the final map and the initial sample plots.

On the second level, the modeling of association groups is performed. The number of sampling points varies significantly from 9 to 154: average 48.96, median 35, standard deviation 39.95. For this reason, no transformation and no spatial rarefication of the number of sample plots is applied. For every association group, the formation mask is applied. The confusion matrix is also calculated.

#### **3. Results**

### *3.1. Pre-Processing of Samples*

Table 1 shows the hierarchical structure of formations (in columns) and association groups (in rows) of the Moscow Region, in accordance with the previously published results of the ecological-phytocenotic classification [43]. The cells indicate the number of field sample plots for the association groups.

Formation A: The composition of spruce forests is complex (combinations of spruce with birch, aspen, pine and broad-leaved species) and it is similar to the composition of the zonal primary coniferous broad-leaved forests. The proportion of silviculture is high (mainly monodominant spruce forests) [44]. The species composition of the subordinate stories (grass-dwarf shrub and moss stories) is represented by the full spectrum of transitions from boreal to nemoral types. A relatively small number of community types is noted in the composition of boreal spruce forests (small herb and small herb-green moss groups). While, subnemoral (small herb-broad herb) and nemoral (broad herd) groups of spruce forests have a higher coenotic diversity, due to the higher participation of other tree species, and the variety of combinations of dominant land cover species.


*Forests* **2020**, *11*, 1088 1 DShG: dwarf shrubs-small herb-green moss, Sh: small herb, ShBh: small herb-broad herb, Bh: broad herb, MhBh: moist herb-broad herb, Gm: grass-marsh, H: herb,DHS:dwarfshrubs-herb-sphagnum.2 Mergedformations:oakandbroadleaf(EandF).3 Mergedformations:grayalderandblackalder(JandK).

Formation B: spruce–aspen/birch forests. This group of formations includes communities where spruce and birch are represented in equal proportions (with a small admixture of aspen). Spruce-small-leaved forests are interpreted by many researchers [47,48] as a short-term stage of spruce forests that form at the felling site as a result of both spontaneous succession and the development of spruce silviculture. The composition of the vegetation of the ground stories is close to that of spruce forests (Formation A).

Formations C and D: pine-spruce forests, and pine forests with spruce and birch, locally with linden, oak and hazel. Pine and pine-spruce forests on uplands are not completely indigenous communities and represent a successional stage in transition to mature forest communities. The absence of pine regeneration in automorphic habitats indicates a derivative origin of pine forests after fires and fellings, as well as in the composition of silvicultural forests. In one case, succession is accompanied by active recovery of spruce forests (Formation C). In the other case, in habitats with nutrient-rich soils, the broad-leaved succession is observed (Formation D). The broad-leaved species displace the pine and pine-spruce communities after a few decades. Since the ecological range of the pine is wide, the pine forests can be found on soils with different texture and moisture regime. Respectively, the pine forests' typological diversity is higher compared to spruce forests (Table 1).

Formation E: oak forests with linden, spruce and birch. This group characterizes broad-leaved forests with a predominance of oak in the first story of the stand, the participation of spruce in the first and second sub-stories and species of the nemoral group in the lower stories of communities.

Formation F: broad-leaved-spruce forests. Forests with a mixed composition of spruce and broad-leaved species (oak, linden, maple) and species of the nemoral group in the lower stories of communities. This is primary (indigenous) communities, which are usually replaced by linden or spruce forests during felling.

Formation G: linden with oak, locally with spruce and birch. Broad-leaved forests with a predominance of linden in the first or second sub-stories and nemoral species in the lower stories. Spruce is occasionally represented in the upper stories and in undergrowth. Oxalis and some other boreal species are involved in the grass story. These communities are primarily only on the slopes of ravines and river valleys and on the uplands in the central and northern parts of the region. Here, they are derivatives of coniferous-deciduous forests [49,50]. Indigenous linden nemoral grass forests are represented only in the southern part of the region.

Formation H: birch forests with spruce and aspen, locally with oak and linden, and Formation I: aspen forests with spruce and oak. The predominance of small-leaved birch and aspen forests in the region is associated with the formation of young forests in fellings. In recent decades, it is also associated with a regenerative succession on massively abandoned tillage. Most of the community types are secondary and can develop in a wide range of habitat conditions. The typological diversity of small-leaved forest communities of birch and aspen is associated with a wide ecological tolerance and the ability to grow on soils of different texture, moisture and nutrient richness. However this forest formation creates conditions for the restoration of conditionally primary communities.

Formations J and K: gray and black alder. These communities are more often referred to as primary forests, preferring either waterlogged or rich brook habitats. However, there is another point of view, according to which gray alder forests are considered derivative forests, developing on the site of broad-leaved spruce communities [51,52] (Table 1). Each forest formation is represented by one or multiple association groups:


Spatial rarefication for formations is performed in accordance with the methodology: A (10 km), B (10 km), C (1 km), D (5 km) and H (10 km). Two pairs of ecologically similar formations are merged: broadleaf formation (E and F) and alder formation (J and K). This made it possible to group the sample plots in accordance with the ecological similarity of syntaxa and achieve a relatively equal sample size: (82–112 sample plots in each sample).

The following non-forest habitats' field data were collected during field surveys: small leaf scrub (L), meadows (N), open marshy habitats (O) and agricultural fields (P). The following non-forest habitats are taken from the Global Forest Change dataset: cuts (M) and water objects (Q) [46]. Settlements (R) are taken from Openstreetmap (OSM) layers [53]. These habitats are mapped and merged with the final map through non-forest mask (Table 2).

**Table 2.** Number of non-forest field data points and data sources.

