Biological invasion leaves traceable genetic tracks on each stage of invasion, which offer indispensable information for making proper management plans and effective control practices (Arredondo et al., 2018; Shirsekar et al., 2021). Colonization of a species takes place in a complex manner associated with four consecutive stages: an initial introduction with a small-sized population, colonization followed by dramatic population growth, and population expansion (Richardson et al., 2000; Pysek et al., 2004; Theoharides and Dukes, 2007; Blackburn et al., 2011). Although the primary ecological and evolutionary forces driving each invasion process may greatly differ (Theoharides and Dukes, 2007; Reed et al., 2020), genetic changes are inevitable due to the demographic changes and the accompanying altered selective pressure. Accordingly, invasion processes can be reconstructed by trailing genetic marks through a population genetics study, which is particularly important for invasive species with lack of well-documented records.

In the initial stages, introduced species likely experience both genetic and stochastic challenges because species colonization is often associated with low propagule pressure, i.e., the combined effect of the propagule size and number (Baker, 1955; Lande, 1988; Simberloff, 2009; Blackburn et al., 2015; Dlugosch et al., 2015; Blackburn et al., 2020). The small founding populations derived from the limited propagule pressure accompany several genetic challenges, e.g., deprivation of genetic diversity, loss of advantageous alleles, elevated genetic load, and altered selection pressure (Frankham, 1995; Blackburn et al., 2015; Dlugosch et al., 2015). However, despite the complications, some introduced species successfully colonize and widely spread in the new area. The genetic paradox (Allendorf and Lundquist, 2003) can partly be explained by the continued propagule pressure (Simberloff, 2009; Blackburn et al., 2020). The colonizing alien species might continuously receive propagules from genetically divergent sources, ameliorating the problems caused by the small founding populations. The flow of continued propagule pressure could increase the existing genetic variation, and it could also bring up new genotypes that are better adapted to novel environments (Novak and Mack, 1993; Mccauley, 1997; Sax et al., 2007; Dlugosch et al., 2015; Barker et al., 2019). In fact, multiple introductions of divergent lineages and/or congeneric species are rather commonly observed in widespread invasive plants, e.g., saltcedar, Tamarix spp., and Japanese knotweed, Reynoutria japonica (Lee, 2016; VanWallendael et al., 2021).

During the last phase of biological invasion, several factors including population characteristics, e.g., growth rates and landscape heterogeneity, affect the rates and direction of the spread (Theoharides and Dukes, 2007). Environmental gradients across the landscape, dispersal vectors, and demographic constraints are the main determinants of the range expansion, which largely influence the genetic architecture of the invasive species (von der Lippe and Kowarik, 2007; Hulme, 2009; Wilson et al., 2009; Sherpa and Després, 2021). Several studies showed the importance of dispersal vectors, e.g., water courses and roads, as potential corridors (Christen and Matlack, 2009; Aronson et al., 2017; Lee et al., 2018; Ward et al., 2020). Also, several studies tested the barrier effects of environmental factors such as temperature and precipitation (Koncki and Aronson, 2015; Lee et al., 2018). In plants, among-population genetic variability is strongly influenced by barriers and corridors of dispersal. Landscape genetics, the amalgamation of molecular population genetics and landscape ecology (Manel et al., 2003), grants a powerful tool for understanding the dynamics of range expansion, particularly for invasive plants. It can provide ecological and evolutionary insights into the landscape factors controlling the dispersal of plants (Ouborg et al., 1999). Anthropogenic aspects should also be taken into consideration because potential effects of strong and complex human interventions such as multiple introductions and human-mediated long-distance dispersal may rescind the landscape effects.

