Impact of Invasive Tree Species on Natural Regeneration Species Composition, Diversity, and Density

Invasive tree species decrease ecosystem resilience with negative impacts on natural regeneration. The influence of alien tree species on ecosystems is unevenly recognized and does not always account for different habitat specificity. We assessed the impacts of the three most frequent invasive tree species in European forests: Prunus serotina Ehrh., Quercus rubra L., and Robinia pseudoacacia L. on natural regeneration diversity, species composition, and density. We hypothesized that invaded forest types, in comparison with non-invaded, will differ in terms of species composition, will have lower taxonomic, functional, and phylogenetic diversity of natural regeneration, and will have lower densities of native tree species. We used a set of 189 study plots (200 m2) in a systematic design, established in various forest types in Wielkopolski National Park (West Poland). We analyzed impacts of forest type, accounting for soil C:N ratio, soil pH, and light availability on natural regeneration (woody species up to 0.5 m height) species composition, diversity, and density. We found an overlap of species composition among invaded and non-invaded forests and low impacts of invasive species on taxonomic diversity and functional richness. We found no impacts on phylogenetic diversity and other functional diversity components. In contrast, we found that the natural regeneration of forest-forming tree species reached lower densities in invaded than non-invaded forest types. However, sub-canopy and shrub species reached higher densities in invaded than non-invaded forest types. We confirmed that invasive tree species affect natural regeneration by decreasing the regeneration density of native tree species (in eight of nine tree species studied), species composition homogenization, and supporting natural regeneration of sub-canopy and shrub species. Therefore, the restoration of invaded forests requires eradication of invasive tree species to decrease propagule pressure and to stop decreases in the abundance of native tree species’ natural regeneration.


Study Design
We investigated the effects of invasive tree species within a set of 189 systematically distributed study plots, with nine plots in each of 21 blocks (Figure 1) [68]. The blocks were chosen according to the following criteria: presence of fruiting population of studied invasive tree species, homogeneity of stand and understory species composition, and their being surrounded by other forest types. Prior to the establishment of study plots, we examined WNP forest management plans to locate potential study sites, and made field visits to ensure they met these criteria before investigation. We established blocks in monocultures of invasive species (or in the case of P. serotina in P. sylvestris monocultures dominated by P. serotina). We stablished plots to assess the effects of decreasing propagule pressure on invasive tree species' natural regeneration [64]. The center of each block was one study plot, in the center of a monoculture of invasive species. We established four additional plots systematically, off of the North, South, East, and West sides of the central plot, nearly outside the alien species monoculture. Then, we established another set of four plots, located 30 m away from each of the four additional plots. We decided to establish study plots using a systematic design to avoid subjectivity. Each study plot was a 200 m 2 rectangle, resulting from pooling two 10 × 10 m plots to avoid potential relocation errors and pseudoreplications, as described in our previous study [68]. Due to the systematic design three study plots occurred in non-forest ecosystems, and were excluded from analyses (thus final n = 186; Figures A1 and A2). Although our study plots were arranged in a block design, due to a high level of spatial heterogeneity, they usually represented various soil types and vegetation types. As resource availability gradients were more important drivers of clustering among plots than block design (see details in [64]), we assumed that our study plots are independent observations. This system was developed to assess impacts of distance from seed source on abundance of alien tree species' natural regeneration. Therefore, edge plots could be affected by the proximity of invasive species-dominated stands. However, due to spatial heterogeneity and high abundance of alien tree species in WPN, both plot groups (near invaded stands and 30 m from invaded stands) can be invaded by other alien tree species (Figures A1 and A2).
We divided study plots into categories representing the main types of forest ecosystems in WNP, based on dominant tree species (i.e., with the highest basal area), soil type, plant association, and invasion level (Table 1). Due to different ecological requirements and management history, we divided Pinus sylvestris forests into two groups: poor and plantation. Poor P. sylvestris forests refer to mesic habitats, mostly growing on brunic arenosols, where this species naturally occurs. In contrast, P. sylvestris plantation refers to forests occupying more fertile soils, usually cambisols and luvisols, which naturally host Quercus-Acer-Tilia forests. We classified forests as invaded by P. serotina when its density in the shrub layer (individuals with height over 1.3 m) exceeded 500 ind. ha-1. We also classified forests as invaded by Q. rubra or R. pseudoacacia when its contribution to basal area exceeded 25% (however, it mostly exceeded 75%). Both Q. rubra and R. pseudoacacia were also intentionally planted in WNP [66,67]. We decided to divide forest types according to dominant species rather than into invaded/non-invaded due to dominant-specific effects, which are not included in measured variables, e.g., soil biota [32], decomposition rate [59] or annual litterfall [69]. The age (of dominant tree species) of the forests studied ranged from 23 to 139 years [70], however only ten plots had tree stands younger than 60 years old.

