For the population genomic analysis, samples were collected in six sampling sites located in the Central and Eastern Andes, within the Northern Andes of Colombia, where the species has established naturalized, non-planted populations and with individuals in reproductive stage (Fig. 1; Suppl. material 1, Table SS1).

In the Central Andes, three sampling sites were selected. One site is located in the Lagunilla sector of Murillo, Tolima Department (Central-TOL1), within the buffer zone of Los Nevados National Natural Park and characterized by a subpáramo ecosystem. The other two sites are situated in high Andean forests: the Alaska sector of Murillo, Tolima (Central-TOL2), and San Felix, Caldas Department (Central-CAL). Both sites are subject to disturbance from livestock activity.

In the Eastern Andes, the selected sampling sites include Bogotá D.C. (East-DC), the capital of Colombia and the area closest to the Eastern Hills, where U. europaeus was first recorded in the Northern Andes; Chía, Cundinamarca Department (East-CUN); and Ventaquemada, Boyacá Department (East-BOY). All three are located within high Andean Forest ecosystems and are embedded in a landscape matrix that has been modified and fragmented by human activities.

At each sampling site, mature leaves were collected from 30 adult individuals (N = 180), each separated by at least five meters to reduce the likelihood of sampling clonal individuals. Samples were labelled and preserved in silica gel. Subsequently, leaf tissue was macerated with liquid nitrogen and stored at -20°C until DNA extraction.

Genomic DNA (gDNA) was extracted using the CTAB II protocol (Doyle and Doyle 1987), following a standardized procedure for the collected plant material. DNA integrity was assessed by electrophoresis on 0.8% agarose gels. DNA concentration was quantified using a spectrophotometer and purity was evaluated by the A260/A280 and A260/A230 ratios. High-quality gDNA from 180 individuals was sent to LGC Genomics (Berlin, Germany) for genotyping by sequencing (GBS). The ddRAD-seq libraries were prepared by digesting gDNA with PstI and ApeKI restriction enzymes, and sequencing was performed using 150 bp paired-end reads on Illumina NextSeq 500/550 v.2 and NovaSeq 6000 platforms.

Sequencing data were demultiplexed using Illumina bcl2fastq v2.20, allowing one or two mismatches (or N) in barcode recognition, depending on barcode uniqueness. Adapter sequences were removed and reads not matching with restriction sites at the 5’ end were discarded. Reads were trimmed with Trimmomatic v.0.39 (Bolger et al. 2014) using a sliding window of 10 bases, retaining only those with an average Phred quality score ≥ 20. Reads containing Ns or shorter than 20 bases were also removed. Quality control of the final reads was performed using FastQC v.0.11.9 (Andrews 2010).

Due to the absence of a reference genome for U. europaeus, a de novo clustering approach was used. Combined reads from all individuals were clustered using CD-HIT-EST v.4.6.1 (Fu et al. 2012), allowing a maximum sequence divergence of 5%. Singletons and clusters supported by fewer than 20 reads were discarded. These clusters, representing consensus ddRAD loci assembled from multiple individuals, served as an artificial reference that enabled consistent alignment and comparison across samples in the absence of a complete genome. Trimmed reads were then aligned to the reference clusters using Bowtie2 v.2.2.3 (Langmead and Salzberg 2012). Although Bowtie2 is commonly used for mapping against long contigs, it is also efficient and accurate for aligning short Illumina reads to fragmented or reduced representations of the genome and is therefore appropriate for ddRAD-seq consensus clusters (see Bowtie2 manual; Chambers et al. 2023). Finally, all alignments were merged with SAMtools v.1.17 (Danecek et al. 2021) and represent ddRAD loci shared among samples.

Given that U. europaeus may exhibit predominantly disomic inheritance with potential segmental components, parallel analyses were performed under both inheritance assumptions to capture its genomic complexity. Because many population genomic tools have limited compatibility with polyploid data, analyses were conducted using diploid and hexaploid datasets when possible (ploidy = 2 and ploidy = 6, respectively). In cases where polyploid data could not be processed, only the diploid dataset was used. Comparing both approaches allowed us to assess the robustness of genetic patterns inferred under different inheritance models.

