Description of the study area and history

The study area is public land, located in northern Greece at the Greek–Bulgarian borders and it extends to an area of 6934.51 ha (Fig. S1 [suppl.]). This site is of great importance for biodiversity conservation as it belongs to Natura 2000 network (site code: GR1140002), according to the Habitats Directive 92/43/EEC. The vegetation of the area is dominated by the tree species: Pinus sylvestris L., Betula pendula Roth and oaks (Quercus frainetto Ten, Q. petraea (Matt.) Liebl.), which form pure and mixed stands. Other tree species present are: Picea abies (L.) H.Karst., Fagus sylvatica L. and Pinus nigra Aiton (the last one is found in small reforested areas). Pure or mixed stands of Scots pine and Silver birch cover about 80% of the area (approximate 5600 ha). Many rare and endangered species of fauna (listed in the Standard Data Form of the site) and flora (Eleftheriadou, 1990) of the region are directly dependent on the intra-forest environment of the region. The importance of these rare species has been recognized as they are included in lists II and / or IV of the Habitats Directive 92/43/EEC or in List I of Directive 79/409/EEC rare birds.

The area is also part of the European Green Belt (EGB). The EGB, running 12,500 km throughout Europe along the former Iron Curtain is the Europe’s backbone of an ecological network, and a living monument and global symbol for transboundary cooperation in nature conservation and sustainable development (Schindler et al., 2011). The heterogeneity and the intense dynamics of the studied area can allow a wide range of factors affecting forest dynamics (Teixeira et al., 2009).

The history of the site and the use of the land have historically played an important role in the vegetation formation of the area. Before 1922, there were some small settlements nearby, inhabited by nomads. As evidence, a traditional settlement is visible at aerial photographs taken in 1945. Extensive meadows and sparse remnants of forest stands existed at that time in the area (Zagas, 1990; Oikonomakis & Ganatsas, 2012). After 1946, the nomads abandoned the region and almost the entire Central Rhodope was designated as a Forbidden Zone. Thus, the existing scattered parental forest stands gradually colonized the whole area. After 1960, the forests were put under the management of Forest Service.

 

Methods applied to detect land cover changes, data processing and analysis

Remote Sensing (RS) accompanied by Geographic Information Systems (GIS) was used to investigate land cover changes, which is considered an inexpensive and practical solution and produces easily land cover maps and solves various environmental problems (Fassnacht et al., 2006). Given the above statement, the following geographical data were used to analyze the land cover changes in forest cover and the expansion of the different forest types in the area:

 • Aerial orthoimages of year 1945 and recent (2009) RGB aerial orthoimages of 0.5 m spatial resolution provided by the Greek cadastral service.

• Orthophotomaps of the year 1975 (forest vegetation maps) from the Hellenic Forestry Service. 

• Forest maps produced by the local Forest Service for the preparation of forest management plans with field measurements, which are updated along with the forest management plans every 10 years. The scale used was 1:20000. These analog maps were corrected using the high resolution Orthophotographs of the Greek cadastral service and a digital up-to-date forest cover map was produced. 

 • Digital Elevation Model (DEM) (labelled as ASTER 30, downloaded from https://search.earthdata.nasa.gov/search. in 30 m spatial resolution) used to create class maps of aspect, slope and elevation of the area.

 • Historical climate data (19 bioclimatic variables) were downloaded from the world climatic data website: http://www.worldclim.org. (accessed in 25-3-2020). These variables are used often in ecological niche modeling (Hijmans et al., 2005; O’Donnell & Ignizio, 2012; Fick & Hijmans, 2017).

• Auxiliary data such as: CORINE 2000 (http://geodata.gov.gr/dataset/corine-2000) Natura 2000 (https://www.eea.europa.eu/data-and-maps/data/natura-11/natura-2000-spatial-data/natura-2000-shapefile-1) outlines, Google maps (https://www.google.gr/maps). and ground observations of the area.

 

All raster and vector dataset projected into the Hellenic Geodetic Reference System (HGRS ’87).

Initially, landcover maps for each time period were created from the available orthoimages. Detailed photointerpretation (minimum mapping unit: 0.1 ha) was conducted to classify forest cover using basic principles (Lillesand et al., 2015). As a result, a three-class cover map was created, categorized as: Forest1 – sparse forests - canopy cover 10-40%, Forest2 – dense forests – canopy cover 40100%, Grasslands (non-forested areas)

Forest types were distinguished using the latest available forest maps of the Forest Service. The dataset was corrected by any convenient way including photointerpretation and in-situ observations as well as other auxiliary data, such as Land Use/Land Cover (LULC) maps of CORINE 2000 dataset and google maps. The final vegetation cover maps included the following categories: 1) forests dominated by Silver birch, 2) forests dominated by Scots pine. The more interesting categories with the greatest expansion were 3) mixed birch-pine forest dominated by Silver birch, 4) mixed pine-birch forest, 5) grasslands – barelands.