Data Collection
We investigated natural regeneration within study plots in July 2015. Within each study plot we counted seedlings (trees germinated in the current year) and saplings (one-year-old and older trees up to 0.5 m height) of each tree and shrub species, including both native and non-native species. Taxonomic nomenclature follows GBIF [71]. In total, we counted 54,593 plants, representing 53 species (Table A1). As it was impossible to assess whether a particular individual originated from seed or clonal growth without excavating the root system, we treated each shoot as an individual. This did not allow us to account for the clonal growth of some species.
Forests 2020, 11, 456 4 of 20 vegetation types. As resource availability gradients were more important drivers of clustering among plots than block design (see details in [64]), we assumed that our study plots are independent observations. This system was developed to assess impacts of distance from seed source on abundance of alien tree species' natural regeneration. Therefore, edge plots could be affected by the proximity of invasive species-dominated stands. However, due to spatial heterogeneity and high abundance of alien tree species in WPN, both plot groups (near invaded stands and 30 m from invaded stands) can be invaded by other alien tree species ( Figure A1,A2).   We characterized potential predictors influencing seedling survival using light availability, soil pH and soil C:N. We measured light availability in August 2016 using a LAI-2200 plant canopy analyzer (Li-Cor Inc., Lincoln, NE, USA). We measured diffuse non-interceptance (DIFN), the fraction of open canopy above the understory. Within each study plot, we took eight series of ten measurements. We decided to measure DIFN in summer, during the maximum canopy development, to account for maximum light availability limitation. We measured soil pH from six soil subsamples per study plot.
Subsamples were collected systematically along plot borders, from 0-10 cm depth. Then, we measured pH of soil solution in distilled water after 24 h using an electronic pH-meter. The same samples were analyzed for total carbon content and total nitrogen content, using the dry combustion method (PN-ISO 10694:2002) and Kjeldahl's method (PN-EN 16169:2012), respectively, in an accredited laboratory (Jars Sp. z.o.o., Legionowo, Poland). We calculated C:N ratio as a synthetic index of soil fertility. We also used natural regeneration density of the invasive species studied to show their effects on natural regeneration, at the level of competition between saplings.

Natural Regeneration Species Composition
All analyses were conducted using R software (The R Foundation for Statistical Computing Platform; Vienna, Austria, version 3.5.3) [72]. We used detrended correspondence analysis (DCA), as a method of indirect ordination, developed for analyzes across long environmental gradients [73]. Due to the long fertility gradient (from nutrient poor to rich soils) in our dataset and higher stability of DCA than multidimensional scaling, we decided to use DCA. We also conducted preliminary multidimensional scaling analysis, which had a poorer fit to the data than DCA. DCA is a modification of correspondence analysis, a common method of multidimensionality reduction in vegetation analyses. In DCA gradients are standardized by standard deviation, to avoid artifacts connected with gradient folding (so-called arch-effects or horseshoe-effects). We performed DCA using the vegan package [74], using study plots as analysis units and all density of all woody species (expressed in ind. 200 m −2 ) recorded in the plot as variables. To stabilize variance and overcome wide ranges of particular species densities, we log-transformed them prior to analysis. We also checked correlation of vectors representing soil pH, soil C:N ratio, and DIFN with ordination results using permutation tests (n = 999) to check how measured variables were correlated with gradients revealed by DCA.

