Drivers of invasive tree and shrub natural regeneration in temperate forests

We assessed drivers of ecological success along resource availability gradients for three invasive woody species: Prunus serotina Ehrh., Quercus rubra L. and Robinia pseudoacacia L. We aimed to check how much of invasion success, measured by invader biomass, is explained by propagule pressure and plant community invasibility. Using 3 years of observations from 372 study plots (100 m2 each) in temperate forests of Wielkopolski National Park (Poland) we investigated the hierarchy of predictors and partial dependencies using the random forest method. Our study indicated that propagule pressure explained more variance in success of invaders than invasibility—describing availability of resources and competitors in understory vegetation. We also found different responses of seedlings and saplings, connected with dependence on stored carbohydrates, which decreased seedling responses to resource availability gradients. However, resource availability (light and leaf litter predictors) had greater influence than predictors describing understory vegetation. Based on importance and response strength the species studied may be arranged by decreasing requirements for soil fertility and acidity: P. serotina < Q. rubra < R. pseudoacacia, whereas for light requirements and competition vulnerability the order is: P. serotina > Q. rubra > R. pseudoacacia. However, low light requirements of R. pseudoacacia may be biased by high proportion of sprouts supplied by parental trees. Results provide guidelines for effective management of invasive woody species in forest ecosystems and describe complex interactions between factors studied on ecological success of invaders.