Variant calling and SNP genotyping for both datasets were performed using Freebayes v.1.3.6 (Garrison and Marth 2012) on the high-performance computing cluster at the Bioinformatics and Computational Biology Center (BIOS) in Manizales, Colombia. Resulting VCF files were initially filtered to retain only bi-allelic SNPs using BCFtools v.1 (Danecek et al. 2021). Subsequently, a two-step quality filtering process was applied using the filter function. First, SNPs with quality score (QUAL) < 20 were removed. For the diploid dataset, genotypes with sequencing depth (DP) < 5 or genotype quality (GQ) < 15 were marked as missing, whereas for the hexaploid dataset the DP threshold was set to < 20 according to Onoue et al. (2022). Then, SNPs with more than 10% missing data and those with a minor allele frequency (MAF) ≤ 0.05 were excluded. Finally, for the diploid dataset, SNPs showing suspiciously high heterozygosity (i.e. > 90% of individuals heterozygous) were identified and excluded using custom R scripts. A conservative threshold was applied to avoid removing potentially informative SNPs for population structure, as can occur when all loci deviating from Hardy–Weinberg equilibrium across the entire dataset are excluded (Pearman et al. 2022). The filtering process was applied separately to diploid and hexaploid datasets, resulting in two final VCF files (Suppl. material 2 and 3). All filtering steps, except the heterozygosity filter, were implemented using the following BCFtools commands: bcftools view -m2 -M2 -v snps input.vcf | bcftools filter -S. -e ‘QUAL<20 || FMT/DP<5 || FMT/GQ<15’ (for diploid dataset) or ‘QUAL<20 || FMT/DP<20 || FMT/GQ<15’ (for the hexaploid dataset) | bcftools filter -e ‘F_MISSING > 0.1 || MAF <= 0.05’ -O v -o filtered_snps.vcf.

Genetic analyses were conducted using a multi-locus, genome-wide approach that enabled a comprehensive assessment of genetic diversity within sampling sites and of the population genetic structure. This approach provides robust insights into the evolutionary history of populations by including loci with different genealogies, thereby increasing statistical power, reducing false discovery rates, and improving the precision of population genetic parameter estimates (Carling and Brumfield 2007).

Observed (H_o) and expected (H_e) heterozygosity were estimated at the overall level for both diploid and hexaploid datasets using the poppr v.2.9.4 and vcfR v.1.14.0 packages in R v.4.2.1 (Kamvar et al. 2014; Knaus and Grünwald 2017; R Core Team 2021). For the diploid dataset, H_o, H_e, and the number of private (i.e. unique) alleles (A_p) were estimated per sampling site using Arlequin v.3.5.1 (Excoffier and Lischer 2010). In addition, the fixation index F_IS, Hardy-Weinberg equilibrium test (P_HWE), and heterozygote excess test (P_E-Het) were calculated per sampling site using genepop v.1.2.3 in R (Rousset 2008). Overall, F_IS values for each dataset were estimated using hierfstat v.0.5-11 in R (Goudet 2005).

To evaluate genetic differentiation, Wright’s fixation indices (F_ST) were estimated pairwise among sampling sites for both diploid and hexaploid datasets using 10,000 permutations in the StAMPP v.1.6.3 package in R (Pembleton et al. 2013). Overall genetic differentiation for both datasets was assessed using the unweighted mean of pairwise F_ST values (F_STp), since the Weir and Cockerham (1984) estimators are not directly applicable to polyploid data. A false discovery rate correction (Benjamini and Hochberg 1995) was applied to adjust for multiple testing. To determine the percentage of variance explained within and among sampling sites, an analysis of molecular variance (AMOVA) was performed on the diploid dataset using Arlequin, with 10,000 permutations. A total of 201 grouping scenarios (hypotheses) were evaluated, representing all possible combinations of sampling sites into 2 to 5 groups, to explore alternative patterns of population structure. The models with significant among-groups variation components (F_CT) were retained for interpretation.