Burcucumber, Sicyos angulatus L. (Cucurbitaceae; 2n = 24), is an annual vine widely distributed in the east side of North America (Walker, 1973; Waminal and Kim, 2015). The plant commonly occupies open spaces and moist areas such as river banks, flood plains, agricultural land, fences along roads, and rarely woods (Walker, 1973). In the last few decades, burcucumber has become one of the most troublesome invasive plants in the northern hemisphere (Smeda and Weller, 2001; Kurokawa et al., 2009; Lee et al., 2015). According to the Global Biodiversity Information Facility (GBIF), the species occurs from over 4,000 localities in 10 countries, of which eight countries are out of the native range. Given the violent vining habit over small native plants and crops, the plant poses enormous economic—particularly for corn and soybean production—and ecological threats (Walker, 1973; Smeda and Weller, 2001; Farooq et al., 2017). In one growing season, about 18,000 g of biomass can be expected from a single burcucumber plant (Esbenshade et al., 2001). The adult plants easily outgrow the natives replacing the natural vegetation, which poses great threats to biodiversity and ecosystem integrity.

Burcucumber exhibits a series of life history traits, contributing to be a notorious invasive plant. The species is primarily an outcrossing plant (monoecious, maybe selfing—see European and Mediterranean Plant Protection Organization, 2010) flowering from summer to fall with three to 15 fertile flowers on each capitate head (Arifin and Okamoto, 2020). Diverse groups of generalist insects including bees and flies pollinate burcucumber, which have been proposed as a key component of invasion success for outcrossing plants (Ward and Johnson, 2013; Montesinos et al., 2016; Arifin and Okamoto, 2020). The plant primarily reproduces with seeds, and the seeds’ annual production is fairly high (Esbenshade et al., 2001; Farooq et al., 2017; Arifin and Okamoto, 2020). The seeds germinate throughout the whole growing season and grow fast reaching 30 cm tall in a day under favorable moist and temperature conditions (Walker, 1973; Smeda and Weller, 2001). When fully mature, the fruits are covered by spines that can be attached to passing animals efficiently and assist dispersal to rather distant areas (Walker, 1973). Burcucumber also has an ability to colonize a wide range of habitats due to its physiological plasticity. For example, the plant can germinate and grow in areas with 5°C–40°C temperatures and survive from wet to semiarid soil moisture conditions. It can even tolerate mildly saline soil (Farooq et al., 2017).

In Korea, there is no well-documented information on the introduction time and the invasion pathways of burcucumber. However, based on the scattered records, it is presumed that the species was transferred from North America as a grafting stock for cucumbers and watermelons (National Institute of Environmental Research, 2005). The first observation of burcucumber was made by a local citizen in 1989, but there was a complete lack of supporting information such as a proper publication and/or herbarium sheets (National Institute of Environmental Research, 2005). The first herbarium record of burcucumber was found in 1990, and it was collected from a central region (National Institute of Environmental Research, 2005; Korea National Arboretum, 2019). In Japan, one of the neighboring countries of Korea, the first occurrence of burcucumber was reported much earlier than the one in Korea (in 1950; Kurokawa et al., 2009). Given the frequent commercial trades and travels between the two neighboring countries, the introduction date of burcucumber in Korea was probably before the 1990s. Over the past few decades, burcucumber has rapidly expanded its range throughout South Korea replacing native plants (National Institute of Environmental Research, 2005; Moon et al., 2007). The plant is currently one of the 10 most abundant invasive plants in Korea, particularly in the riparian system (Park and Lee, 2018). On average, ca. 40% of the riparian vegetation is composed of burcucumber (Park and Lee, 2018). Unfortunately, the causal mechanisms of the successful burcucumber invasion in Korea have never been explored, thus remaining elusive.

The aims of our study are twofold: 1) to assess the genetic diversity pattern that might have been associated with multiple introductions and the demographic history and 2) to determine the contributing landscape factors during the rapid range expansion. Former genetic studies of burcucumber in Japan suggested the possibility of multiple introductions, but the sample size and the number of molecular markers used were limited (Kurokawa et al., 2009; Kobayashi et al., 2012). Despite the aggressive infestations in many places across the globe, the determinants of the successful burcucumber invasion have not been empirically explored. We employed a large number of molecular markers sampled from the whole genome and collected many populations throughout the entire country to provide a detailed and precise assessment of population structure. We also inferred the invasion history by competing likely invasion scenarios with the approximate Bayesian computation (ABC) approach. Additionally, through a species distribution model (SDM) and correlation analyses, we developed the isolation-by-resistance (IBR) model and compared it with the conventional isolation-by-distance (IBD) model to identify the contributing landscape factors during the rapid range expansion of burcucumber invasion in Korea.

