Drivers of Solidago species invasion in Central Europe-Case study in the landscape of the Carpathian Mountains and their foreground.

AIM: The invasion process is a complex, context-dependent phenomenon; nevertheless, it can be described using the PAB framework. This framework encompasses the joint effect of propagule pressure (P), abiotic characteristics of the environment (A), and biotic characteristics of both the invader and recipient vegetation (B). We analyzed the effectiveness of proxies of PAB factors to explain the spatial pattern of Solidago canadensis and S. gigantea invasion using invasive species distribution models. LOCATION: Carpathian Mountains and their foreground, Central Europe.
METHODS: The data on species presence or absence were from an atlas of neophyte distribution based on a 2 × 2 km grid, covering approximately 31,200 km2 (7,752 grid cells). Proxies of PAB factors, along with data on historical distribution of invaders, were used as explanatory variables in Boosted Regression Trees models to explain the distribution of invasive Solidago. The areas with potentially lower sampling effort were excluded from analysis based on a target species approach.
RESULTS: Proxies of the PAB factors helped to explain the distribution of both S. canadensis and S. gigantea. Distributions of both species were limited climatically because a mountain climate is not conducive to their growth; however, the S. canadensis distribution pattern was correlated with proxies of human pressure, whereas S. gigantea distribution was connected with environmental characteristics. The varied responses of species with regard to distance from their historical distribution sites indicated differences in their invasion drivers. MAIN
CONCLUSIONS: Proxies of PAB are helpful in the choice of explanatory variables as well as the ecological interpretation of species distribution models. The results underline that human activity can cause variation in the invasion of ecologically similar species.

INTRODUCTION

Biodiversity and the function of ecosystems are threatened by global change drivers such as changes in land use and climate, as well as biological invasions (Linders et al., 2019; Sala et al., 2000). Invasive species alter a wide range of ecosystem services, including provisioning, regulation, and cultural and supporting functions, and they are particularly hazardous for biodiversity maintenance, human welfare, and the economy (Charles & Dukes, 2007; Chytrý et al., 2009; Hejda et al., 2009; Pejchar & Mooney, 2009; Vilà & Ibáñez, 2011). Globalization (e.g., international trade and travel) and climate change (e.g., global warming, droughts, and floods) can interact, which can in turn increase the level of biological invasions (Catford et al., 2009; Le Maitre et al., 2004; Pino et al., 2005; Seebens et al., 2015). As the total number of invasive species increases, some sites may host several alien species (Kuebbing & Nuñez, 2015).

The invasion process is a complex phenomenon, driven by numerous interacting processes, and the effect of this interaction is highly contingent on the context (Chamberlain et al., 2014; Frost et al., 2019). Consequently, drivers of plant invasion can vary depending on the specific region and habitat (Taylor et al., 2016). Nevertheless, invasions have a common pattern, which can be summarized as the joint effect of propagule pressure, abiotic characteristics of the environment, and biotic characteristics of both the invader and recipient vegetation (Catford et al., 2009), the so‐called PAB framework. Propagule pressure (P) includes dispersal and geographical constraints, while abiotic characteristics (A) comprise environmental and habitat constraints and biotic characteristics (B) describe the internal dynamics of the vegetation and community interactions (Catford et al., 2009). All these factors operate at different spatial scales (Czarniecka‐Wiera et al., 2020; Milbau et al., 2009) and are influenced by human activity (Essl et al., 2011). In practice, different indices can be applied as proxies of propagule pressure and abiotic and biotic conditions in modeling plant invasion process (Bazzichetto et al., 2018; Beaury et al., 2020; Chytrý et al., 2008; Szymura et al., 2018).

