A total of 133 populations of P. distans agg. were sampled during 2018–2020. The main study area where the majority of the populations were sampled spanned Central and partly Eastern Europe. A few additional populations for comparison purposes were sampled in Southern Europe (see Supplementary Table 1 for a complete list of localities). Usually, 15 individuals per population were collected. Plants were collected at least 1 m apart to minimize repeated sampling of the same genotype. Herbarium vouchers are deposited in the herbarium WU.

For the purpose of statistical evaluation, the localities were categorized into (semi)natural and anthropogenic. As the distinction between purely natural and seminatural habitats was sometimes blurry, and most of the sites are indeed under long-term human management (such as pastures), we took a pragmatic approach and classified all habitats occurring in the regions with known historical occurrences of halophytes and having an extensive management as natural. Thus, the status of natural localities was assigned not only to clearly natural sites like shores of saline lakes, but also to saline meadows, pastures and similar habitats. Anthropogenic habitats included mainly roadsides (including agricultural roads) but also salt mines, salt processing sites, and dump sites.

Genome size of all collected plants was determined using flow cytometry. Only fresh tissues were used. The simplified two-step procedure of nuclear isolation and staining (Otto, 1990) modified for plant tissues (Doležel et al., 2007) was employed. Fluorescence intensity of 3500 particles was analyzed using a Partec CyFlow ML flow cytometer (Sysmex, Münster, Germany) equipped with a 532 nm (green) diode-pumped solid-state laser (100 mW output), employing propidium iodide (PI) as a fluorescent stain. Bellis perennis (2C = 3.38 pg; Temsch et al., 2021) was used as an internal standard. Each plant was measured individually.

The resulting histograms were evaluated using FloMax 2.52 (Sysmex, Münster, Germany), recording mean fluorescence and coefficient of variation for all fluorescent peaks. Only analyses with coefficients of variation below 2.5% were accepted for further statistical analyses as intra-cytotype variation was detected (see Results) and any suboptimal measurements would easily obscure the results. Genome size (GS, ratio of the mean fluorescence of the sample to the mean fluorescence of the internal standard multiplied by the genome size of the internal standard) was then calculated for each sample. Samples were clustered based on GS using the finite mixture Gaussian distribution modelling in R 4.1.2 (R Development Core Team, 2022). The FitGMM function from the package MGMM (McCaw, 2021) was used. Two values of the number of components (number of groups modelled, i.e., parameter k) were tested, namely eight and nine, and the maximum number of iterations (parameter maxit) was set to 10⁶. All other parameters were left at their default values.

To calibrate the FCM results, chromosomes were counted for 39 individuals representatively covering all the genome size groups discovered. Protocols of Pijnacker and Ferwerda (1984) and Belyayev et al. (2018) were used, with minor modifications. Fresh root tips of pot-cultivated plants were pretreated in ice cold water for 24 h and fixed in fresh fixative solution (pure ethanol and glacial acetic acid 3:1, v/v). The fixed material was stored in the fixative solution at -24°C. The fixed root tips were rinsed in double distilled water (ddH₂O; 2 × 5 min.) and citric buffer (10 mM sodium citrate, pH 4.8, 1 × 5 min.). Roots were subsequently treated in 0.3% (w/v) enzymatic solution [0.3% (w/v) cellulase, 0.3% (w/v) cytohelicase and 0.3% (w/v) pectolyase, Sigma St. Louis, MO, USA, in 10 mM citric buffer] and incubated in humid chamber for 1 h at 37°C. Digested roots were transferred into ddH₂O and kept at 4°C. Chromosome preparations were made by smear method (Pijnacker and Ferwerda, 1984). The root meristems were carefully transferred on clean slides and stirred by needle in 40 µL of 75% acetic acid on a warm plate (49°C) for 3 min., fixed in 300 µl fixative solution, washed in pure ethanol and air-dried. Multiple cells per specimen (at least three) were counted. Metaphase plates and chromosome counts were observed by Zeiss Axio Imager.Z2 microscope system at a magnification of 1000x. A few representative plates were photographed, and final figures were prepared using Adobe Photoshop version 21.1.3.