The number of genetic clusters (K), representing the number of inferred populations to which individuals could be assigned, was estimated using the Structure v.2.3.4 tool (Pritchard et al. 2010) implemented in ipyrad v.0.9.92 (Eaton and Overcast 2020), which uses Bayesian clustering methods. Analyses were conducted separately for diploid and hexaploid datasets, including all individuals. VCF files were converted to hdf5 format using the vcf_to_hdf5 tool. Each run included a burn-in of 10,000 iterations followed by 100,000 MCMC, testing K values from 2 to 6. The most likely number of clusters was determined using both the average log-likelihood (estLnProbMean) and the ΔK method (Evanno et al. 2005). Genetic assignment barplots for K = 2–6 were generated using the toyplot v.0.18.0 tool.

To explore patterns of genetic variation among individuals and populations, visualize population structure, and assess the degree of differentiation between sampling sites in the Northern Andes, principal component analysis (PCA) and discriminant analysis of principal components (DAPC) were also performed for both datasets using the adegenet v.2.0.1 package in R (Jombart 2008; Jombart and Ahmed 2011). For DAPC, 300 principal components (PCs) were initially retained and validated through cross-validation, dividing the data into 90% training and 10% testing sets using the xvalDapc function with. Thirty repetitions were performed to ensure robustness. One hundred PCs were retained and a value of K = 6 was set based on the Bayesian information criterion (BIC) (Jombart and Collins 2015).

The relationship between geographic distance and genetic differentiation among sampling sites was evaluated for both datasets (diploid and hexaploid) using a Mantel test (Mantel 1967), based on Spearman’s correlation and 10,000 permutations with a free permutation approach (Jackson and Somers 1989). Geographic distances among sampling sites were calculated using Google Earth Pro. Genetic differentiation matrices (pairwise F_ST) were transformed to F_ST/(1-F_ST) to linearize the relationship with geographic distance (Rousset 1997; Jiao et al. 2024). The Mantel test was conducted using the vegan v.2.6-4 package in R (Oksanen et al. 2022). In addition, clustering trees based on Kosman/Leonard genetic distance were constructed for each dataset using the stats v.4.2.1, PopGenReport v.3.0.7 (Adamack and Gruber 2014), and dendextend v.1.19.1 (Galili 2015) packages in R.

To infer the introduction history of U. europaeus from the Eastern to the Central Andes, an Approximate Bayesian Computation analysis based on the diploid dataset was performed using the Random Forest algorithm (ABC-RF), implemented in DIYABC-RF v.1.2.1 (Collin et al. 2021). The main objective was to reconstruct the direction and pattern of introduction into the Central Andes, originating from previously established populations in the Eastern Andes, which historical records identify as the first introduction sites in Colombia (Ríos 2005; Baptiste et al. 2010). Although the initial introduction from the native range was not explicitly modelled, an unsampled “ghost ancestral population” was included to represent the unknown genetic source introduced into the country (Lee and Son 2022). To reduce model complexity and enhance statistical power of scenario selection, the sequential approach proposed by Fraimout et al. (2017) was adopted. This approach involves the evaluation of two consecutive sets of scenarios, which were defined based on historical records and previous genetic structure analyses that revealed clear clustering among the sampling sites.

In a first step, three introduction scenarios were tested for the Eastern Andes sites East-CUN and East-DC, located near the Eastern Hills of Bogotá, where the earliest records of the species in Colombia are concentrated. The scenarios evaluated were: (a) independent introductions of East-CUN and East-DC from a shared ancestral population; (b) introduction to East-DC from the ancestral population, followed by introduction to East-CUN from East-DC; and (c) introduction to East-CUN from the ancestral population, and a subsequent introduction to East-DC from East-CUN (Suppl. material 1, Fig. S1). This step aimed to reconstruct the initial spread of the species within the Eastern Andes before reaching the Central Andes. The East-BOY site was excluded from this step because it did not provide relevant information regarding the likely invasion route into the Central Andes. Genetic structure analyses indicated that East-BOY represents a highly structured and genetically distinct population, exhibiting even greater divergence from Central Andes sites.