Natural Regeneration Species Diversity
For each study plot we calculated taxonomic species richness and Shannon's diversity index, Faith's phylogenetic diversity (PD), mean nearest taxon distance (MNTD), and mean pairwise distance (MPD), and several indices of functional diversity: functional richness (FRic), divergence (FDiv), dispersion (FDis) and evenness (FEve). All these indices were calculated using both native and non-native species. To compute phylogenetic diversity indices, we extracted the phylogenetic tree covering all taxa noted in study plots from the V.PhyloMaker package [75]. We measured Faith's PD as a sum of all branch lengths in the plant community, MPD as mean pairwise distance between all taxa in the community and MNTD as mean of minimal distances between all taxa in the community using the PhyloMeasures package [76]. To avoid the influence of species richness on these variables, we standardized them according to the null model, i.e., mean and standard deviation (SD) of phylogenetic diversity metrics for a given species richness. For each set of phylogenetic tree tips, the observed metric has been standardized by subtracting the mean and dividing by SD. Null model mean and SD were calculated among all sets of tips that have the same number of elements as the tips set, with uniform probability among all possible tips samples of a given species richness [76]. To calculate functional diversity, we used a set of traits extracted from LEDA [77], BiolFlor [78], BIEN [79] and global wood density database [80], as well as Ellenberg ecological indicator values [81], as proxies of species requirements. In total we used 12 traits (ecological indicator values for light, temperature, moisture, soil reaction and nitrogen, flowering start and duration, pollination mode, specific leaf area, seed mass, height and wood density), with completeness from 39% (soil nitrogen and reaction indicators) to 100% (pollination mode, flowering, seed mass and height).
We assessed differences in species diversity between stand types using linear models accounting for soil pH, soil C:N ratio, DIFN, natural regeneration density of the invasive species studied, and forest type. We included soil and light variables in models to account for their impacts on natural regeneration.
We evaluated influence of forest type using Tukey's posteriori test, implemented in the multcomp package [82]. For species richness we used generalized linear models (GLM), assuming Poisson distribution of the dependent variable. Although models could reveal differences among forest types with p-values < 0.05, in Tukey's tests we applied single-step adjustment of p-values, to account for multiple hypothesis testing. Single-step adjustment decreases the probability of committing Type I error (i.e., rejection of the true null hypothesis), and also accounts for correlations among variables tested [82].

Natural Regeneration Density
We assessed impacts of forest type on natural regeneration density of the species studied, using zero-inflated Poisson GLM, implemented in the pcsl package [83]. This is a hurdle model, estimating both probability of occurrence (assuming binomial distribution) and abundance (assuming Poisson distribution). In our case, we assumed that species presence (binomial part of the model) depends on soil C:N ratio, soil pH, natural regeneration density of the invasive species studied and DIFN, while the count part of the model (information about species density)-also depends on forest type. Therefore, we assumed that tree species effect is more important for a particular species abundance than presence, accounting for impact of seed source, which is one of the most important drivers of natural regeneration [7,64]. As zero-inflated Poisson GLM provides estimates for both count and binomial parts of the distribution, standard posteriori tests are not available. Therefore, we used ceteris paribus comparisons, i.e., model predictions assuming all variables at constant (mean) levels with an exception of the considered variable [84,85]. In our case, we compared mean predictions for each forest type, assuming mean values of DIFN, natural regeneration density of the invasive species studied, soil pH and C:N ratio, to assess how forest type influences a particular species density in similar soil and light availability conditions. This analysis allows us to describe the dissimilarity among groups based on likelihood ratio tests [84]. For that reason, sometimes, GLM revealed p-values >0.05, while likelihood ratio tests revealed differences between forest types.

Impact of Alien Tree Species on Natural Regeneration Species Composition
The DCA ordination grouped study plots along two gradients, representing variability of soil pH and light availability ( Figure 2). Forest types overlapped each other in the ordination space, revealing continuity of vegetation. Along the DCA1 axis, forests were ordered from sites with high light availability, dominated by P. sylvestris, to the most shaded Quercus-Acer-Tilia and R. pseudoacacia forests. Axis DCA2 revealed shifts along soil reaction and fertility gradients. Composition of natural regeneration in particular forest types differed along the DCA1 axis and division into forest types was correlated with DCA results (r 2 = 0.359, p < 0.001). Points representing 50% of the species composition variability of P. sylvestris forests invaded by P. serotina overlapped between areas representing divided types of P. sylvestris forests-poor stands and plantations. Points representing Q. rubra forests partially overlapped the area of Q. petraea forests, which were close to poor P. sylvestris stands. Similarly, the area representing Q. rubra stands only partially overlapped the area of F. sylvatica forests, shifted more into the lower part of ordination space. Forests with R. pseudoacacia mostly overlapped the area represented by Quercus-Acer-Tilia forests, as well as P. sylvestris plantations, which were similar to almost all forest types.