A spatial analysis was performed in order to estimate land cover changes and their forest type distribution. The spatial and quantitative distribution of birch and other species was investigated in the newly established forests in the later periods (after 1945). In the area of newly established forests, each factor’s (elevation, slope, aspect, bioclimatic variables, distance from parental clusters) spatial output was examined separately within the spatial database created with the appropriate GIS analysis, such as spatial overlay, reclassification, spatial query functions with either vector or raster layers (Oikonomakis & Ganatsas, 2012).

To facilitate statistical analysis, a random network of points was used with the limit being at least 50 meters apart. Following application of that limit, the number of random points was limited to about 2300 points, which were finally used. All the environmental and autecological factors were calculated for each point, and, thus, a random network of points (n=2267) with all the information (forest type, elevation, aspect, slope, bioclimatic variables, distance from parental clusters from the pre-existing forest in 1945) for statistical analysis was created. The Generalized Linear Model (GLM) (Guisan et al., 2002) was applied, in order to model the available data and to estimate the impact of each factor on the patterns of birch and pine forest changes. The data were separated in two time periods (1945-1975 and 1975-2009) and two GLMs were, thus, applied, to enable comparison of results between the first and the second period of forest expansion. The dependent variable is by nature categorical, the aspect variable was converted to categorical and the other variables were used as continuous.

Bioclimatic variables were tested for collinearity and correlations. Strong correlations among them can result in over-fitting of models in species distribution modelling (Pradhan, 2016). The variables were tested through multiple regression models and Variance Inflation Factors (VIF) were calculated. The cutoff value >10 was selected as in other studies (Shekede et al. 2018; Dong et al. 2020; Ncube et al. 2020). As a result, only two of them were selected: Mean Temperature of Driest Quarter (3-month interval) - BIO9; Precipitation of Driest Quarter (BIO17). These two variables delimit the temperature and drought tolerances of forests species.

Boxplots were used as a data exploration method for further insight on spatial distribution of forest species. The InterQuartile Range (IQR, Q1 – Q3), which repre-sents 50% of the distribution was used as a robust scale measure (Tarr et al., 2012) to better understand the dis-tribution of forest species in relation with each factor. All the variables were also converted to categorical ones, so as to cross-tabulate them with forest formation variables and to estimate differentiations from expected values, per-forming the chi-square χ² test.

The following categories were selected: 1) Depen-dent variable: 1: pure birch forest, 2: mixed birch-pine, 3: mixed pine-birch and 4: pure pine forests, 2) Aspect: Nor-th, NE, East, SE, South, SW, West, NW, 3) Slope: <10, 10-20, 20-30, 30-40, 40-50, >=60%, 4) Elevation: <800, 800-900, 900-1000, 1100-1200, 1200-1300, 1300-1400, 1400-1500, >=1500 m., 5) Distance from parent clusters: <100, 100-200, 200-300, 300-500, 500-800, >=800, 6) BIO9: equal intervals of 1ºC (temperature) were crea-ted and 7) BIO17: equal intervals of 3mm (rainfall) were created. SPSS statistical package was used for the statis-tical analysis.

ResultsTop

Tree colonization pattern in the non-forest areas

