RETRACTED ARTICLE: Major environmental factors and traits of invasive alien plants determining their spatial distribution

Abstract

Background

As trade increases, the influx of various alien species and their spread to new regions are prevalent and no longer a special problem. Anthropogenic activities and climate changes have made the distribution of alien species out of their native range common. As a result, alien species can be easily found anywhere, and they have nothing but only a few differences in intensity. The prevalent distribution of alien species adversely affects the ecosystem, and a strategic management plan must be established to control them effectively. To this end, hot spots and cold spots were analyzed according to the degree of distribution of invasive alien plants, and major environmental factors related to hot spots were found. We analyzed the 10,287 distribution points of 126 species of alien plants collected through the national survey of alien species by the hierarchical model of species communities (HMSC) framework.

Results

The explanatory and fourfold cross-validation predictive power of the model were 0.91 and 0.75 as AUC values, respectively. The hot spots of invasive plants were found in the Seoul metropolitan area, Daegu metropolitan city, Chungcheongbuk-do Province, southwest shore, and Jeju island. Generally, the hot spots were found where the higher maximum temperature of summer, precipitation of winter, and road density are observed, but temperature seasonality, annual temperature range, precipitation of the summer, and distance to river and sea were negatively related to the hot spots. According to the model, the functional traits accounted for 55% of the variance explained by the environmental factors. The species with higher specific leaf areas were more found where temperature seasonality was low. Taller species preferred the bigger annual temperature range. The heavier seed mass was only preferred when the max temperature of summer exceeded 29 °C.

Conclusions

In this study, hot spots were places where 2.1 times more alien plants were distributed on average than non-hot spots (33.5 vs 15.7 species). The hot spots of invasive plants were expected to appear in less stressful climate conditions, such as low fluctuation of temperature and precipitation. Also, the disturbance by anthropogenic factors or water flow had positive influences on the hot spots. These results were consistent with the previous reports about the ruderal or competitive strategies of invasive plants instead of the stress-tolerant strategy. The functional traits are closely related to the ecological strategies of plants by shaping the response of species to various environmental filters, and our result confirmed this. Therefore, in order to effectively control alien plants, it is judged that the occurrence of disturbed sites in which alien plants can grow in large quantities is minimized, and the river management of waterfronts is required.

Background

As the economic costs to control and manage biological invasions have been getting severe in recent decades (Diagne et al. 2021), the government in many countries request more efficient and practical management plans for invasive species. For this reason, the demand for systematic analysis of the nationwide distribution of alien plant species kept increasing. In this regard of view, annually collected nationwide survey data for invasive species in South Korea provide an excellent opportunity to meet those demands.

Invasive species distribution has been usually modeled by single species distribution modeling, but it is less suitable to analyze the community data comprised of many different invasive species. Joint species distribution modeling confers a great opportunity to analyze this kind of community data (Warton et al. 2015; Abrego et al. 2017), and it can be used to find the location and the distribution of the hot spots and cold spots of invasive alien plants. Finding the hot spots of invasive alien plants is essential to building a risk map for more efficient management. Also, these local clusters can be related to the combination of relevant abiotic factors, such as climate variables, topographic factors, and anthropogenic factors, and these relationships would deepen our understanding of the pattern of invasive species.

Moreover, understanding the role of functional traits of invasive species when they pass through environmental filters or face disturbances is one of the critical questions in invasion science (Mouillot et al. 2013; Cadotte et al. 2015; Pearson et al. 2018). In community data, species traits should be treated as essential predictors to answer why some taxa are more abundant than others in the common environment. Therefore, including the interaction between the environment and traits could improve the power of the joint species distribution model.

In this paper, we aimed to systematically analyze the nationwide distribution of invasive alien plants (IAPs) to find their hotpots and the environmental factors associated with those places. Also, we included functional traits of the species to find how the traits affect the responses to the environment. We used the spatial joint species distribution model with climatic variables, topographic variables, and disturbance-related variables, and also, we included functional traits and phylogenetic relatedness among the species as predictors.

Materials and methods

Survey