Related to the propagule pressure, the biological invasion correlates with many anthropogenic factors, such as density of the communication network, percentage of urban areas, gardening, and the fragmentation of natural habitats. Such factors can serve as a proxy of propagule pressure (Foxcroft et al., 2011; Pollnac et al., 2012; Štajerová et al., 2017; Szymura et al., 2018; Vilà & Ibáñez, 2011). In addition, economic and demographic variables reflect the intensity of human activities; therefore, socioeconomic factors such as gross domestic production and human population density can be important in predicting the invasion level (Essl et al., 2011; Hulme, 2017; Pino et al., 2005; Pyšek & Richardson, 2010) because they correlate with trade intensity and communication network density (Hulme, 2009). Among the abiotic interactions with the greatest impact on a large spatial scale (continental, regional), climate is considered the most critical in limiting the geographic distribution of species (Hulme, 2017; Thuiller et al., 2007). In terms of resource availability, invasive species usually prefer productive habitats where they are able to achieve competitive dominance (Czarniecka‐Wiera et al., 2020; Peltzer et al., 2016; Perkins et al., 2011). In addition, environments with high variability in resource availability, resulting from periodic external supply (e.g., surface runoff) or destruction of local vegetation that previously used the resources (e.g., human disturbances, abandonment of agricultural crops), are more susceptible to invasions than habitats with stable availability of resources (Davis et al., 2000; Kulmatiski et al., 2006; Rejmánek, 1989). Given the biotic characteristics of the invader and receipt communities, the limiting similarity hypothesis proposes that the invasion by alien species will be successful if the native species of the recipient community differ from the invader in terms of functional traits and resource requirements (MacArthur & Levins, 1967), which decreases competition for resources (Funk et al., 2008). Thus, the functional traits of the invader should not overlap with traits of native plants occurring in the invaded community, which will allow it to occupy an empty niche and successfully invade the community (Funk et al., 2008; Hejda & de Bello, 2013). Because some sites can be invaded by several species simultaneously, determining the interaction between invaders is critical for understanding their distribution (Kuebbing & Nuñez, 2015). For example, the local species assemblage can be driven by a priority effect, and the effect is particularly strong when interacting species have similar use of resources (Vannette & Fukami, 2014). In practice, the abundance and composition of invasive species are also related to landscape characteristics (e.g., habitat fragmentation, patch size, shape, and connections), habitat type, land use, and the composition of the surrounding landscape because these factors correlate with propagule pressure and habitat quality and availability (Basnou et al., 2015; Chytrý et al., 2009; González‐Moreno et al., 2013; Štajerová et al., 2017; Szymura et al., 2016).

Because of the complexity of biological invasion, better understanding of the underlying factors and their management is challenging. As tools for obtaining reliable and repeatable information for biological analyses as well as nature conservation and management of the invaders, invasive species distribution models (iSDMs) are considered useful (Lozano et al., 2020; Zurell et al., 2020). Modeling species’ environmental requirements and mapping their distributions through space and time help to identify the main introduction pathways and secondary spread and the areas and land use types that are more prone to invasion. These various threads could be woven into a strategy of prevention and elimination of invasive plant species on a regional scale (Lozano et al., 2020). Despite their deficiencies (e.g., problematic species‐environment equilibrium; Gallien et al., 2012; Hattab et al., 2017), iSDMs are still useful in the face of accelerating global changes and data deficiencies, as well as limited research funding (Yates et al., 2018). The PAB approach, despite its obvious advantages for selection of explanatory variables and model results interpretation, has rarely been used within an invasive species distribution modeling framework (but see Bazzichetto et al., 2018; Czarniecka‐Wiera et al., 2020; Lozano et al., 2020).

Goldenrod species from North America represent successful invaders in Europe, Asia, Australia, and New Zealand (Gusev, 2015; Szymura & Szymura, 2013; Ye et al., 2019; Zhang & Wan, 2017). In Central Europe, two invasive Solidago species occur, S. gigantea Aiton (giant goldenrod) and S. canadensis L. (Canadian goldenrod). Due to their high environmental impact, wide range of distribution, and locally high abundance, invasive Solidago species have to be controlled in Europe (Fenesi et al., 2015; Sheppard et al., 2006; Skórka et al., 2010). They have been proposed for addition to the list of hazardous alien species that threaten ecosystems, habitats, or other species in European Union countries (CABI, 2018; EPPO, 2020; Tokarska‐Guzik et al., 2015). Unfortunately, the eradication of widely established invasive plant species, such as Solidago, is not feasible. The management strategies need to integrate different options that account for the distribution and abundance of the invader, its environmental niche, and the areas that are likely to experience high impacts (Nagy et al., 2020; Shiferaw et al., 2019; Woodford et al., 2016). Management needs to consider intrinsic factors related to the biology and ecology of the invader, as well as extrinsic environmental factors, such as dispersal vectors and invasion pathways (Shiferaw et al., 2019).

Solidago canadensis and S. gigantea differ with regard to ecological niche in their native range (Johnson, 1995; Werner et al., 1980) and the time of introduction into Europe (Tokarska‐Guzik, 2005). However, previous studies suggest that these two species do not differ regarding their habitat preferences in Central Europe, and observed differences in their spatial distribution patterns emerge from historical contingency and limitation in long‐range dispersal (Szymura & Szymura, 2016). The two Solidago species occupy different areas and rarely form mixed‐species stands (Szymura & Szymura, 2016). In this study, we aimed to find the main drivers of Solidago species’ invasion at a regional scale, using a species distribution model and applying the PAB framework for selection of adequate explanatory variables and for ecological interpretation of the models. The distribution models can be used for mapping of invasion probability at a regional level to facilitate invasion control at a macroecological scale.

MATERIALS AND METHODS

Studied species

Goldenrod species are hemicryptophytes (shoots are annual and newly sprout each spring) with rhizomes; they are insect pollinated and self‐incompatible, with inflorescences forming at the top of each shoot which can produce up to 10,000–20,000 wind‐dispersed seeds per one ramet (Bielecka et al., 2017; Guzikowa & Maycock, 1986; Moran et al., 2017; Schmid et al., 1988). The seeds of S. canadensis and S. gigantea have a high germination percentage (Weber, 2000; Weber & Jakobs, 2005), but in dense, well‐established Solidago stands, seed germination and seedling emergence are exceptional. The clone size increases via horizontal rhizomes, and the death of an established genet is a rare event (Meyer & Schmid, 1999a, 1999b).