During the studied 64-year period, the forest expanded dramatically, occupying the existing non-forest areas. This resulted in almost all the area being covered by forests. Specifically, in 1945, the non-forest area was extended to 61.7% of the total studied area, and the forests covered only 38.3%. In 2009, after the gradual colonization of forest species, the area was found to be covered by forest by 97.8%, while only a very low percentage (2.23%) thereof remained as forest openings (Table 1). These very few areas that remained uncovered by forests are either remote areas (grasslands) or rocky areas which do not favor tree establishment. Photointerpretation of the available images showed that some areas have not been reached yet and remain as forest openings, as the direction of the expansion indicates. Rocky areas are also photo-interpreted in small areas. Forest species colonized almost all the existing non-forested areas (grasslands), which were found to cover 3,993.0 ha or 59.5% of the total area (Table 1). The area covered by forests and typical landcover examples of the three years 1945, 1975 and 2009 are displayed geographically and exemplified in Fig. 1.

The rate of forest expansion was lower during the first period (1945-1975), and significantly higher during the next period, 1975-2009. In 1975, the grasslands diminished to a percentage 37.13% of the total area. Specifically, the annual rate of this forest expansion was 55.06 ha/year in the first post-war period, and 68.86 ha/year during the second period (Table 1). An analysis of the expansion rates and behavior for each of the dominant forest species revealed that:

Birch increased (almost doubled) its occupied area during the studied period. In 1945, the species covered an area of 616.8 ha (mean altitude: 1124 m, mean BIO9: 14.04ºC, mean BIO17: 113.32 mm, main aspects: NorthNW-West slopes), and by gradual colonization of the open non-forest land, increased by another 536.4 ha (mean altitude: 1064 m, mean BIO9: 14.10ºC, mean BIO17: 114.44 mm, main aspects colonized: South – South-West slopes) reaching to a total area of 1153.2 ha pure birch forests in 2009 (Table 2, Fig. 1). Simultaneously, the mixed stands dominated by birch almost tripled, from 299.0 ha (mean altitude: 1141 m, mean BIO9: 13.97ºC, mean BIO17: 114.31 mm, main aspects: NW-West slopes) to 824.5 ha (mean altitude: 1175 m, mean BIO9: 13.92ºC, mean BIO17: 113.97 mm, main aspects: South-SW slopes) and eventually covered 1123.5 ha in the study area. A similar trend was also observed for Scots pine. P. sylvestris colonized another part of the bare land, resulting in the whole area being covered by forests. Similarly to birch, the mixed stands dominated by P. sylvestris increased considerably, from 406.6 ha (mean altitude: 1334 m, mean BIO9: 13.56ºC, mean BIO17: 115.73 mm, almost equal aspect slopes) to 1076.7 ha (mean altitude: 1211 m, mean BIO9: 13.79oC, mean BIO17: 115.34 mm, main aspects: SE-South-SW slopes).

These species trends and the combination thereof resulted in the formation of much more mixed stands than pure stands of both dominant tree species, leading to a coexistence of the two pioneer forest species B. pendula and P. sylvestris, in the greatest part of the afforested area.

As a consequence, mixed forests (especially the type: P. sylvestris– B. pendula) are the most abundant among the newly established forests, covering a large part of the newly forested area (1901.1 ha - 46.74%), while pureforests (especially B. pendula forests) can primarily be found at the old forest areas (1119.4 ha – 43.29% total pure forests). The distribution of forest cover in the time periods examined is shown geographically in Fig. 2 a) and Fig. 2 b). Fig. 2 c) also shows the percentages of the occupied area in three time periods. The percentage of the landcover is greater for P. sylvestris in new forests. Especially in the first period of expansion, pure and mixed P. sylvestris forests cover a large area (54.45% of the area forested in the period 1945-1975) which indicates that it has been spreading faster than B. pendula in pure or mixed formations. In the second period of expansion (1975-2009), the percentages for pure pine or birch forests are significantly lower than those for mixed forests.

 

Table 1. Land cover area in studied years, land cover changes and annual rates of change

* Land cover classes:

• Forest1 – sparse forests - canopy cover 10-40%

• Forest2 – dense forests – canopy cover 40-100%

• Grasslands (non-forested areas)

 

Figure 1. a) Land cover maps of the examined years and examples of landcover in each year. b) Examples of forest colonization in grasslands-barelands. Delineation in the examples is of the year 1975 and shows a typical mixed pine-birch forest installation in grassland

 

Table 2. Land cover area in studied years, land cover changes and annual rates of change

 

Factors affecting the expansion of forest species