Based on the most likely scenario identified in the first step and on genetic clustering patterns revealed in previous analyses, four alternative introduction scenarios were subsequently modelled for the Central Andes sampling sites (Central-TOL2, Central-TOL1, and Central-CAL). These scenarios evaluated whether the sites were established through three independent introduction events (scenarios 1 and 2), or through two introduction events: one involving the Tolima sites (Central-TOL2 and Central-TOL1), and another involving the Caldas site (Central-CAL; scenarios 3 and 4; Suppl. material 1, Fig. S2).

In each modelling step, 30,000 and 80,000 simulated datasets were generated under a uniform prior, respectively. Effective population size priors ranged from 10 to 10,000 for sampling sites and from 10 to 100,000 for the unsampled ghost ancestral population. Divergence time priors ranged between 10 and 1,000 generations, with the constraint that the most recent divergence time was smaller than the older ones (t1 < t2 for the first step, and t1 < t2 < t3 < t4 for the second step). Population size changes following the introductions were also included. The data were summarized using summary statistics for SNPs and linear discriminant analysis (LDA). The best scenario was selected using the Random Forest algorithm with 10 replicates of 1,000 trees. Posterior probabilities and prior error rates were estimated to evaluate model performance (Fraimout et al. 2017).

In summary, both diploid and hexaploid datasets were used to estimate overall heterozygosity indices (H_o and H_e), F_ST, Bayesian clustering, PCA, DAPC, Mantel test, and clustering trees. In contrast, only the diploid dataset was used for AMOVA, per-site estimates of H_o, H_e, A_p, F_IS, P_HWE, and P_E-Het, and ABC-RF analyses (Suppl. material 1, Fig. S3), due to software limitations with polyploid data.

A total of 472.5 million raw paired-end reads were obtained from the 180 individuals, averaging approximately 2.5 million reads per sample. The reference cluster assembly consisted of 359,271 contigs. Variant calling generated 394,021 SNPs for the diploid dataset and 82,388 SNPs for the hexaploid dataset. After quality control and filtering, 3,891 SNPs were retained in the diploid dataset (with 5.6% missing data) and 2,994 SNPs in the hexaploid dataset (with 6.1% missing data). These variants were distributed across 1,867 contigs in the diploid dataset and 1,582 contigs in the hexaploid dataset. Both datasets included all 180 individuals, and quality control confirmed high sequencing quality across all samples (Suppl. material 1, Table SS2).

At the overall level, observed heterozygosity (H_o) exceeded expected heterozygosity (H_e) in both datasets: H_o = 0.302 and H_e = 0.252 for the diploid dataset, and H_o = 0.483 and H_e = 0.222 for the hexaploid dataset (Table 1). In the diploid dataset, this pattern was consistent across all sampling sites (Table 1). As mentioned in the methodology, all other diversity estimates (except for overall heterozygosity) were based solely on the diploid dataset due to software limitations. The Hardy-Weinberg equilibrium test showed highly significant deviations (P_EHW < 0.0001) both at the overall level and per sampling site, likely due to an excess of heterozygotes (P_E-Het < 0.0001; Table 1). These results were further supported by negative inbreeding coefficients (F_IS), which persisted even after filtering loci with fixed heterozygotes, suggesting that the observed heterozygosity excess underlies these patterns.

Among sampling sites, East-BOY exhibited the highest observed heterozygosity (H_o = 0.460) and the lowest F_IS (-0.493), suggesting high genetic diversity and excess of heterozygotes. In contrast, East-DC showed no private alleles, while Central-TOL1 had the highest number (A_p = 5), highlighting differences in allelic richness across sampling sites (Table 1).

The evaluation of genetic differentiation revealed moderate population structuring in the diploid dataset and weaker structuring in the hexaploid dataset. In the diploid dataset, the overall fixation index estimated using Weir and Cockerham’s method (F_ST = 0.100) was nearly identical to the unweighted mean of pairwise values (F_STp = 0.118), indicating consistent differentiation estimates regardless of the approach used. In contrast, the hexaploid dataset showed a lower overall differentiation (F_STp = 0.074), reflecting weaker population structure. Significant pairwise F_ST values were observed among all sampling sites in both datasets (Fig. 2), with the strongest differentiation between East-BOY/Central-CAL and East-BOY/Central-TOL1, and the lowest between East-DC/East-CUN.

