Genomic differentiation and conservation relevance of isolated marginal Primula farinosa populations in Central-Eastern Europe

Quaternary climatic oscillations produced heterogeneous range dynamics in European cold-adapted plants. Primula farinosa L. has been discussed in the context of the interglacial expansion model for moisture-demanding cold-adapted species. However, it remains unclear whether isolated low-elevation populations in Central European calcareous fens and wet meadows represent long-term relicts or recently established populations. We used RADseq-derived SNPs from four sampling localities in the Balaton Region of Hungary and one low-elevation comparative locality at Lipce in the Julian Alps, Slovenia, to assess population-genomic structure, differentiation, allelic diversity and demographic history. Population structure was analyzed using SNP-based maximum-likelihood tree inference, PCA, model-based ancestry clustering and pairwise F ST estimates. We further quantified rarefied allelic richness and private allelic variation and performed exploratory population-specific historical effective population-size trajectories from folded site-frequency spectra. Population structure analyses revealed pronounced differentiation between Lipce and the Hungarian samples. Within Hungary, the two neighbouring localities of Lesencetomaj and Tapolca were nearly indistinguishable, whereas Köveskál and Nemeshany populations showed stronger differentiation and reduced allelic diversity. Rarefied allelic richness was highest in Lipce and Lesence and lowest in Nemeshany, with Köveskál only slightly richer. Collectively, the Balaton Region carried a distinct allelic pattern and comparable rarified richness at the broader regional scale. Demographic reconstructions suggested population-specific effective population-size changes over recent to Late-Glacial timescales, although older intervals were weakly resolved and exact timing was uncertain. Low-elevation Primula farinosa populations in Hungary and Slovenia are strongly differentiated, and this isolation is unlikely to be explained solely by recent anthropogenic wetland fragmentation. The higher rarefied allelic richness of Lipce and Lesence suggests that these populations retain less-eroded genetic variation, whereas the depauperate allelic profiles of Köveskál and especially Nemeshany are consistent with recent founder effects, bottlenecks, strong drift or long-term small effective population size. These findings indicate that the Balaton Region populations represent a conservation-relevant and regionally differentiated component of P. farinosa genomic diversity. Although the exact timing and historical origin of their isolation remain unresolved, the observed genomic pattern is compatible with, but does not demonstrate, long-term regional persistence, including a possible Late Glacial refugial scenario.

The effects of Quaternary glacial-interglacial cycles on European plant distributions have often been interpreted through the classical temperate-refugia paradigm (Hewitt, 1999). Under this model, many temperate taxa persisted during glacial maxima in southern refugia, especially in the Iberian, Italian and Balkan peninsulas, and subsequently expanded northwards after climatic amelioration (Hewitt, 2004; Tzedakis et al., 2013). This paradigm has been highly influential, but it does not describe the histories of cold-adapted taxa. Species associated with arctic, boreal, montane or cold wetland habitats, often termed as ‘cold-adapted’, may respond idiosyncratically to climatic oscillations because their distributions are constrained not only by temperature, but also by growing-season length, moisture, soil chemistry, competition and the availability of suitable habitats (Birks, 2008; Crawford, 2008; Stewart et al., 2010).

Two main hypotheses have been proposed for cold-adapted species during Late Quaternary climatic cycles. The better-established interglacial contraction hypothesis assumes that cold-adapted taxa expanded during glacial periods into extensive ice-free peripheral, low-elevation or lowland tundra-steppe habitats, and subsequently contracted during interglacials into areas where suitable cool or moist conditions persisted, including northern regions, mountain systems and locally favourable refugial habitats (Hewitt, 2004; Stewart et al., 2010). This scenario is supported by several studies of alpine and arctic organisms, including European high-mountain plants, arctic-alpine invertebrates and cold-adapted vertebrates, which show patterns consistent with larger glacial distributions followed by postglacial or interglacial retreat and fragmentation (Kropf et al., 2003; Muster and Berendonk, 2006; Dalén et al., 2007; Westergaard et al., 2011; Espíndola et al., 2012). Under this model, present-day disjunctions are interpreted mainly as remnants of formerly broader cold-stage ranges.

The alternative interglacial expansion hypothesis proposes a different demographic mechanism. Some cold-adapted species may not have expanded widely during full glacial conditions, despite the colder climate, because suitable habitats were not determined by temperature alone. For moisture-demanding, weakly competitive plants, the availability of persistently wet, open and calcareous habitats may have been more limiting than regional temperature, especially where full-glacial landscapes were dominated by open, continental and often dry vegetation mosaics (Schmitt, 2007; Harrison and Prentice, 2003; Fletcher et al., 2010; Magyari et al., 2010). According to this hypothesis, such species persisted during glacial maxima in relatively moist low-elevation, peripheral or near-glacial refugia, and expanded during warmer interglacials when retreating glaciers exposed new open habitats and when suitable wet habitats became available at higher altitudes and latitudes (Birks, 2008; Stewart et al., 2010). This mechanism is expected particularly for cold-adapted species with broad altitudinal or ecological amplitude, rather than strict high-alpine specialists, and for taxa whose present ranges combine Central European mountain occurrences with southern-boreal or lowland northern occurrences (Ehrich et al., 2008; Lienert et al., 2002; Oliver et al., 2006; Ronikier et al., 2008; Šmídová et al., 2011; Urbaniak and Kwiatkowski, 2023).

The biological logic of the interglacial expansion model is especially relevant for wetland and fen specialists. During the Holocene, warming and the spread of temperate forests and tall-herb vegetation likely increased competition in many low-elevation landscapes. Cold-adapted plants with low competitive ability could therefore persist outside mountains only where edaphic or hydrological conditions maintained open vegetation, for example in calcareous fens, river gravels, cliffs or other treeless refugial habitats (Birks, 2008; Hájek et al., 2006, Hájek et al., 2011). The key point is that low-elevation populations of cold-adapted plants are not automatically ecological anomalies or recent sink populations. In some cases, they may represent remnants of formerly broader distributions, preserved by local habitat conditions while surrounding landscape became climatically or competitively unsuitable.

Primula farinosa L. is key empirical system in which the interglacial expansion hypothesis has been proposed to explain present-day distribution and genomic diversity (Theodoridis et al., 2017). It is a cold-adapted, insect-pollinated, heterostylous and mostly outcrossing perennial associated with moist calcareous habitats, including calcareous fens, wet meadows, lake and stream margins and, at higher elevations, calcareous grasslands (Hambler and Dixon, 2003; Richards, 2003). In southern Europe, including the Alps, Carpathians and Pyrenees, it is primarily a montane to alpine species, whereas in northern Europe it also occurs in low-elevation flatlands (Hambler and Dixon, 2003; Richards, 2003; Theodoridis et al., 2013). Its broad, fragmented European distribution and its inferred Late Quaternary evolutionary history make it a particularly informative system for testing how cold-adapted plants responded to climatic cycles (de Vos et al., 2014; Stewart et al., 2010).