The statistical analysis revealed significant differences in some of the environmental parameters between the areas colonized by birch and the areas colonized by Scots pine, probably due to differences in their niche requirements, or due to environmental barriers ( distance from forest edges).

 

Environmental drivers

Especially, during the first post-war period (1945-1975), birch was found to colonize lower altitudes (Q1:1159 - Q3:1282 m) compared to P. sylvestris, which was found to colonize higher altitudes (Q1:1262 - Q3:1416 m). In the second period of expansion, both species would spread to lower altitudes (Fig. 3). Also, birch was found to show preference to North-NW-West slopes, compared to P. sylvestris, which showed preference to south faced slopes. (Fig. 4). GLM models, which were implemented separately for the two time periods of forest expansion (1945-1975 and 1975-2009), showed that all environmental factors played an important role, except slope and BIO9 (Table 3), probably due to the fact that they overlapped with the other factors. According to Wald chi-square statistic, the greatest influence among the environmental factors comes from the elevation. Elevation inevitably correlates with bioclimatic variables because at higher altitudes the temperature is lower and rainfall is greater.

According to the distribution frequencies of the sample, in general (for the whole period 1945-2009) birch was found to colonize lower altitudes (Q1:1063 m, Q3:1253 m), compared to P. sylvestris, which was found to be established in higher altitudes (Q1:1219 m, Q3:1398 m) (Fig. S2 [suppl.]). Elevation and bioclimatic variables also diverged significantly from the expected values (Table 4) of possible species distribution (expected values), which entails that they are major factors with important influence to the final species geographical distribution.

The aspect in study area is distributed unevenly (Fig. S3 [suppl.]). The forested area in 1945 is dominated by North-NW-West aspect slopes. On the contrary, newly shaped forests are dominated by SW-South-SE aspect slopes. For this reason, new forests had to colonize southern slopes for the largest part. The geographical distribution analysis focused on the preference of species to certain aspect slopes, applying a chi-test to compare observed values with expected values.

A GIS analysis in combination with a statistical analy-sis (χ² test) showed that pure B. pendula forests showed preference to Northern aspects (mainly West-NW-Nor-th) (Fig. 4), while mixed forests with B. pendula as dominant species do not deviate from expected values in terms of their distribution (Fig. S4 [suppl.]). Conversely, P. sylvestris pure forests show preference to southern aspects (mainly SW-South), and mixed forests (P. sylvestris-B. pendula) preferred East, NE and SE aspects (Fig. S4 [suppl.]). In the second period of expansion, P. sylvestris forests still preferred North and South aspects, but the differences are not so obvious, as the χ² test values are lower (Table 4). Interestingly, the remained bare lands and grasslands differ in their environmental parameters, compared to the areas occupied by forests, being more abundant in S and SW aspects (Fig. S4 [suppl.]). This suggests that S and SW aspects were the least preferable (or less feasible) areas for colonization of tree species.

Generally, all new forest types were distributed in lower altitudes because it is easier for them to expand downhills. As a consequence, new forests colonized drier and hotter sites as shown in the climatic variable boxplots (Fig. 3). The direction followed by all forest types was from colder and rainier sites to hotter and drier sites. This is shown in Fig. 3 gradually from mature forests (<1945) to later-established forests (1945-1975) and finally to younger forests (1975-2009). The chi-test revealed higher values than expected for B. pendula in lower temperatures and lower rainfalls (BIO9: 13-14ºC, BIO17<113mm) in the first period of species colonization (1945-1975). It also revealed higher temperatures and lower rainfalls than expected (BIO9: 15-17ºC, BIO17<113mm) for the second period of birch colonization (1975-2009). On the contrary, the χ² revealed higher values than expected for P. sylvestris, in the coldest and wettest areas of its expansion, in both periods of its colonization (for period 1945-1975: BIO9<13ºC, BIO17>119mm and for period 1975-2009: BIO9<14ºC, BIO17>119mm).

Finally, slope inclination does not seem to play an important role in the manner of expansion of the examined tree species, as the two species and their mixed forest types demonstrate a similar behavior, especially in the first period of forest expansion. Generally, slope is the environmental factor with the lowest and least signif icant χ2 test values (Table 4). Generally, goodness of fit for GLMs shows that in the second period of expansion (1975-2009) the model fitted better, since the AIC value for the second model is lower (2539.9) than in the first period (2830.6).

 