The nationwide survey for the alien plant species was conducted by the National Institute of Ecology, South Korea from 2015 to 2019. For this period, 20 scientists conducted convenient sampling in most provinces (165 districts) as much as possible. At all the survey points, they placed temporary plots with variable sizes and recorded all the alien plants found within the plots. They conducted a pilot survey in Jeju island in 2015 and divided the mainland of Korea into three regions, and they surveyed one region per year.

We had 10,287 sample points in total after data cleaning. As for the environmental factors, we collected 19 bioclimatic variables at the resolution of arc 30 seconds from WorldClim version 2 (Fick and Hijmans 2017), ASTER GDEM version 3 at the resolution of 1 arcsecond (NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team 2019), distance to the river, distance to the sea, and road density within a 1-km circle. All the environmental variables were resampled to the resolution of 10 × 10 km².

All the sample points were aggregated to the 10 × 10 km² grid by pooling the recorded species to remove duplicates within a pixel of an environmental variable. The aggregated sample points had their coordinates at the centroid of the original sample points. The values of environmental variables were extracted by the aggregated sample points and multicollinearity was removed by a stepwise procedure using the threshold of VIF ≤ 5 by usdm R package version 1.1.18 (Naimi et al. 2014). The selected variables were isothermality (bio 03), temperature seasonality (bio 04), max temperature of the warmest month (bio 05), precipitation of the driest month (bio 14), precipitation of the warmest quarter (bio 18), distance to the river, distance to the sea, and road density. We only focused on species with a frequency of greater than 1% to remove rare species to get robust results.

Functional traits

According to the Leaf-Height-Seed model proposed by Westoby (1998), functional traits related to leaf characteristics, plant height, and seed are essential factors to define plant ecology strategy schemes. These are related to the efficiency of resource capture and utilization or competitive ability in various habitats (Díaz et al. 2016). We selected specific leaf area (SLA), leaf dry matter content (LDMC), plant height, and seed dry mass. These characteristics were relatively easily and largely obtainable variables from the functional trait database. Although morpho-anatomical (soft) traits, such as the variables we used, have lower predictive power than physiological (hard) traits, their combination data can provide a better explanation on species response along environmental gradients (Belluau and Shipley 2018).

We downloaded the data for seed dry mass, SLA, and LDMC traits from the TRY database (Kattge et al. 2020). To summarize trait values among multiple measurements from multiple references, we first averaged trait values per reference and species and calculated median values per species. We imputed missing values by the median values of its congeners to reduce data bias if possible. Finally, we had 126 species left after listwise deletion for the missing trait values.

Statistical modeling

We adopted the hierarchical model of species communities (HMSC) framework as a spatial joint species distribution model to explain what environmental factors are related to the occurrence of IAPs and their joint species richness (Fig. 1) (Ovaskainen et al. 2017; Tikhonov et al. 2020). We included the functional traits per species to analyze the contribution of the functional traits to the level of species response to the environmental factors. Focal IAPs were phylogenetically related to each other; therefore, we incorporated phylogenies to the model to account for the non-independence of traits among taxa. A phylogenetic tree was generated through the V.PhyloMaker R package version 0.1.0, in which a mega-tree of 74,533 vascular plant species is implemented (Jin and Qian 2019). The data were too spatially extensive to compute Bayesian JSDM; therefore, we used nearest neighbor Gaussian process (NNGP) approaches in the latent factor structure of HMSC (Tikhonov et al. 2019).

Fig. 1

The schematic summary of the statistical modeling. The responses of IAPs to the environmental variables were modeled using their occurrence data, traits data, phylogenetic tree, and environmental variables. The predicted distribution map of IAPs can be made by projecting the model onto environmental variable layers

The HMSC framework uses these matrices to model in spatial context: spatial coordinates of sampling unit, species occurrence at each sampling unit, environmental variables at each sampling unit, a phylogenetic covariance matrix of focal species, and trait values of each species.