In the Central Andes, Central-TOL1 and Central-TOL2, the two geographically closest sites, exhibited relatively low genetic differentiation (diploid: F_ST = 0.084; hexaploid: F_ST = 0.054) compared to their higher differentiation with Central-CAL (diploid: F_ST = 0.152 and 0.136; hexaploid: F_ST = 0.081 and 0.079). Similarly, in the Eastern Andes, East-BOY showed high genetic differentiation from all other sites, while East-DC and East-CUN, being geographically closer, presented lower levels of differentiation (diploid: F_ST = 0.066; hexaploid: F_ST = 0,051; Fig. 2).

The Mantel test revealed a moderate but significant positive correlation between genetic and geographic distance, suggesting a pattern of isolation-by-distance (diploid: r = 0.479, p value = 0.022; hexaploid: r = 0.507, p value = 0.012; Suppl. material 1, Fig. S4). AMOVA results showed that most of the genetic variation was found within sampling sites (86–87%), with lower proportions of variation among sites within groups (7–9%) and among groups (3–5%). A total of six structuring scenarios showing significant F_CT values were retained. The sampling site East-BOY clustered separately from the other populations in five scenarios. In contrast, Central-TOL1 and Central-TOL2 were consistently grouped together in four scenarios (Suppl. material 1, Table S3).

In the structure analysis, although the ΔK method identified K = 3 as the most likely number of genetic clusters for both datasets, the average log-likelihood (estLnProbMean) supported K = 6, which reflects the finer-scale genetic structure observed among sampling sites (Suppl. material 1, Figs S5, S6). In this context, Q represents the proportion of ancestry (membership coefficient) of each individual assigned to each genetic cluster inferred by Structure (Pritchard et al. 2000). In the diploid dataset, results showed a clearly defined assignment pattern in most populations (Fig. 3A). Central-TOL1 and Central-TOL2 exhibited complete assignment of individuals to a single cluster (100% with Q > 0.8), followed by East-BOY with 93.3% of individuals predominantly assigned. East-DC (76.7%) and Central-CAL (70%) showed intermediate levels of assignment, while East-CUN presented the most diffuse pattern, with only 60% of individuals strongly assigned. These results suggest a partially defined genetic structure, with well-differentiated populations in the Central Andes and higher levels of admixture among populations from the Eastern Andes. In the hexaploid dataset, overall patterns were similar (Fig. 3B). Central-TOL1 again showed complete assignment (100%), while East-BOY and Central-TOL2 had high assignment levels (93.3% and 86.7%, respectively). In contrast, Central-CAL, East-CUN, and East-DC showed lower proportions of strongly assigned individuals (66.7% in all three cases), indicating a higher degree of admixture compared to the diploid dataset. Overall, both datasets revealed a general correspondence between genetic structure and sampling sites. Central-TOL1, Central-TOL2, and East-BOY displayed greater differentiation, whereas Central-CAL, East-CUN, and East-DC showed higher genetic admixture. The detected admixture reveals genetic connections among some Central and Eastern populations, indicating recently shared ancestry or restricted gene flow, consistent with the significant F_ST values observed.

Principal component analyses (PCA) confirmed the genetic differentiation among sampling sites, revealing six distinguishable groups in both datasets. The first two principal components explained 15.92% of the total variation in the diploid dataset and 12.08% in the hexaploid dataset (Fig. 4A, B). In the Central Andes, Central-TOL1 and Central-TOL2, both geographically close, showed greater genetic similarity, while Central-CAL appeared more differentiated. In the Eastern Andes, East-CUN and East-DC, also geographically close, were genetically similar, while East-BOY was more distinct. However, the PCA did not reveal a clear differentiation between the two regions (Eastern vs Central Andes). Discriminant analysis of principal components (DAPC), supported by the BIC, also indicated K = 6 as the most likely number of genetic groups (Suppl. material 1, Figs S7, S8).