Theodoridis et al. (2017) explicitly tested the interglacial expansion hypothesis in P. farinosa by integrating species distribution models projected across the last 72 thousand years (ky) with genome-wide reduced-representation sequence data and approximate Bayesian computation. Their analyses supported a distribution-wide interglacial expansion rather than the classical expectation of broader glacial distributions followed by interglacial contraction in cold-adapted taxa. They also argued that this result is biologically plausible because P. farinosa has broader climatic preferences than several closely related alpine/arctic primroses, prefers wet open microhabitats, and is a weak competitor; consequently, it may have been constrained by dry glacial steppes but favoured by the availability of open, moist habitats after glacial retreat (Lindborg and Ehrlén, 2002; Theodoridis et al., 2013, Theodoridis et al., 2017).

The Hungarian low-elevation occurrences of P. farinosa are especially relevant in this framework. In the Carpathian Basin, including Hungary, vegetation reconstructions indicate dry, Artemisia -rich steppes at low altitudes during the Last Glacial Maximum (LGM) (Magyari et al., 2010; Feurdean et al., 2014; Magyari et al., 2014). Yet, the species has long-known natural occurrences (Borbás, 1900) from fen meadows in the Balaton Region (incl. the Balaton Uplands), where persistence depends on very specific wet, open, calcareous microhabitats (Salamon-Albert and Morschhauser, 2003). These occurrences are commonly reported as ‘glacial relicts’ in the local literature (Borbás, 1900; Borhidi et al., 2012). However, such terminology requires careful evaluation, because present-day occurrence in isolated fen habitats does not by itself establish uninterrupted glacial persistence. Central European calcareous fens are groundwater-fed, nutrient-limited habitats with distinctive plant and bryophyte assemblages, and several fen specialists show distributions consistent with long-term regional persistence or Holocene-scale relict dynamics (Hájek et al., 2006, Hájek et al., 2011).

This relict interpretation raises a testable population-genomic expectation and defines the contribution of the present study. If Balaton Region populations represent long-term relicts (i.e., preglacial or glacial) or remnants of older regional persistence, they may be expected to retain genomic signatures often associated with ‘rear edge’ populations (Hampe and Petit, 2005), such as strong differentiation, reduced contemporary connectivity, and population-specific components of allelic diversity. Conversely, if they are recent Holocene colonists or introductions, they may show patterns more consistent with recent founder events, reduced allelic richness and limited private variation. Because our sampling includes only one comparative low-elevation Alpine-region population, Lipce is treated as a regional comparative reference for evaluating the distinctiveness of the Balaton Region populations.

By examining these alternatives with RADseq-derived SNPs from five Central European populations, the study extends the species-wide interglacial-expansion framework for Primula farinosa to a finer regional scale and evaluates whether isolated low-elevation fen populations preserve distinct components of genomic diversity relevant for conservation. Specifically, we asked: (i) how genomic variation is structured among the Balaton Region populations and the comparative Lipce population; (ii) whether allelic variability and genomic differentiation are compatible with long-term regional persistence, recent colonization or strong drift in isolated populations; and (iii) whether site-frequency-spectrum-based demographic reconstruction reveals population-specific changes in effective population size, while acknowledging the limited ability of these data to resolve the exact timing of population origin or isolation.

Sampling and study populations

RADseq data were generated for 47 individuals assigned to five sampling localities: Lipce, Köveskál, Nemeshany, Tapolca: Körtvélyes and Lesencetomaj: Sövényháza (Figure 1; Table 1). These localities were initially treated as separate population units in the bioinformatic analyses, with the labels Lip, Kov, Nem, Kor and Sor, respectively. Lipce was included as a low-elevation comparative locality from the Alpine region, occurring in marshland habitat in the Julian Alps, Slovenia. All other localities represent isolated low-elevation occurrences of P. farinosa in the Balaton Region, Hungary, traditionally interpreted as glacial relicts in the local floristic literature. Tapolca: Körtvélyes and Lesencetomaj: Sövényháza are located in the same headwater basin of the Lesence Stream in close geographic proximity, and they are separated by an approximately 0.5 km of unsuitable habitat, including a paved road. They were therefore initially analysed as separate sampling localities to test whether genetic differentiation was detectable at this fine spatial scale. Because population genomic analyses showed little or no differentiation between them, they were subsequently combined and treated as a single Lesence biological population in downstream analyses. Analyses aimed at describing fine-scale structure retained the original Körtvélyes and Sövényháza labels.

Primula farinosa is a strictly protected species in Hungary, and the field sampling, which involved the removal of two leaves from each sampled individual, was conducted in the presence of the local nature conservation authority (Balaton-felvidéki National Park Directorate) under the permit issued by the Pest County Government Office (‘Pest Megyei Kormányhivatal’; permit no. PE-KTFO/1568-18/2022). The leaf tissue samples were immediately placed into a resealable sachet filled with approximately 10 g silica gel, and individual coordinates (latitude, longitude) were collected. We aimed to keep at least 5 metres between the sampled individuals, and we also endeavoured to spatially cover the full extent of the population’s area to represent the variability of the populations in the analyses.

DNA extraction and RADseq library preparation

Total genomic DNA was extracted from silica-dried leaf material using a modified CTAB protocol based on Sramkó et al. (2014). Leaf fragments were homogenized with sterile sand and polyvinylpolypyrrolidone (PVPP), and the powder was incubated in 2% CTAB extraction buffer containing β-mercaptoethanol and polyvinylpyrrolidone (PVP). Proteins and other contaminants were removed by ammonium acetate treatment, followed by chloroform:isoamyl alcohol purification. DNA was then recovered by precipitation with 96% ethanol, the pellets were rinsed twice in 70% ethanol, and the final DNA was dissolved in 10 mM Tris-HCl (pH 8.0). Isolate quality was first checked on 1% agarose gels using standard electrophoresis. DNA concentration and purity were then estimated by Qubit v3 fluorometry (Thermo Fisher Scientific Ltd., Waltham, MA, USA) and NanoDrop Lite v3 spectrophotometry. Only high-molecular-weight, visually intact, DNA extracts with a fluorometric/spectrophotometric concentration ratio between 0.9 and 1.4 were retained for library construction.