The sum of the predicted probability of all the species represents the predicted species richness. For spatial cluster analysis, local Moran’s I of the predicted species richness was calculated for each pixel using queen contiguity-based weights (Anselin 1995) using spdep R package version 1.1.5 (Bivand and Wong 2018). All the analyses were performed in R 4.0.2 (R Core Team 2020).

Results

The explanatory power of the model was evaluated by AUC and Tjur’s R squared (Tjur 2009). AUC was evaluated for each species (Table S1), and the average AUC was 0.910 ± 0.058 (mean ± SD). Tjur’s R squared was 0.27 ± 0.16 (mean ± SD). Meanwhile, the predictive power of the model was calculated by fourfold cross-validation, and the average AUC was 0.753 ± 0.101 (mean ± SD). Tjur’s R squared was 0.12 ± 0.11. There are no strict guidelines of an acceptable range of AUC, but greater than 0.7 is generally acceptable according to the rule of thumb (Hosmer Jr et al. 2013). Tjur’s R squared can be similarly interpreted as the R squared of linear regression, but it is generally much lower than R squared by the nature of presence-absence data.

The richness of IAPs was remarkably high in the Seoul metropolitan area (suburban), Chungcheongbuk-do Province, southwest shore, Daegu Metropolitan City, and Jeju island (Figs. 2 and 3). These hot spot areas for IAPs have 33.5 IAPs on average, which was 2.1 times greater than non-hot spot areas (15.7), and they were explained by the combination of environmental factors (Table 1). Among the variables, the climate-related variables and the spatial random variable were the major variables.

Fig. 2

The predicted richness of IAPs in all the districts in South Korea. The sum of all the predicted probabilities of occurrences of IAPs indicates the predicted richness of IAPs

Fig. 3

Spatial clusters and outliers of the predicted richness of IAPs. High-high clusters represent hot spot areas, and the low-low clusters represent cold spot areas. The others indicate spatial outliers

Table 1 Importance of the environmental and anthropogenic variables and the variance explained by traits for each variable. The importance of each variable is described as relative importance. The responses of IAPs to the variables can be partially explained by their traits. This was quantified by the portion of the variance explained by traits (%Var explained by traits)

We included functional traits in the model to increase the explanatory power and predictive power according to the assumption that functional traits of species can explain the magnitudes and the signs of the response of the species to the environmental variables. Of the total variance explained by the environmental variables, 55% was accounted for by traits. Therefore, the functional traits were the important predictor for the response of species to the environmental factors as we expected.

We also checked whether the variable importance changed by family (Fig. S1). We found that the major environmental variables were slightly different among families. The majority of these differences were seemingly rooted in the different traits among families; therefore, we further analyzed the relationship between traits and environmental factors.

The species richness of invasive alien plants was positively correlated with isothermality (bio 3), max temperature of the warmest month (bio 5), precipitation of the driest month (bio 14), and road density (Fig. 4). On the other hand, it was negatively correlated with temperature seasonality (bio 4), precipitation of the warmest quarter (bio 18), and distance to the river and sea.

Fig. 4

The predicted richness of IAPs over the environmental gradients. These were calculated by the sum of all the response curves of predicted probabilities of occurrences of IAPs. The environmental variables used in the analysis were, in order, isothermality (a), temperature seasonality (b), max temperature of the warmest month (c), precipitation of the driest month (d), precipitation of the warmest quarter (e), distance to the river (f) and sea (g), and road density (h). Gray dots indicate raw data. The shaded areas show 95% credible intervals

Only the influence of traits on the responses of IAPs to the environmental variables of which a 95% credible interval does not contain zero were shown in Fig. 5a (the absolute value of support level ≥ 0.95 or ≤ − 0.95). The species with higher SLA preferred the lower temperature seasonality (Fig. 5b). The species with higher height preferred the environment of low isothermality (Fig. 5c). On the other hand, the relationship between seed mass trait and the high max temperature of the warmest month was positive, but it showed a threshold-like pattern (Fig. 5d). The community weighted mean of seed mass was monotonic until 29 °C of the max temperature of the warmest month, but it increased at over that temperature.

Fig. 5