For the three main cytotypes (diploids, tetraploids, and hexaploids; see Results), we predicted their potential occurrences in unsampled regions of the main study area and identified any climate-driven limits of their distributional ranges using climatic niche modelling. The maximum entropy modelling approach implemented in Maxent 3.4.1 was employed (Phillips and Dudik, 2008) using the R package wallace (Kass et al., 2018). Nineteen bioclimatic variables from WorldClim (Fick and Hijmans, 2017) were downloaded in the highest available resolution (30 arc seconds ≈ 1 km²). Spearman correlation coefficients between all pairs of the bioclimatic variables were calculated, and only one variable was retained from each group of highly correlated (|r|>0.7) variables (i.e., bio01, bio02, bio03, bio04, bio06, and bio12). Duplicate data were removed by spatial thinning at a 1 km distance, and 10,000 background points were randomly sampled without replacement from the study area (the study area was defined as a minimum convex polygon around all the sampling sites in the main study area, expanded by a buffer zone of 2 degrees). The occurrences were partitioned into testing and training bins using the jackknife method. The maxent models were run with the regularization multiplier (RM) values ranging from 0.5 to 10 (incremented by 0.5) and five alternative settings for feature classes (i.e., types of functions used to fit the response curves; Merow et al., 2013) and their combinations (i.e., L, LQ, H, LQH, LQHP, where L = linear, Q = quadratic, H = hinge, and P = product; see Phillips et al., 2006; Phillips and Dudik, 2008). The Akaike information criterion (AIC; Bozdogan, 1987) was used to select the model with the highest predictive power.

In total, 1414 FCM samples were measured, 1153 of which met the CV cut-off criterion (see Supplementary Table 2 for an overview of all the measured samples). Raw flow cytometric data are available at the Zenodo digital repository (Kúr et al., 2023). Chromosome counting confirmed the presence of five distinct ploidy levels, i.e. diploid (2n = 14), tetraploid (2n = 28), pentaploid (2n = 35), hexaploid (2n = 42), and heptaploid (2n = 49) ( Figure 1 ). The finite mixture Gaussian distribution modelling of the GS for k = 8 found eight clear clusters, delimiting two clusters within the tetraploids and three clusters within the hexaploids ( Figure 2 , Table 1 ; see Supplementary Figure 1 for a chart of the raw C-values). Running the same procedure for k = 9 did not produce interpretable results as it distinguished one more group among the tetraploids which, however, was very small and entirely overlapped with the other fitted group (see Supplementary Figure 2 for the results of the finite mixture Gaussian modelling for k = 9). Hence, we used the clusters fitted for k = 8 for the definition of the GS groups used for further analyses. The GS groups were delimited based on strict criteria: (1) boundaries between overlapping heteroploid groups were defined by the individuals with highest/lowest GS that had counted chromosomes; (2) boundaries between overlapping homoploid groups were defined as the 1^st and 99^th percentiles of the corresponding Gauss curves. Individuals not falling within the defined ranges were left as unassigned. The result was one diploid (2x), two tetraploid (4xA, 4xB), one pentaploid (5x), two hexaploid (6xA, 6xB, 6xC), and one heptaploid (7x) groups ( Figure 2 ).

Thirty-eight percent of the populations were mixed, i.e., containing individuals from different GS groups ( Figure 3 ). Diploids (cytotype 2x) in the main study area occurred in a confined region in Central Germany. Additionally, two diploid populations were also found in N Italy ( Figure 4 ). Suitable climatic conditions for the diploids were predicted mostly in the area of their current distribution (in the main study area), with some minor additions in South-Western Germany and North-Eastern France ( Figure 5 ; the best maxent model L, RM = 1, AIC = 185.8, omission rate = 0.125 ± 0.354 [mean ± standard deviation]).

Tetraploids in the main study area were found nearly exclusively in the Pannonian Basin (with one outlying population in the Transylvanian Plateau in Romania) ( Figure 4 ). The 4xA cytotype clearly prevailed; the 4xB occurred only sporadically in mixed populations with the 4xA cytotype. Outside the main study area, plants with the GS fitting the 4xB cytotype were also found in N Italy. The predicted climatic suitability for the tetraploids highlights mostly the area of their current distribution in the main study area as having the highest occurrence probability ( Figure 5 , the best maxent model LQHP, RM = 1.5, AIC = 647.4, omission rate = 0.038 ± 0.196 [mean ± standard deviation]).