The native habitats of S. canadensis are tall‐grass prairies, infrequently grazed pastures, abandoned farmlands, roadsides, and waste areas in North America (Johnson, 1995; Werner et al., 1980). Solidago gigantea prefers moist habitats, such as woods, stream edges, and woodland borders (Johnson, 1995). In Europe, S. gigantea and S. canadensis occupy similar habitats and prefer fallow lands and ruderal habitats on moist to mesic sites, such as abandoned farmlands, scrub, roadsides, forest edges, grasslands, wetlands, and riversides (Szymura & Szymura, 2013, 2016). Invasive goldenrods are highly competitive for nutrients, water, and space, and they release allelopathic compounds that inhibit growth of other plants (Gusev, 2015; Ledger et al., 2015; Werner et al., 1980; Zhang & Wan, 2017). Due to prolific vegetative propagation, they form dense stands and decrease the biodiversity of plants (Chmura et al., 2016; Ye et al., 2019; Zhang & Wan, 2017); arthropods (de Groot et al., 2007), including pollinators (e.g., wild bees, hoverflies, and butterflies) (Lenda et al., 2020; Moroń et al., 2009, 2021) and ants (Kajzer‐Bonk et al., 2016; Lenda et al., 2013); and birds (Skórka et al., 2010).

Solidago canadensis was the first alien Solidago species recorded in Europe, in 1648, while S. gigantea was first recorded in 1758. The species were found in the territory of Poland about 100 years later, S. gigantea in 1853 and S. canadensis in 1872 (Tokarska‐Guzik, 2005). After S. canadensis and S. gigantea were introduced into botanical gardens, they were distributed among gardeners. The plants were attractive and easy to grow as ornamental plants, and they were useful for beekeepers (Guzikowa & Maycock, 1986; Roháčová & Drozd, 2009; Weber, 1997; Zihare & Blumberga, 2017). Recently, Solidago species have become widely distributed throughout Poland. According to the stages of invasion (Blackburn et al., 2011), S. canadensis and S. gigantea are now fully invasive species, with individuals dispersing, surviving, and reproducing at multiple sites in a wide variation of habitats over an extensive spatial area (E category).

Study area and species distribution data

The study area comprises approximately 31,200 km2 in the southeast part of Poland, which extends from latitude 50.2° to 49°N and longitude from 19° to 23°E (Figure 1). This area is diversified due to environmental conditions mostly shaped by the altitude ranging from 160 to 2,503 m a.s.l. Additional factors underlying diversity are correlated with climate, land use systems, land relief, and human population density. In the northern part, the lowland areas are used for agriculture and the foothills are dominated by forests, and the southern part has high mountains with alpine vegetation. In addition to the north–south altitudinal gradient, there is also a climatic gradient of continentality, with higher temperature range in the eastern part of the study region (Szabo‐Takacs et al., 2015) which, in the studied region, correlated strongly with decreasing eastward precipitation (Appendix S4, Table S3). The study area includes a densely populated industrial landscape (Silesia), urban agglomerations (largest city Kraków), and moderately populated agricultural areas, as well as sparsely populated areas in the mountains. The detailed characteristics of the study area (climate, topography, land use structure, and human population density) were previously described by Szymura et al. (2018).

FIGURE 1

The study region location (green) on a background of land relief (a), and distribution of communication network and settlements on the background of altitude within the study region (b)

The data on distribution of the studied Solidago species were obtained from the atlas Distribution of Kenophytes in the Polish Carpathians and their Foreland (Zając & Zając, 2015), which shows maps of species presence or absence in a 2 × 2 km grid in the Polish part of the Carpathian Mountains and their foreland, Central Europe. The fieldwork designed for the purpose of compiling the atlas was based on a survey of flora in particular regions (e.g., mountain ranges, particular towns, and surrounding areas) and exploration focused exclusively on neophytes in given regions. These observations were supplemented with additional data from phytosociological relevés, herbarium records, and published materials. The fieldwork was carried out by several dozen professional botanists as well as graduate students, focusing on a predefined 2 × 2 km grid for sampling (Zając A., personal information). This work represents a “survey” type of data, according to Elith et al. (2020) nomenclature. Such data, with true absence records, enable species distribution models to be less biased and to perform better, compared with presence‐only records, the “collection” data type (Barbet‐Massin et al., 2012; Elith et al., 2020). This distinction is of particular importance for examination of wide‐ranging and tolerant species (Brotons et al., 2004). To reduce the possible effect of lower sampling effort in some regions (Bailey et al., 2017; Yang et al., 2013), the potentially undersampled squares were excluded from modeling. For this purpose, we used a “target group approach” (Chapman et al., 2019; Phillips et al., 2009) and a previously established model which explains neophyte richness (the “target group” in this case) as a function of environmental and socioeconomic variables in the studied region (Szymura et al., 2018). We assumed that the squares with the highest negative model residuals (i.e., squares where recorded neophyte richness was much lower than predicted by the model) indicated potentially undersampled regions. After preliminary testing, we decided to exclude from modeling 25% of squares (1950 squares) with the highest negative residual values and simultaneously without any invasive Solidago records (for details of this calculation see Appendix S1).

