* 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 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, Cameron, Parra, Jones, & Jarvis (2005). Before the analyses, the Pearson correlations between each pair of explanatory variables were checked. If the coefficient exceeded 0.7, the variable with the weaker ecological meaning was eliminated to avoid collinearity (Dormann et al., 2013). For details see Appendix 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 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 occurrence in the 1950s (distance, for details see Appendix, 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 modelled using a boosted regression trees (BRT) technique (De’Ath, 2007; De’Ath & Fabricius, 2000) employing packages gbm, dismo, and Biomod2 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 3900 forS . canadensis and 3850 for S . gigantea . Models for each species were constructed using all explanatory variables, and then simplified to obtain the parsimonious model. The BRT modelling and simplification of models were done based on Elith, Leathwick, & Hastie (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, Elith, Lahoz‐Monfort, & Guillera‐Arroita, 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 forS 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. The performance of the models was evaluated using area under the receiver-operating characteristic curve (AUC).The ecological interpretation of the model relied on the relative influences of explanatory variables and drawing response curves for each explanatory variable (Elith, Ferrier, Huettmann, & Leathwick, 2005). Eventually, maps of projected S. canadensis and S. giganteaprobability of occurrence were drawn (Figure 6). The probability of species presence in a given 2 × 2 km square was modelled for particular spatially blocked cross-validation runs and averaged.