The sequenced 3RAD library was 211 Gbp and produced ~14 million reads with an average GC content of 39.5%. We called over 150,000 SNP loci from the initial SNP identification. After a series of screening processes, 2,696 SNPs remained for the downstream analyses. On average, the genetic diversity of burcucumber varied across populations although the differences were not significant ( Table 1 ). The adjusted allele numbers ranged from 1.21 (GO; see Table 1 for the population acronyms) to 1.73 (HF; Table 1 ). We found the lowest mean He (= 0.06) in the UJ population at which the number of collected samples was the smallest, whereas the highest mean He was observed in HF (mean He = 0.24; Table 1 ). Pairwise genetic divergence among populations estimated as F _ST greatly varied across the population pairs (0.10, NI/NO – 0.77, SG/UJ; Table 2 ). The mean F _ST was 0.37, and all F _ST values were statistically robust (P < 0.01). Based on the sign and Wilcoxon’s rank tests, most populations sampled in our research experienced recent population bottlenecks ( Table 3 ). In five populations (YT, GF, GJ, DH, and YD), the P values were high (P > 0.05) showing that the populations are likely in the equilibrium state of population demography ( Table 3 ). However, among the five populations, three populations showed a significant allele frequency mode shift compared with the equilibrium state, indicating the possibility of recent population bottlenecks ( Table 3 ).

The results of clustering analyses differ between PCoA and fastSTRUCTURE (Figs. 1 and 3). In our PCoA analysis, the first three axes only accounted for ~15% variation ( Figure 3 ). UJ and YT were clearly separated from the remaining populations along the first axis, yet the variation explained by PC1 was about 6% (plots associated with second and third provided in Figure S1 ). The second axis (PC2 = 4.9%) separates SG populations from the rest, whereas SG and GJ were clustered and departed from the rest of the populations along the third axis (PC3 = 4.3%; Figure 3 ). The optimal K, the number of clusters at HWE, inferred from FastSTRUCTURE ranged from 2 to 7 (K = 3 presented in Figure 1 ; plots for K = 2–7 provided in Figure S2 ). The overall pattern of clustering reflected the originated river basin, and it was consistent throughout the varying cluster numbers with slight differences as K increases ( Figure S2 ). Here, we summarized the result of the clustering pattern for K3 since it offers a simple and spatially more intuitive visual presentation. Overall in K3, populations within the river basins shared similar assignment patterns forming a deme ( Figure 1 ). The HG, GG, and ND river basins (see Table 1 for the river basin acronyms) were primarily assigned to a cluster represented in red, but a small fraction of the genomic composition was also affiliated with the remaining two clusters ( Figure 1 ). The two populations in the YS river basin showed a contrasting assignment pattern. YF shared a rather similar pattern with the populations within the HG, GG, and ND river basins, whereas YT was much alike with populations within the DR, the east coastal streams ( Figure 1 ). Populations in DR streams exhibited a little more complex pattern than populations located in the rest river basins, which are composed of more admixed genotypes with all the three genetic clusters ( Figure 1 ). In pPost-hoc AMOVA analyses associated with the FastSTRUCTURE clustering pattern, the average population differentiation across the river basins (the percentage variation explained = 4.5%; F _CT = 0.05, P < 0.01) was much less than the population differentiation within each river basin (the percentage variation explained = 32%; F _SC = 0.34, P < 0.01; Table 4 ).