The Kosman/Leonard dendrograms showed six well-defined genetic groups corresponding to the sampling sites, which were clustered into two main clades. In the diploid dataset, Central-TOL1 and Central-TOL2 (Central Andes) showed a close genetic similarity with East-CUN (Eastern Andes), while East-BOY and East-DC (Eastern Andes), grouped together and exhibited genetic similarity with Central-CAL (Central Andes) (Fig. 5A). A similar clustering pattern was observed in the hexaploid dataset, although Central-TOL2 appeared more genetically similar to East-CUN than Central-TOL1 (Fig. 5B). These patterns indicate closer genetic affinities, i.e. more recently shared ancestry, between some Eastern and Central Andean sampling sites than among sites within each region.

In the ABC-RF analysis, the first step indicated that the most likely scenario for the introduction of U. europaeus to the Eastern Hills involved at least two independent introduction events (scenario 1): one into East-CUN and another into East-DC. In both cases, propagules originated from an ancestral, unsampled source population introduced into the country (Fig. 6; Table 2; Suppl. material 1, Fig. S1A).

The second step of the analysis suggests that the most likely origin of the sampling sites in the Central Andes involved two independent introduction events: one into the Tolima sites (Central-TOL1 and Central-TOL2), originating from East-CUN; and another into Central-CAL, originating from East-DC (Fig. 6; Table 2). The most supported scenario (scenario 4) further proposes that the introduction in Tolima occurred first at Central-TOL1 from East-CUN, followed by dispersal to Central-TOL2 from Central-TOL1 (Fig. 6). However, the distinction between this and the alternative scenario (scenario 3), in which Central-TOL2 was colonized first, was only moderate, as both received a similar proportion of RF votes (Table 2; Suppl. material 1, Fig. S2C–D). These results suggest that the introduction of U. europaeus to the Central Andes occurred through at least two pathways, each tracing back to distinct Eastern Andes sources, consistent with a scenario of multiple introductions followed by localized dispersal events.

Genetic diversity and population structure analyses performed with both diploid and hexaploid datasets revealed similar patterns. Diversity estimates showed consistent trends across assumed ploidy models, and the genetic structure analyses indicated similar clustering among sampling sites. These results suggest that U. europaeus populations in the Northern Andes exhibit a predominantly disomic behaviour, such that the diploid model adequately describes the patterns observed in the hexaploid dataset. Any potential segmental behaviour would likely be restricted to specific multicopy genomic regions, such as the ETS and ITS (Bellot et al. 2023), rather than representing a genome-wide pattern.

Hornoy et al. (2013) proposed a two-stage invasion process for U. europaeus: initial dispersal from Spain to Northern Europe, followed by anthropogenic introductions elsewhere. During this process, despite the loss of rare alleles, overall genetic diversity was not significantly reduced. Although native populations were not included in the present study, the results may be consistent with Hornoy et al. (2013) in that the loss of rare alleles does not necessarily lead to a substantial reduction in total genetic diversity. This comparison should, however, be interpreted with caution, as SNP and microsatellite markers differ in mutation rates and evolutionary dynamics, even though both marker types can reveal similar patterns of genetic structure (Coates et al. 2009).

The populations analysed showed no evidence of heterozygosity loss. On the contrary, a consistent excess of heterozygotes was detected, which may be explained by several factors. On the one hand, the founding individuals may have been predominantly heterozygous, providing an initial evolutionary advantage (Bock et al. 2015). Likewise, hexaploidy can facilitate the accumulation and maintenance of genetic variation by allowing multiple chromosomal sets (Mounger et al. 2021). This excess may also be associated with hybrid vigour, as heterozygous individuals often exhibit higher fitness and adaptive capacity. Another possibility is the presence of ongoing positive selection favouring heterozygous genotypes, thereby maintaining high levels of genetic diversity even after introduction events (Mounger et al. 2021).

Beyond these intrinsic mechanisms, demographic processes commonly associated with recent invasions could contribute to the high levels of heterozygosity observed. In particular, introduced populations can exhibit elevated heterozygosity due to the recent admixture of multiple introduced lineages, which rapidly increases genetic variability by combining alleles from different source regions (Keller et al. 2014). However, it remains unknown whether Northern Andes populations descend from multiple European genetic lineages, as native-range populations were not included. Nevertheless, given that the bioclimatic conditions of the Northern Andes differ both from the native range and from other invaded regions, this excess heterozygosity may have facilitated climatic niche expansion (Ángel-Vallejo et al. 2025) and, consequently, enhanced the invasive success (Keller et al. 2014) of the species in this region.