Introduction
Among functional groups of invasive plant species, woody plants are distinguished by their long lifespan (Richardson 1998;Richardson et al. 2000;Richardson and Rejmánek 2011). For that reason their establishment and first year of life is crucial for further development and dispersal (Baraloto et al. 2005). Thus, as dispersal is the key factor that allows naturalized species to become invasive (Richardson et al. 2000;Chytrý et al. 2008;Dyderski and Jagodziński 2016), ecological success of young regeneration determines the ultimate success of woody invaders. In comparison with short-lived (herbaceous) taxa, assessment of woody species requires more complex studies, also accounting for more predictors and potential interactions (e.g. Richardson 1998;Pyšek et al. 2014;Brundu and Richardson 2016).
Natural regeneration of both native and invasive species faces numerous limiting factors. The abiotic conditions most limiting to growth of young trees are usually water, light and nutrient availability (Niinemets and Valladares 2006;Ellenberg and Leuschner 2010). However, levels of these factors are also influenced by other species co-occurring within the plant community (Tilman 1986). Competitive potential of a plant species is a derivate of its life strategy (Westoby 1998;Grime 2006), which results from its functional traits and their level of fitness to environmental conditions, particularly expression of resource acquisition strategies, e.g. specific leaf area, leaf photosynthetic capacity or specific root length (Cornelissen et al. 2003;Pierce et al. 2013;Díaz et al. 2016;Kunstler et al. 2016). Natural regeneration of woody species competes with different species and may ''win'' and reach the next level of development only in particular conditions. Thus, ecological success of invasive species depends on characteristics of the recipient plant community-habitat invasibility, which is determined by level of resources (Davis et al. 2000;Godefroid et al. 2005;Funk 2008;Paquette et al. 2012;González-Muñoz et al. 2014), and potential competitors (Chmura 2004;Sierka 2005, 2007;Godefroid et al. 2005;Knight et al. 2008).
Biological invasions are driven by three connected factors: habitat invasibility, alien species invasiveness and propagule pressure-quantity and quality of propagules able to establish in new sites (Richardson et al. 2000;Davis et al. 2005;Jeschke 2014). Among them, propagule pressure is the obligatory factor determining invasion success (Lonsdale 1999;Lockwood et al. 2005;Vanhellemont et al. 2009). Numerous studies have confirmed its importance at both large (Křivánek et al. 2006;Pyšek et al. 2009Pyšek et al. , 2015Woziwoda et al. 2014) and small spatial scales (Pairon et al. 2006;Deckers et al. 2008;Vanhellemont et al. 2009;Jagodziński et al. 2015;Bonilla and Pringle 2015). Propagule pressure seems to be a crucial driver of relationships between habitat invasibility and ecological success of invasive species. Although most of these interactions were studied within relatively short resource gradients, covering few habitat types (Vanhellemont et al. 2009;Aslan et al. 2012;Terwei et al. 2013;Jagodziński et al. 2015), Davis et al. (2005) conceptualized the relationship between invasibility and propagule pressure as the concept of invasion pressure, which provided a framework to merge both groups of factors. Our previous study (Jagodziński et al. 2018a) also revealed the interactions between propagule pressure and invasibility, but only for one species and in an experimental system.
Despite existing knowledge of factors determining success of biological invasions, it remains unknown whether interactions between, and relative importance of, propagule pressure and invasibility across a wide range of forest types are similar to the limited types reported on by previous studies. Comprehensive review of the community ecology of invasive species (Gallien and Carboni 2017) indicated a lack of clarity on whether processes filtering species colonizing new communities vary along environmental gradients. Moreover, there is an insufficient number of studies on invasive species spread conducted in protected areas, which are especially threatened by biological invasions (Hulme et al. 2014).
We aimed to assess drivers of biomass along resource availability gradients for the three most frequent invasive woody species in temperate Europe (Wagner et al. 2017)-Prunus serotina Ehrh., Quercus rubra L. and Robinia pseudoacacia L.-to assess how much of invasion success (moving along an introduction-naturalization-invasion continuum), measured by invader biomass, is explained by propagule pressure and plant community invasibility. We hypothesized that: (1) propagule pressure explains more variance in ecological success of invaders than invasibility, according to previous studies (e.g. Lonsdale 1999;Jagodziński et al. 2015Jagodziński et al. , 2018a, (2) total biomass of seedlings depends more on propagule pressure than for saplings, as seedlings are more dependent on dispersal than older plants, which require suitable conditions for survival (e.g. Beckage et al. 2005;Knight et al. 2008;Rodríguez et al. 2017), and (3) understory competition, expressed by functional diversity components, have larger effects on ecological success of the invasive species studied than surrogates of resource availability, as light competition limits occurrence of the species studied (Cierjacks et al. 2013;Jagodziński et al. 2015Jagodziński et al. , 2018a. Because the three species studied differ in biology and ecology, we decided to characterize their invasiveness by comparing their responses to propagule pressure and habitat predictors on the background of their life history traits.

Studied species
In European woodlands three alien woody species are most frequent: Prunus serotina Ehrh., Quercus rubra L. and Robinia pseudoacacia L. (Wagner et al. 2017). These species have been broadly recorded as invasive in Europe and all of them were widely introduced via forestry, in the eighteenth, nineteenth and seventeenth century, respectively (Muys et al. 1992;Woziwoda et al. 2014;Vítková et al. 2017). These species came from eastern North America and are widely distributed in Western and Central Europe, but their frequency decreases eastwards. P. serotina in its native range is a tree with valuable timber, but in Europe its introduction was not successful either in terms of timber production or soil improvement (Muys et al. 1992;Starfinger et al. 2003;Aerts et al. 2017). Its presence strongly influences nutrient cycling (Aerts et al. 2017;. P. serotina is dispersed mainly by birds, up to 600 m from the seed source (Pairon et al. 2006;Jagodziński et al. 2015), however ca. 80% of fruits fall beneath the crown of the mother tree (Pairon et al. 2006). Q. rubra produces medium-value timber in Europe and due to its long lag-time is considered a weakly invasive species. Its dispersal is limited, as acorns mainly fall close to the parents and are not preferred by birds (Myczko et al. 2014;Bieberich et al. 2016), which are the main long-distance vectors of this species. R. pseudoacacia is a wind-dispersed pioneer tree species, important for wood production and providing nectar sources for pollinators, and is also widely used as an ornamental tree (Cierjacks et al. 2013;Vítková et al. 2017). This species transforms soils, due to symbiosis with nitrogen-fixing organisms (Rice et al. 2004). All these species invade forest ecosystems of temperate Europe, ranging from less fertile coniferous sites to the most fertile riparian sites.

Study area
We conducted our study in Wielkopolski National Park (WNP; W Poland; 52816'N, 16848'E). WNP covers 7584 ha and conserves mostly forest ecosystems and very diverse geomorphology connected with the last glaciation. According to the nearest meteorological station in Poznań (c.a. 15 km from WNP) mean annual temperature in 1951-2010 was 8.4°C and annual precipitation was 521 mm. Forests of WNP were heavily transformed by human activity, especially by former forest management, replacing mixed and broadleaved forests with monocultures of Scots pine. Moreover, before WNP establishment in 1957, the area was a place of numerous introductions of alien trees and shrubs, thus WNP is the national park in Poland with the highest number of alien woody species (Purcel 2009;Gazda and Szwagrzyk 2016). For that reason WNP provides a wide range of soil fertility and tree species composition, as well as variable propagule pressure of alien tree species, which makes WNP a good area for invasion ecology studies. All three species studied are abundant and are present both as an admixture in the tree stands and as monoculture tree stands (Purcel 2009). These advantages, together with known history of forest management and alien species introductions, make WNP a valuable area, especially for studies on biological invasions in forest ecosystems.

Study design
We used a set of 378 plots (squares 100 m 2 in area) arranged in 21 blocks: nine for Q. rubra and six for P. serotina and R. pseudoacacia. We established more blocks for Q. rubra due to lower densities of this species, resulting in lower detection rates. Moreover, within blocks only the central part of each block is located in a parental (monoculture) stand of invasive species (for P. serotina, which occurs only as an admixture in tree stands, these are tree stands with high density of fruiting trees), and the remaining plots are also invaded by the other two species. One pair of plots was located within the invasive species monoculture (maternal stand). A second pair of experimental plots was located along each of four sides of the Drivers of invasive tree and shrub natural regeneration in temperate forests 2365 maternal stand (N, S, E, and W), nearly outside the stand, at the invasion edge (Rodríguez et al. 2017), and a third pair of experimental plots was located 30 m further out from the second set of plot pairs (Fig. 1). This design produced 18 experimental plots with different distances from propagule sources-within the maternal stands, outside maternal stand borders and at a distance of 30 m from maternal stands. Due to the high level of spatial heterogeneity in tree stands, longer distances may result in biases in estimation of the distance effects. Therefore, four classes of distance from propagule source were distinguished: in the propagule source, near the propagule source, 30 m from the propagule source and further than 30 m (plots without a propagule source of the considered species in the nearest neighborhood). Despite concerns, this did not generate pseudoreplications, due to high spatial heterogeneity manifested in high variability of natural regeneration within plots. Due to systematic distribution of the plots, six of them were located in non-forested vegetation paths and these plots were excluded from analyses, thus the final number of plots was n = 372.  , with a center in the propagule source-an invasive species monoculture or for P. serotina a tree stand with a dense fruiting P. serotina layer light availability ranged from 0.007 to 0.251 with an average of 0.044 ± 0.001.

Data collection
In July of 2015, 2016 and 2017 we counted the number of seedlings (defined as individuals germinated in a particular year) and saplings (defined as individuals at least one year old and with height \ 0.5 m) within each plot, and in 2015, 2016 and 2017 we also measured root collar diameters (RCD) and heights (H) of the seedlings and saplings. This data was used for biomass calculation for each year. We treated sprouts (i.e. specimens generated through vegetative reproduction from root suckers) as saplings, due to their usually larger dimensions and connectivity with their parental organism. After measuring 28,703 seedlings, we found that height measurements of up to 30 seedlings is enough to cover the range of variance in seedling heights, and thus in 2017 for that age class we measured up to 30 heights, due to high seedling densities in some study plots (in cases of 50 plots with P. serotina and 23 with R. pseudoacacia).
For these plots we used mean dimensions of seedlings to calculate individual biomass and multiplied by seedling density to produce biomass of seedlings. This did not affect calculated biomass, due to low SE of the mean for seedling total biomass (TB). In total we measured 39,664 plants. During these inventories we also investigated vegetation species composition and abundance using the Braun-Blanquet method: for each plant species in the understory (forest layer up to 50 cm) we assigned one of nine cover degrees in an alphanumeric scale. Moreover, in 2015 we assessed tree stand structures for study plots within larger (0.02-0.20 ha) plots, where diameter at breast height was measured for all trees and basal area (BA) was computed for each tree species. In August 2016 we measured canopy openness index (diffuse non-interceptance; DIFN) in each plot, using an LAI-2200 plant canopy analyzer (Li-Cor Inc., Lincoln, NE, USA). For each plot we recorded four series of ten samples, using the methodology of Machado and Reich (1999). We chose August for the measurement time, as our previous studies showed that canopy openness was lowest during this month due to maximal canopy foliage development (Knight et al. 2008); thus, we accounted for the minimal light availability within our model. Light was an important predictor in previous invasion ecology studies (e.g. Knight et al. 2008;Jagodziński et al. 2018a). To account for rates of nutrient cycling and acidity we used litter mass and pH. We decided to use litter predictors instead of soil, as litter is the main source of nutrients for the uppermost soil horizons. However, these predictors are only proxies for direct measurements of soil predictors. Litter was sampled in March 2017, when its amount was stabilized after winter. For each pair of plots we collected four samples from circular plots (0.16 m 2 ). Woody debris with diameter [ 1 cm was excluded from litter samples. Litter samples were dried in an oven at 65°C to constant mass and visually divided into two parts: recognizable and decomposed (unrecognizable) parts of litter. Next we weighed both samples using a balance with an accuracy of 1 g and determined proportion of the decomposed part of litter, which is usually higher in tree stands with low rates of organic matter cycling. Litter pH was assessed using a pH-meter in distilled water solution after 24 h.
To assess biomass of natural regeneration we destructively sampled 647 trees in July 2017 (Table S1). We used biomass as a measure of ecological success, because this predictor increases with both increasing density and dimensions, reflecting space filled by the species, and is also claimed to be a good measure of plant fitness (Younginger et al. 2017). For each pair of plots we surveyed an area with radius of 5 m around the plot borders and randomly selected up to five specimens for each species, according to available number of specimens and species densities within plots. We excluded heavily damaged and browsed plants, unless there were no alternative specimens in the area examined. This accounted for the joint effects of lower growth and resistance to herbivory in suboptimal sites. This approach also resulted in unequal numbers of sample trees per species: 356 sample trees of P. serotina (195 seedlings and 161 saplings), 133 of Q. rubra (72 seedlings and 61 saplings) and 158 of R. pseudoacacia (94 seedlings and 64 saplings). These unequal numbers resulted from unequal distributions of natural regeneration of the species studied within the study plots. Each sample tree was dug up, cleaned and divided into roots, stems with branches and leaves. We excluded acorns which were still attached to Q. rubra and Q. petraea seedlings from the total biomass (TB). For each sample tree we also measured RCD and H.
Within the sampled tree dataset we found that 14 of 161 (8.7%) P. serotina and 23 of 64 (35.9%) R. pseudoacacia saplings had root suckers, and for these sample trees we did not include belowground biomass, because their root system was part of the parental plant. However, despite this belowground biomass exclusion, we included them in the dataset because their aboveground biomass also contributed to the understory. We did not assess whether each of the plants that were not destructively sampled were vegetative or generative reproduction, to not influence their survival probability, which would undermine future usage of permanent plots.

Data analysis
All analyses were conducted using R software (R Core Team 2017). Prior to modelling datasets were centered, scaled and processed using Yeo-Johnson power transformations (Yeo and Johnson 2000) to stabilize variance, increase normality of distributions and overcome problems with different magnitudes of variables. This preprocessing was carried out using the caret::preProcess() function (Kuhn 2008). As our data has hierarchical structure, which may influence model outcomes (Roberts et al. 2017), we decided to assess potential blocking effects of the plot layout in the field. We checked whether the most important predictors describing resource availability-DIFN, litter mass and pH-showed clustering related to the blocks. We performed k-means clustering using six clusters (the optimal cluster number was estimated by the elbow method) using the stats::kmeans() function. Then we compared clusters with principal components analysis of centered and scaled values of the analyzed predictors and we visually inspected how the blocks were arranged within the ordination space ( Figure S1). This analysis showed that potential blocking effects were related to the availability of resources rather than spatial proximity, and therefore for further analyses we decided to use modelling methods that did not account for the arrangement of blocks.
To predict TB of the species studied we used allometric equations, following Jagodziński et al. (2018a, b). Because our sample trees were collected within plots and blocks, they differed in dimensions and plot design expressed different resource availabilities rather than spatial structures ( Figure S1). We quantified how much variance in TB was explained by dimensions and by resource availability predictors (litter mass, pH and DIFN) using boosted regression tree models (Elith et al. 2008). Influence of the latter predictor explained 10.2, 9.0 and 0.0% for Q. rubra, P. serotina and R. pseudoacacia, respectively, thus we decided to use the simplest approach-linear and non-linear allometric models. From ten formulas of allometric relationships used by Jagodziński et al. (2018a, b) we chose the model with the lowest AIC (Table S2). For P. serotina TB = 0.03917 9 DRC 2.13334 9 H 0.49044 (with mean error of 0.059 and R 2 = 0.880), for Q. rubra TB = 0.004935 9 DRC 1.672448 9 H 1.371159 (with mean error of 0.067 and R 2 = 0.852) and for R. pseudoacacia TB = 0.004002 9 DRC 1.078374 9 H 1.203052 (with mean error of 0.013 and R 2 = 0.955).
We analyzed vegetation patterns using functional traits provided by BiolFlor (Klotz et al. 2002), BryoAtt (Hill et al. 2007) and LEDA (Kleyer et al. 2008) databases, as well as ecological indicator values (EIV; Ellenberg and Leuschner 2010). For each plant species we extracted EIV for light, fertility, soil reaction, moisture, continentality and temperature, as well as canopy height, leaf mass, size, dry matter content, specific leaf area (SLA), growth form, seed mass and number per shoot, reproduction mode and Grime's (Grime 2006) life strategy (Table S3). These traits were used for calculation of functional diversity components: functional dispersion (FDis), functional divergence (FDiv), functional evenness (FEve) and functional richness (FRic). These components describe distribution of species' functional traits for each sampled plant community within a hypervolume of traits (Mason et al. 2005;Laliberté and Legendre 2010;Pla et al. 2011), indicating prevalence of processes shaping species composition of the community. High values of FDiv indicate numerous functional ways of resource acquisition, which is associated with higher competition, similar to high FEve, indicating lack of one dominant type of resource acquisition. High values of FRic and FDis indicate numerous functional plant types and low level of habitat filtering in the plant community (Kotowski et al. 2010;Hedberg et al. 2014;; but see Kraft et al. 2015). Moreover, to characterize dominant plant strategies, we used three communityweighted mean (CWM) values of functional traits: SLA, height and seed mass, according to the LHS concept (Westoby 1998) and global spectrum of plant functions (Díaz et al. 2016). All functional diversity components and CWMs used for analyses of particular invasive species studied were calculated excluding that species, to avoid circular reasoning and biased variable importance, connected with different functional traits (Thomsen et al. 2016), for example, overestimated importance of seed mass CWM for Q. rubra, with its high value of this trait. Usage of functional diversity components allows conclusions about impact of plant community on the species of interest (Kurokawa et al. 2010;D'Astous et al. 2013;Jagodziński et al. 2017;Czortek et al. 2018).
In our dataset some predictors were correlated, because they reflect different aspects of phenomena studied. For example, propagule source presence/ absence and basal area correlation coefficient r was 0.56 for P. serotina, but the former predictor reflected mere presence of seed source and the latter predictor reflected its quantity. For that reason and for better ecological interpretation we decided not to exclude them from analyses. Only three pairs of predictors had correlation coefficients (r) [ 0.6-species richness and Shannon's index (0.73), FRich and species richness (0.65) and litter mass and proportion of decomposed part of litter (0.6). To ensure that our results were not biased by year-to-year dynamics, in modelling we used averaged vegetation predictors and biomasses for each plot. This averaging was also necessary as some predictors, e.g. tree stand and litter predictors as well as DIFN, were measured only once and combining them with more variable predictors would bias their importance. To assess influence of factors describing resource availability, propagule pressure and competitors we used the random forest method (RF; Breiman 2001), trained using the caret::train() function (Kuhn 2008). We chose RF due to its good performance in cases of collinear variables, non-normal distributions of predictors and high predictive power. As an alternative, we also tested boosted regression trees (Elith et al. 2008), which showed lower fitness-lower coefficients of determination and higher RMSE (Fig. S2). Due to differences in survival, we assessed biomass of seedlings, saplings and total natural regeneration of each species separately. For each RF model we produced two results/outputs -importance of variables and partial dependence plots. Variable importance is expressed by mean decrease in accuracy (%IncMSE), which is a percentage increase of mean squared error of the result when the considered variable was permuted. Partial dependence plots express impact of a single variable on a response when all remaining predictors within a given model are constant.

Biomass of natural regeneration
Natural regeneration of R. pseudoacacia, Q. rubra and P. serotina, was found in 165, 194 and 239 plots, respectively. Biomass of natural regeneration for these three species within study plots ranged from 0.00 to 257.19 kg ha -1 , with an average of 1.52 ± 0.36 kg ha -1 (Fig. 2). However, medians for each species seedlings and saplings did not exceed 1 kg ha -1 . For P. serotina total biomass ranged from 0 to 9.61 kg ha -1 , with an average of 0.81 ± 0.09 kg ha -1 ; for Q. rubra from 0.00 to 257.19 kg ha -1 , with an average of 1.04 ± 0.69 kg ha -1 and for R. pseudoacacia from 0.00 to 1.81 kg ha -1 , with an average of 0.07 ± 0.01 kg ha -1 . Mean proportions of total seedling biomass within plots were 30.7 ± 2.3, 29.8 ± 2.7 and 48.7 ± 3.2% for P. serotina, Q. rubra and R. pseudoacacia, respectively.
Hierarchy of predictors RF models explained from 56.5% (R. pseudoacacia seedlings) to 78.3% of variance (P. serotina total) in the biomass of the species studied and included from 3 to 10 predictors ( Table 1). Two of the most important predictors for each species and biomass type were presence of parental trees in the plot or in the nearest neighborhood and basal area of parental trees in the tree stand (Table 1). Exceptions were saplings and total biomass of Q. rubra, for which the second most important predictor was litter mass. Nevertheless, %IncMSE connected with nearest presence of propagule source ranged from 30.2 to 50.2%. Basal area of parental trees had the highest importance for P. serotina and the lowest for Q. rubra. This predictor also was more important for seedlings than saplings. The second group of predictors was connected with leaf litter predictors. The most important of them was litter mass, however its importance for R. pseudoacacia was lower than for the other species. Leaf litter predictors were the most important for P. serotina, especially for seedlings. Litter pH usually had lower importance than litter mass, with the exception of total R. pseudoacacia biomass. Fraction of decomposed material in litter was included only in the model of total P. serotina biomass. Light availability, expressed by DIFN, had the highest importance for R. pseudoacacia and the lowest for Q. rubra. Among predictors expressing tree stand features, number of tree species had higher importance than tree stand BA. Tree stand species richness was most important for R. pseudoacacia, and least for Q. rubra. Tree stand BA had the highest importance for R. pseudoacacia. Among understory predictors the highest importance was species richness, which had the highest importance for P. serotina and Q. rubra. Its importance was higher for saplings and total biomass than for seedlings. Understory species diversity, expressed by Shannon's index had the highest importance for R. pseudoacacia. CWMs describing functional traits of understory had the lowest importance for Q. rubra and the highest for P. serotina. For an average, the most important CWM was plant height. For Q. rubra the most important was CWM of SLA. Functional diversity indices was the group of predictors with the lowest importance. However, among species FRic, FDis and FEve were most important for P. serotina and FDiv for Q. rubra. In most cases values of importance of functional diversity components were higher for saplings than seedlings.
Impact of predictors on studied species biomass Partial dependence plots (Fig. 3) revealed that response of saplings to presence of propagule source in the vicinity was lower than response of seedlings. All species biomasses increased with proportion of  parent trees in basal area, up to thresholds at which response curves reached plateaus. This threshold was reached first by Q. rubra, then by R. pseudoacacia and finally by P. serotina. However, the response curve was steepest for R. pseudoacacia total biomass. Increasing litter mass was associated with increasing biomass of P. serotina and Q. rubra and decreasing biomass of R. pseudoacacia. Only seedlings of P. serotina showed a steeper response curve at the higher range of litter mass than saplings and total biomass. Responses to DIFN differed among species studied: biomasses of P. serotina and Q. rubra increased with increasing DIFN, however biomass of P. serotina seedlings in the highest fraction of DIFN range decreased. Biomass of Q. rubra decreased with increasing tree stand species richness and biomasses of P. serotina and R. pseudoacacia increased. The positive response of R. pseudoacacia was steeper in the highest range of tree stand species richness, especially for seedlings. This species showed a similar response to tree stand BA. The biomasses of all species studied increased with increasing understory species richness, however in the case of seedlings the size effect was minimal.

Model limitations
Despite the high number of predictors studied, some of them are only proxies for direct measurements. For example, litter traits are surrogates for soil chemistry, which is crucial for plant growth (e.g. Mueller et al. 2012;Aerts et al. 2017). We also did not divide regeneration into generative and vegetative specimens. Another potential drawback may result from low variability of climatic conditions, which drive the distributions of species studied at large spatial scales . Resource availability gradients do not account for water availability, as we did not sample dry and wet habitat types, due to their relatively low abundance, thus soil moisture effects were not studied. Nevertheless, the large number of study plots and temporal variation increased robustness of the conclusions.

Impact of propagule pressure on species invasion success
Almost all models showed that propagule pressure was the most important influence on invasion success, expressed by biomass of natural regeneration. Propagule source presence and propagule quantity both contribute to increasing propagule pressure, however, they each account for a different aspect of propagule availability. Different species reacted in a different way to presence and quantity: Q. rubra depended more on propagule vicinity than R. pseudoacacia and P. serotina was the least dependent among species studied. However, P. serotina had the highest importance of propagule quantity, similar to other studies (e.g. Chabrerie et al. 2008;Vanhellemont et al. 2009;Terwei et al. 2013). The important role of propagule pressure has been confirmed for numerous woody species (e.g. Pyšek et al. 2009;Sinclair and Arnott 2015;Rodríguez et al. 2017), and also for the species studied (e.g. Vanhellemont et al. 2009;Woziwoda et al. 2014;Jagodziński et al. 2015). Ecological success of P. serotina and Q. rubra, expressed either as seedling density (Pairon et al. 2006;Riepšas and Straigyté 2008;Jagodziński et al. 2015) or biomass (Jagodziński et al. 2018a) decreased with distance from propagule source. Although Q. rubra acorns are not preferred by birds (Myczko et al. 2014;Bieberich et al. 2016), its maximum dispersal distance referred to in the literature is higher than that of P. serotina (1500 vs. 600 m), which is eaten by numerous bird species (Bartkowiak 1970;Deckers et al. 2008). The high importance of fruiting tree BA may also reflect the fact that ca. 80% of fruits fall beneath the crowns of P. serotina and birds may disperse only 20% of them (Pairon et al. 2006). Lower importance of propagule source BA for Q. rubra may suggest better dispersal of this species, e.g. by rodents (Bieberich et al. 2016). Another explanation may result from with high seed mass and stored carbohydrates allowing growth in unsuitable conditions (Ziegenhagen and Kausch 1995). Under this assumption, even a single tree in the stand may produce an effective number of seedlings, and due to their low densities, may result b Fig. 3 Partial dependence plots showing impacts of particular variables (abbreviations and predictors importance in Table 1) on response functions of biomass for biomass types and species studied, when all remaining predictors were constant. Response and variables were centered, scaled and processed by Yeo-Johnson power transformation Drivers of invasive tree and shrub natural regeneration in temperate forests 2373 in biomass similar to the other species despite the low contribution of Q. rubra to tree stand BA. Propagule pressure increases the probability of alien species establishment in the plant community (Lonsdale 1999;Lockwood et al. 2005;Vanhellemont et al. 2009). For that reason plots with lower availability of seeds or sprouts are less invasible than those with higher availability. Due to probability of vegetative reproduction, effects of these variables may be biased, as the species studied are able to sprout (Closset-Kopp et al. 2007;Vítková et al. 2017). Although 10-50% of P. serotina natural regeneration may result from vegetative reproduction (Closset-Kopp et al. 2007), our study reports a fraction (8.7%) similar to our previous study in Rogów Arboretum . For R. pseudoacacia it was 35.9%, but due to young age of fruiting (Burns and Honkala 1990), there were no differences between presence of fruiting and sprouting trees in the vicinity of the plots. The high advantage of R. pseudoacacia in plots with presence of propagule source may also result from its capacity for nitrogen fixation, which fertilizes the soil (Rice et al. 2004;Cierjacks et al. 2013).

Impact of resource availability on studied species biomass
The most important factors related to availability of resources were those connected with leaf litter. These predictors were most important for P. serotina. However, the highest %IncMSE in this study was for litter mass as a predictor of Q. rubra saplings, and the highest effect size was for seedlings. The presence of Q. rubra in the overstory usually maintains a thick layer of leaf litter, due to its low decomposition rate (Horodecki and Jagodziński 2017). Thus, in this particular case its effect cannot be clearly separated. Nevertheless, both species showed similar tendencies for litter pH-their biomass decreased with increasing litter pH, which is an indicator of higher fertility. Both P. serotina and Q. rubra have been reported from fertile deciduous forests (Chmura 2004(Chmura , 2013Woziwoda et al. 2014;Dyderski et al. 2015;Jagodziński et al. 2015). P. serotina was claimed to prefer the poor and medium fertility sites of coniferous forests, and the species was planted there as a soil improver (Muys et al. 1992;Starfinger et al. 2003), but also reached high densities on habitats of fertile deciduous forests . Our study confirmed earlier observations that P. serotina reaches higher ecological success in less fertile and more acidic plots (Zerbe and Wirth 2006;Chabrerie et al. 2008;Knight et al. 2008;Halarewicz 2012;Terwei et al. 2016). In contrast, R. pseudoacacia reached the highest biomass in plots with low litter mass, low proportion of decomposed material and high pH. This indicated that although R. pseudoacacia may occur in less fertile sites, it has the highest growth potential in the most fertile habitats. For that reason this species effectively invades riparian forests (e.g. Terwei et al. 2013;Dyderski et al. 2015;Marozas et al. 2015).
One of the most surprising results was low importance of light availability for the majority of species studied, in contrast with other studies indicating higher importance of light availability on invasion success (e.g. Knight et al. 2008;Paquette et al. 2012;Rodríguez et al. 2017). The relatively low light requirements shown by R. pseudoacacia were opposite to other studies (Groninger et al. 1997;Cierjacks et al. 2013). However, due to relatively high proportion of seedlings in the biomass, which was highest in the stands with lowest light availability, it may be an effect of low competition within the understory. Nevertheless, seedling survival is low (Dyderski unpubl.). Another important factor is high proportion of sprouts, which may be supplied with nutrients by parental trees, and partially independent of light availability. Groninger et al. (1997) used only trees obtained from generative reproduction, thus our results are not comparable. However, in plots with high light availability R. pseudoacacia also reached high biomass. Dominance of propagule pressure over invasibility may be connected with their interaction, conceptualized by Davis et al. (2005), who elaborated a model showing how propagule pressure may modify invasibility, with the two variables working together to produce 'invasion pressure'.

Impact of interactions with overstory and understory on biomass
Interactions with other members of plant communities have been claimed to be drivers of invasibility from the beginning of invasion ecology, when Elton (1958) formulated the hypothesis of biotic resistance. This concept assumed that higher species richness in plant communities may decrease invasibility. This concept was recently confirmed for small plots (Brown and Peet 2003;Parker et al. 2010;Iannone et al. 2016), however at larger spatial scales recent studies confirmed that relationships between species richness of native and alien species or invader ecological success (expressed as abundance or density) are positive (Lonsdale 1999;Stohlgren et al. 1999Stohlgren et al. , 2006Knight et al. 2008;Dyderski et al. 2015, but see Parker et al. 2010). Increasing ecological success of the invasive species studied here in plots with higher species richness confirms this theory. Understory species richness was most important for P. serotina and Q. rubra saplings and total biomass, but not for seedlings. This may result from the effects of seed storage, allowing supplies of carbohydrates to seedlings which are independent of competition. The species with the lowest seed mass-R. pseudoacacia-did not show a similar response.
CWM of seed mass had the highest significance for P. serotina and R. pseudoacacia saplings. The biomass of P. serotina increased with increasing seed mass CWM. This may be driven by frequently cooccurring Quercus petraea and Q. rubra with high seed mass at a similar ecological scale, but also may reflect bird-mediated dispersal (Pairon et al. 2006;Deckers et al. 2008;Kurek et al. 2015;Dylewski et al. 2017). Terwei et al. (2016) interpreted a similar relationship obtained in floodplain forest as evidence for P. serotina shade tolerance, which has not been confirmed in our study. Biomass of R. pseudoacacia decreased with increasing CWM of seed mass, which may reflect its ruderal character and preference for disturbed sites (Cierjacks et al. 2013;Vítková et al. 2017), as low seed mass is connected with disturbance tolerance (Westoby 1998). P. serotina responded negatively to height CWM, which shows its negative response to understory competition. Similarly Terwei et al. (2013) found a negative relationship between density of P. serotina and herb layer cover.
Differences between studied species and management implications According to our results, management of species studied should primarily focus on removing propagule sources. If these species are planted according to their commercial value (Woziwoda et al. 2014) or ecosystem services (Vítková et al. 2017), invasion management should be focused on plant communities in the neighborhoods of plantations by shaping conditions not supporting the spread of invaders, e.g. by limiting light availability, planting competing native species or removing emerging saplings or juvenile specimens. Introduction of these best management practices for invasive woody species into forestry is crucial for limiting further invasion of these species (Brundu and Richardson 2016).