In contrast, hexaploids were ubiquitous through the main study area ( Figure 4 ). The most common cytotype was 6xC, which typically formed uniform populations ( Figure 3 ). The cytotype 6xA was relatively rare, confirmed in just 10 localities, and it mostly occurred in mixed populations with 6xB and 6xC. The cytotype 6xB was relatively abundant, again mostly in mixed populations with 6xC. Outside the main study area, plants with the GS fitting the 6xB and 6xC cytotypes were also found in Italy and France. Climatic modelling for the hexaploids showed most of the main study area as suitable, with lower occurrence probabilities only in mountainous regions ( Figure 5 , the best maxent model H, RM = 2.5, AIC = 2200.1, omission rate = 0.025 ± 0.158 [mean ± standard deviation]).

Pentaploids (cytotype 5x) were found exclusively in mixtures with other cytotypes, mostly tetraploids, in the Pannonian Basin. Similarly, heptaploids (cytotype 7x) were found exclusively in mixed populations, predominantly with hexaploid cytotypes, and they were distributed randomly across the whole study area ( Figures 3 , 4 ).

There was a clear difference in the cytotype frequency between populations from natural and anthropogenic habitats. A majority of the populations from anthropogenic habitats were hexaploid (with the cytotype 6xC clearly prevailing). Tetraploids were also rarely found in anthropogenic habitats (roadsides) in the areas of their natural occurrence (Eastern Hungary, Eastern Austria), as were heptaploids. Diploids and pentaploids were not found in anthropogenic habitats at all ( Figure 3 ).

Multiple homoploid groups differing in the genome size were found among the tetraploids and the hexaploids. Although the boundaries between the groups were not sharp, it is unlikely that the resulting pattern is an artifact as stringent filtering of the FCM analyses was used (only samples with coefficient of variation below 2.5% were kept). Clearly, two groups among the tetraploids (4xA and 4xB) and three groups among the hexaploids (6xA, 6xB, and 6xC) were present.

The minor (i.e., less frequent) cytotypes 4xB, 6xA, and 6xB were observed nearly exclusively in mixed populations with the respective major (i.e., dominant) cytotypes 4xA and 6xC, suggesting the minor cytotypes are evolutionarily derived from the major ones (the same seems to be true for the heptaploids which were observed nearly exclusively in the presence of the 6xC cytotype). The observed intraploidal GS variation could be caused by the presence of aneuploidy. Polyploid lineages, especially young ones, tend to produce unbalanced gametes, which may be viable and produce aneuploid offspring (Ramsey and Schemske, 2002). However, our data does not support this scenario in the case of the tetraploid and hexaploid P. distans agg. as we did not find any aneuploids in the chromosome counts of any of the 39 individuals counted, including the 15 individuals of the minor cytotypes ( Supplementary Table 2 ). This finding is congruent with the general lack of aneuploidy in the genus Puccinellia, with only P. maritima being the exception (Hughes, 1976; Rice et al., 2015). Thus, a more probable scenario is that the minor cytotypes are generated via the segregation of homeologous chromosomes of different sizes in the parental cytotypes. This would require the 4xA and 6xC cytotypes to be of allopolyploid origin, which may well be the case. This mechanism explaining intraspecific genome size variation was suggested for example in Festuca pallens Host and is most likely to occur in young rapidly radiating groups (Šmarda and Bureš, 2010).

The most probable pathway of the formation of the pentaploids is a cross between tetraploids and hexaploids as no aneuploids or triploids, which could otherwise enable their formation, were found. The pentaploids were indeed found exclusively in mixed populations with tetraploids and/or hexaploids, and they were discovered exclusively in regions with a common co-occurrence of both tetraploids and hexaploids (i.e., the Pannonian Basin), which further corroborates this scenario.

Contrary to the conclusions of Moravcová et al. (2001), who reported tetraploids of P. distans agg. only from Eastern Hungary, we discovered tetraploid cytotypes in an extensive area ranging from Eastern Austria and Southern Slovakia through Hungary to Romania. The most likely reason of the discrepancies with the literature is an unintended bias in the sampling of Moravcová et al. (2001) who missed the tetraploids’ localities in the western part of their study area.

