Plant response to habitat amount and configuration in Swedish forests

There is an intense debate about whether habitat fragmentation has a negative or positive effect on biodiversity. We examined whether species richness and incidence of forest plants were negatively or positively associated with fragmented forest configuration. We also analysed whether the results support the fragmentation threshold hypothesis with fragmentation effects only in landscapes with small habitat amount.


| INTRODUC TI ON
While there is a consensus that habitat amount has a positive effect on species distribution and abundance, the role of habitat fragmentation has raised plenty of debate (Andrén & Andren, 1994;Fahrig, 2003Fahrig, , 2017Haddad et al., 2017;Hadley & Betts, 2016).
According to theory, fragmentation, partitioning of habitat into smaller units, may have a negative effect on populations through increased edge effect, small population sizes and Allee effects, decreased migration, changes in resource flow, and disconnection of positive biotic interactions and demographically interacting units (Haddad et al., 2015(Haddad et al., , 2017Ibáñez, Katz, Peltier, Wolf, & Connor Barrie, 2014). For instance, Haddad et al. (2015), in their synthesis of fragmentation experiments, state that fragmentation negatively affected population abundance, species richness and ecosystem functions in forests. Contrary to these ideas, Fahrig (2003Fahrig ( , 2017Fahrig ( , 2018 recently challenged the paradigm of the negative effect of fragmentation. She argues that most studies have shown that effects of fragmentation per se on species richness and occurrence and abundance of species are positive rather than negative (Fahrig, 2017). Fahrig (2018)  The revived debate on the importance of habitat fragmentation covers several important aspects concerning landscape effects both on species distributions and on community composition. In largescale systems covering entire landscapes, it is usually unrealistic to perform manipulative experiments disentangling the two main effects to estimate the interaction between area and fragmentation.
One possibility is to design an observational study containing replicate landscapes with the same amount of habitat but in different configurations (Fahrig, 2003). However, even in such comparisons, fragmentation effects may depend on habitat amount (Didham, Kapos, & Ewers, 2012). The fragmentation threshold hypothesis (Andrén & Andren, 1994;Villard & Metzger, 2014) describes an interaction between habitat amount and configuration where there are no effects of fragmentation at high levels of habitat amount.
Negative fragmentation effects start when habitat loss reaches a threshold, and there is also an extinction threshold at a very low level of habitat amount where persistence becomes impossible.
Although the literature on fragmentation is huge, empirical studies testing fragmentation thresholds are still rare and evidence limited (see Swift & Hannon, 2010 for review). For example, Rueda et al. (2013Rueda et al. ( , 2015 observed fragmentation effects in forest birds and plants that were negative for some species and positive for other species. Single-species analyses of the interaction between habitat amount and configuration have usually been performed on a small number of species (Villard & Metzger, 2014). Plants are an understudied taxonomical group in this respect: despite the ecological importance of plants as putative drivers for biodiversity in whole ecosystems through cascade effects (Schleuning et al., 2016), out of 118 studies and 232 single-species responses included in Fahrig's (2017) fragmentation review, only one study dealt with single-species responses of plants.
In this study, we examine the effects of forest amount and configuration on incidence of forest plant taxa in the province of Södermanland, Sweden. We define forest configuration as spatial patterning of forest, measured with different landscape pattern metrics (Hadley & Betts, 2016;Wang, Blanchet, & Koper, 2014).
We include the whole species pool of forest plants in the province of Södermanland in our analysis. Forest is suitable as a focal habitat in this study because our study landscapes include a broad range of forest covers, enabling an efficient test of the effect of forest configuration at different levels of forest amount. For taxa with a positive relationship to forest cover, we test for the effect of forest configuration and its interaction with forest amount. Finally, we inspect the form of the forest amount-configuration interaction in different taxa to explore how commonly the results follow the fragmentation threshold model in relation to other interaction forms.

| ME THODS
Our study area, the province of Södermanland, is located in the hemiboreal zone in south-eastern Sweden. The forested area of the province decreased several centuries ago due to fuel wood consumption and agricultural expansion. Forest loss was followed by reforestation, predominantly with production forests, starting in the 19th century and accelerating during the 20th century (Cousins, Auffret, Lindgren, & Trank, 2015;Peterken, 1996).
Sörmlands flora (Rydberg & Wanntorp, 2001), an atlas for the flora of the province of Södermanland, was used as the source for species distribution data. Distribution maps in the atlas cover the occurrence of vascular plant species in 2.5 km × 2.5 km quadrats in the province. The data were gathered by the Botanical Society of Stockholm (Botaniska Sällskapet i Stockholm) through inventories carried out by professional and amateur botanists in 1980-1999. Due to the design of the inventory, the sampling effort varied among the quadrats which increased the variability of the data but, in our opinion, has not caused any substantial bias in our analyses, clumpiness, critical thresholds, forest landscapes, fragmentation, habitat availability, habitat configuration, species incidence, species richness mainly because re-inventories were systematically carried out to minimize the number of false absentees (Rydberg & Wanntorp, 2001). Good road access could possibly affect the sampling effort (Rydberg & Wanntorp, 2001). This would favour quadrats with fragmented forests, since there were more roads in quadrats with fragmented forest configuration, especially when forest area was low (autocovariate regression, results not shown). However, our results showed a general tendency that such quadrats had a low species richness and, in many cases, low incidence of forest plant taxa (see Results).
The full data for Sörmlands flora, consisting of 1,687 plant taxon distributions, were downloaded from Artportalen, the Swedish Species Observation System (www.artpo rtalen.se). In line with Artportalen, the nomenclature follows the Swedish taxonomic database Dyntaxa (2015). Sörmlands flora primarily uses species as taxonomic units, but the data also include some taxa at within-species level when they have a distinctive distribution and ecology. Taxa that predominantly occur as cultivated were excluded from the analyses, similarly as aquatic plants. For forest cover data, we used the data- The total forest cover of Södermanland was 4,633 km 2 (55.2%). We included in the analysis quadrats having a land area of more than 50%. In quadrats with more than 50% aquatic habitats, terrestrial land mostly consists of small islands and shorelines habitats. We excluded this subset of quadrats (9% of all quadrats), because they have specific ecological conditions and a characteristic flora due to exposure to wind, water and ice erosion, water level fluctuations, and in seashores salinity (Jerling, 1999). To estimate the total land area in each quadrat, we used the Land Cover data of the Swedish version of the European CORINE Land Cover database (Naturvårdsverket, 2014), with 25 m × 25 m pixel size, and year 2000 as the reference year. In addition, we excluded the quadrats that were partially outside Södermanland, and four quadrats at small, isolated islands that did not have any neighbours, impairing the autocovariate regression analysis.
We measured forest area, as a proportion of total area, and forest configuration in each 2.5 km × 2.5 km quadrat. We use the term habitat configuration, that can be fragmented or aggregated, instead of the term habitat fragmentation when referring to our own data, because fragmentation is appropriately defined as a configuration change rather than a static landscape pattern (Hadley & Betts, 2016), and we did not include the historical development of the landscapes in our analysis. We carried out separate sets of models using two measures of habitat configuration: the clumpiness index (CLUMPY in Fragstats 4;McGarigal, Cushman, & Ene, 2012) and edge density (total length of forest edge in km per hectare). To minimize the confounding effects of correlation between forest area and configuration, the quadrats with extreme values of forest area were excluded from the analyses. In the models including clumpiness, the forest area of the quadrats ranged from 15% to 75% (Figure 1a), whereas the quadrats having 25% to 80% forest area were used in the models including edge density ( Figure 1b). As a result, 1,096 and 1,103 quadrats out of 1,375 quadrats that fulfilled the other criteria were included in the clumpiness and edge density models, respectively. Figure  nearest-neighbour distance or the largest patch index were not appropriate in our case because a large part of our study area has a high forest cover without a true patchy structure (see Figure 1a, rightmost forest maps). The same forest crossing the border and entering the same quadrat in several places is difficult to deal with when using patch-based indicators but does not confound the clumpiness index and edge density.
The statistical analyses were performed in R 3.6.0 (R Core Team, 2019). Single-species analyses were carried out with autocovariate logistic regression by generalized linear models (GLM) using binomial error distribution and logit link. Quasibinomial error distribution was used for models with overdispersed data. Species richness was analysed with ordinary multiple autocovariate regression. Multiple regression yields, according to Smith, Koper, Francis, and Fahrig (2009), unbiased estimates of the relative importance of habitat amount and configuration. One alternative could have been to use generalized additive models (GAM) enabling more fine-scale modelling of responses for different levels of predictors. In single-species analyses, GAM models can potentially show both extinction threshold and fragmentation threshold in the same model, which is not technically possible with two-way logistic regression. However, GAM models require a subjective visual interpretation of surface plots. We chose to use logistic regression, because it was straightforward to carry out for a large number of taxa, and the results were comparable and easy to interpret.
Tests with Moran's I showed that the GLM residuals were spatially autocorrelated for most of the taxa of the study (p < .05, moran.test function in the spdep package of R; Bivand & Piras, 2015). Similarly, the residuals of the multiple regression tests of species richness were spatially autocorrelated (p < .05). To account for spatial autocorrelation, we included an autocovariate, calculated by spdep, in GLM and multiple regression (Bardos, Guillera-Arroita, & Wintle, 2015;Dormann et al., 2007). Correlogram analysis of forest-related taxa with ncf package (Bjornstad, 2019) showed that, on average, significant autocorrelation vanishes after 15 km, so this distance was used as the neighbourhood size.
Forest area and the measures of configuration were standardized before the analyses and back-transformed to their original ranges in figures. The interaction between forest area and configuration was included in the model when it was significant. Complete separation in logistic regression that inflates the coefficient estimates was present in 79 and 74 of the clumpiness and edge density models, respectively. We could not fit these models despite trials with penalized maximum likelihood estimation (Heinze & Schemper, 2002). Complete separation occurred in taxa with low incidence, at highest 0.8% (11 or fewer quadrats occupied out of approximately 1,100 quadrats). When these taxa were omitted from further analyses, the clumpiness and edge density models both included 1,002 taxa.
Taxa with significant positive association with forest (p < .05; two-sided test) were regarded as forest taxa, in total 163 and 119 taxa (16.2% and 11.9% of all tested taxa) in the clumpiness and edge density models, respectively. Species richness was measured as the number of forest taxa in each quadrat using this criterion. This method to select forest taxa makes direct comparisons of effects sizes of forest area and configuration difficult since configuration effect can vary freely, while forest area effect is always positive and significant. An alternative option would have been to use independent data to find the species that have forest as their primary habitat, but we did not find such data for our study area with sufficiently high quality.
In the single-species analyses, we first inspected the group of forest taxa without significant forest area × configuration interaction. We were interested in the direction of the configuration effect and its strength in comparison to the forest area effect. We used the binomial test to analyse whether the configuration effect was significantly positive or negative for so many taxa that it is unlikely that it occurred by random. When using the critical value p = .05, the null hypothesis is rejected when significantly larger proportion of taxa than 5% × 1/2 = 2.5% has a significant positive or negative configuration effect. We also used the binomial test to analyse whether forest area usually had a stronger effect than forest configuration or vice versa.
When forest area × configuration interaction is significant, the direction and strength of configuration effect depends on forest area. We thus inspected the shape of interaction, that is how the slope of configuration changes between positive, zero and negative within the range of observed values of forest area. The model allows six possible interaction shapes. Denoting negative, zero and positive slopes with −, 0, and +, respectively (Appendix S1), the possible shapes, categorized according to how the slope of configuration changes from low to high forest area, are: −0, −0+, 0−, 0+, +0 and +0−. For example, +0− indicates that the slope of configuration is positive at low, zero at intermediate and negative at high habitat area. Of all taxa with significant positive relationship to forest area, the expected proportion of each shape by random is the critical value of significant interaction times the probability of each shape, that is 5% × 1/6, which we used as H 0 in the binomial tests. To be more precise, under random effect the shapes −0+ and +0− have a lower expected proportion than the other shapes, because random variation tends to push either the lower or the upper turning point outside the analysed range of habitat area. To minimize the risk for type II error, we carried out an additional binomial test with 5% × 1/9 as the expected proportion for these shapes, and to avoid type I error in testing the other shapes, we used a conservative value of 5% × 1/4 as the expected proportions. These additional tests yielded qualitatively similar results as the first set of tests, except in one case (see Results).
We used the Johnson-Neyman technique to examine the interaction shape by calculating the simple slopes of configuration along the range of the observed forest area (Preacher, Curran, & Bauer, 2006). We estimated the forest area range where the 95% confidence intervals of the simple slope of configuration included zero, implying that the effect of configuration was non-significant, and the regions of significance where the confidence interval of the slope was entirely negative or positive (Preacher, Rucker, & Hayes, 2007). Using the R package boot (Canty & Ripley, 2016), we calculated bootstrap estimates of the regions of significance for each taxon. The bootstrap sample size was 1,000. The confidence intervals were calculated by using the adjusted bootstrap percentile interval (type bca in the boot.ci function). We present results of individual taxa for descriptive purposes, but we are aware that our statistical approach is not designed to identify the statistical significance of responses to forest configuration in specific taxa, because we carried out many simultaneous tests without a correction for family-wise error rate.

| Single-species responses to clumpiness
In the single-species tests of responses to forest area and clumpiness, 163 taxa were positively associated with forests, that is the incidence of these taxa had a significant positive relationship with forest area in the 2.5 km × 2.5 km quadrats (p < .05, autocovariate logistic regression; Table 1; see Appendix S1 in supporting information for information on single species). Sixty-eight of the forest-associated taxa were also sensitive to clumpiness, that is there was a significant main effect of clumpiness on incidence, a significant forest area × clumpiness interaction, or both (Table 1).
There were 132 forest-associated taxa without significant forest area × clumpiness interaction, out of which 37 taxa had a significant clumpiness main effect (Table 1). The effect of clumpiness on incidence was significantly positive for 33 taxa (Table 1, Figure 2a). This number of taxa was higher than expected by random (33 positives out of 132 forest-associated taxa without significant area × clumpiness interaction, H 0 = 2.5% of taxa, p < .001, one-sided binomial test). These taxa thus had a higher incidence when forests were aggregated, indicating that they showed a negative response to fragmented forest configuration. An example of one such taxon, Maianthemum bifolium, is shown in Figure 3a. There were four forest-associated taxa with a negative effect of clumpiness (i.e., a positive response to fragmented configuration), which is not different from what is expected by random (4 negatives out of 132 forest-associated taxa without area × clumpiness interaction, H 0 = 2.5% of taxa, p = .420, one-sided binomial test).
Forest area had a greater impact on incidence than clumpiness. For the 33 taxa with significant positive association with both forest area and clumpiness, but without significant interaction, forest area had a larger standardized regression coefficient than clumpiness in 28 taxa and vice versa only in five taxa (significant difference from 1:1 ratio, p < .001, two-sided binomial test). Comparing average values of the standardized regression coefficients, the size of the clumpiness effect was 59% of the size of the forest area effect (Figure 2c).

A significant interaction between forest area and clumpiness
was observed in 31 of the 163 forest-associated taxa (19.0%), which is a higher proportion than expected by random (H 0 = 5% of taxa; p < .001, one-sided binomial test). The interactions had three different interaction shapes (Table 1). In the first type of interaction shape, including 19 taxa, the slope of clumpiness was significantly positive when forest area was small and non-significant when forest area was large (Figure 3b). In other words, the incidence increased with clumpiness, but only in landscapes with a small forest area. The proportion of this interaction type was higher than expected by random (observed proportion 11.7%, H 0 = 0.83%, p < .001, one-sided binomial test). The second type, including eight taxa, had a positive slope of clumpiness when forest area was small, non-significant slope in intermediate forest area and a negative slope when forest area was large (Figure 3c; observed proportion 3.7%, H 0 = 0.83%, p = .003, one-sided binomial test). In addition, the slope of clumpiness was non-significant in small forest area and negative in large forest area in four taxa (observed proportion 2.5%). This interaction type was also more common than expected by random when using H 0 = 0.83% (p = .048, one-sided binomial test) but not with the more conservative null hypothesis H 0 = 1.25% (see Methods; p = .148).

| Single-species responses to edge density
In the analyses of single-species responses to edge density, we found 119 taxa with significant association with forest area (Table 2). Of these 119 taxa, 50 taxa had a significant association with edge density, a significant interaction between forest area and edge density, or both (Table 2). Twenty-six taxa had only a significant main effect of edge density and no significant interaction (Table 2). Edge density had a TA B L E 1 Summary of relationship between forest area, clumpiness and incidence of taxa in 2.5 km × 2.5 km quadrats in Södermanland, Sweden Out of 163 taxa with significant positive association with forest area:  Twenty-four of the 119 forest-associated taxa (20.2%) had a significant area × edge density interaction effect (Table 2), which is a higher proportion than expected by random (H 0 = 5% of taxa, p < .001, one-sided binomial test). Of these 24 taxa, 10 taxa showed negative response at low levels of forest area, no response at intermediate levels and a positive response at large forest area (Table 2, Figure 3d). Nine taxa showed a negative response to edge density at low levels, but no response at high levels of forest area (Table 2, Figure 3e). Both these interaction shapes were more common than expected by chance (H 0 = 0.83% of taxa, p < .001, one-sided binomial test). There were three additional types of interaction shapes, but they were not more common than expected by random (Table 2, H 0 = 0.83%, p > .05, one-sided binomial test).

| D ISCUSS I ON
We found that fragmented forest configuration often had a negative effect on the incidence of the forest-associated plant taxa of Södermanland. In some taxa, negative response to fragmented configuration was independent of the amount of forest in the landscape.
Some other taxa showed a forest area threshold, where fragmented configuration only had a negative effect when forest area was small.
Such a forest area threshold was also present for species richness of forest-associated plants. It was sensitive to fragmented configuration in landscapes with small, but not large, forest area. There were also taxa with a positive response to fragmented configuration in landscapes with large forest area. When forest area was small, these species responded negatively to fragmented configuration.
Even when considering that some taxa showed both negative and positive responses to fragmentation, negative responses dominated, which contrasts with the findings of Fahrig's (2017) review that ecological responses to habitat fragmentation are more often positive than negative. Another important point that she makes is that even when there are significant responses to fragmentation, they are generally weak. Thus, the question is whether the negative responses to fragmented habitat configuration in our data are strong enough to have ecological importance. Out of all forest taxa, fragmented configuration had a direct negative effect, irrespective of the level of forest area, on 20% and 18% of taxa as measured by clumpiness and edge density, respectively. The proportion of taxa with a negative effect in low forest area and no effect in high forest area was 12% and 8% by using clumpiness and edge density, respectively. A negative effect in low forest area turning to positive effect in high forest area was present in 5% and 8% of taxa by using the two configuration measures.
Summing up, 34%-37% of forest taxa were negatively affected by fragmented configuration at least when forest area was low. It should be noted that direct comparisons of effect sizes of forest area and configuration are not meaningful because the criterion to select the forest taxa was that they had a significant positive association with forest area, whereas configuration effects could vary without restrictions. This applies to effect sizes in the species richness analyses as well as in the single-species analyses. One solution is to add a restriction based on configuration effects, that F I G U R E 3 Incidence and species richness of forest taxa in 2.5 km × 2.5 km quadrats in relation to forest area and configuration in Södermanland, Sweden. (a) Maianthemum bifolium, significant positive response for forest area and clumpiness (p < .05) and non-significant interaction. (b) Melampyrum sylvaticum, significant interaction between forest area and clumpiness. Simple regression shows a significant positive response to clumpiness at low values of forest area (black grid) and no response at high values of forest area (grey grid). (c) Response of Salix pentandra to clumpiness, positive response at low values and negative response at high values of forest area. (d) Response of Lysimachia europaea (Trientalis europaea) to edge density, negative response at low values and positive response at high values of forest area. (e) Lycopodium annotinum, negative response to edge density at low values and no response at high values of forest area is inspect only taxa having both a significant positive response to forest area and a significant negative response to fragmented configuration. In this group of taxa, forest area usually had a stronger effect than fragmented configuration, and the configuration effect was, on average, 59%-71% of the forest area effect, depending on the configuration measure. To conclude, forest area effects were stronger than configuration effects, but up to one-third of forest taxa were negatively affected by fragmented habitat configuration at least in low forest area, and when the effect sizes could be compared, the configuration effect was approximately two-thirds of the forest area effect. From a biodiversity management point of view, priority should be given for large forest area in our study system but avoiding fragmentation can be considered as a management goal, at least when it is not in conflict with management for large forest area.
The question is why our study deviates from findings of Fahrig's (2017) review. Response to fragmentation is inevitably specific to taxon and habitat. One possibility is that our study system, forest plants in hemiboreal forests, is an exceptional case, having negative responses to fragmented habitat configuration that are rare in other taxa and ecosystems. Undeniably, certain features of our taxa and ecological system differ from the ones in the studies included in Fahrig's (2017) review. However, it is difficult to find obvious reasons why these distinctive features would cause an opposite response to landscape configuration in our study compared to the majority of the reviewed studies. Fahrig's review contains 37 significant fragmentation effects on plants. However, only one of the studies included analyses at single-species level, as in our study, and it showed that fragmentation had a negative effect on the occurrence of the invasive bunchgrass Bothriochloa ischaemum (Alofs & Fowler, 2010). All the other plant studies had species richness as the response variable and indicated a positive response to fragmentation, which is also opposite to the findings of our study. In other taxa than plants, single-species studies commonly showed positive responses to fragmentation (Fahrig, 2017). De Camargo et al. (2018 also found that single-species responses of birds of southern Ontario, Canada, were more often positive than negative. It is not possible to assess from Fahrig's (2017) review whether studies carried out in the same biome (hemiboreal zone) and habitat type (forest) as our study showed positive or negative responses, since Fahrig (2017) only presents comparisons of biomes in broad categories of subtropical/tropical and temperate/boreal/tundra. In both categories, positive responses to fragmentation dominated.
Methodological factors, such as the measure of habitat configuration, may affect the outcome of fragmentation studies. The measures of habitat configuration that we used, clumpiness and edge density, yielded similar results for both incidence of single TA B L E 2 Summary of relationship between forest area, edge density and incidence of taxa in 2.5 km × 2.5 km quadrats in Södermanland, Sweden Out of 119 taxa with significant positive association with forest area: 2012) that we used does not have this property. We used edge density and the clumpiness index because they were not correlated with forest area and they worked more consistently than patch-based metrics in landscapes with high amount of forest. Wang et al. (2014) also showed that edge density and clumpiness were efficient metrics in differentiating between landscapes with different habitat configurations.
There are several studies that are not compatible with the criteria to be included in Fahrig's (2017) review that have reported negative fragmentation effects (Fletcher et al., 2018 For single species, between-patch heterogeneity may mean that although the average incidence in patches is low, there is a high probability that somewhere in the landscape there exist favourable conditions for the species. Fragmentation may also stabilize predator-prey interactions, if fragmented landscapes provide hiding places to prey populations (Fahrig et al., 2019). More studies are needed about mechanisms behind positive and negative fragmentation effects and how they combine to produce the total effect of fragmentation. The hypothesis that landscape-scale mechanisms more often produce positive fragmentation effects than patch-scale mechanisms also needs be tested, because there are positive fragmentation effects even at a patch scale. For instance, edge effects may be positive instead of negative. In our study system, Vinter, Dinnétz, Danzer, and Lehtilä (2016) observed that species richness was higher at forest edge than at forest core, both for all understorey plants and for forest specialist species.
Our approach was to statistically partition species responses to habitat amount, habitat configuration and their interaction. We found several types of interactions between habitat amount and configuration. The most common interactive response showed a fragmentation threshold (Andrén & Andren, 1994;Villard & Metzger, 2014), where fragmented forest configuration had a negative effect on the incidence when forest area was low but did not affect the incidence when forest area was high. Species richness also showed this type of response. In another common response type, the response to fragmented forest configuration changed from negative to non-significant and to positive with increasing forest area. In landscapes with a lot of forest, fragmented forest configuration entails several small gaps instead of a few larger areas of matrix (Figure 1). A possible reason for this result is that many forest understorey plants benefit from increased light availability in small gaps (Reader & Bricker, 1992;Svenning, 2000), whereas larger openings favour non-forest species (Gálhidy, Mihók, Hagyó, Rajkai, & Standovár, 2006). Plant dispersal may also be easier in a landscape with small gaps. These results support forest management models that recommend small harvest units instead of large contiguous clear-cuts in landscapes with high forest cover (Fahrig, 2017;Lindenmayer & Franklin, 2002). However, the effect of harvest unit size is species-dependent. A rain forest keystone species, the army ant Eciton burchelli, showed higher levels of local extinction when the forest had many small tree harvest sites compared to harvest of larger continuous areas (Boswell, Britton, & Franks, 1998). In that case, migratory behaviour and habitat specificity of army ants were important factors behind their sensitivity to fragmentation. Whitmore (1989) described that the difference in gap size can result in different species combinations due to the variation among the tree species establishment success in gaps of varying sizes. This idea is implemented in gap-based silvicultural systems that aim to enhance biodiversity by varying gap size (Kern et al., 2017;Millar, Stephenson, & Stephens, 2007;Mori, Spies, Sudmeier-Rieux, & Andrade, 2013). However, the effect of variable gap size will also be dependent on environmental conditions (Kern et al., 2017). This ongoing theoretical development of forest configuration effects will clearly have important implications for future sustainable forest management.
In summary, the incidence of many plant taxa in our study system was affected both by forest area and configuration. The effect of forest area was stronger than the effect of habitat configuration, but our results also show that habitat configuration is of importance explaining the patterns of species incidence. Our results showed that fragmented habitat configuration of landscapes had negative effects for approximately one-third of the forest species and for species richness in our study system. There is still relatively little empirical research about how single species are affected by the amount and configuration of habitat or by the combination of both factors (Villard & Metzger, 2014). We gathered empirical evidence for a large amount of species from readily available sources, a method also used, for example, by De Camargo et al. (2018). A large sample size enabled a powerful test of the effect of forest configuration at different levels of forest area. The asset of the method is that remote sensing data and species distribution maps for many taxa are available over large spatial scales. The existence of fragmentation thresholds also shows that, in our study system, habitat area is one of the mediating factors affecting the fragmentation effect.

ACK N OWLED G EM ENTS
The study was financed by The Foundation for Baltic and East European Studies. We are grateful to Sofia Dinnetz for her work with plant distribution maps.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data will be made publicly available before publication at Dryad Digital Repository.

B I OS K E TCH E S
Kari Lehtilä is interested in plant ecology and evolution in fragmented habitats. He also studies environmental and societal effects of management of biodiversity.
Tiina Vinter is a PhD student in plant ecology with interest in examining how spatial and temporal distribution of species is associated with landscape structure. She works with forest landscapes in local scales as well as in scales embracing several countries in northern Europe.
Patrik Dinnetz is interested in habitat characteristics and dynamics and its relationship with distribution of different taxa.
He also studies host, parasite and disease dynamics in periurban landscapes.
Author contributions: All authors conceived the ideas; K.L. analysed the data and wrote the first draft of the manuscript; all authors contributed to the writing of the final version.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.