The effects of traits on the response of IAPs to the environmental variables. Only the significant relationships between traits and environmental variables were shown (a). The positive support level indicates positive relationships, and the negative support level indicates negative relationships. For each relationship, the predicted community trait mean over environmental gradient was shown: SLA to the temperature seasonality gradient (b), plant height to the isothermality gradient (c), and seed dry mass to the max temperature of the warmest month (d). Gray dots indicate raw data. The shaded areas show 95% credible intervals

Discussion

Environmental conditions with high species richness of IAPs were determined by plotting the estimated species richness with changes in each environmental variable. As a result, it is confirmed that it has a linear relationship with selected environmental variables.

In general, hot spots of invasive species were observed in the less stressful climate conditions due to the low fluctuation of annual temperature (bio 3 and bio 4). Also, the higher species richness of IAPs is accompanied by the higher the temperature in summer (bio 5), during which the plant growth peaks, and the higher precipitation in winter (bio 14), during which water stress was severe. On the other hand, the richness of IAPs tended to slightly increase when the precipitation in summer (bio 18) decreased, which is presumably due to the fact that rainfall is concentrated in summer in Korea, and this heavy rainfall might adversely affect the growth of IAPs.

To sum up the responses of invasive alien plants to the climate variables, the hot spots of invasive plants were observed in less stressful climate conditions characterized by the low variability of both temperature and precipitation. This result is consistent with the ruderal or competitive strategies of IAPs (Guo et al. 2018). Also, the increasing pattern of the richness of IAPs when they are coming closer to the river corresponds to the previously reported phenomenon that the riparian wetlands are more susceptible to the invasion of alien species than other ecosystems (Pysek and Prach 1994; Hood and Naiman 2000). The vulnerability of the riparian ecosystem to plant invasions is often explained by the periodic disturbance providing an opportunity for seedling establishment and the role of water flow as a dispersal agent (Pysek and Prach 1994).

The disturbance by anthropogenic factors or water flows had positive influences on the hot spots. These results were consistent with the previous reports about the ruderal or competitive strategies of invasive plants. According to Dawson et al. (2017), coastal regions tend to have higher species richness on a global scale. Our study also showed that the closer the distance to the sea, the higher the richness of IAPs appeared. The presence of ports, which are typical pathways of invasion, can be a plausible reason for the high richness of IAPs in coastal regions (Hulme 2009; Kaluza et al. 2010). The result of Benedetti and Morelli (2017) also supported the high richness of IAPs in the areas having high road density. This positive relationship between roads and IAPs indicates that the roads may facilitate the spread of invasive species (Joly et al. 2011; Meunier and Lavoie 2012; Brisson et al. 2010).

Functional traits are relevant and important predictors for the response of each taxon to the various environment. In a single species distribution model, the functional traits of the focal species add no surplus information, but in the community data, the coefficient of environmental variables can be fine-tuned by the functional traits of each taxon. Therefore, the functional traits are essential variables to provide additional information on the distribution of all invasive species in a single unified framework. The efforts to predict the response of IAPs to the environment or their invasiveness by their traits have been made continuously, but the context-dependence characteristic of invasion hampered its prediction (Moravcová et al. 2015; Pearson et al. 2018; Novoa et al. 2020). There were still few attempts to include functional traits in a species distribution model (Regos et al. 2019; Vesk et al. 2021), but a more systematic study should be accumulated to generalize the role of the functional traits.

In our results, the plants living in a stressed environment, such as high standard deviation of temperature, tend to have a stress-tolerant strategy, which is represented by a low SLA value and a long leaf lifespan (Reich et al. 1992; Reich et al. 1997). On the other hand, the plants living in a more benign climate tend to have ruderal/competitive strategies, which is represented by a high SLA value (Lambers and Poorter 1992; Reich et al. 1997) and short leaf lifespan (Grime 1994). In this regard of view, SLA values indicate the strategy of IAPs towards temperature seasonality.