RADseq libraries were prepared according to a classical single-digest RAD protocol (Baird et al., 2008) using PstI (CTGCA/G) (New England Biolabs, Ipswich, MA, USA) as the restriction enzyme. Custom-modified Illumina TruSeq P1 adapters were ligated to the restriction overhangs, after which the DNA was fragmented using a Bioruptor Pico instrument (Diagenode SA, Seraing, Belgium) with five cycles of 30 s sonication at high power followed by 30 s pauses. Fragment pools were subjected to double-sided size selection with the SPRI Select Beads Kit (Beckman Coulter Inc., Brea, CA, USA) using SPRI beads in ratios of 0.65× and 0.7×. Custom-modified P2 adapters were subsequently ligated with T4 DNA ligase (New England Biolabs). Library preparation was carried out in sub-libraries of 12 samples, which were then pooled in equimolar proportions. Adapter-ligated fragments were amplified with Phusion High-Fidelity PCR Master Mix with HF Buffer (New England Biolabs); PCR artifacts were reduced by limiting amplification to 14 cycles and by carrying out eight replicate reactions that were pooled afterwards. The pooled library was size-selected again with SPRI beads using a bead ratio of 0.7×, targeting an insert distribution around 400 bp. Final library size profile and concentration were assessed on a Bioanalyzer system (Agilent Technologies Inc., Santa Clara, CA, USA). Sequencing was performed at the Genomic Medicine and Bioinformatics Core Facility of the University of Debrecen, Hungary, on an Illumina NextSeq 500 platform, generating 150-bp single-end reads.

RADseq processing, assembly optimization and SNP discovery

Raw RADseq reads were demultiplexed and quality filtered with process_radtags in Stacks v2.2 (Catchen et al., 2013). Reads were assigned to samples using inline barcodes and PstI as the restriction enzyme. Data were processed as Phred+33 FASTQ files with barcode rescue, standard cleaning and quality filtering enabled. Reads were truncated to 140 bp, and reads exceeding the length limit were discarded. Initial read quality was assessed with the FastQC software tool, followed by adapter removal and base-level quality filtering (Q ≤ 5) using fastp v1.0.1 (Chen et al., 2018).

In the lack of a reference genome for Primula farinosa, we applied a de novo RAD locus assembly and SNP calling approach using the denovo_map.pl module of Stacks. Before final SNP discovery, assembly parameters were optimized following the parameter-space exploration strategy proposed for RADseq datasets (Paris et al., 2017). Optimization was performed on ten representative samples by varying the minimum depth required to create a stack within individuals (m = 3–6) and the maximum number of mismatches allowed within individuals and during catalogue construction (M = n = 2–6). For each parameter combination, denovo_map.pl was run and the populations module was used to export loci present in at least 80% of individuals. Parameter combinations were compared using the number of retained loci and SNPs, sample coverage, total retained loci and individual missingness. This screening identified m = 3, M = 5 and n = 5 as the most appropriate setting for the final de novo assembly.

The final de novo assembly was therefore performed in Stacks using denovo_map.pl with m = 3, M = 5 and n = 5. SNPs were exported with the populations module, requiring loci to be present in at least four populations (-p 4). The raw VCF file generated by Stacks was filtered using vcftools v0.1.16 (Danecek et al., 2011). Sites with more than 20% missing data were removed, and the dataset was restricted to biallelic SNPs. To reduce the effect of physical linkage, one SNP was randomly selected from each retained RAD locus after filtering using a custom bash script. In the final filtering step, SNPs were retained only if mean read depth across individuals was between 5 and 60. The resulting one-SNP-per-RAD-locus dataset contained biallelic SNPs present in at least four populations, with no more than 20% missing data, mean depth between 5 and 60, and at most one randomly selected SNP per RAD locus. This dataset was used for population structure, pairwise F ST, allelic-richness and private-allele analyses.

Maximum-likelihood SNP-based tree inference

Maximum-likelihood phylogenetic tree reconstruction was performed in IQ-TREE v2.4.0 using the final filtered one-SNP-per-RAD-locus dataset (Minh et al., 2020). Because the alignment contained only variable positions, ascertainment-bias correction was applied during model selection and tree inference following the ASC approach for variable-only data (Lewis, 2001). The best-fitting model was selected using ModelFinder under MFP+ASC according to the Bayesian information criterion (Kalyaanamoorthy et al., 2017) as implemented in IQ-TREE. Nodal support was evaluated with 1,000 ultrafast bootstrap replicates and 1,000 SH-like approximate likelihood-ratio test replicates (Anisimova and Gascuel, 2006; Hoang et al., 2018).

Population structure analyses

Principal component analysis was performed in PLINK v1.9.0-b.8 (Chang et al., 2015) using the final filtered one-SNP-per-RAD-locus dataset (Purcell et al., 2007). To ensure compatibility with PLINK, all SNPs were assigned to a single artificial chromosome and given unique identifiers based on their original scaffold and position. The VCF was converted to binary PLINK format using --make-bed and --double-id, after which the first 20 principal components were calculated. Eigenvectors and eigenvalues were imported into R, where the percentage of variance explained by each principal component was calculated. The first two principal components were visualized in R using ggplot2 and the ggrepel package for non-overlapping text labels, with samples colored by population and labeled by individual number (Wickham, 2016).

Ancestry clustering was performed with ADMIXTURE v1.3.0 (Alexander et al., 2009) using the final filtered one-SNP-per-RAD-locus dataset. Analyses were run for K values from 1 to 5, corresponding to the range from a single inferred ancestry component to the number of sampled populations. The preferred K was identified using cross-validation error, while higher K values were inspected as exploratory summaries of hierarchical structure rather than as independently confirmed demographic units.

A NeighborNet network was generated in SplitsTree v6.5.1 (Huson and Bryant, 2024) to provide a supplementary visualization of genetic relationships and potential conflicting signal among individuals. Pairwise genetic distances were calculated in SplitsTree using uncorrected p -distances, and the NeighborNet algorithm (Bryant and Moulton, 2004) was used to construct the network. The resulting network was interpreted as a distance-based visual summary of population structure.

Pairwise F ST and clustering based on genetic differentiation

Pairwise genetic differentiation among populations was quantified using Weir and Cockerham’s estimator of F ST (Weir and Cockerham, 1984) as implemented in vcftools. Calculations were based on the final filtered one-SNP-per-RAD-locus dataset. Population-specific sample lists were supplied to vcftools using the --weir-fst-pop option, and all pairwise comparisons were performed among the five analyzed populations. For each comparison, the weighted Weir and Cockerham F ST estimate was extracted from the output and arranged into a symmetric pairwise F ST matrix. To assess uncertainty around pairwise differentiation estimates, we additionally used the site-wise Weir and Cockerham F ST values reported by VCFtools and performed non-parametric bootstrap resampling across SNPs/loci. Because the final dataset contained one randomly selected SNP per RAD locus, SNP-based resampling approximated locus-level resampling. For each population pair, 10,000 bootstrap replicates were generated by sampling SNPs with replacement, and 95% confidence intervals were calculated from the 2.5% and 97.5% percentiles of the bootstrap distribution. To visualize relationships among populations based on genetic differentiation, the resulting pairwise F ST matrix was imported into PAST v4.17 and used as a user-supplied distance matrix for UPGMA clustering (Hammer et al., 2001). The resulting dendrogram was interpreted as a phenetic summary of relative genetic differentiation among populations, rather than as an explicit demographic or phylogenetic model.

Allelic richness and private alleles