Impact of Alien Tree Species on Natural Regeneration Diversity
We found that forest type had impacts on species richness, MNTD, FRic, FDis and FEve (Supplementary Table S1, Figure 3). Invaded and non-invaded stands differed most strongly in species richness in P. sylvestris plantations as invaded plots had 6.6 ± 0.7 species while non-invaded plots had 9.8 ± 0.5. Less distinct differences occurred in R. pseudoacacia forests, which had 9.6 ± 0.8 species, while Quercus-Acer-Tilia had 11.6 ± 0.8. However, these forest types did not differ according to the Tukey's test. Similarly, forests with Q. rubra had 7.4 ± 0.8 species, while Q. petraea had 9.4 ± 0.7 (but not significantly different at p = 0.05). Forests with F. sylvatica had the lowest species richness -4.2 ± 0.4 species. Similar trends (but not included in final models) occurred for the Shannon index. No differences in phylogenetic diversity were revealed by Tukey tests. However, lower values were found in R. pseudoacacia than in Quercus-Acer-Tilia, and higher values in invaded than non-invaded P. sylvestris plantations. For functional diversity we found that only invaded P. sylvestris plantations had lower FRic than other types of forests (with the exception of F. sylvatica forests). For FDis we found the lowest values in Quercus-Acer-Tilia forest, and the highest in non-invaded P. sylvestris plantations. In contrast, we found the highest FEve in the invaded P. sylvestris plantations while the lowest was in R. pseudoacacia forests and non-invaded P. sylvestris plantations. We did not find any

Impact of Alien Tree Species on Natural Regeneration Diversity
We found that forest type had impacts on species richness, MNTD, FRic, FDis and FEve (Supplementary Table S1, Figure 3). Invaded and non-invaded stands differed most strongly in species richness in P. sylvestris plantations as invaded plots had 6.6 ± 0.7 species while non-invaded plots had 9.8 ± 0.5. Less distinct differences occurred in R. pseudoacacia forests, which had 9.6 ± 0.8 species, while Quercus-Acer-Tilia had 11.6 ± 0.8. However, these forest types did not differ according to the Tukey's test. Similarly, forests with Q. rubra had 7.4 ± 0.8 species, while Q. petraea had 9.4 ± 0.7 (but not significantly different at p = 0.05). Forests with F. sylvatica had the lowest species richness-4.2 ± 0.4 species. Similar trends (but not included in final models) occurred for the Shannon index. No differences in phylogenetic diversity were revealed by Tukey tests. However, lower values were found in R. pseudoacacia than in Quercus-Acer-Tilia, and higher values in invaded than non-invaded P. sylvestris plantations. For functional diversity we found that only invaded P. sylvestris plantations had lower FRic than other types of forests (with the exception of F. sylvatica forests). For FDis we found the lowest values in Quercus-Acer-Tilia forest, and the highest in non-invaded P. sylvestris plantations. In contrast, we found the highest FEve in the invaded P. sylvestris plantations while the lowest was in R. pseudoacacia forests and non-invaded P. sylvestris plantations. We did not find any impacts of natural regeneration density of the invasive species studied on any biodiversity index included in the study. impacts of natural regeneration density of the invasive species studied on any biodiversity index included in the study. Letters denote groups which did not differ at confidence level α = 0.05 after multiple hypotheses adjustment (see Table S1), n.a. indicates that forest type has been not included in the final model for a particular variable.