Based on the LDA plot for the prior model checking, the location of the observed data was within the range of the simulated data, which suggests that ca. 3,000 simulated data sets per model were suitable for inferring the model choice ( Figure S3A ). Likewise, the out-of-bag error rates were stabilized with an increasing number of trees indicating that the number of trees produced for each model was sufficient ( Figure S3B ). The classification vote was the highest in scenario 7 (195) with a posterior probability of 0.67, which was followed by scenario 4 (127). The global and local errors were 0.093 and 0.156, respectively. The DIYABC-RF analysis strongly suggested that there were multiple introductions from two genetically divergent sources during the Korean cucumber invasion (see scenario 7 in Figure 2 ). The time of introduction (t1) was estimated approximately 70 years, whereas the divergence time in origin for the two genetically distinct sources dated back a little less than 300 years ( Table 5 ).

With DEM and bioclimatic variables, we developed the best-fitting species distribution model and prediction using the maximum entropy method for the burcucumber (MaxEnt; AUC = 0.83 ± sd 0.012; Figure S4 ). In the model, the most contributed variable was B19 (precipitation of the coldest quarter; 33.7% of total contribution), followed by DEM (21.8%), B17 (precipitation of the driest quarter; 6.4%), and B13 (precipitation of the wettest month; 5.8%, Table S1 ). The contributions of other variables for the model were below 5%. The optimal habitats inferred from the SDM were distributed around big cities with a river, where high disturbances are expected ( Figure S5 ). The resistance distance estimated varied greatly across population pairs (ranging from 7, TF/GJ to 105, YT/GO; Table S2 ). In our MLPE models, models with F _ST outperform the linearized F _ST whereas the three different distances computed showed a similar model fit ( Table S3 ). The model fits of the linearized F _ST dropped about 20% compared with the F _ST , and the confidence intervals for the correlation coefficient values include zeros for all three distance models ( Table S3 ). The best predictor of the genetic discontinuity was the log-transformed Euclidean distance ( Table 6 ). Models associated with resistance distance showed low support (ΔcAIC ≥ 10; Burnham and Anderson, 2002).

The major components of successful colonization consist of high propagule pressure, multiple introductions, release from natural enemies, and competitive life history traits (Sakai et al., 2001; Theoharides and Dukes, 2007; Hirsch et al., 2017). Despite the short colonization period (less than a half-century) and supposedly low propagule pressure during the initial introduction, burcucumber invasion in Korea has been a great success (National Institute of Environmental Research, 2005). However, given the lack of well-documented records, the invasion history of burcucumber has been largely unknown. Through molecular and landscape data, we partly reconstructed invasion paths of burcucumber in Korea. The ABC and Bayesian assignment results revealed that the burcucumber invasion has been tightly associated with multiple introductions from more than two genetically divergent sources from the origin. In contrast, the landscape genetics results indicated that landscape factors did not primarily contribute to the burcucumber invasion in Korea.

Most populations investigated in the study harbored rather low genetic diversity ( Table 1 ). The genetic diversity observed in our study was a little lower than or comparable with the ones measured from other invasive plants. For instance, the expected heterozygosity observed (mean He = 0.17) in the study was lower than the long-lived invasive shrubs’ [Tamarix ramosissima, mean He based on 1,996 SNPs = 0.30 (Lee et al., 2018); Frangula alnus, mean He based on 133 SNPs = 0.29 (de Kort et al., 2016)], an invasive perennial herb’s [Reynoutria japonica, mean He based on 12,912 SNPs = 0.30 (VanWallendael et al., 2021)], and an annual exotic’s [Mimulus guttatus, mean He based on 69 SNPs = 0.32 (Pantoja et al., 2017)]. The reduced genetic variability was a little surprising as the mating system of burcucumber is mostly outcrossing (Arifin and Okamoto, 2020). Outcrossing plants, in general, are genetically more diverse than self-compatible plants (Hamrick and Godt, 1996). Due to lack of native samples, the direct comparison between native and the introduced populations of burcucumber could not be carried out in this study. Accordingly, extreme care must be taken when interpreting the lowered genetic diversity. Nevertheless, considering the recent population bottlenecks observed from most populations examined in South Korea ( Table 4 ), the low genetic diversity likely resulted from the founder effect that is commonly expected for a colonizing species.