To complement the population-structure and demographic analyses, we quantified rarefied allelic richness and private allelic variation from the filtered SNP dataset. The VCF file was imported into R v4.5.1 using the vcfR v1.13.0 package (Knaus and Grünwald, 2017) and converted to a genind object with adegenet v2.1.10 (Jombart, 2008). Rarefied allelic richness was estimated using the allelic.richness function of hierfstat v. 0.5-11 (Goudet, 2005) after conversion of the genind object to hierfstat format. Allelic richness was calculated for each locus and population, and mean allelic richness per population was obtained by averaging locus-specific estimates across loci while excluding missing values. This rarefaction-based approach allowed comparison of allelic richness among populations while accounting for differences in sample size (El Mousadik and Petit, 1996; Kalinowski, 2004). Uncertainty around mean rarefied allelic richness was assessed by non-parametric bootstrap resampling across loci. For each population and regional group, loci were sampled with replacement 10,000 times, the mean rarefied allelic richness was recalculated for each replicate, and 95% confidence intervals were calculated as the 2.5% and 97.5% percentiles of the bootstrap distribution.

Private alleles were quantified using the private_alleles function of poppr v. 2.9.6 (Kamvar et al., 2014). For each population, we calculated the total number of private alleles, the number of loci carrying at least one private allele and the mean number of private alleles per locus. Because raw private-allele counts are sensitive to sample size and missing data, interpretation emphasized rarefied allelic richness and mean private alleles per locus rather than total private-allele counts alone. Mean private alleles per locus were bootstrapped by resampling loci with replacement 10,000 times, and 95% confidence intervals were calculated from the 2.5% and 97.5% percentiles of the bootstrap distribution.

In addition to population-level analyses, we repeated the same rarefied allelic-richness and private-allele calculations at a broader, regional level. For this analysis, individuals were reassigned according to a regional population map, allowing comparison between the Balaton region and Lipce. The same workflow was then applied.

Historical effective population size reconstruction

Historical changes in effective population size were inferred using Stairway Plot v2.2 (Liu and Fu, 2015, Liu and Fu, 2020) from population-specific folded site-frequency spectra (SFS). Analyses were based on independently generated population-specific SFS datasets. Demographic inference was treated as exploratory because the SFSs were derived from reduced-representation RADseq data, were folded due to the absence of reliable ancestral-state information, and were based on SNPs retained after population-specific locus callability and SNP filtering. Specifically, we used the Stacks populations module v2.2 to generate population-specific VCF files from the final de novo assembly (denovo_map.pl run with m = 3, M = 5 and n = 5). For each demographic unit, the populations module was run to restrict the analysis only to non-missing variants (flags: -p 1 and -r 1.0) in all successfully genotyped individuals, using only individuals assigned to the focal population, thereby avoiding the requirement that loci be shared across multiple populations. Folded spectra were generated separately for four populations: Köveskál, Lesence, Nemeshany and Lipce. For each population, the number of chromosomes was set according to the number of successfully genotyped diploid individuals included in the corresponding SFS: 16 chromosomes for Köveskál, 36 for Lesence, 26 for Nemeshany and 16 for Lipce. Population-specific callable sequence lengths from the corresponding population-specific Stacks populations output were used as the L parameter in the Stairway Plot blueprint files: 3,692,520 bp for Köveskál, 2,367,401 bp for Lesence, 3,502,068 bp for Nemeshany and 3,686,762 bp for Lipce. Folded SFSs were supplied as integer counts of segregating sites, and whether_folded was set to true.

Folded site-frequency spectra were generated from these population-specific VCF files produced separately for each demographic unit. At each biallelic SNP, REF and ALT alleles were counted across diploid individuals. Sites with missing genotypes within the focal population were excluded so that each retained SNP contributed the full population-specific number of chromosomes to the SFS. Because ancestral states were unavailable, spectra were folded by assigning each SNP to the minor-allele-count class. Monomorphic sites within a population were not included in the polymorphic SFS line, whereas the total number of callable sites was supplied separately as the L parameter.

Stairway Plot analyses were run independently for each population. The proportion of sites used for training was set to 0.67, and 200 bootstrap input files were generated for each population. The number of random breakpoints was set to 4, 7, 11 and 14, following recommendations for reduced-representation-scale datasets (Liu and Fu, 2020). Demographic histories were scaled using a mutation rate of 5.28 × 10 –8 mutations per site per generation assessed for our focal species (Theodoridis et al., 2017) and a generation time of 2.5 years. This latter parameter was given as five years for P. farinosa by Theodoridis et al. (2017), but our extensive observations of the Hungarian populations suggest a shorter generation time: the first flowering usually takes place in the second, sometimes in the third year. Therefore, we set a generation time of 2.5 years for the analysis. Nevertheless, to evaluate the effect of generation-time uncertainty on the temporal scale of the inferred trajectories, Stairway Plot was also run with a generation time of 5 years as a sensitivity analysis. This analysis was used to assess whether the main qualitative population-specific trajectories were robust to generation-time assumptions, while recognizing that both calendar timing and absolute Ne estimates may vary with scaling parameters. The time axis was plotted in thousands of years before present. Because demographic inference was based on RADseq data and folded spectra, the resulting trajectories were interpreted as broad historical trends in effective population size rather than high-resolution reconstructions of demographic change.

RADseq assembly and filtered SNP dataset

Demultiplexing yielded 63.3 million reads across 48 samples, with a mean of 1.32 million reads per sample. Sample-level read numbers ranged from 206,467 to 2,150,312 reads, indicating generally high sequencing representation across individuals, although one sample (Kov07) had markedly lower read recovery (16% of the mean read count) and was removed from the dataset. The optimized Stacks assembly retained a large set of RAD loci suitable for population-genomic analysis. The final m3M5n5 assembly retained 66,614 loci from 121,578 assembled loci after applying the requirement that loci be present in at least four populations. Across these loci, 9,741,500 sites were represented, including 194,907 variable sites before downstream filtering. Population-level callable sequence lengths were broadly comparable among populations, ranging from 8,778,221 sites in Lipce to 9,612,655 sites in Nemeshany. After removing sites with more than 20% missing data, restricting the dataset to biallelic SNPs, filtering by mean read depth and retaining a single SNP per RAD locus, the final filtered one-SNP-per-RAD-locus dataset contained 32,628 high-quality variants that were used for all exploratory analyses of population structure, allelic diversity, but not for demographic history. For historical effective population size reconstruction, the starting final m3M5n5 dataset was filtered for 24,673 (Köveskál), 15,012 (Lesence), 20,228 (Lipce), and 14,249 (Nemeshany) biallelic SNPs that were present in all individuals within the population from 3,692,520 bp (Köveskál), 3,686,762 bp (Lipce), 2,367,401 bp (Lesence), and 3,502,068 bp (Nemeshany) callable sites.

Maximum-likelihood tree reconstruction