Pairwise F_ST values from both datasets were significant among all sampling sites, indicating moderate population structuring and restricted gene flow (non-random mating) even among geographically close sites. Although the Central and Eastern Andes are separated by the Magdalena River Valley and lie above 2600 m a.s.l., a potential biogeographic barrier, no clear genetic clustering by mountain range (Eastern vs Central Andes) was detected. This suggests that some level of gene flow does occur among sites, but that it is not sufficient to fully homogenize the populations.

Natural dispersal in U. europaeus occurs mainly through explosive dehiscence, zoochory, and hydrochory (Roberts and Florentine 2021), mechanisms that primarily facilitate short-distance movement. Consequently, long-distance gene flow would require additional processes capable of overcoming the species’ inherent dispersal limitations. In the Northern Andes, the limited availability of effective local agents for its spread suggests that the contribution of natural dispersal to gene flow at larger spatial scales is likely modest.

Limited gene flow within and between mountain ranges may be facilitated by human activities. Intentional dispersal, such as the use of U. europaeus in live fences (Ríos 2005), could connect sites across the Central and Eastern Andes, while accidental dispersal through the transport of seeds on machinery, vehicles, or agricultural materials (Gouldthorpe et al. 2006) likely promotes gene flow within each mountain range. Additionally, road edges provide favourable conditions for U. europaeus, enhancing colonization, establishment, and expansion, thereby connecting populations that would otherwise remain isolated (León-Cordero et al. 2016). In the Northern Andes, populations are commonly observed along roadsides, which further facilitates local gene flow.

Clustering analyses revealed that the East-DC and East-CUN populations, located near the Eastern Hills of Bogotá, exhibited higher genetic admixture and lower differentiation than the other populations. This pattern suggests that propagules introduced into other sampling sites likely originated from areas near Bogotá (Ríos 2005; Baptiste et al. 2010). This inference is supported by ABC-RF results, which identified at least two independent introduction events of U. europaeus into the Eastern Andes, followed by introductions into the Central Andes. In addition, multiple independent introduction events within the Central Andes were detected, including the spread of propagules from East-CUN into Tolima, and a separate introduction from East-DC into Caldas.

Regarding the colonization of U. europaeus in the Central Andes, the ABC-RF analysis showed some uncertainty about the sequence of establishment of the Tolima populations. However, the most likely scenario indicates that Central-TOL1 was first established from East-CUN, followed by the colonization of Central-TOL2 from Central-TOL1, whereas Central-CAL originated independently from East-DC more recently. This sequence aligns with the introduction times reported by Ángel-Vallejo et al. (2025), who estimated the time since introduction through interviews with local residents regarding how long the plant had been present on their land. Based on these accounts, Central-TOL1 is the oldest population (~40 years), followed by Central-TOL2 (~33 years), and Central-CAL (~18 years). Together, this evidence supports the inference of independent introduction events in the Central Andes. Although the East-BOY population was not included in the ABC-RF analysis, clustering and genetic structure results indicate that it is genetically distinct, suggesting a separate introduction event, possibly from a different propagule source or colonization pathway.

In conclusion, the results of this study indicate that populations of U. europaeus in the Northern Andes maintain consistently high levels of heterozygosity. This pattern supports the hypothesis that polyploidy may contribute to preserving genetic variation and enhancing the species’ adaptability to a variety of environmental conditions, thereby facilitating its persistence and expansion in invaded ecosystems. Genetic structure analyses revealed recent differentiation among populations and restricted gene flow. While biogeographic features such as the inter-Andean valley may limit dispersal in the region, our results did not detect an effect of this barrier on the clustering pattern and instead suggests that other factors, including human-mediated dispersal, play an important role in shaping genetic connectivity. Finally, the inferred introduction history points to at least two independent introduction events in the Eastern Andes, followed by multiple introductions into the Central Andes, highlighting the complexity of U. europaeus invasion in high mountain ecosystems. Future studies incorporating environmental data and genome-wide selection analyses could provide insights into the role of local adaptation in the success of this species.