LDMC was not significant in our result, but it might be due to the relatively high Pearson correlation coefficient (r=−0.46) between LDMC and SLA; therefore, the contribution of LDMC could be masked by that of SLA. Both SLA and LDMC reflect the leaf economics, but they do it in the opposite way. The plants with high SLA and low LDMC values tend to have acquisitive economics, but those with low SLA and high LDMC tend to have conservative economics (Wright et al. 2004; Pierce et al. 2013).

According to the study of Moles et al. (2009), plant height and isothermality showed a positive correlation and relatively high explanatory power (R² = 0.222). This means that when the annual range of temperature becomes smaller relative to the mean diurnal range, then the dominance of the taller species increases in the plant community. The shorter species are preferred in the stressed condition by large temperature fluctuations. On the other hand, the taller species, which are more competitive for light resources (Westoby 1998; Aan et al. 2006; Vojtech et al. 2008), are preferred in benign conditions with small temperature fluctuation. This trade-off feature of the investment in height was discussed in Falster and Westoby (2003). In other words, shorter IAPs are preferred in a stressed condition due to environmental fluctuation, but taller IAPs are preferred in a less stressed condition.

In a more stressed condition, plants tend to have larger and heavier seeds than in a benign condition. This is partially due to the positive influence of seed weight on the establishment success (Harper et al. 1970; Smith and Fretwell 1974; Pluess et al. 2005). Temperature stress is one of the causes of the increasing seed weight; for example, Pluess et al. (2005) found the pattern of increase in seed weight of Alpine plant species with increasing elevation (Pluess et al. 2005). Not only the cold stress but also heat stress can result in inhibition of development, growth, and yield of plants (Lobell and Asner 2003; Lobell and Field 2007). Our study showed that the plants with a heavier seed trait were preferred when the maximum temperature in summer exceeded 29 °C.

We used the spatially explicit species distribution model because species occurrences are usually spatially autocorrelated. The 60% of the explained variation was contributed by spatial random effect, and it indicated that the spatial autocorrelation was strong, even after the environmental and anthropogenic variables were considered. This random effect is not just noise but modeled by latent factors to quantify the spatially structured unknown variation. According to this result, a species distribution model should consider spatial autocorrelation to increase the predictive power of the model. We used a broad range of environmental and anthropogenic variables, but there might still be important variables accounting for the portion of the random effect variation. These variables will be further explored in future studies.

This study located the hot spots of IAPs nationwide and found relevant environmental or anthropogenic factors. Also, it confirmed that the functional traits are relevant and important factors to determine the responses of IAPs to the environment. In the future, this fundamental research can be used to build a risk map by considering the expansion rate of IAPs and the socio-environmental impact. The risk map will support building more efficient and practical management plans for the IAPs.

The original survey database we used also contains other alien species in different taxa, such as mammals, fishes, amphibians, reptiles, and insects. The coordinates of their occurrence points can be obtained from the National Institute of Ecology, Korea, by submitting an official document of request at no charge. Combining IAPs data with other alien taxa, considering biotic interactions, would provide a new insight to our understanding of invasion science.

Conclusions

In this study, we systematically analyzed nationwide distribution patterns of IAPs in Korea to find where the hot spot areas of IAPs are located and which abiotic or disturbance-related factors are related to the hot spots. Moreover, we further analyzed the interaction between environment and functional traits. The predicted species richness of IAPs was high in Seoul metropolitan area (suburbs), Daegu metropolitan city, Chungcheongbuk-do Province, southwest shore, and Jeju island. This distribution map can be used to build a risk map. The hot spots of invasive plants were expected to appear in benign climates, such as low fluctuation of temperature and precipitation. Also, hot spots are related to the high road density and proximity to the river or sea. The functional traits are closely related to the ecological strategies of plants by shaping the response of species to various environmental filters, and our result confirmed this. In less stressed conditions, the IAPs having higher SLA and plant height were preferred. In heat stress environments, the IAPs having heavier seed mass increasingly appeared. These results were consistent with the previous reports about the ruderal or competitive strategies of invasive plants instead of the stress-tolerant strategy. Also, the found relationship between the traits and environmental factors can be helpful to predict invasion success based on functional traits.