The maximum-likelihood tree revealed strong population-level structuring (Figure 2). When rooted with Lipce as the external comparative population, it was separated from all Hungarian populations by a long internal branch with maximal nodal support. Within the Hungarian samples, Nemeshany and Köveskál formed well-supported, coherent clusters, with only one Köveskál individual (#06) clustering with the Tapolca: Körtvélyes and Lesencetomaj: Sövényháza samples. These latter two populations showed close affinity and extensive intermixing, consistent with their geographic proximity in the Lesence Basin. The Lesence samples were sister to Nemeshany and Köveskál within the Hungarian populations.

Population structure

PCA ordination of allelic variation (Figure 3) supported the main relationships inferred from the phylogenetic tree. PC1 explained 19.3% of total variation and clearly separated Lipce from all other populations. PC2 explained 6.3% of total variation and differentiated Nemeshany from the remaining Hungarian populations. Subsequent axes each explained substantially less variation: PC3 explained 5.4%, PC4 5.1% and PC5 4.8%, and the remaining PCs each explained less than 5% of the variation. Köveskál occupied an intermediate but coherent position. In contrast, Tapolca and Lesencetomaj clustered closely together and showed limited intermixing. Again, individual #06 from Köveskál showed affinity to the Lesence group (i.e. Tapolca: Körtvélyes and Lesencetomaj: Sövényháza).

ADMIXTURE cross-validation identified K = 2 as the best-supported model, with the lowest cross-validation error (CV = 0.38879; Figure 4A). Under K = 2, individuals from Lipce were assigned almost entirely to one ancestry component, whereas the Hungarian populations were dominated by the alternative component (Figure 4B). Higher values of K revealed additional structure among the Hungarian populations, particularly the separation of Nemeshany and Köveskál, while Körtvélyes and Sövényháza remained closely related. However, these values were inspected as exploratory summaries of hierarchical population structure, but were not interpreted as independently confirmed historical ancestral populations or discrete demographic units.

A NeighborNet analysis provided a concordant visual summary of genetic relationships with a high fit value (99.6). This high fit indicates that the displayed split network captured nearly all variation in the pairwise distances among individuals. Lipnice individuals are clearly separated from all Hungarian samples by a long split, while the Hungarian samples formed a more compact, reticulate cluster, with Köveskál and Nemeshany showing differentiation, while the Lesence samples — Körtvélyes and Sövényháza — are more intermingled around the central part of the network (Figure 5).

Pairwise F ST

Pairwise F ST estimates confirmed the strong differentiation of Lipce and the close relationship between Tapolca and Lesencetomaj (Figure 6). The highest values (Table 2) involved Lipce: Nemeshany–Lipce (weighted F ST = 0.2617; mean site F ST = 0.159, 95% CI = 0.155–0.162), Köveskál–Lipce (weighted F ST = 0.2387; mean site F ST = 0.137, 95% CI = 0.134–0.140), Tapolca–Lipce (weighted F ST = 0.2200; mean site F ST = 0.124, 95% CI = 0.121–0.126) and Lesencetomaj–Lipce (weighted F ST = 0.2122; mean site FST = 0.115, 95% CI = 0.113–0.118). Differentiation among Hungarian populations was lower, with the smallest estimate observed between Tapolca and Lesencetomaj (weighted F ST = 0.0015; mean site F ST = -0.00144, 95% CI = -0.00240 to -0.000486). In contrast, Köveskál was weakly but consistently differentiated from both Lesence localities (mean site F ST = 0.0172, 95% CI = 0.0158–0.0185; mean site F ST = 0.0166, 95% CI = 0.0154–0.0178), whereas Nemeshany showed stronger differentiation from all other Hungarian populations. Bootstrap resampling of site-wise F ST estimates indicated narrow confidence intervals for all pairwise comparisons, showing that the relative differentiation patterns were robust to marker resampling (Table 2). As a phenetic summary of relative F ST differentiation, UPGMA dendrogram based on the pairwise F ST matrix grouped Tapolca and Lesencetomaj together, placed Köveskál close to this pair, and separated Nemeshany from the other Hungarian populations.

Allelic richness and private alleles

Rarefied allelic richness and private-allele analyses revealed clear differences in population-level allelic diversity (Table 3). Mean rarefied allelic richness was highest in Lipce (1.321; 95% CI: 1.32–1.33) and Lesence (1.302; 95% CI: 1.30–1.31), intermediate in Balaton (1.266; 95% CI: 1.26–1.27) and Köveskál (1.248; 95% CI: 1.24–1.25), and lowest in Nemeshany (1.178; 95% CI: 1.17–1.18). These results indicate that Lesence retains comparatively high allelic diversity among the Hungarian populations, while Köveskál and especially Nemeshany are genetically depauperate.

Private-allele summaries showed a similar but sampling-sensitive pattern. Nemeshany and Köveskál had few private alleles and low mean numbers of private alleles per locus (0.149 and 0.170, respectively). Lipce and Lesence showed higher rarefied allelic richness and more private variation than Nemeshany and Köveskál. At the regional level, Balaton had the highest total number of private alleles and the largest number of loci carrying private alleles; however, raw total private-allele counts were interpreted cautiously because they are affected by unequal sample size and missing-data structure.

To account for the evident dependence of private allele numbers on sample size, we calculated mean private alleles per population and regional groups. In this respect, Lipce showed the highest value (1.033; 95% CI: 1.005–1.062), followed by Lesence (0.795; 95% CI: 0.778–0.812). Köveskál and Nemeshany had substantially lower values, with overlapping or nearly overlapping confidence intervals (Köveskál: 0.170; 95% CI: 0.162–0.178; Nemeshany: 0.149; 95% CI: 0.140–0.159). Thus, private variation was concentrated mainly in Lipce and Lesence, whereas the two more depauperate Hungarian populations carried relatively little private variation per locus. At the regional level, the pooled Balaton Region group showed high private variation, but this comparison should be interpreted cautiously because the Balaton group contained substantially more individuals than Lipce (Table 3).

Historical effective population size trajectories

Stairway Plot analyses on independently generated population-specific folded SFSs suggested population-specific differences in historical effective population-size trajectories (Figure 7). Lesence showed the highest inferred Ne over much of the reconstructed time range, followed by a marked decline towards the present. Köveskál showed a broadly similar but less pronounced pattern, whereas Nemeshany displayed a strong reduction from formerly higher inferred Ne. Lipnice showed a comparatively lower and more stable trajectory. To assess the sensitivity of temporal interpretation to generation time, the analysis was repeated using a five-year generation time (Supplementary Figure 1). This sensitivity analysis recovered the same broad population-specific patterns, including the comparatively high inferred Ne of Lesence, the relative stability of Lipce, and the pronounced decline in Nemeshany, but shifted the inferred timing of demographic changes further into the past and altered the absolute Ne estimates. These differences confirm that the demographic reconstructions are sensitive to scaling assumptions. Nevertheless, all reconstructions showed reduced temporal resolution in older parts of the curves, approximately beyond the Late Glacial interval. These results should therefore be interpreted as broad, population-specific demographic signals rather than precise estimates of the timing or magnitude of historical events. The apparent loss of signal beyond approximately 15 kya was treated as a limitation of the information contained in the folded RADseq site-frequency spectra, not as evidence for a specific demographic event. We therefore interpret these reconstructions as exploratory and do not use them to infer the exact timing of population origin or isolation.

This study shows that the low-elevation populations of Primula farinosa sampled in Hungary and Slovenia contain geographically structured genomic diversity. The combined evidence from SNP-based tree inference, PCA, ancestry clustering, allelic richness and private allelic variation points to a hierarchically structured and regionally distinctive genomic pattern. Our analyses provide strong evidence for the regional-level population-genomic differentiation between the sampled P. farinosa populations. This result has direct conservation relevance because the Balaton Region populations are not genetically interchangeable with Lipce, and several Hungarian populations are also differentiated from one another. However, the results do not allow the exact timing of this differentiation to be resolved. In particular, the present RADseq dataset cannot distinguish unequivocally among uninterrupted local glacial persistence, older regional persistence with habitat shifts, remnants of an early-Holocene distribution, or more recent founder and drift processes. We therefore interpret the Balaton Region populations as a conservation-relevant component of regional genomic diversity, while treating long-term relict persistence as a plausible but not directly demonstrated historical scenario.

Evolutionary history of the Hungarian low-elevation populations

The Hungarian populations are best interpreted as components of the local Balaton Region gene pool of P. farinosa, whose present genomic structure has been shaped by geographic isolation, restricted gene flow, drift and relatively recent local demographic processes. Population structure analyses (Figures 2–5) place all Hungarian samples apart from Lipce, indicating that the Balaton Region populations are genetically differentiated from the Slovenian comparative population. This first-order split is important because Lipce is also a low-elevation population close to the Alps; therefore, the Hungarian populations are not simply genetically similar low-elevation satellites of the Slovenian Alpine-margin population. The Hungarian stands instead appear to preserve a separate regional genomic signature.

Within Hungary, the Lesence Basin populations, Körtvélyes and Sövényháza, are nearly indistinguishable in all structure analyses (Figures 2–5) and show almost no pairwise differentiation (Figure 6). Their treatment as a single biological population is therefore justified for broad population-genomic interpretation, even though management should still consider local habitat conditions separately. Lesence also has one of the highest rarefied allelic-richness values in the dataset, suggesting that it has avoided stronger recent bottlenecks affecting the smaller and more isolated Hungarian stands. This contrast between field observations of repeated apparent disappearance (Salamon-Albert and Morschhauser, 2003) and genomic diversity suggests that demographic census size and retained genomic variation may be partly decoupled. Lesence may have experienced recent demographic fluctuation or local disappearance–recolonization dynamics while still retaining a relatively diverse regional gene pool. By contrast, Nemeshany and Köveskál form recognizable genetic clusters (Figures 2–5) and show reduced allelic richness (Table 3). Together with the demographic reconstructions, which indicate declines in inferred effective population size in these populations (Figure 7; Supplementary Figure 1), this pattern is most consistent with demographic bottlenecks followed by strong genetic drift in small, isolated stands. This interpretation is especially plausible for Nemeshany, which combines the lowest rarefied allelic richness with a pronounced reduction in inferred Ne. Köveskál shows a less extreme but broadly similar pattern.

The position of Lesence within the Hungarian clade is especially relevant, but should be interpreted cautiously. When rooted with Lipce as the external comparative population, the SNP-based tree (Figure 2) placed Lesence (Körtvélyes and Sövényháza) as sister to the Nemeshany and Köveskál group within the Hungarian samples. This topology indicates that Lesence occupies a distinct position in the sampled regional genomic structure, whereas Nemeshany and Köveskál form a more closely related pair that have subsequently accumulated drift and lost allelic variation. While this should not be interpreted as direct evidence that Lesence gave rise to the other Hungarian populations, because a SNP tree is not an explicit demographic model, the interpretation is reinforced by the allelic-richness pattern. The allelic-richness pattern is nevertheless consistent with Lesence having experienced less severe genetic erosion than Nemeshany and Köveskál: Lesence was almost as rich as Lipce, whereas Nemeshany was the most depauperate population and Köveskál was only slightly richer (Table 3). Thus, Lesence may represent a relatively less eroded component of the Balaton Region gene pool, although alternative explanations, including differences in effective population size, drift, bottleneck severity or historical connectivity, cannot be excluded.

It is also notable that one Köveskál individual (#06) consistently groups with individuals from Lesence (Figures 2–4). That individual was sampled in the western part of the Köveskál population and consistently shows affinity to Lesence in independent analyses suggesting it is not an analytical noise. Although alternative scenarios, including natural gene flow, shared ancestry, incomplete lineage sorting, local genetic heterogeneity, or sample-label/genotyping issues, cannot be fully excluded, the documented sowing of Lesence propagules for conservation reasons at Köveskál provides an independent management context that makes anthropogenic introduction a plausible hypothesis for this pattern. Indeed, records from the Balaton-felvidéki National Park Directorate indicate repeated sowing of Lesence propagules at Köveskál since 2006, when the habitat began to suffer from drying. Despite these repeated introductions, the Lesence-like genomic assignment was detected in only a single individual from the western part of the population. This may indicate limited success of the sowing efforts and/or the persistence of the prevailing local genotype, although confirmation would require genetic analysis of the actual material used for sowing or of individuals known to derive from the introduced propagules. Given the short-term persistent seed bank reported for the species (Lindborg and Ehrlén, 2002), the limited genomic signal of sowing may also reflect restricted recruitment from introduced seeds, although this was not directly tested here.

Thus, the Hungarian populations combine evidence of regional distinctiveness, local demographic erosion and recent management-mediated genetic input, emphasizing that their present genomic structure likely reflects multiple historical and contemporary processes rather than a single recent event.

Relationship to the interglacial-expansion model in Primula farinosa

At the species-wide scale, P. farinosa has already been presented as a strong example of the interglacial-expansion model for cold-adapted plants (Theodoridis et al., 2017). This contrasts with many European high-mountain plants, for which molecular syntheses and comparative phylogeographic studies often point to survival in peripheral, nunatak or lowland refugia followed by postglacial range rearrangement (Schönswetter et al., 2005; Holderegger and Thiel-Egenter, 2009; Stewart et al., 2010). This scenario is biologically plausible because P. farinosa is not a strict high-alpine specialist: it is a moisture-demanding, weakly competitive plant of open calcareous habitats. During cold stages, many Central European lowland landscapes were characterized by cold steppe or steppe-tundra vegetation mosaics under strongly continental climatic conditions. Such landscapes may have contained extensive open habitats, but they were not necessarily suitable for a fen specialist whose persistence depends on continuously wet, calcareous and competition-limited microhabitats (Birks, 2008; Fletcher et al., 2010; Theodoridis et al., 2013, Theodoridis et al., 2017). Therefore, potential refugia for P. farinosa should be sought primarily in locally buffered wet habitats, such as calcareous fens, spring-fed wetlands, lake margins or other groundwater-influenced open sites, rather than in generalized steppe landscapes.

Our results add a finer regional layer to that picture. Pairwise F ST provides the most direct measure of the extent of genetic differentiation among the sampled populations. The narrow bootstrap confidence intervals of our analysis indicate that the observed hierarchy of genetic differentiation is not driven by stochastic marker sampling. Viewed against the broad comparative patterns summarized by Gamba and Muchhala (2020)—where outcrossing and temperate herb species generally show lower population differentiation, whereas small-insect pollination is associated with higher differentiation—the Lipce–Balaton Region comparisons indicate substantial regional differentiation within a conspecific, predominantly outcrossing temperate herb, with weighted F ST values of 0.2122–0.2617 (Figure 6). These values are far above the almost complete lack of differentiation between Körtvélyes and Sövényháza (weighted F ST = 0.0015), and they are also higher than the differentiation among most Hungarian populations. The contrast is informative. Very low F ST within Lesence is consistent with recent connectivity or a shared local gene pool, whereas high F ST between Lipce and the Balaton Region populations is compatible with prolonged restriction of gene flow (including drift, founder effects, bottlenecks and long-term small effective population size) and/or independent regional demographic histories.

While several microevolutionary processes can contribute to the observed pattern of strong regional differentiation, these results do not contradict the interglacial-expansion model. Rather, they suggest that the present-day Balaton Region occurrences are unlikely to be the product of very recent connectivity with the Alpine-margin low-elevation population. If such P. farinosa had spread widely across the region during the Holocene and only very recently contracted into scattered fen remnants, we would expect weaker differentiation between the Slovenian and Hungarian low-elevation populations, and probably more homogeneous ancestry among the Hungarian populations. Instead, the strong Slovenian–Hungarian split and the internal structure among Hungarian stands imply that low-elevation persistence and fragmentation operated over a longer or more complex timescale. Comparable cases in European wetland and subalpine plants show that strongly separated lineages may persist in geographically scattered refugial or peripheral habitats and can remain historically informative even when demographic timing is uncertain (Lienert et al., 2002; Urbaniak and Kwiatkowski, 2023). Thus, although the present results do not prove uninterrupted in situ survival at each sampled Hungarian site since the Late Glacial, they strongly support the conservation-genetic conclusion that the Balaton Region populations represent an isolated and conservation-relevant regional component of P. farinosa diversity. Their possible relictual origin remains a plausible historical interpretation, but it should be treated as a hypothesis requiring further testing rather than as a demonstrated result of the present RADseq dataset.

This interpretation fits the ecology of P. farinosa as a moisture-demanding, weakly competitive species of open calcareous habitats. In the Carpathian Basin, dry LGM vegetation and later Holocene forest expansion would both have constrained the species outside specialized wet calcareous microhabitats (Magyari et al., 2010, Magyari et al., 2014; Feurdean et al., 2014). Central European relict-rich fens provide a more relevant context than comparisons with strict high-mountain species. Hájek et al. (2011) showed that several fen plants and terrestrial snails display Holocene-scale relict distribution patterns, while Dítě et al. (2018) argued that many Western Carpathian glacial-relict candidates are linked to habitats that were regionally more widespread under colder or more open past landscapes, including peatlands and mineral-rich fens. Thus, any refugial interpretation of the Balaton Region populations should be understood in terms of locally buffered wetland refugia, such as calcareous fens, spring-fed wetlands or lake-margin habitats, rather than generalized high-mountain or steppe-refugial models (Stewart et al., 2010; Hájek et al., 2006, Hájek et al., 2011). The Balaton populations therefore make biogeographic sense as isolated low-elevation fen populations maintained by local edaphic and hydrological conditions, although their exact historical origin remains unresolved. Seen against this background, the Hungarian populations of P. farinosa fit a recognized Eurasian pattern: southern or Central-Eastern European low-elevation populations of cold-adapted taxa survive only where local habitat conditions buffer them from the surrounding warmer, more competitive landscape.

Demographic reconstruction and unresolved timing of isolation

One aim of this study was to test whether demographic reconstruction could provide an independent evidence for the relict status of the Balaton Region populations. However, the Stairway Plot analyses did not allow a reliable reconstruction of the timing of population origin or isolation. The generation-time sensitivity analysis further supports a cautious interpretation of the Stairway Plot results: the 2.5-year and 5-year analyses produced broadly similar relative patterns among populations, but differed in the calendar timing of inferred changes and in absolute Ne estimates. Although the population-specific curves suggest demographic heterogeneity among the sampled populations, the older segments of the trajectories, particularly beyond approximately 15 kya, were weakly resolved, probably reflecting the limited information content of folded site-frequency spectra derived from RADseq data. Demographic reconstruction methods based on linkage disequilibrium could provide a powerful alternative to SFS-based inference, but they require reliable genomic coordinates of markers. Because our de novo RADseq dataset lacked absolute genomic positions, such analyses were not feasible in the present study.

Therefore, these analyses should be interpreted only as broad indicators of population-specific demographic differences, not as precise reconstructions of Late Glacial or Holocene history. The present data cannot distinguish whether the Hungarian populations represent continuous local survivors, older regional relicts that shifted among suitable fen habitats, or remnants of an early Holocene lowland distribution that became isolated before modern habitat fragmentation. Nevertheless, when considered together with the strong regional differentiation, the distinct Balaton Region genomic signature and the reduced allelic richness of some Hungarian populations, the demographic results are compatible with potential long-term isolation, while not providing direct evidence for it.

Comparison with other low-elevation relicts and fragmented wetland specialists

The general phylogeographic pattern observed here resembles several other European cases in which lowland, peripheral or isolated populations of cold-adapted plants retain strong regional structure but differ in genetic diversity depending on isolation history and local demographic stability (Ægisdóttir et al., 2009; Reisch and Rosbakh, 2021). Comparisons with wet-habitat species at southern or fragmented parts of their ranges provide a more appropriate context than analogies with strict high-mountain specialists. These comparisons are relevant because, like P. farinosa in the Balaton Region, such species often persist in isolated wetland habitats where local hydrology, openness and low nutrient availability maintain suitable conditions despite an otherwise unsuitable regional landscape.

In Saxifraga hirculus L., a circumpolar species of wet habitats, European populations were strongly differentiated and had lower haplotype diversity than Alaskan populations, yet, chloroplast data did not identify clear European migration routes or refugia (Oliver et al., 2006). This is a useful cautionary parallel: strong geographic genetic structure can be well established even when the location and age of historical refugia remain unresolved. The cloudberry, Rubus chamaemorus L., provides another wet-habitat comparison. Phylogeographic work combined with pollen evidence showed that this peatland species could persist through dry glacial climates in restricted regions (Ehrich et al., 2008), while later analyses found strong differentiation between Czech and Norwegian populations and unique regional allele compositions in marginal populations (Leišová-Svobodová et al., 2018). In Ligularia sibirica Cass., isolated Central European wetland populations retained substantial within-population diversity and showed geographically structured differentiation, demonstrating that remnant wetland populations need not share a single simple diversity signature (Šmídová et al., 2011). For Swertia perennis L., another cold-adapted wetland plant with disjunct European lowland and montane occurrences, Lienert et al. (2002) reported that fen size and isolation affected isozyme variability. A later phylogeographic study of the species recovered several major Central European lineages and interpreted some north-eastern Polish and German populations as relics of formerly wider periglacial genetic groups (Urbaniak and Kwiatkowski, 2023).

Species-specific and congeneric studies provide useful context. Reisch et al. (2005) found significant altitudinal genetic differentiation among ten northern Alpine P. farinosa populations, with higher-altitude populations tending to retain greater within-population diversity; the authors attributed this pattern to phenological differences and vegetation-mediated restrictions on gene flow. A RAPD survey of four Lithuanian P. farinosa populations likewise found each sampled population genetically distinct (Augėnaitė et al., 2010). In north-eastern Switzerland, habitat fragmentation in P. farinosa itself reduced population size and altered performance-related ecological conditions in fen fragments (Lienert and Fischer, 2003). These studies make the present allelic-richness pattern more interpretable: the low diversity of Nemeshany and Köveskál is not unexpected for small, isolated fen populations, especially in an outcrossing and heterostylous species where demographic decline can quickly interact with mate limitation, pollinator limitation and genetic drift. More generally, they show that habitat configuration, isolation and restricted gene flow can be associated with substantial population structure in wetland-associated Primula populations, without necessarily providing a direct estimate of the age of that structure.

Allelic richness and private variation: conservation implications

Rarefied allelic richness is central to the conservation interpretation because it allows comparison among populations with different sample sizes. Lipce and Lesence have the highest values, suggesting that these populations retain more standing variation and may have experienced less severe recent genetic erosion. In contrast, Nemeshany and Köveskál showed lower allelic richness, with Nemeshany showing the lowest rarefied allelic richness and low private-allele values. However, given the modest population-level sample sizes, small differences among populations should not be overinterpreted. Low diversity can indicate recent founder events or bottlenecks, but in rear-edge and relict populations it can also reflect long-term persistence at small effective population size (Hampe and Petit, 2005). The elevated differentiation among several Balaton Region populations is compatible with this latter possibility, but, admittedly, it is not diagnostic of it, because similar patterns may also result from recent drift, founder events, bottlenecks or restricted gene flow. Broad comparative work on European plants also shows that isolated peripheral populations may combine reduced within-population diversity with stronger among-population differentiation, so low allelic richness should not be equated automatically with low historical or conservation value (Reisch and Rosbakh, 2021). Such populations may still contain unique allele combinations, locally adapted genotypes or population-specific components of genetic variation, even if their genome-wide diversity is reduced.

The private-allele results should be read cautiously. At the regional scale, Balaton has many private alleles and many loci carrying private alleles, but this total is affected by the larger number of samples pooled into the regional group. This limitation is especially relevant for the regional comparison, where the pooled Balaton Region sample contains substantially more individuals than Lipce. For this reason, we do not interpret raw private-allele totals as directly comparable estimates of regional uniqueness. Instead, we emphasize rarefied allelic richness as the main standardized diversity measure and use mean private alleles per locus, with bootstrap confidence intervals, as a supplementary description of private variation under the observed sampling design. At the population level, this analysis showed higher mean private alleles per locus in Lipce and Lesence than in Köveskál and Nemeshany, whereas raw regional private-allele totals were retained only as descriptive summaries. The stronger conservation point is therefore not the absolute number of private alleles, but the concordant evidence that the Balaton Region populations are genomically distinctive from Lipce, and several Hungarian populations are also distinguishable from one another (Table 3). For conservation, this argues against replacing one population with another or treating all low-elevation populations as a single exchangeable unit. Lesence may be valuable as a relatively allele-rich Hungarian reservoir, whereas Nemeshany and Köveskál may be valuable because they represent differentiated and potentially historically isolated fragments of the Balaton Region gene pool. Both types of value matter, because conservation based only on high diversity can overlook historically distinctive but genetically eroded populations.

Broader significance of southern low-elevation refugia

The Balaton Region populations belong to a broader biogeographic phenomenon: low-elevation occurrences of cold-adapted or high-mountain plants at southern latitudes can persist where local habitat conditions reproduce some of the ecological circumstances otherwise associated with colder regions or higher elevations. In Europe, such habitats include calcareous fens, spring systems, wet meadows, river gravels, cliffs and other open, competition-limited refugia (Birks, 2008; Hájek et al., 2006, Hájek et al., 2011; Horsák et al., 2015; Holderegger and Thiel-Egenter, 2009). Similar logic also applies elsewhere in temperate regions. Circumpolar and East Asian arctic-alpine studies show that disjunct mountain or peripheral populations can retain complex histories of isolation, connectivity and refugial persistence (Eidesen et al., 2013; Ikeda et al., 2016), while Lake Superior shoreline populations illustrate how cool, wet, exposed microhabitats can maintain southern disjunct arctic-alpine plant assemblages far from their main ranges (Hillman and Nielsen, 2024). These populations are often small, fragmented and genetically idiosyncratic, but that is precisely why they are biogeographically valuable.

In this wider context, Hungarian P. farinosa should not be viewed as a botanical curiosity sitting awkwardly outside the “real” mountain range of the species. Rather, its presence is compatible with a relict-refugial system in which low-elevation wetland habitats can preserve a genetic structure long after regional climate and vegetation have changed. The conservation implication is clear: maintaining multiple isolated populations is more important than simply protecting the largest or most allele-rich site. Lesence, Köveskál and Nemeshany each contribute differently to the regional genomic pattern. Lesence retains higher allelic richness, Nemeshany and Köveskál preserve differentiated but genetically eroded components, and the Balaton Region as a whole is strongly distinct from the Slovenian Alpine-margin low-elevation population. In light of the potential rear-edge nature of these populations (see Hampe and Petit, 2005), future conservation should therefore combine genetic management with practical conservation. In practical terms, this should be combined with continued maintenance of suitable habitat conditions, especially hydrology, open fen vegetation and low nutrient status. Additional non-genomic monitoring, including population size, recruitment and floral morph ratios in this distylous species, would be valuable future management considerations, although these parameters were not directly tested in the present study.

Taken together, our results indicate that the Balaton Region populations of P. farinosa represent a distinct and conservation-relevant component of regional genomic diversity. The strongest conclusion of this study is that these populations contain geographically structured and non-interchangeable genomic variation: they are strongly differentiated from Lipce, show internal structure within Hungary, and differ in the extent of allelic erosion. This conservation-genetic conclusion does not depend on demonstrating uninterrupted Glacial persistence. Although the present data do not resolve the exact timing of their isolation, they argue against treating them as recent, interchangeable remnants of a single homogeneous low-elevation population. Their conservation should therefore maintain the distinctiveness of individual populations while also addressing the local habitat processes that determine their persistence.