Explanatory variables and statistical analysis

We prepared a data set of environmental variables that can be considered as proxies of propagule pressure, abiotic environment, and biotic characteristics, based on the PAB framework (Catford et al., 2009; Table 1). These proxies were chosen based on the results of previous study on Solidago (Szymura et al., 2016) and the most influential drivers of neophytes in the region (Szymura et al., 2018).

TABLE 1

Explanatory variables selected for modeling invasive Solidago distribution. Variables in bold type were used in the final model, and the remaining variables were excluded from further analysis due to collinearity

Explanatory variable	Abbreviation	Probable sphere of PAB framework
Communication routes (railways and roads) density	communication	P
Shannon's diversity index of landscape	SHDI	B
Urban area percentage	urban	P
Cropland area percentage	cropland	B
Forest area percentage	forest	B
Human population density	density	P
Income per capita	income	P
Topographic roughness index	TRI	A
Topographic position index	TPI	A
Average annual temperature	temperature	A
Topographic wetness index	TWI	A
Temperature seasonality	Ts	A
Annual sum of precipitation	precipitation	A
CaCO3 content	Ca	A
K content	K	A
N content	N	A
P content	P	A
pH in H2O	pH	A
Distance to nearest introduction site S. canadensis	distance_S.can	P
Distance to nearest introduction site S. gigantea	distance_S.gig	P
Presence of competing Solidago species a	competitor	B

Presence of one invasive Solidago species in the same 2 × 2 km square was considered as an explanatory variable for the other; that is, in model for S. canadensis, its presence explained the presence of S. gigantea and vice versa.

The anthropogenic variables were derived from CORINE 2012 database (urban), the Central Statistical Office of Poland (income), and Statistics Poland (density). The income (as an estimator of wealth) is directly correlated with trade intensity and thus reflects the potential to alien species propagule transportation by trade or accidentally (Hulme, 2009; Pyšek et al., 2010). The length of communication routes (communication) was obtained from the Polish Geographical Objects Database (BDOO). The other data were calculated from the CORINE 2012 database (cropland, forest, SHDI). A Digital Elevation Model for Europe (EU‐DEM) was used to calculate the topographic metrics (TPI and TWI). Maps prepared by Ballabio et al. (2019) using data from Land Use and Cover Area frame Survey (LUCAS) were used to calculate soil characteristics (content of N, P, K, and soil pH). The climate data (precipitation, temperature) were derived from a climatic model developed by Hijmans et al. (2005). Before the analyses, the Pearson correlations between each pair of explanatory variables were checked. If the coefficient exceeded 0.7, one of the correlated variables was eliminated to avoid collinearity (Dormann et al., 2013). For details, see Appendix S3 Table S2. The average values of the variables were calculated for each 2 × 2 km grid cell acquired from Zając and Zając (2015), and the landscape diversity (SHDI) was expressed by Shannon's diversity index.

Maps showing the historical distribution of goldenrods before their spreading phase (Tokarska‐Guzik, 2005) were used to calculate the distances from a focal 2 × 2 km square to the nearest site of goldenrod introduction in the 1950s (distance, for details see Appendix S3, Map S2). To check whether the presence of one Solidago species in a 2 × 2 km square explained the presence of the second species (possible priority effect), the data on distribution of the potential competitor were used as an explanatory variable (competitor). All the calculations and map handlings were done using QGIS, SAGA GIS, and FRAGSTAT software.