On the other hand, our study did not discover any tetraploids in Central Germany as reported in Paule et al. (2017). There are three possible explanations of this result. First, there might be insufficient sampling in our study. A second possible explanation is that the tetraploid records from Central Germany are wrong. Third, a cytotype turnover (i.e., a replacement of one cytotype by another due to selective pressure), as observed for example by Trávníček et al. (2010) in Vicia cracca L., could also explain the observed discrepancy. While an unintentional oversight of tetraploids by us is possible, it is unlikely as our sampling in the area of alleged tetraploids’ occurrence was fairly thorough (18 localities), and we even had material from one exact locality where tetraploids were previously reported (Bad Salzungen). Paule et al. (2017) already suggested the tetraploid counts of P. distans agg. from Germany, all reported by Link (1992), to be incorrect. These counts were recorded in scope of a diploma thesis, and based on the provided photographic documentation it is difficult to interpret them as tetraploid. So, although we cannot fully exclude a historical presence of a tetraploid cytotype of P. distans agg. in Central Germany (a critical inspection of the published records, possibly including revisions of any existing herbarium material, would be useful), it seems currently more likely, given the absence of any other sources corroborating the findings of Link (1992), that the tetraploid P. distans agg. cytotypes have never been present in the German inland.

Diploids, in accord with the already known occurrences (Amarell and Gregor, 2021), were confirmed in a small area in Central Germany. Climatic niche modelling supported a general habitat unsuitability in other regions. Outside the main study area, diploid populations were also found in N Italy. Although in the latest edition of Flora d’Italia (Pignatti et al., 2017), P. distans is described as having diploid, tetraploid, pentaploid, and hexaploid cytotypes, we were unable to find any original publications reporting chromosome counts of P. distans agg. from Italy, and it is not clear whether the mentioned chromosome counts come from Italian populations at all. There are only a few published records of diploid Puccinellia from the European inland: Avdulow (1931) – European part of Russia, Montserrat and Montserrat (1988) – Spain, and Amarell and Gregor (2021) – Germany. Therefore, the diploids from Italy, found so far only in two close localities of mud volcanos, probably represent the first verifiable report of diploids from this country. The great cytotype diversity observed in the eight measured Italian localities alone also suggests that the patterns of diversity in P. distans agg. in Southern Europe may be as complex as in Central Europe, thus deserving a separate study.

The prevalence of the hexaploids (especially the 6xC cytotype) in anthropogenic habitats suggests their high propensity to adventive spread. Judging from their observed occurrences in localities far away from any natural habitats, it seems that the hexaploids are capable of an efficient long-range dispersal. It is therefore likely that a vast majority of literary reports of P. distans from road verges in Europe (Hetzel, 2006; Šerá, 2008; Šerá, 2010; Šerá, 2011) is attributable to the hexaploids. It is also likely that the hexaploid cytotypes have pervaded many natural habitats. We have detected three mixed-ploidy hexaploid-tetraploid populations and one hexaploid-diploid population in the main study area, 75% of which were in natural habitats. Although the hexaploids in the natural habitats might potentially represent local native genotypes, it is more probable that, at least in some cases, the widespread hexaploids represent an invading lineage in the early stage of its incursion into natural habitats. Alarmingly, indirect indications of the hybridization between the tetraploids and the hexaploids in the form of the presumably hybridogenous pentaploid cytotype exist.

The low expansion potential of the non-hexaploid cytotypes can most probably be attributed to their different life histories. Their lower vigor compared to the hexaploids could most probably manifest in lower drought and temperature resistance and/or lower seed set (Moura et al., 2021). Polyploidy has been indeed demonstrated as an important determinant of invasiveness/expansiveness in plants (te Beest et al., 2012; Walczyk, 2018; Moura et al., 2021). Our climatic models show a low habitat suitability in areas beyond the current distributions of the diploids and tetraploids. Although climatic models have their limitations (Oliveira et al., 2021) and other than climatic factors certainly play a role in the cytotypes’ distributions, it appears unlikely that these cytotypes can spread much outside their current distribution areas and therefore pose any threat to other native populations.