The inflated genetic divergences among the regional populations [S. angulatus mean F _ST = 0.37] might also be driven by the drastic decline in population sizes. F _ST , the among-population differentiation measure used in the study, is inversely related to the effective population size as suggested by Wright (1943); Wright (1950); F _ST ≈ 1/4N_e m +1), which was empirically verified by a number of studies, e.g., Jordan and Snell (2008) and Whiteley et al. (2010). The small population size likely resulted from the founder effect might also be the cause of the increased F _ST in our study. Alternatively, the mating system of burcucumber (mostly outcrossing by small insects, Arifin and Okamoto, 2020) might partly contribute to the high F _ST values. Life history traits largely impact the genetic diversity pattern in both within- and among-population levels (Hamrick and Godt, 1996). A wind-pollinated long-lived tree, in general, exhibits high genetic diversity with reduced population differentiation, whereas a herbaceous plant with insect pollination often harbors relatively small genetic variation and high population divergence (Hamrick et al., 1981; Olson et al., 2016). Particularly, the herbs pollinated by small insects likely show greater among-population divergence, which may also be the causal mechanism of the increased population divergence observed in our study (Gamba and Muchhala, 2020).

The spatial clustering pattern suggests multiple introductions of divergent lineages during burcucumber invasion in Korea. In K = 3, four populations (GF, GJ, YT, and NT) were assigned nearly to a single cluster, suggesting that those populations might be the sources of an early establishment ( Figure 1B ). Except for YS, populations within the same river basin shared a similar clustering pattern ( Figure 1 ). Contrary to the assignment pattern of the populations within a river basin, we found an abrupt shift in population allele frequencies between YF and YT along the YS river basin, which is repeated between TF and GJ, neighboring populations in DR streams ( Figure 1 ). The burcucumber fruits carrying water-impermeable seeds (Walker, 1973) can travel by floating on water stream, thus contributing to gene flow among populations within a river basin. Considering the geographic proximity, the abrupt change of allele frequency observed between the two neighboring populations along the YS river basin and DR streams was a rather surprising result. The steep allele frequency change likely resulted from human interventions such as repeated independent introductions from multiple sources that are genetically diverse. The pattern of the multiple introductions coincides with the historical records, although those records are incomplete and insufficient to reconstruct the whole burcucumber invasion in Korea. Notably, in the occurrence records, we found that after just 1 year since the first report in 1989, the burcucumber was recorded in several locations throughout South Korea (National Institute of Environmental Research, 2005). The multiple occurrences at such large spatial scale within a year may only be possible by multiple introductions. Alternatively, varying selection regimes responding to sudden environmental changes between the neighboring areas might result in similar genetic discontinuities. However, the effects of selection tend to affect relatively small regions of the genome (e.g., tb1 in maize; Kaplan et al., 1989; Andolfatto, 2001; Charlesworth et al., 2003; Baucom, 2016). Therefore, genetic discontinuity assessed from genome-wide molecular markers in our study is more likely derived from demographic factors such as multiple introductions (Black et al., 2001).

Consistent with the spatial clustering pattern ( Figure 1 ) that we observed, the best model inferred from DIYABC results strongly supported the multiple introductions. The scenario of choice demonstrated that there were at least two genetically divergent sources introduced from the origin ( Figure 2 scenario 7). According to the selected scenario, the two distinct sources diverged into the clustering groups represented by red and green after the initial colonization (see Figure 1 for the cluster colors). Subsequently, the dominant red group has likely spread into the middle west (GS and GF) and diverged into the blue group along the east coastal area while passing through the middle area (NS, NF; Figs. 1 and 2 scenario 7). Although the selected scenario shed some light on the path of the burcucumber invasion in Korea such as a strong association with multiple introductions, some of the detailed assessments should be interpreted with caution. For example, the number of genetically distinct sources may not be precise as we constructed the scenarios with summarized clustering groups based on K = 3. In fact, scenario 4 ( Figure 2 scenario 4) which assumes the three independent introductions and subsequent divergences to the clustering groups represented by red, blue, and green showed the second highest classification votes (votes for scenario 7 = 195; votes for scenario 4 = 127). Considering the lack of genotype data from the origin, the detailed demographic history inferred in our study was not without caveats. However, by employing the ABC approach, we were able to determine the tight association of the multiple introductions with the burcucumber invasion in Korea.