Figure 2. Space-time composition of spread of the examined forest species – a) Silver birch, b) Scots pine and mixtures of them. c) Percentages of landcover of the main forest species of the area.

 

Figure 3. Boxplots showing the expansion of main forest types in relation to slope, elevation and bioclimatic factors over two periods of forest expansion.

Figure 4. Graphs showing the distribution of points of two main species pure forests in relation to the expected values of χ² test.

DiscussionTop

B. pendula is a pioneer and photophilous tree species (Rebele 1992; Suominen K. et. al 2003) which co-lonizes bare lands (Kinnaird, 1974). The same characteristics also appear at the species P. sylvestris (Gaudio et al., 2011; Durrant et al., 2016). B. pendula competes with P. sylvestris in space occupation in the existed non-forested areas, as revealed by the results. This competitive process resulted in the co-existence of the two pioneer forest species in the same area, and the formation of mixed forests. The high open space colonization capacity is attributed to the autecological behavior of the two species. Therefore, according to the results, there was a remarkable forest expansion in the study area, which was greater in the second period of the expansion studied (1975-2009).

The forested area in the period 1945-2009 was predominated by southern slopes (Fig. S3 [suppl.]). As a consequence, the secondary ecological succession occurred mainly in warmer and drier soils. In Northern Europe, B. pendula prefers similar sites to P. sylvestris i.e., dry soils with low solute concentration (Hynynen et al., 2009). Therefore, both species easily expanded in those areas, which adequately fulfilled their niche requirements. However, at late successional stages, the upper story controls competitive interactions between pine and birch, and birch decreases its presence and forms more compact stands, because pine is less sensitive to the competition of the upper story than birch (Paluch & Bartkowicz, 2004). This explains why older forests in the study area are occupied by pure forests to a greater extent. Scots pine in this study colonized bare areas faster than Silver birch. This can be explained by its niche requirements and its initially higher-altitude position, which enables it to spread fast downhills. Also, Scots pine often produces large quantities of seeds, which are small and lightweight and are easily dispersed by wind over relatively long distances. The reproductive cycle of the species starts early, from the 30th year in close stands, and much earlier, from the 7th to 10th year for individuals grown in free space in the studied area, which allows the species to quickly colonize non-forested areas (Zagas, 1990).

Combined geographical and statistical analysis showed that elevation is the most important environmental driver which triggers the differentiation of the two species in terms of colonization. Scots pine trees expanded to higher altitudes (Q1:1219 m, Q3:1398 m) than Silver birch (Q1:1063 m, Q3:1253 m). However, the mixtures of the two species make the establishment more complicated, especially due to the fact that both species have an obvious tendency to spread in lower altitudes, which are easier to colonize. GLM’s higher Wald chi-square values and χ² test showed that elevation is the most important environmental factor because it influences the distribution of the final species, as it differs from the expected distribution.

The GIS analysis showed that pure Silver birch forests showed a preference to N, NW and W aspects, while pure Scots pine forests preferred S and SW aspects (Fig. 4) following the distribution of barelands-grasslands. P. sylvestris probably meets the competence of B. pendula in northern slopes and has greater adaptability to southern slopes which are drier and hotter than northern slopes (Kinnaird, 1974; Kutiel & Lavee, 1999; Sternberg & Shoshany, 2001).

The χ² of bioclimatic variables revealed that Scots pine forests in colder and wetter locations are more abundant than expected. Silver birch did not show clear patterns of preferences, except for the first period of its expansion, in which it showed preference to colder areas. This suggests that it was easier for Scots pine to distribute in better climatic condition sites because of the advantage of the initial higher altitudes of their parental clusters. This contributed to a large concentration of pines in more climatically preferable sites. In contrast, this was probably not feasible for Silver birch, which did not have the opportunity to distribute in wetter and colder sites. Therefore, this indicates that Silver birch is probably located in more unfavorable sites compared to Scots pine.

Slope inclination proved to play a less important role in the expansion of the examined forest species. An older site condition study (Kinnaird, 1974) with findings about the colonization conditions of Silver birch has also showed that slope did not affect their seedlings and saplings growth. The study also suggested that there was a low preference of southern slopes, which is consistent with the present study findings. Finally, the study verifies that gaps in woodland and bare soils had the highest densities of seedlings that survived and developed better at these sites, which somehow explains the rapid expansion of birch in the present study area.