Impacts on Natural Regeneration Densities of the Most Frequent Native Tree Species
Analyses of regeneration densities for particular tree species revealed that eight of nine dominant or co-dominant tree species reached higher densities in their native optimum forest types than in invaded forests (Figure 4, Supplementary Table S2). The exception was Ulmus minor, which reached the highest densities in R. pseudoacacia forests. However, some species reached similar densities in invaded forest types, not paired with their optimum native types, for example Carpinus betulus had similar densities in invaded P. sylvestris plantations and Tilia cordata in invaded P. sylvestris poor forests. Three of the six shrub species studied reached their highest densities in invaded forests, two (Corylus avellana, Frangula alnus) in R. pseudoacacia forests and one (S. aucuparia) in P. sylvestris plantations invaded by P. serotina. In contrast, Sambucus nigra reached the highest densities in Quercus-Acer-Tilia forest and Prunus avium and Euonymus europaea in P. sylvestris plantations. We found that sapling densities of the invasive species studied negatively influenced the binomial part of density models for three species: A. platanoides, P. sylvestris and S. aucuparia, and negatively influenced the count part of models for A. platanoides and S. nigra. For A. pseudoplatanus, C. betulus, E. europaea, F. sylvatica, F. excelsior, P. sylvestris, Q. petraea, S. aucuparia, and T. cordata, we found positive but low influences of sapling densities of the invasive species studied on natural regeneration density (estimate ± SE varied from 0.0016 ± 0.0008 to 0.0078 ± 0.0004, Supplementary  Table S3). FEve-functional evenness, FDiv-functional divergence) within studied forest types. Letters denote groups which did not differ at confidence level α = 0.05 after multiple hypotheses adjustment (see Table S1), n.a. indicates that forest type has been not included in the final model for a particular variable.

Impacts on Natural Regeneration Densities of the Most Frequent Native Tree Species
Analyses of regeneration densities for particular tree species revealed that eight of nine dominant or co-dominant tree species reached higher densities in their native optimum forest types than in invaded forests (Figure 4, Supplementary Table S2). The exception was Ulmus minor, which reached the highest densities in R. pseudoacacia forests. However, some species reached similar densities in invaded forest types, not paired with their optimum native types, for example Carpinus betulus had similar densities in invaded P. sylvestris plantations and Tilia cordata in invaded P. sylvestris poor forests. Three of the six shrub species studied reached their highest densities in invaded forests, two (Corylus avellana, Frangula alnus) in R. pseudoacacia forests and one (S. aucuparia) in P. sylvestris plantations invaded by P. serotina. In contrast, Sambucus nigra reached the highest densities in Quercus-Acer-Tilia forest and Prunus avium and Euonymus europaea in P. sylvestris plantations. We found that sapling densities of the invasive species studied negatively influenced the binomial part of density models for three species: A. platanoides, P. sylvestris and S. aucuparia, and negatively influenced the count part of models for A. platanoides and S. nigra. For A. pseudoplatanus, C. betulus, E. europaea, F. sylvatica, F. excelsior, P. sylvestris, Q. petraea, S. aucuparia, and T. cordata, we found positive but low influences of sapling densities of the invasive species studied on natural regeneration density (estimate ± SE varied from 0.0016 ± 0.0008 to 0.0078 ± 0.0004, Supplementary Table S3).  . Mean (± SE) predicted density (ind. 200 m −2 ) of the most frequent native tree (first three rows) and shrub (last two rows) species within studied forest types. Letters denote groups for which predicted responses did not differ at confidence level α = 0.05 after likelihood ratio tests (for full models see Table S2). . Mean (± SE) predicted density (ind. 200 m −2 ) of the most frequent native tree (first three rows) and shrub (last two rows) species within studied forest types. Letters denote groups for which predicted responses did not differ at confidence level α = 0.05 after likelihood ratio tests (for full models see Table S2).