Goldenrod species spatial pattern of distribution was modeled using a boosted regression trees (BRT) technique (De’Ath, 2007; De’Ath & Fabricius, 2000) employing packages gbm (Greenwell et al., 2020), dismo (Hijmans et al., 2020), and Biomod2 (Thuiller et al., 2020) in the R environment. After initial examinations, the BRT settings were applied: tree complexity, 5; bag fraction, 0.5; learning rate, 0.001; and cross‐validation, 10 fold. The optimal number of trees was 3,900 for S. canadensis and 3,850 for S. gigantea. Models for each species were constructed using all explanatory variables and then simplified to obtain the parsimonious model. The BRT modeling and simplification of models were done based on Elith et al. (2008) suggestions. Then, the modeling, using the tuned model parameters and a minimal set of explanatory variables, was performed in Biomod2 package with spatially blocked cross‐validation (Valavi et al., 2019). We applied 5‐fold cross‐validation, using spatial blocks constructed based on 10 × 10 km squares for S. canadensis and 20 × 20 km squares for S. gigantea. The sizes of the squares were chosen based on spatial autocorrelation of raw distribution data (Roberts et al., 2017), and the blocks were constructed using BlockCV package within the R environment (Valavi et al., 2019). For details of this approach, see the Appendix S1. The performance of the models was evaluated using area under the receiver operating characteristic curve (AUC). Following qualitative descriptions, an AUC value in the range of 0.7–0.8 can be considered a good prediction, 0.8–0.9 as a very good prediction, and above 0.9 as an excellent prediction (Šimundić, 2009). The ecological interpretation of the model relies on the drawing response curves for each explanatory variable (Elith et al., 2005) and calculation of relative influences of explanatory variables. The relative influence is calculated after model training (calibration) and prediction making. Then, values of one of the variables are randomized, and a new prediction is made. The correlation between this new prediction and the standard prediction is calculated and is considered as an estimation of the variable importance in the model: A high correlation score between the two predictions shows that the randomized variable has little influence on making a prediction and is not considered important for the model in its prediction. In contrast, a low correlation means a significant difference in the prediction making, showing that the variable is important for the model. The variable importance is calculated as 1‐correlation, repeated for each variable independently and for each spatially blocked cross‐validation run (Thuiller et al., 2012). Eventually, maps of projected S. canadensis and S. gigantea probability of occurrence were drawn (Figure 5). The probability of species presence in a given 2 × 2 km square was modeled for particular spatially blocked cross‐validation runs and averaged, employing the “projection” function in the Biomod2 package. Additionally, maps of squares projected to be suitable for invaders and not colonized yet were produced for conservation purposes (Appendix S7, Map S4).

FIGURE 5

The projected probability of presence of the invasive Solidago species. The optimal cutoff value was 0.205 for S. canadensis and 0.539 for S. gigantea

RESULTS

Goldenrod species were observed in 60.5% of the squares (in 3,544 out of 5,850 finally examined squares). Solidago gigantea was the most frequent species (53.1%, 3,107 squares) followed by S. canadensis (21.4%, 1,255 squares).

Solidago gigantea localities were widespread throughout almost the entire area, aside from the higher altitudes in the southern part of the study region. The S. canadensis was concentrated in the western part of the study area, while being sporadically dispersed in the eastern part and also avoiding the southern fragment with higher altitudes (Figure 2).

FIGURE 2

Distribution of invasive Solidago species (orange color) in studied region. The light gray color show distribution of squares with confirmed Solidago absence. Squares excluded from analysis, are not shown (left blank)

The average value of AUC was 0.836 for S. candensis and 0.786 for S. gigantea. Despite some differences in model evaluations of particular spatially blocked folds, the models for S. canadensis generally performed better than those for S. gigantea (Figure 3a). The parsimonious (simplified) model for S. canadensis relied on a higher number of explanatory variables than those for S. gigantea.

FIGURE 3

The values of area under curve (AUC) for simplified models of S. canadensis and S. gigantea distributions (a), and variable importance for each variable involved in the simplified models of (b) S. canadensis and (c) S. gigantea. The boxplot shows the results of runs in the spatially blocked, 5‐fold cross‐validation. The values of variable importance close to 1 indicate high variable importance to the model, while those close to 0 have low importance. Abbreviations: cropland, cropland area percentage; density, human population density; temperature, average annual temperature; Ts, temperature seasonality; distance_S.can, distance to nearest introduction site S. canadensis; distance_S.gig, distance to nearest introduction site S. gigantea

Both species reacted to climatic conditions, expressed by the annual average temperature (temperature) and temperature seasonality (Ts), as well as the distance from the initial introduction sites (distance) (Figure 3b,c). Moreover, the spatial pattern of distribution of S. canadensis was also explained by anthropogenic factors, such as human population density (density) as well as the percentage of agricultural lands (cropland). The full list of all variables included in the final models, along with their relative influence, is shown in Figure 3b,c.

The modeled response of species on particular variables is shown in Figure 4. The distribution of both species was climatically limited, with the species being unlikely to occur in regions with an average annual temperature below 5.5°C. The probability of S. canadensis occurrence increased with human population density (density) (Figure 4), as well as distance from its introduction site (distance_S.can), with squares placed 100 km distant from the initial sites of introduction having the highest probability. The distribution of S. gigantea was also correlated with the pattern of its initial introduction (distance S. gig), and the probability of its occurrence generally decreased with the distance (Figure 4), reaching the lowest value at about 40 km and fluctuating above it.

FIGURE 4