The geographical analysis of the spread of the two main species of the area showed that P. sylvestris expansion is faster than B. pendula’s expansion. The rate of expansion depends on the seed and pollen dispersal distances and the mortality of the seedlings. Kuparinen et al. (2010) suggests that the maturation age is 13 years for B. pendula and 20 years for P. sylvestris, and that birch performs better than pine in seed and pollen distance dispersal. If this is the normal condition of the two species expansion, birch should have established better in the study area. Nevertheless, Scots pine occupied more non-forested space and faster than Silver birch in this area. One possible explanation is that Scots pine individuals start their reproductive cycle very early in open space. Another possible explanation is that the mortality of birch seedlings may be higher than that of pine, due to lack of water, and drought. Besides, Scots pine is more resistant to drought and poor-nutrient soils (Gaudio et al., 2011). These characteristics give an advantage to P. sylvestris over B. pendula in hotter and drier conditions, and, thus, it can expand better during summer periods or prolonged summers. Another advantage of P. sylvestris in the area is that its parental clusters are located at a higher altitude, which allows easier spreading at downhills and faster and more distant colonization. Finally, it was proved that the expansion of Silver birch is more dependent on distance to parental clusters than the expansion of Scots pine.

The newly-formed forests differ from the pre-existing (in 1945) forests, which were mainly pure forests of the two main species. Species interspecific competition is probably the main reason why at a later successional stage the two species tend to form pure (non-mixed) forests. According to the bibliography, the two pioneer tree species tend to form pure stands as the time passes from their early life until they become mature stands. Valkonen & Ruuska (2003) showed that as the number of birches increases in mixed stands, there is a negative effect on the diameter growth and the maximum branch diameter in pine, while the height of P. sylvestris is not affected. The latter means that pines compete with birch for the upper story and suffer an impact in growth. Paluch & Bartkowicz (2004) make reference to the tendency of young pines to concentrate around the old ones; this was attributed to pine elimination by birches in gaps, which were more rapidly growing and strongly competing for water. This is the possible explanation for the extended areas of pure Silver birch forests in old forests of the area. It is also reported that there is a negative influence of young and old pines on birches, probably because of the worse light conditions, which explains the successive isolation of the two species. Other studies have concluded that B. pendula is able to modify its crown architecture and alter its strategy to compete with different neighbors (Lintunen & Kaitaniemi, 2010). B. pendula and P. sylvestris, as early-successional species, respond to shade with an increase in stem height and changes to biomass, and B. pendula fiercely competes with conifer species such as P. sylvestris, especially in fertile humus soils (Dehlin et al., 2004).

 

Table 3.Tests of Model Effects for variables in GLM models for main forest formations against environmental and autecological factors. Df=degree of freedom

Dependent Variable:

FOREST_CAT = Forest categories (1=B. pendula pure, 2=B. pendula mixed, 3= P. sylvestris mixed, 4=P. sylvestris pure)

PIN_DIST_45 = Distance from P. sylvestris forest existing in 1945,

BET_DIST_45 = Distance from B. pendula forest existing in 1945,

PIN_BET_DIST_45 = Distance from mixed P. sylvestris - B. pendula forest existing in 1945,

BET_PIN_DIST_45 = Distance from mixed B. pendula - P. sylvestris forest existing in 1945,

PIN_DIST_45 = Distance from P. sylvestris forest existing in 1975,

BET_DIST_45 = Distance from B. pendula forest existing in 1975,

PIN_BET_DIST_45 = Distance from mixed P. sylvestris - B. pendula forest existing in 1975,

BET_PIN_DIST_45 = Distance from mixed B. pendula - P. sylvestris forest existing in 1975

 

Table 4. χ² test values for established forests after 1945. Observed vs. Expected va-lues for pure and mixed B. pendula and P. sylvestris forests. Ν=degrees of freedom.

BIO9 and BIO17 are explained in methodology

 

Figure 5. Boxplots of main forest types, showing the distribution of main forest types in relation to distance from parental clusters for two periods of forest expansion.