* 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.