The modeled responses of Solidago species for particular environmental variables. The shape of the response was modeled using the evaluation strips method (Elith et al., 2005), with spatially blocked, 5‐fold cross‐validation. The graphs are sorted according to decreasing value of variables’ importance, upper panel for S. canadensis, lower for S. gigantea. Abbreviations: cropland, cropland area percentage; density, human population density; temperature, average annual temperature; Ts, temperature seasonality; distance_S.can, distance to nearest introduction site S. canadensis; distance_S.gig, distance to nearest introduction site S. gigantea

The results of the projections are presented in Figure 5. The average cutoff values, calculated based on the AUC values, were 0.205 for S. canadensis and 0.539 for S. gigantea. In a comparison of the observed distribution with the model's prediction, the number of squares suitable for the invaders and not colonized yet (including the undersampled squares excluded from the analyses) increased by 45% (1,255 squares with presence versus. 2,293 predicted) for S. canadensis and 36% (3,107 squares with presence vs. 4,897 predicted) in the case of S. gigantea. For detailed maps see Appendix S7, Map S4.

DISCUSSION

The model's performance in interpreting the AUC values (Šimundić, 2009) should be considered as good for S. gigantea and very good for S. canadensis, despite the relatively limited number of explanatory variables retained after the model's simplification. Moreover, in the case of species with broad environmental tolerance, such as the studied Solidago, the model's performance is usually lower than it is in comparing with specialist species, both plants and animals (Guisan et al., 2007; Regos et al., 2019). The model's performance is improved by variables that can be interpreted as proxies of P, A, and B factors (see discussion below); however, the importance of the variables differed considerably between particular P, A, and B factors, as well as species studied.

Ecological interpretation of the models

Propagule pressure

The recent distributions of examined species were correlated with initial patterns of their introductions in the 1950s. Quite surprisingly, the two species revealed an opposite relationship to these historical patterns. In the case of S. gigantea, the pattern was rather simple and intuitive: The probability was highest in squares closer to the sites of initial distribution. However, S. canadensis quite surprisingly was the most likely to occur in squares 100 km from the initial sites of introduction. These results suggest different mechanisms of long‐range dispersals, not related to biological issues, since their seeds, dispersal mechanism, and flowering time are similar (Weber, 2000; Weber & Jakobs, 2005).

Recently, S. canadensis was considered to have a higher ornamental value (because of larger size, bigger inflorescences, and clump occurrence) than S. gigantea. As a result, it is offered by garden shops, but S. gigantea is not (Szymura M. personal observations, data from internet shops offering ornamental plants). A similar pattern of trade has been described in Estonia, Central Europe, where only S. canadensis is offered in markets (Ööpik et al., 2013). Moreover, the honey from S. canadensis has recently been promoted on social media, without supporting scientific data, as a “superfood” with healing properties. This claim could encourage beekeepers to produce goldenrod honey, which would lead to further spread of S. canadensis and exacerbate its existing negative environmental impact (Lenda et al., 2020). Thus, it could be assumed that the long‐range dispersal of S. canadensis is recently enhanced by humans.

The distribution of S. canadensis is positively correlated with human population density. This straightforward correlation breaks if the population density exceeds 5,000 ind km−2. This happened in a few of the most densely inhabited squares, representing strict city centers. It was generally found that the plant species richness in areas with moderate levels of urbanization (e.g., suburban areas) exceeded the richness recorded in nonurbanized areas as well as in central, urban core areas (McKinney, 2008). The lack of a further increase in alien species richness in strict city centers, despite the high propagule pressure, was explained by the loss of suitable areas for plants (McKinney, 2008). Such generally limited neophytes’ richness caused by population density has previously been shown for this region (Szymura et al., 2018).

The results of the modeling support the assumption that recent S. gigantea dispersal has occurred mostly spontaneously without any human aid, while S. canadensis dispersal is still related to human presence and, additionally, intentional transport over longer distances via, for example, Internet commerce (Lenda et al., 2014). This pattern is partially related to longer invasion history of S. gigantea comparing to S. canadensis, with caused that the studied region is exposed to S. gigantea seed for a longer time.

Abiotic factors

The variables representing abiotic environment (A) are the most important for model performance for both species; however, the impact of these variables was more pronounced in the case of S. gigantea, compared with S. canadensis.

The distribution of both species was restricted climatically, and their presence was unlikely in areas with an average yearly temperature below approximately 5.5°C. The temperature corresponds with the altitudinal zonation of vegetation in the studied region and relates to a lower limit of the montane zone, starting from an altitude of approximately 600–850 m a.s.l. in the studied region. The negative effect of cold climate on the distribution of both Solidago species studied is in accordance with studies examining their potential distribution in Europe, which indicated that northern Europe as a region is outside their climatic requirements (Weber, 2001). Although both species can be observed sporadically at higher altitudes, their typical upper limit is 1,200 m a.s.l. (Moran et al., 2017; Weber & Jakobs, 2005). In the case of S. gigantea, positive correlations have been found between the mean temperature and growth parameters, and high spring temperatures (above 24°C) are advantageous for germination (for review, see Weber & Jakobs, 2005). Solidago canadensis plants are taller at lower attitudes, and at higher altitudes, they are not able to develop seeds because of the limited length of the vegetation period (Moran et al., 2017). It should be noted that the data referred to here regarding altitude come from the central Alps, while the climate in the Carpathian Mountains is more severe; therefore, the upper limits of the vegetation zones are at lower altitudes in the Carpathian Mountains compared with the Alps (Ellenberg, 1988; Pawłowski, 1972).