The result of our landscape genetics models was not consistent with what was expected from the previous studies. Because circuit theory-based IBR accounts for spatial heterogeneity, IBR, in general, outperforms the conventional IBD (McRae, 2006; McRae and Beier, 2007; Emel et al., 2021). However, in our landscape genetics model, circuit resistance was not the best predictor of genetic discontinuity, suggesting a general lack of environmental effects on population divergence. The resistance distance was computed based on the SDM; thus, the resistance is primarily attributed to the temperature-related variables and DEM, which is also correlated with the temperature. Given the temperature-dependent germination and growth habits of burcucumber (Walker, 1973; Smeda and Weller, 2001; Onen et al., 2018), it was unexpected that the resistance distance was not the key predictor of the population divergences. The repeated introductions from genetically divergent sources might have driven the weak environmental effects. Smith et al. (2020) provided a compelling example of the reduced effect of environmental gradient on genetic structure in the non-native range. Consistent with our results, the most significant predictor of genetic discontinuity for a widespread weed, Plantago lanceolata, was the geographic distance in the colonizing area, whereas native ones were influenced more by environmental factors (Smith et al., 2020). Coupled with the population bottlenecks found in our study, multiple introductions and the genetic admixture likely lowered the environmental effects on the spatial structure among populations. With the aid of demographic events, burcucumber populations might have rapidly expanded their range without a need for adaptation to the local environments. Some might argue that our landscape genetics approach does not consider the underlying assumption, the gene flow-drift equilibrium among populations (Marko and Hart, 2011; Epps and Keyghobadi, 2015). However, the equilibrium is an unlikely scenario for an invasive plant with the rapid and recent history of colonization such as the burcucumber that invaded in Korea.

Burcucumber invasion in Korea has been tightly linked to a suite of factors including the continued propagule pressure from multiple and deliberate and/or accidental introductions from unknown seed sources, unplanned human aids in dispersal, founder effects, river network effects, and maybe countless unaccounted factors. As shown in a number of empirical studies (Hirsch et al., 2017; van Boheemen et al., 2017; Bertelsmeier and Keller, 2018), multiple introductions vectored by anthropogenic causes enhance the propagule pressure and strongly influence the adaptive potential for a widespread invasive species (Simberloff, 2009; Dlugosch et al., 2015; Bertelsmeier and Keller, 2018; Blackburn et al., 2020). Not all invasive species, therefore, suffer from environmental constraints in the introduced range; instead, a merger of divergent genetic lineages from multiple populations can ameliorate the adaptive challenges. Our work thus addresses the importance of preventing secondary introductions, particularly for aggressive weedy plants like burcucumber. A flow of continued propagule pressure and a chance of admixture through the human-mediated secondary introductions may increase the chance of further range expansion in Korea. It is highly recommended to build a periodic and well-planned monitoring system for the new introductions particularly to the ports along the East and West Seas. Preventing further introductions of burcucumber must be the number one priority in management strategy in Korea. We also hope that our study motivates to further investigate the genomic and population-level diversity pattern of the native region and fully reconstruct the invasion paths.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, PRJNA822832.

S-RL and DCS conceived ideas and designed the study. DCS prepared research grant. Field sampling was planned and performed by S-RL. S-RL designed the laboratory work and performed the genetic and statistical analyses. S-RL wrote the manuscript. DCS edited the manuscript. All authors contributed to the article and approved the submitted version.

This work was supported by a research project of the Korea National Arboretum (KNA1-2-39, 21-2).

We are grateful to a postdoctoral fellow, Tae-Young Choi, for the laboratory and analytical assistance. We also thank the graduate assistant, Eun-Su Kang, with the undergraduate assistants, Tea-Young Eom, Hyeon-Su Kim, and Geon-Ho Kim, for their help with field sampling and laboratory work.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.