Continuity in Species Composition of Natural Regeneration
Our study revealed low differences in species composition between similar types of invaded and non-invaded forests. This might be caused by two issues: vegetation continuity and the spread of invasive tree species. Due to systematic design of our study, we sampled vegetation patches which did not always resemble textbook examples of plant communities [86,87]. Therefore, our study design covered a continuum of plots from invaded to non-invaded, as well as light availability and soil fertility gradients [88]. A continuous transition between invaded and uninvaded plots could be regarded in conditions of only one invasive species. However, in the study area, spatial heterogeneity of geomorphology and dominant tree species led to overstory dominant-specific impacts of particular native tree species. For example, F. sylvatica canopies, despite functional similarity to Q. petraea, transmit less light, and therefore reduce species richness and abundance. This did not allow us to apply a simple invaded-uninvaded approach [89] to assess effects of the invasive species studied, due to diverse types of uninvaded ecosystems. A more important factor is spread of invasive tree species' natural regeneration to the non-invaded plots. High densities and biomass of invasive tree species [64,90], recorded also in plots with non-invaded canopies, contributed to the biotic homogenization of regeneration layers [91,92]. DCA ordination scores of invasive tree species studied showed that Q. rubra and Q. petraea, as well as R. pseudoacacia and Acer platanoides, occupied similar niches, regardless of forest type. Due to high dispersal capacity (over 500 m) of wind-dispersed R. pseudoacacia [55] and animal dispersed P. serotina [93,94] and Q. rubra [61,95], the species studied were able to effectively colonize all study plots. Similarly, Šibíková et al. [40] confirmed the homogenizing effect of R. pseudoacacia on forest vegetation. However, models did not reveal influences of sapling densities of the invasive species studied on natural regeneration diversity indices, and found limited impacts on regeneration densities of particular species. This would indicate a low fulfilment of ecological niches within study sites, providing opportunity for growth of numerous species. Low impacts of invasive species saplings densities suggests that native tree species are limited by overstory-understory interactions with invasive species rather than by competition in the regeneration layer.

Low and Context-Dependent Impacts of Particular Invasive Species on Natural Regeneration Diversities
Our study confirmed the widely-known context-dependence of invasive species impacts [96][97][98]. For example, P. serotina decreased functional richness in P. sylvestris plantations, while we did not observe this in poor stands of P. sylvestris. This might be connected with higher likelihoods of P. sylvestris plantations occupying rich soils to be colonized by multiple species [99,100]. This might be connected with both higher light transmittance [26,101], as well as higher rates of bird visitation, and resulting defecation of bird-dispersed seeds [94,102]. Pinus sylvestris plantations provide suitable conditions for P. serotina growth which, due to ability for quick growth, is able to overgrow native tree species [103,104]. In poor stands, where soil fertility is a limiting factor for tree species, light competition is lower and despite the limited pool of species able to occur, functional richness is similar in both invaded and non-invaded stands. Context-dependence is also clear when comparing plantations of native (P. sylvestris) and non-native (Q. rubra and R. pseudoacacia) species, as plantations of P. sylvestris, due to higher light transmittance, host more species than Q. rubra. Dissimilarity in natural regeneration species composition between P. sylvestris and R. pseudoacacia plantations may result from higher proportion of shade-tolerant and nitrogen-demanding species beneath R. pseudoacacia canopies.
We found lower species richness and higher functional richness and dispersion in R. pseudoacacia and Quercus-Acer-Tilia stands (although with p > 0.05, but biologically relevant). Due to nitrogen fixation, forests with R. pseudoacacia have been reported with similar understory species richness as native forests [30]. In our case, we found lower species richness in R. pseudoacacia forests, which might be connected with our focus on woody species rather than understory herbaceous species. Differences between Q. rubra and Q. petraea forests also were low. This might be connected with functional similarity of both dominant tree species. Although Q. rubra is known to transmit less light than Q. petraea or Q. robur [25,58], we accounted for light availability in our models, to exclude this effect from comparisons among forest types.
Low effects on community metrics of biodiversity resulted from homogenization and replacement of native tree species by alien species. Moreover, this might be connected with lower survival probability of both alien and native tree species in invaded stands [50,68]. Our results also revealed a lack of differences in phylogenetic diversity among forest types. This is in line with Piwczyński et al. [105], explaining this pattern of shaping phylogenetic structure of the plant community by higher importance of soil pH, fertility, and light availability, rather than species identity. Due to the higher differences in particular species densities than community structure among forest types studied, phylogenetic diversity based on a species presence-absence matrix might not be sensitive enough to detect the influence of invasive tree species.