The species distributions were also correlated with temperature seasonality, which in the studied region is also related to the precipitation pattern (Appendix S4, Table S3). Solidago canadensis is more abundant in the western part of the study region, which has lower temperature seasonality and higher precipitation, while S. gigantea avoids the southern part of the region with higher precipitation and also lower temperature seasonality. Previous studies examining the potential range of this species in Europe (Weber, 2001) suggested that these aspects (continentality gradient and precipitation) did not restrict their distribution in this part of Europe. Therefore, the extent to which the observed relation is causal is not clear, and the possibility exists that it reflects a peculiarity of the distribution in the studied region.

The models did not indicate that soil properties and land relief features are among the crucial factors explaining the distributions of the invaders. Both species are known to have rather broad tolerance to soils (Szymura & Szymura, 2016; Weber & Jakobs, 2005; Werner et al., 1980), which could explain why soil properties were not relevant in studied region. Observations from early phase of invasion on studied region, up to 1989s, underlined the role of river valleys, as a main route of invasion (Tokarska‐Guzik, 2005). The results obtained here show that the species are broadly widespread and their invasion is no longer related to watercourses.

Biotic factors

Because of the character of the data (observation for 2 × 2 km grid), we had no detailed information regarding invaded habitats. However, the data still allowed testing the hypothesis regarding species co‐occurrence at landscape scale and the effect of dominant land cover/land use forms. Results from other region of Central Europe revealed existence of large areas dominated by a single invasive Solidago species, where the presence of another was unlikely. This spatial pattern results, most likely, from priority effect (Szymura & Szymura, 2016). In the studied region, we had no evidences for such phenomenon: The presence of one species did not explain the absence of the other. The species rarely formed mixed stands (Szymura & Szymura, 2016), but considering grain size used in this examination (square 2 × 2 km) it can be assumed that they could co‐occur in the same landscape. We also found that the presence of S. canadensis is rather unlikely in a landscape dominated by agricultural areas. It could be linked to high use of herbicides and a small amount of available area for invasive goldenrod habitats (e.g., abandoned fields, meadow, pastures) in lands with intense, large‐scale agriculture (Szymura & Szymura, 2016; Szymura, Szymura, & Wolski, 2016).

The relatively low importance of variables that can be related to biotic interactions does not necessarily mean that biotic interactions did not shape invasion pattern. It is more likely related to the grid size in this study (2 × 2 km), while the biotic interactions occur mostly in the closest vicinity of the studied individuals. Such data can potentially be derived from other sources of information, namely phytosociological relevés, which document species composition and abundance in small plots (~25 m2 for herbal vegetation).

Conservation implications

The two species differed regarding prominent constraints: both were limited climatically, avoiding cold, mountain climate, but S. canadensis with a still limited range was also related to proxies of human pressure. Based on the results, it can be hypothesized that recent dispersal of S. gigantea in the studied region has happened mostly spontaneously, while the spread of S. canadensis could be related to trade and intentional introductions. Given the wide range of distribution of both species, their successful eradication in the region seems unlikely. The eradication of Solidago is not easy and must include the establishment of native vegetation to prevent reinvasion of the Solidago on the site. It needs a long time and financial effort (Szymura et al., 2019; Szymura, Szymura, & Wolski, 2016). However, local eradication in mountains, above 600–850 m a.s.l. where the species occur infrequently may still be feasible and could be considered as a management option. In the case of S. canadensis, proscription of its sale could restrict its further spread. Assuming the successful restriction of the trade, eradication in the eastern and central parts of the region, where the species is still uncommon, will be achievable. Similarly, the control of invasive plant species populations in human settlements and their surrounding area seems to be a reasonable method. In contrast, the management of S. gigantea should focus on areas with a high value for nature conservation that are close to already existing populations of this species. Among management of invasive Solidago stands, the mowing, grazing, flooding, and combination of these methods are considered (Nagy et al., 2020, 2021). Herbicide use should be banned because of its environmental impact, including the effect on native vegetation (Schulz et al., 2021; Weidlich et al., 2020); moreover, its long‐term effect is not better than mechanic methods (Szymura et al., 2019). Nonetheless, the model's prediction suggests a considerable increase of invaded areas by both species. The location of suitable squares that are not yet colonized suggests that the expansion of the invaders will take place by range filing rather than increasing the range (Appendix S7, Map S4). In Central Europe especially prone to invasion are abandoned agricultural lands (Bartha et al., 2014; Fenesi et al., 2015), therefore a policy preventing agricultural land abandonment is desirable to counteract the further increase of goldenrods invasion level. The model outputs seem to be transferable into other areas with similar climate, land use history, economy, and invasion history, including the Carpathian Mountains and the surrounding regions in Slovakia, Ukraine, Hungary, and Romania. However, we did not have enough data to directly test the possible model application, using the limited number of explanatory variables, to maximize the iSDM transferability (Petitpierre et al., 2017). In this context, the procedure of model simplification (Elith et al., 2008), which reduced the number of explanatory variables, seems to be a great advantage of the BRT modeling technique.

Model limitations, and methodological problems

Recently, numerous iSDMs have been based on presence‐only data and employ so‐called background points (pseudo‐absences). Nonetheless, data not only on species presence but also their true (i.e., confirmed) absence are considered more relevant for modeling (Barbet‐Massin et al., 2012; Brotons et al., 2004; Elith et al., 2020). Unfortunately, confirmed absence data are problematic because they need a high sampling effort (Barbet‐Massin et al., 2012; MacKenzie & Royle, 2005) to be realistic. Our results show that exclusion of squares with low sampling effort improves the model's performance. This suggests an issue of sampling bias, which can be ameliorated by appropriate procedures. Our approach seems to be promising, but it needs further study in order to better understand its operation. The typical assumption, such as higher sampling effort in densely populated areas and near roads, is not adequate for invasive species because they typically occur in urban areas and along communication routes (Niinemets & Peñuelas, 2008; Szymura, Szymura, & Wolski, 2016).

Another problem consists of causality in our model: The approach applied represents a correlative type of model that is unable to directly capture the underlying processes driving the observed patterns of distribution. Contrary to the correlative approach, the mechanistic (or process‐based) models, which are built using explicit descriptions of biological mechanisms, are free from this disadvantage (Yates et al., 2018). In result, mechanistic and hybrid models have recently been recommended for modeling species distribution (Zurell et al., 2016). Studies in simulated systems reveal great potential of mechanistic models as BioGEEM in examination of eco‐ecological questions (Cabral et al., 2019). However, they still meet a numerous obstacles in practical implementation, especially in macroecological and biogeographical applications (Cabral et al., 2017), since simulation of large, species‐rich ecosystems is challenging (Cabral et al., 2019). They need appropriate formulation including detailed data on species response to environment, preferably coming from experiments, which are typically unavailable (Yates et al., 2018; Zurell et al., 2016). In practice, the models rely to a considerable degree on parametrization based on observational data, and as a result, the difference between correlative and mechanistic models is often fuzzy (Yates et al., 2018). Similarly, there is a problem with incorporating the effect of long‐range seed dispersal by wind. It needs detailed data regarding wind direction and velocity, seed dispersal kernel, and local population demography (Neubert & Caswell, 2000). As result, models are developed for restricted areas (e.g., Baker, 2017; Williams et al., 2008), which caused particularly useful for modeling dispersal of newly established populations (Gallien et al., 2010). To conclude, the recent state of knowledge regarding processes driving Solidago invasion restricts application of mechanistic or hybrid models of their invasion to a regional spatial extent.

CONCLUSIONS

The PAB framework enhanced the iSDM by helping in the selection of explanatory variables, as well as the ecological interpretation of the models. Nonetheless, in practice it needs high‐quality data that are typically unavailable to fulfill this approach, especially regarding biotic interactions. In case of plant invasion, adequate data on the biotic component could be delivered by phytosociological relevés. The employment of maps showing the historical distribution of invasive species enhanced the modeling by revealing differences in patterns of species spread into a region. In result, the model reveals that two alien species with similar ecology and biology can vary considerably in their invasion pattern due to direct human interference. Therefore, the conservation options, derived from iSDM, should be focused on a particular species, not groups of species, even if they have similar ecology and are closely related taxonomically.

The presence/absence data, in addition to their pre‐eminence compared with opportunistic, presence‐only data for species distribution modeling purposes, are still prone to some bias. Results of this study suggest that the bias is correlated with mistakenly reported species absence. Exclusion of the potentially undersampled plots increased the model performance; however, additional data are needed (e.g., richness of target species group).

CONFLICT OF INTEREST

The authors declare no competing interest.

AUTHOR CONTRIBUTIONS

Peliyagodage Chathura Dineth Perera: Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (equal); Visualization (supporting); Writing‐original draft (lead); Writing‐review & editing (lead). Tomasz H. Szymura: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (equal); Software (equal); Supervision (equal); Visualization (equal); Writing‐original draft (equal); Writing‐review & editing (equal). Adam Zając: Data curation (equal); Resources (lead). Dominika Chmolowska: Software (equal); Visualization (equal); Writing‐review & editing (equal). Magdalena Szymura: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Supervision (lead); Visualization (equal); Writing‐original draft (equal); Writing‐review & editing (equal).

Supplementary Material

Click here for additional data file.