Decreasing Ability of Native Trees to Regenerate
Our study revealed that native forest-forming species usually reached lower densities in invaded forest types than in non-invaded. This is in line with Kowarik et al. [14], who found that forest-forming tree species occurred with higher frequency in plots with native tree species than with alien. Similarly, the invasive shrub Ligustrum lucidus decreases abundance of native tree species in subtropical Argentina [48]. Our results contradict Terwei et al. [45], who found a negative impact of P. serotina on U. minor natural regeneration and no impact of R. pseuodacacia. Also, Terwei et al. [45] did not find either a positive impact of P. serotina or a negative impact of R. pseudoacacia on C. betulus. The only similar result is the impact of R. pseudoacacia on Q. petraea natural regeneration, in line with the positive impact of R. pseudoacacia on Q. robur [45]. The differences might be connected with the different type of ecosystem studied, i.e., riparian forest in northern Italy.

Study Limitations
Despite accounting for the main gradients of species diversity in the vegetation-soil fertility, pH and light availability [90], our results might depend on other factors not studied. The most important is basing the study on one year of observation. Our earlier study within the same study plots revealed that seedling survival rates have significant seasonal variation [68]. However, in this study, we accounted for seedlings and saplings (i.e., natural regeneration which survived one year), which are more persistent [5,45]. We could also miss the effect of higher light availability in spring, which could affect germination, as we did not measure light availability before canopy development. Some minor bias might be caused by mismatch between year of natural regeneration investigation and light availability measurement. However, we did not report significant changes in canopy cover in both study years. Another factor is not accounting for deer browsing, which is an important limiting factor of natural regeneration [11,106]. Although we did not account for dispersal limitation [44,102], due to proximity of study plots to heterogenous stands, we assumed lack of seed dispersal limitation caused by the distance to the propagule sources. Another potential limitation might be connected with study design-all studied stands were located near invaded stands, thus some of them might be affected by edge effects. On the other hand, spatial heterogeneity and proximity of other invaded stands did not allow us to distinguish clear edge zones, and we may treat them as different stands. Bearing in mind all acknowledged drawbacks, we believe that results of our study allow for inference about the impacts of the invasive tree species studied on natural regeneration in temperate European forests.

Conclusions
Our study revealed context-dependent influences of invasive tree species on natural regeneration. We found that similar invaded and non-invaded forest types have similar seedling and sapling banks, due to homogenization caused by natural regeneration of invasive trees. We found differences in natural regeneration in terms of taxonomic and functional diversity among forest types, connected with environmental filtering, by both alien and native tree species. Invasive tree species decreased the extent of natural regeneration of native forest-forming tree species, while they supported the natural regeneration of some shrub or sub-canopy species. The latter occupy regeneration niches of forest-forming species. This led to the persistence of invasive tree species by decreasing ecosystem resilience. Therefore, the restoration of invaded forest patches requires the eradication of invasive tree species, not only as propagule sources, but also as factors decreasing the ability for the spontaneous regeneration of native tree species.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/4/456/s1, Table S1: Models of species diversity as functions of forest type, DIFN, invasive sapling density, soil C:N ratio and soil pH (best fit model, according to AIC). Species richness was estimated using Poisson generalized linear models (GLM) while other metrics of biodiversity -by linear regression; Table S2: Models of particular species' natural regeneration density (ind. 200 m −2 ) as functions of forest type, DIFN, invasive sapling density, soil C:N ratio and soil pH, estimated using zero-inflated Poisson generalized linear models.

Acknowledgments:
We are grateful to the two anonymous Reviewers for helpful comments on the earlier draft of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Spatial Distribution of Forest Types within Blocks and Species Frequency within Studied Forest Types
Forests 2020, 11, x FOR PEER REVIEW 13 of 21 Figure A1. Schematic distribution of forest types within study blocks 1-12. Schemes do not account for distances between plots. For details of study design see Figure 1. Figure A1. Schematic distribution of forest types within study blocks 1-12. Schemes do not account for distances between plots. For details of study design see Figure 1.  Figure A2. Schematic distribution of forest types within study blocks 13-21. Schemes do not account for distances between plots. For details of study design see Figure 1. Table A1. Frequency (%) of species recorded during the study within forest types (see Table 1 for explanations).