Hybridization may endanger the rare North Apennine endemic Cirsium bertolonii

We examined populations of North Apennine stenoendemics Cirsium bertolonii in the Apuan Alps and Tuscan-Emilian Apennines and found individuals morphologically shifted to co-occurring C. acaulon or C. erisithales. Hybrid status of these intermediates was confirmed by flow cytometry, morphometrics and amplified fragment length polymorphism (AFLP). We interpreted these hybrids taxonomically as Cirsium ×sagrense (C. acaulon × C. bertolonii) and C. ×abetonense (C. bertolonii × C. erisithales). Estimated genome size (2C) was 2244 ± 31 Mbp for C. ×sagrense and 2152 ± 99 Mbp for C. ×abetonense. Their genomic GC content was 38.95 ± 0.35% and 38.77 ± 0.26%, respectively. Diploid chromosome number 2n = 34 was counted for C. bertolonii, and the previously reported 2n = 12 needs to be considered erroneous. We found C. bertolonii to be a gynodioecious species (like many other Cirsium species) that is not reproductively isolated by ploidy level or chromosome number from co-occurring congeners. The relatively frequent occurrence of C. ×sagrense in Monte Sagro (locus classicus of C. bertolonii) and the prevalence of C. ×abetonense in Alpe Tre Potenze suggest that hybridization occurs repeatedly in C. bertolonii, as also confirmed by older herbarium specimens. C. ×abetonense produces ripe achenes and F2 hybrids or backcrosses, as some other Cirsium hybrids do. Genetic erosion/swamping via interspecific hybridization can therefore pose a risk to the genetic integrity of C. bertolonii, as it does for some other narrowly endemic high-mountain Cirsium species in Europe.


Introduction
Interspecific hybridization, i.e., reproduction between two or more species, is a common and ongoing process in populations of plants (Arnold 1992). The percentage of species that hybridize varies among families, but around 25% of plant species are known to create interspecific hybrids with at least one other species (Mallet 2007). Although the majority of hybridization events are evolutionary "dead ends" (Mayr 1942), as hybrids in most cases perish due to their reduced fitness in parental habitat (Barton and Hewitt 1985;Barton 2001;Burke and Arnold 2001), interspecific hybridization has been shown to have important evolutionary consequences (Arnold 1992;Mallet 2007;Soltis and Soltis 2009;Whitney et al. 2010;Abbott et al. 2013;Ellstrand and Rieseberg 2016;Alix et al. 2017;Goulet et al. 2017;Marques et al. 2019) and may even result in the origin of a new species (Rieseberg and Wendel 1993;Gross and Rieseberg 2005;Brochmann et al. 2000;James and Abbott 2005;White et al. 2018) or lead to species extinction (Levin et al. 1996;Rhymer and Simberloff 1996;Todesco et al. 2016).
Interspecific hybridization is particularly dangerous for rare stenoendemic species whose populations can easily get overwhelmed by large populations of widespread hybridizing congeners; disequilibrium in pollen may result in production of hybrid seeds at the expense of conspecific seeds of the rarer species (Arnold et al. 1993;Levin et al. 1996;Todesco et al. 2016). This is often the case for species that Handling Editor: Karol Marhold. are bound to small islands (Brochmann 1984;Rieseberg and Gerber 1995;Plume et al. 2013), to specific rarely occurring substrates or habitats (Balao et al. 2014;Vít et al. 2014) or to isolated mountain ranges (Krahulcová et al. 1996;Čertner et al. 2015;Gómez et al. 2015;Ma et al. 2019).
One of most narrowly distributed stenoendemics of the genus is Cirsium bertolonii Sprengel, a species occurring scarcely, in only two adjoining mountain ranges of the North Apennines, i.e., in the Apuan Alps and the Tuscan-Emilian Apennines separated by the valley of the Serchio River in the North Italian province of Tuscany (Caruel 1865; Archbald 1874; Gibelli and Pirotta 1882;Baroni 1897Baroni -1908Ferrarini 1984;Zanotti and Cristofolini 1985;Ferrarini 2000;Romagnoli and Foggi 2005;Pierini and Peruzzi 2014;Tomaselli et al. 2019). A remote population from Sasso Fratino (Bottacci et al. 2003) is considered an error (Viciani and Gonnelli 2014) and also rather doubtful is the record from Monti Sibillini (Zanotti and Cristofolini 1985;Conti et al. 2005;Pignatti 2018). Populations from Alpe della Luna, previously included in C. bertolonii, are now considered a different species, C. alpis-lunae Brilli-Catt. & Gubellini (Brilli-Cattarini and Gubellini 1991). Due to the similarity in spiny habitus and pale-yellow flowers, as well as similarity in ecological preferences, some authors consider C. bertolonii a subspecies or even a variety of the Alpine C. spinosissimum (L.) Scop. (cf. Fiori 1925(cf. Fiori -1929Werner ap. Heywood 1975;Werner 1976;Ferrarini 2000). Based on morphological and biochemical analyses, C. bertolonii has been recognized as a well-delimited species and was typified by Zanotti and Cristofolini (1985). This conclusion was also confirmed by our molecular analyses in which C. spinosissimum was shown to be less related to C. bertolonii than some other accepted European Cirsium species (Michálková et al. in prep.), so their taxonomical integration into one species is thus groundless. Despite the presence of over ten other Cirsium species in the Apuan Alps and the Tuscan-Emilian Apennines and their surroundings (Ferrarini 2000;Pierini and Peruzzi 2014), there are no documented records of interspecific hybrids of C. bertolonii, and its hybridization has only been speculated (Caruel 1865;Pierini and Peruzzi 2014). The chromosome number 2n = 12 detected for C. bertolonii by Garbari (1970) differs substantially from the number 2n = 34 prevailing in other European Cirsium species (Werner 1976; Bureš et al. 2004Bureš et al. , 2018. Although the count was later corrected to 2n = 34 (by Miceli and Garbari 1976), it could be speculated that C. bertolonii might be reproductively isolated from the rest of the genus, perhaps by chromosome arrangement or ploidy level, as is the case for two other rare high-mountain endemic thistles-C. appendiculatum Griseb. in the Balkanids and C. waldsteinii Rouy in the Carpathians-which are both tetraploid (2n = 68; Bureš et al. 2018).
Morphological identification of hybrids can be more reliable if corroborated by genetic methods. One such method is the whole-genome, dominant marker AFLP (amplified fragment length polymorphism), because it reveals the additive presence of nuclear genomic profiles from both parents in the examined hybrid individual. Nuclear profiles of potential parental species should, however, be based on samples from numerous distinct "pure" populations in which the second species is absent, in order to prevent including genetically eroded individuals (Balao et al. 2014;Michálková et al. 2018). In Cirsium, AFLP has been successfully applied to the identification of hybrids with high mountain stenoendemics C. rufescens (Segarra-Moragues et al. 2007), C. carniolicum, and C. greimleri . Flow cytometry is an additional "whole genome" marker for detection of interspecific hybrids and estimation of their ploidy level, genome size, and GC content (Bureš et al. 2003;Suda et al. 2007;Vít et al. 2014;Prančl et al. 2018;Popelka et al. 2019) that has been applied in Cirsium hybrid s' detection Michálková et al. 2018).
Using AFLP, morphometrics, and flow cytometry, we examined natural populations of Cirsium bertolonii with an emphasis on morphologically intermediate individuals, and confirmed or disproved their hybrid nature using reliable molecular markers, and thus considered the potential threat of the genetic erosion of this rare endemic mountain species. The potential reproductive isolation of C. bertolonii from co-occurring congeners mediated by different chromosome number and/or ploidy level was also analyzed in our study, using flow cytometry and chromosome counting.

Plant sampling strategy and sample processing
Populations of Cirsium bertolonii were examined in 2017 and 2018 in the surroundings of Monte Sagro in the Apuan Alps and Alpe Tre Potenze in the Tuscan-Emilian Apennines (Online Resource 1). For the purpose of AFLP and flow-cytometric analyses, we sampled two pieces of young upper undamaged leaf/plant from 18 individuals of C. bertolonii and 17 of putative hybrids. For the purpose of morphometric analyses, a lateral-view photo of the terminal capitulum was taken, 3-8 flowers with pappus were extirpated from the terminal capitulum using tweezers, and one fully developed middle leaf was sampled for each of the above-mentioned individuals (of C. bertolonii and putative hybrids); some other parameters (see below) were measured directly in the localities; mature achenes were also collected from some of the same individuals later in the season. Sampled plants were at least 10 m apart in order to avoid clones. Because C. bertolonii is very rare and we did not want to negatively influence its natural populations, the number of individuals kept as herbarium specimens was limited. Additionally, we only collected leaves and shoots, leaving the rhizome for the plant to regenerate. The specimens are kept in BRNU (acronyms for herbaria mentioned in this study follow Thiers 2018).
To confirm the hybrid nature of putative hybrids and to determine their parental species, the AFLP sample-set of C. was not included because it occurs in lowlands. Distantly related species of the genus Lophiolepis Cass. (previously included in Cirsium; cf. Ackerfield et al. 2020;Del Guacchio et al. 2022) occurring also in Tuscany were represented by L. eriophora (L.) Del Guacchio, Bureš, Iamonico & P.Caputo in our sample-set. Intergeneric hybridization between Lophiolepis and Cirsium is negligible (Wagenitz 1987;Bureš et al. 2004;Stace et al. 2015), so we did not expect Lophiolepis taxa to participate in hybridization with C. bertolonii. Samples of potential parental species were collected (as described above) in populations distant from the Apennines (Online Resource 1) to exclude any gene flow from C. bertolonii. The morphometrics and flow cytometry measurements were not intended to detect potential parents, but to confirm them. The flow cytometric sample-set of C. bertolonii and putative hybrids was therefore complemented by 9 samples of C. acaulon and 6 of C. erisithales; similarly, the morphometric dataset was enhanced by 40 and 40 samples of these species (Online Resource 1).
Samples for AFLP were placed in plastic bags and later stored in a deep freezer at − 80 °C until analyzed. Samples for flow cytometry were placed in plastic bags with a few droplets of water and analyzed after a few days (while still fresh). Morphometric samples-leaves and flowers-were put into paper or cellophane bags, respectively, and herbarized. The total sample-set thus contained 235 individuals of ten species and two hybrids for AFLP, 50 for flow cytometric analyses (3 spp. + 2 hybrids), and 129 for morphometric analyses (3 spp. + 2 hybrids), from which achenes of two individuals were subsequently used for chromosome counting of C. bertolonii; these samples originated from 106 localities in Austria, Bosnia & Herzegovina, Czech Republic, Italy, Slovakia, Slovenia, Spain, and Switzerland (Online Resource 1).

AFLP
DNA from deep-frozen leaves was extracted using the commercial kit NucleoSpin Plant II (MarcheryNagel) with the extraction buffer PL2, according to the manufacturer's instructions. The AFLP fingerprinting analysis follows Michálková et al. (2018) and was performed according to the protocol in the AFLPtm plant mapping kit (Applied Biosystems). The DNA (300 ng) was double-digested at 37 °C for 4 h with EcoRI and MseI enzymes and ligated to EcoRI and MseI adaptors in the same reaction. The preselective amplification was performed with EcoRI + A and MseI + C primers in a GenePro(Bioer) thermocycler. The following four selective primer pair combinations were chosen: 6-FAM-EcoRI + ACT/MseI + CTC, VIC-EcoRI + AGG/ MseI + CAG, NED-EcoRI + AGC/MseI + CAT, PET-EcoRI + ACA/MseI + CTA. The products of selective amplification were mixed with a GS500 LIZ size standard and Hi-DiTM Formamide (Applied Biosystems) for fragment analysis on an ABI 3500 Genetic Analyzer (Life Technologies). The error rate (according to Bonin et al. 2004) was estimated from 23 replicate samples. To perform sizing and scoring of the raw data, we used GeneMarker 2.4.0 (Soft-Genetics). A panel of scorable peaks (loci) for each primer combination was created manually, and only fragments from 50 to 550 bp were scored. Scoring of samples was first done automatically by the software and then checked and corrected manually. Two samples were excluded due to low amplification (Online Resource 1). The final file was then exported as a binary matrix.
Bayesian analysis implemented in STRU CTU RE 2.3.4 (Pritchard et al. 2000) was used to examine the nature ("pure" species or hybrid) of samples from Monte Sagro and Alpe Tre Potenze together with representatives of other potential Cirsium-parents. To confirm that there are 10 different species and that there is no divergence of C. bertolonii between the two regions (Apuan Alps and Tuscan-Emilian Apennines), different K values (K = 7 to K = 13) were tested. For each K, 15 separate runs were performed, each with 100,000 iterations and 100,000 burn-in using an admixture model with independent allele frequencies. The ploidy level was set as "2"-diploid-since all species (except for tetraploid C. vulgare) are diploid. STRU CTU RE HARVESTER v0.6.94 (Earl and von Holdt 2012) was used to summarize the output files and to construct the Evanno plots (Online Resource 2; Evanno et al. 2005). The STRU CTU RE diagram of computed assignment probabilities was generated using the R package pophelper (Francis 2017).
The Bayesian approach was complemented with distancebased analyses. Relationships between species occurring in Monte Sagro and Alpe Tre Potenze and the putative hybrids were visualized as a NeighborNet in SplitsTree software 4.14.4 (Huson and Bryant 2006) and as a PCoA plot in Canoco5 (ter Braak and Šmilauer 1998). Both analyses were based on Jaccard's binary distances.

Flow cytometry
The estimations of genome size and genomic GC content were performed on two CyFlow flow cytometers (Partec GmbH; recently Sysmex) using the internal standard Bellis perennis L. (2C = 3089.89 Mbp, 39.54% GC, Veselý et al. 2012), whose genome size and genomic GC content were derived from comparisons with the completely sequenced Oryza sativa subsp. japonica Shig.Kato 'Nipponbare' (International Rice Genome Sequencing Project 2005). The analyses follow the protocol of Šmarda et al. (2008). For the genome size and ploidy level estimations, we used propidium iodide as a fluorochrome, and for the genomic GC content estimation, the samples were co-processed with DAPI (4',6-diamidino-2-phenylindole) fluorochrome. To improve the signal/background ratio, the original OTTO I solution (Otto 1990) was mixed 1: 1 with 0.1 M hydrochloric acid and two drops of Tween. This "acid" buffer reduces the influence of secondary metabolites on the estimation. The genomic GC content was calculated using an ATGCFlow spreadsheet prepared by P. Šmarda: http:// www. sci. muni. cz/ botany/ syste mgr/ downl oad/ Festu ca/ ATGCF low. xls (Šmarda et al. 2008). 5000 nuclei were analyzed in one estimate. If any individual was measured repeatedly, the genome size and genomic GC content were averaged for subsequent statistical analyses based on individuals. (For original cytometric data, see Online Resources 3, 4).

Morphometrics
The selection of analyzed morphologic features was focused on traits in which the co-occurring, potentially parental species Cirsium acaulon and C. erisithales most differ from C. bertolonii, such as foliage, leaf and involucre spinosity, shape and division of leaves, and shape of involucre. We selected traits that can be measured either: (1) by simple techniques directly in the field; or (2) on photographs taken in the field; or (3) from minimally destructive samples of certain plant parts, so as not to damage the natural population of rare C. bertolonii. The directly measured values were then combined using simple mathematical formulas that best capture the nature of visually apparent differences among above the three above-mentioned species. Habitual features ( Fig. 1a)-lengths of uppermost three stem leaves (1Lf_L, 2Lf_L, and 3Lf_L) and distance between base of terminal capitulum and 3rd leaf-stem connection (3Lf_D)were measured on live plants or herbarium specimens using a plastic or paper ruler. Only foliage leaves were measured, not bracts (= modified leaves subtending the capitulum or its peduncle usually more simple in division than leaves-Br on Fig. 1a). Also, from these samples, the number of welldeveloped capitula (Cap_N) per plant (shoot) was counted (only well-developed capitula were counted, not dwarfed capitular buds) and the sex (female or hermaphrodite) of the plant was determined using a magnifying glass (Online Resource 5). Leaf characters (Fig. 1b, c)-length and width (Lf_L, Lf_W), longest lateral lobe length and width (Lob_L, Lob_W), length of terminal spine in the longest lateral lobe (Spi_L), and sinus depth of the longest lateral lobe (Sin_D)-were measured using a paper or plastic ruler on a herbarized, well-developed middle stem leaf (usually on the 3rd or 4th leaf counted down from the terminal capitulum) by scanning the leaf and ruler with the Department's printer and magnifying the resulting image on a large LCD screen. In the same manner, flower and fruit features (Fig. 1d, e)corolla length (Cor_L), pappus length (Pap_L), and achene length (Ach_L)-were measured on scans of flowers or fruits prepared from the terminal capitulum of live plants or herbarium specimens and put into transparent bags prior to scanning. The capitulum features (Fig. 1f)-involucre length and width (Inv_L, Inv_W), involucre width including bract tips (Inbr_W), capitulum width (Cap_W), and position and width of involucral neck (Nec_P, Nec_W)-were measured on lateral-view photos (with the ruler) of the live, fully flowering terminal capitulum. While measuring the leaf characters (as explained above) we also counted lobes-number of lateral lobes (main leaf lobes, excluding the terminal lobe; Lob_N) and average number of secondary lobes (SLob_N). For example, the leaf in Fig. 1c has 14 lateral lobes; it has an average of 0.29 secondary lobes because four lobes have one secondary lobe and the rest have none. ([4 × 1 + 10 × 0] /14 = 0.29 = SLob_N).
Because many of these features vary depending on environmental conditions, management, plant size, and developmental stage, we converted some of them into their less-variable relative forms, which better discriminate the potentially parental species and hybrids. From the habitual features ( Fig. 1c), we derived ( Fig. 1g) the sum of the lengths of the three uppermost leaves (123Lf_L) and the terminal foliage (Ter_Fol). From the features of leaves and their lobes (Fig. 1b, c), we derived ( Fig. 1h, j) the leaf shape (Lf_Sh), the leaf spinosity (Lf_Sp), the leaf division (Lf_Di), and the leaf lobe shape (Lob_Sh). From the features of the capitulum (Fig. 1f), we derived ( Fig. 1i) the involucre spinosity or the fuzziness (Inv_Sp), the involucre shape (Inv_Sh), and the involucre constriction (Inv_Con). For primary morphometric data of all 28 measured features, see Original Resources 6, 7.
To detect, measure, and eliminate the bias effect of sexual dimorphism (difference between hermaphrodites and females of the same taxon), the role of sex, taxon identity, and their interaction were tested using factorial ANOVA for all morphometric data (Online Resource 8-probabilities of the statistical tests are given after Bonferroni correction for multiple-comparison error-28 tests performed in total). The interspecific and/or intraspecific between-sex differences were visualized using boxplots accompanied by Tukey HSD (post hoc) tests for illustrative purposes (Online Resource 9). F-statistics ("determination value") of each morphological feature has been subsequently calculated using one-way ANOVA (Online Resource 10) in respect to both parental hybrid combinations. All the statistical tests of flow cytometric or morphometric data were calculated and respective graphs prepared in Statistica 13 (www. stats oft. com) and R version 3.6.1.
Canonical discriminant analysis (CDA) was used to visualize the morphological relationships among species and Leaf features (c) measured from the same leaf-length, width and number of lateral lobes (indicated by their serial numbers); 2nd, 4th, 11th, and 13th lobes have one secondary lobe while other ten have none, thus the average number of second-ary lobes is (4 × 1 + 10 × 0)/14 = 0.29 in this example leaf. Achene length (d) and pappus and corolla lengths (e) measured on scanned achenes and flowers prepared from live or herbarized plants and put into zipped plastic bags prior to scanning with the ruler. Capitulum features (f) measured on lateral-view photos (with the ruler) of live fully flowering capitula verify the intermediate position of the hybrids using the R-script Morphotools (Koutecký 2015). The standardized and transformed (to achieve normal distribution) morphological traits were used for the visualization.
Prior to the analysis, 7 characters were excluded. Five characters were excluded due to their strong correlation with other characters (Online Resource 11): pappus length (correlated with corolla length), distance between terminal capitulum and 3rd highest leaf (correlated with terminal foliage), leaf width (correlated with leaf lobe length), leaf spinosity (correlated with spine length), and leaf division (correlated with sinus depth). Capitulum width was excluded because neither of the parental species pairs significantly differed in this character. Achene length was excluded because ripe achenes were not available for some individuals.
DCA including the remaining 21 characters was based on individuals of Cirsium acaulon, C. bertolonii, and C. erisithales, and hybrid individuals were projected passively.

Morphologically intermediate individuals and co-occurring congeners in Cirsium bertolonii populations
We examined five populations of C. bertolonii from the two distinct regions-one large population in the Apuan Alps and four smaller populations in the Tuscan-Emilian Apennines. The population in the Apuan Alps was located on the SW slope of Monte Sagro on marble bedrock in an open, grassy habitat, which is typical for C. bertolonii (Zanotti and Cristofolini 1985). The population was relatively abundant-a few hundreds of individuals/flowering shoots scattered over a relatively small area (approximately 2 ha). All C. bertolonii individuals from this locality were very spiny with only mildly lobed leaves, and their morphological appearance (Fig. 2a,d,g) corresponded to the neotype which comes from this locality (Zanotti and Cristofolini 1985: 35). This population shared its location with another Cirsium species-C. acaulon (Fig. 2c, f, i)-whose abundance was comparable to that of C. bertolonii. Seven morphologically intermediate individuals were found here (Fig. 2b, e, h).
While searching for populations of C. bertolonii in the Tuscan-Emilian Apennines, namely on the slopes of Alpe Tre Potenze, we failed to find the species in open and grassy natural habitats similar to those in Apuan Alps. Nevertheless, we found four small populations in human-influenced habitats-along the roads and hiking trails. Some of these populations were in open forests in sites partly shaded by trees. Unlike Monte Sagro, where we could choose among hundreds of individuals of C. bertolonii, in Alpe Tre Potenze the number of individuals in populations was limited, and therefore, we had to include leaf samples of almost all nonclonal individuals. The individuals resembling C. bertolonii (Fig. 3b, e, h, Online Resource 12b) found in Alpe Tre Potenze were obviously morphologically different from those from Monte Sagro (Fig. 3a, d, g, Online Resource 12a). Namely, their leaves were more deeply lobed and less spiny, which suggested they may be influenced by more shaded habitats or may be genetically eroded. Within the region, there were abundant populations of C. erisithales (Fig. 3c, f, i, Online Resource 12c), and some of them also occurred in proximity to the micro-populations of these individuals.

Molecular verification of the hybrid nature of intermediate individuals and their parental affiliation
There were 604 AFLP polymorphic loci scored in 232 individuals from Cirsium bertolonii, co-occurring C. acaulon and C. erisithales, putative hybrids, and all other 7 potentially parental Cirsium species (occurring in Tuscany). The overall error rate was 0.4% for compared AFLP-replicates.
Within the STRU CTU RE analyses, the highest ΔK for 15 separate runs was yielded for K = 10 (Online Resource 2) confirming there are 10 distinct clusters representing 10 species. The STRU CTU RE analysis (Fig. 4) computed the assignment probabilities of respective individuals for each cluster/species (assignment values range from 0 = no assignment to 1 = 100% assignment) and demonstrated that representatives of morphologically "pure" (non-hybrid) individuals from localities other than Monte Sagro and Alpe Tre Potenze are also genetically "pure," because the assignment probabilities for their respective species clusters were above 0.9 (in accordance with Michálková et al. 2018) and the assignment to other species clusters was negligible. The only exception was one C. oleraceum individual from locali-tyA140 (Online Resource 1), whose assignment probability for the C. oleraceum cluster was slightly below the arbitrary limit of 0.9 (0.891; Fig. 4).
All morphologically "pure" individuals of C. bertolonii collected on Monte Sagro were confirmed to be also genetically "pure," as there was no eminent assignment to other species and the assignment probabilities for C. bertolonii clusters were over 0.9 (Fig. 4). Samples of morphological intermediates between C. acaulon and C. bertolonii from the same locality showed significant admixture of both suspected parents, while the admixture of other analyzed species was negligible (Fig. 4). From all Alpe Tre Potenze individuals, only two individuals found at the first locality-ATP1-had the assignment probability for C. bertolonii cluster over 0.9, and in one of those, there was a visible admixture of the C. erisithales cluster (Fig. 4). In the remaining individuals from the Alpe Tre Potenze, the assignment probabilities for the C. bertolonii cluster were substantially decreased due to considerable admixture of the C. erisithales cluster (Fig. 4), suggesting they are actually hybrids.
Distance-based analyses of the AFLP data-Neigh-borNet and PCoA-largely confirmed the STRU CTU RE results. The parental species-C. acaulon, C. bertolonii, and C. erisithales-formed three separated clusters (Fig. 5, Online Resource 13). Admixed individuals (morphological intermediates) from Monte Sagro were placed between the clusters of C. acaulon and C. bertolonii. Admixed individuals (morphological intermediates) from Alpe Tre Potenze were placed between the clusters of C. bertolonii (from Monte Sagro) and C. erisithales. Summarizing both analyses, we consider morphological intermediates from Monte Sagro to be hybrids between C. acaulon and C. bertolonii and the intermediates from Alpe Tre Potenze to be hybrids between C. bertolonii and C. erisithales. Both "pure" C. bertolonii individuals from locality ATP1 were slightly shifted from the Monte Sagro C. bertolonii toward C. erisithales in NeighborNet and PCoA (Fig. 5, Online Resource 13). This indicates the individuals from the ATP1 locality are introgressed rather than "pure" C. bertolonii. Fig. 2 Habitus of Cirsium acaulon × C. bertolonii hybrid (b, e, h) and of the parental species C. bertolonii (a, d, g) and C. acaulon (c, f, i). Flower head (a-c; bar = 1 cm); stem and its indumentum (d-f; bar = 1 cm); whole plant (g-i; bar = 5 cm)

Genome size, ploidy level, and genomic GC content in hybrids and parental species
Flow cytometric analyses yielded high-resolution histograms with distinct peaks of samples and internal standard in 33 analyzed individuals representing Cirsium acaulon, C. bertolonii, and C. erisithales and in 17 putatively hybrid individuals. The genome sizes of representatives of C. acaulon and C. bertolonii from Monte Sagro were significantly different from C. erisithales, but not from one another (Tukey HSD test; Fig. 6a). The average genome size of C. acaulon × C. bertolonii hybrids (2C = 2244 Mbp; standard deviation = 31 Mbp) was positioned between the parental genome sizes but did not statistically differ from them (Tukey HSD test; Fig. 6a). The DNA contents of C. bertolonii × C. erisithales individuals from Alpe Tre Potenze, excluding the two introgressed individuals, were 2C = 2144; s.d. = 109 Mbp and differed statistically from both parents (Tukey HSD test; Fig. 6a). The genome size of the two introgressed individuals from locality ATP1 was 2C = 2210 Mbp and 2C = 2162 Mbp. Such a low value was not detected among any of the genetically pure C. bertolonii individuals from Monte Sagro (Fig. 6a), which is consistent with Neigh-borNet and PCoA analyses (Fig. 5, Online Resource 13) and Within the measured genomic GC contents, only C. acaulon was statistically different from other analyzed taxa (Tukey HSD test; Fig. 6b), yet the GC contents of hybrids were positioned between the parental GC contents (Fig. 6b). The genomic GC content was 38.20%; s.d. = 0.30% for C. acaulon × C. bertolonii and 38.95%; s.d. = 0.35% for C. bertolonii × C. erisithales hybrids. GC contents of the two individuals from locality ATP1 were higher (mean = 39.3%) compared with other C. bertolonii × C. erisithales individuals (mean = 38.9%).
The chromosome numbers detected for C. bertolonii from the root tips of germinated achenes harvested from plants MS2a and MS2f from the locality Monte Sagro (Apuan Alps) were diploid 2n = 34 (Fig. 7).

Morphometrics
28 vegetative and generative morphological characters were measured in 114 individuals representing Cirsium acaulon, C. bertolonii, and C. erisithales and in 15 putatively hybrid individuals. In the Monte Sagro population, two sexesfemales and hermaphrodites-were identified among the plants of C. bertolonii, so the species should be considered gynodioecious. Therefore, we balanced morphometric sample-sets of C. acaulon and C. erisithales (both also gynodioecious) to contain not only the same number of individuals but also the same proportion of sexes. In three generative features-length of corolla, pappus, and achenes-betweensex differences were detected within the species or hybrids: Hermaphrodites have longer corollas and pappus than females (in the second case, only marginally supported), while females have larger achenes than hermaphrodites (Fig. 8a, b, Online Resources 8, 14a). Actually, the differences in the pappus or corolla length reliably distinguish all three parental species-C. acaulon (longest), C. bertolonii (medium), and C. erisithales (shortest)-regardless of the intraspecific differences between the sexes (Fig. 8a, b, Online Resource 8). However, in the case of hybrids, the differences between hybrids and parental species are not statistically significant unless one takes sex into account, i.e., only females with females or hermaphrodites with hermaphrodites should be compared (Fig. 8a, b, Online Resource 8).
For the determination/delimitation of hybrids the most valuable characters are involucre spinosity, involucre constriction, terminal foliage, distance from capitulum to 3rd highest leaf, number of capitula, and leaf spinosity in the case of C. acaulon × C. bertolonii (Table 1) and corolla showing Bayesian assignment probabilities for 10 clusters derived from AFLP profiles of Cirsium bertolonii and hybrids from Monte Sagro (MS) and Alpe Tre Potenze (ATP) and the remaining 9 species occurring in Tuscany. C. bertolonii individuals sampled at Monte Sagro show assignment to C. bertolonii (cluster assignment probability over 0.9). Hybrids from Monte Sagro show equal assignment to C. acaulon and C. bertolonii. Admixture of other species is negligible. From all individuals sampled in Alpe Tre Potenze, only the first two sampled at locality ATP1 have an assignment probability for the C. bertolonii cluster over 0.9, and the first of them has some admixture of C. erisithales. Remaining individuals from Alpe Tre Potenze are hybrids with equal assignment to both C. bertolonii and C. erisithales length, involucre spinosity, pappus length, length of spines, involucre width (incl. spines of bracts), leaf spinosity, and leaf division in the case of C. bertolonii × C. erisithales (Table 2) because these characters differ significantly not only between the parental species but also between the hybrid and each of its parents (Tables 1, 2).
Canonical discriminant analysis (CDA; Fig. 9) based on 21 of 28 analyzed features-7 characters excluded (see Materials and methods) prior to the analysis-well delimited the parental species. In the "morpho"-space the hybrids of C. acaulon × C. bertolonii were projected intermediately between the parents, as were hybrids of C. bertolonii × C. erisithales, where some individuals were slightly shifted toward C. bertolonii. Such a shift toward C. bertolonii was, however, not statistically significant in the introgressive individuals of C. bertolonii × C. erisithales from ATP1 (marked by red arrow).

Fig. 5
Genetic relationship among individuals collected at Monte Sagro (MS) and Alpe Tre Potenze (ATP) and representatives of cooccurring species visualized using NeighborNet based on Jaccard distances derived from AFLP profiles. Cirsium acaulon = purple circles, C. acaulon × C. bertolonii hybrids = green squares, C. bertolo-nii = yellow circles, C. bertolonii × C. erisithales hybrids = deep blue squares, and C. erisithales = light blue circles. Numbers on labels refer to respective localities (Online Resource 1). While MS1 and MS2 refer to an identical population sampled in two years, ATP1-ATP4 refer to different populations

Applicability of molecular, cytometric, and morphometric approaches in hybrid identification
Although AFLP is currently considered outdated and is slowly being taken over by NGS markers allowing the analyses of tens to hundreds of thousands of markers (Wang 2014), in our study, it undoubtedly demonstrated the hybrid nature and parental affiliation of samples of Cirsium acaulon × C. bertolonii and C. bertolonii × C. erisithaleswithin the STRU CTU RE analysis, hybrids contained admixture from both parental species (Fig. 4) and distance-based analyses positioned them between the respective parents (Fig. 5, Online Resource 13). The STRU CTU RE analysis worked well not only for the three parental species and the hybrids, but also when including 7 more closely related species occurring in Tuscany. The 10 well-delimited STRU CTU RE entities (K = 10, Online Resource 2) correspond to 10 species included in our study, which would indicate that, despite frequent interspecific hybridization in this genus, none of these species is of hybrid origin. Genome size estimation via flow cytometry indicated that the hybrids are diploid, which is the most common ploidy level in the genus (detected in nearly 80% of Cirsium species in Europe and Asia- Bureš et al. 2004Bureš et al. , 2018. Flow cytometry was not a suitable tool for detecting the parental affiliation of either of the two hybrids (when statistical significance > 0.05 was taken into account). This is primarily given by the small differences between the hybrids and their respective parents due to low variation in genome size across the genus Michálková et al. 2018), which is within an order of magnitude of the error of individual cytometric estimation. Although the genome sizes and GC contents of hybrids were positioned more or less between the respective parental values, the difference was not always significant (Tukey HSD test; Fig. 6a, b). This is a consequence of the limited number of analyzed individuals and would likely be overcome by a larger number of samples. Fig. 8 Variability in pappus (a) and corolla length (b) detected in Cirsium acaulon × C. bertolonii and C. bertolonii × C. erisithales hybrids and parental species C. acaulon, C. bertolonii, and C. erisithales with respect to sex (females and hermaphrodites). The length was measured on flowers extirpated from terminal capitulum. Boxplots show median (thick line), interquartile range (boxes), non-outlier range (whiskers), and outliers (empty circles). Boxplots marked with the same letter do not differ significantly at p > 0.05 (Tukey HSD test) Morphometric analysis detected some characters well delimiting respective parental-species pairs and their hybrids (first five/seven characters in Tables 1 or 2, respectively). Although the hybrids were usually positioned between their respective parents, they were usually shifted toward one of the parents (second and third sections in Tables 1 and 2) or occasionally the hybrids were even out of the parental range (e.g., C. acaulon × C. bertolonii has a smaller number of secondary leaf lobes and more spiny leaves than its parents-Online Resource 14n, x-and C. bertolonii × C. erisithales has narrower involucres, a larger number of capitula, and a shorter distance from capitulum to 3rd stem leaf than its parents-Online Resource 14c, h, i). The achenes of C. acaulon × C. bertolonii were also smaller than those of its parents. However, they were not ripe and therefore cannot be considered out of the parental range. This situation well illustrates that hybrid plants are more frequently a mosaicism of parental features rather than their average.

Gynodioecy of Cirsium bertolonii and fertility of its hybrids
The hybrids of C. acaulon × C. bertolonii, as well as two of the nine individuals of C. bertolonii × C. erisithales (BERT × ERIS-ATP2c, BERT × ERIS-ATP4a), did not contain developed pollen in their anther tubes. One could easily interpret this as hybrid sterility. However, these individuals were actually females whose anther sterility has been inherited from the parental species-C. acaulon, C. bertolonii and C. erisithales-which are all gynodioecious, Table 1 Evaluation of morphological characters of hybrid between Cirsium acaulon and C. bertolonii using F-statistics against parental species (= among-groups variation/within-groups variation) and post hoc Tukey in one-way ANOVA *Hybrid plants did not contain ripe achenes. For explanation of the characters see Fig. 1. and the chapter "Material and methods"  i.e., containing hermaphrodites and females (= individuals with hermaphrodite flowers but sterile rudimentary anther tubes (synantherium)-Online Resource 5) regularly in their populations. A similar co-occurrence of females and hermaphrodites (with viable pollen) was previously detected in hybrids of C. acaulon (with C. canum, C. oleraceum, or C. pannonicum) and in hybrids of C. erisithales (with C. heterophyllum, C. oleraceum, or C. palustre) by Bureš et al. (2010). The gynodioecious sexual dimorphism recognized in 275 (~ 2%) of angiosperm genera (Dufay et al. 2014) has so far been found in more than 20 European Cirsium species (Delannay 1978;Bureš et al. 2010). In our study, the gynodioecy of C. bertolonii is reported for the first time.
Besides the sterility/fertility of anther tubes, the conspecific females and hermaphrodites in Cirsium can also differ in other generative features  as also demonstrated in our study for all three analyzed species (C. acaulon, C. bertolonii, and C. erisithales), whose hermaphrodites have longer corollas and pappus but smaller achenes than females ( Fig. 8; Online Resource 14a). Although the larger flowers of hermaphrodites in most gynodioecious species have been well known to evolutionary biologists since the time of Charles Darwin (1877), this fact is often ignored in taxonomic studies and determination keys. In Cirsium, especially among closely related species, the differences between the females and hermaphrodites within the species could even be higher than between the species, as demonstrated for closely related C. greimleri and C. waldsteinii by Bureš et al. (2018). This is, fortunately, not the case in our recent study, where the differences among C. acaulon, C. bertolonii, and C. erisithales are rather larger than between sexes within a species (Fig. 8; Online Resource 14a), but for hybrid delimitation against the parental species, the comparison of females with females or hermaphrodites with hermaphrodites is more reliable (Fig. 8).
The analyzed hybrids of C. bertolonii with C. acaulon at Monte Sagro did not contain ripe achenes which one could interpret as a hybrid sterility. However, in Cirsium, the absence of ripe fruits or low fruit fertility is relatively common even in regular non-hybrid plants (Bureš, pers. observ.). Therefore, we do not consider undeveloped achenes in C. acaulon × C. bertolonii to be conclusive evidence of a reproductive barrier between these species via hybrid sterility. Clarifying this issue would require experimental crossing and subsequent cultivation of hybrids. Fertility of C. acaulon hybrids with other Cirsium species has previously been confirmed, e.g., by Pigott (1968), or by Bureš et al. (2010). In C. bertolonii × C. erisithales hybrids from Alpe Tre Potenze, well-developed achenes were detected, suggesting the hybrids might be fertile, as was demonstrated for some other hybrids of C. erisithales (Bureš et al. 2010).
The hybrid samples from the ATP1 locality in the Alpe Tre Potenze were situated outside the cluster of hybrids C. bertolonii × C. erisithales toward C. bertolonii in STRU CTU RE (Fig. 4) and in distance-based analyses ( Fig. 5; Online Resource 13). A similar shift of these samples toward C. bertolonii can be seen in genome size and genomic GC content (Fig. 6a, b). In morphometric analysis, this shift is more or less present in the majority of characters, except for three characters related to leaf lobes, where the samples from ATP1 were strongly shifted toward C. erisithales or even exceeded it. As a result of these three exceptions, CDA (Fig. 9) did not show the shift to C. bertolonii that was observable in other analyses.
One could speculate that this shift toward C. bertolonii indicates backcrossing (of C. bertolonii × C. erisithales Fig. 9 Spider web projection of canonical discriminant analysis (CDA) of Cirsium acaulon = purple circles, C. acaulon × C. bertolonii hybrids = green squares, C. bertolonii = yellow circles, C. bertolonii × C. erisithales hybrids = deep blue squares, and C. erisithales = light blue circles based on 21 morphological characters. Lines connect respective individuals with the centroids of their taxa. The BC or F2 hybrid (as identified by AFLP) from ATP1 is marked with a red arrow with C. bertolonii), and however, it could also be F2 offspring resulting from random segregation of chromosomal homologs during gametogenesis and the subsequent selfing of a homoploid F1 hybrid (or outcrossing between F1 hybrids prevailing in the analyzed populations in Alpe Tre Potenze). In any case, the occurrence of such genetically shifted individuals supports the premise that primary hybrids (F1) of C. bertolonii × C. erisithales may be at least partially fertile as documented for other Cirsium hybrids (e.g., by Correns 1916;Bureš et al. 2010;Michálková et al. 2018).

The possible threat of genetic erosion of Cirsium bertolonii
Somatic DNA content estimated for C. bertolonii in our study is very similar to DNA contents of other diploid Cirsium species (Bureš et al. , 2023Michálková et al. 2018), and chromosome number 2n = 34 based on achenes sampled from the individuals in the type locality Monte Sagro also undoubtedly confirmed the diploid status of C. bertolonii, as previously suggested by Miceli and Garbari (1976). Therefore, the previously reported 2n = 12 (Garbari 1970) needs to be considered erroneous and neither different ploidy level nor some kind of chromosomal rearrangement protects C. bertolonii against hybridization with other frequently hybridizing and predominantly diploid species of Cirsium, in contrast to the cases of the rare mountain endemic tetraploid (2n = 68) C. appendiculatum in the Balkan Mountains or of C. waldsteinii in the Carpathians . Our study also confirmed that C. bertolonii really can produce (at least partially) fertile hybrids with two relatively more abundant and more widespread diploid congeners: C. acaulon and C. erisithales. C. acaulon occurs almost continually in Western and Middle Europe-from England, France, and Spain to Southern Sweden, the Baltics, Poland, Czech Republic, Austria, Italy, Slovenia, Croatia, Bosnia and Herzegovina, with dispersed occurrences also in Slovakia, Hungary, and Romania; C. erisithales grows in most of the higher European mountains (or their parts) such as in the Easternmost Pyrenees, Massif Central, the Alps, the Northern and Central Apennines, the Dinarides, and the Carpathians (Meusel and Jäger 1992). Across their natural ranges, both species, C. acaulon and C. erisithales, willingly hybridize with most of their sympatric congeners (Pignatti 2018;Talavera 2017;Stace et al. 2015;Bureš 2004;Wagenitz 1987;Nyárády 1964). Therefore, it is not surprising that they hybridize when co-occurring with the rare C. bertolonii, with whom they share an ecological preference for calcareous substrates (Ferrarini 1967;Gabellini et al. 2006;Tomaselli et al. 2019).
The higher prevalence of hybrids between C. bertolonii and C. erisithales may also be due to overlap in their flowering times during late spring-early summer, whereas C. acaulon blossoms later-in summer-early autumn, depending on altitude (Michálková and Bureš pers. observ.). This phenological difference of C. acaulon was notable during our repeated visits to Monte Sagro-in July, C. bertolonii was in pollination optimum, whereas C. acaulon had only developed inflorescence buds and the hybrid already had capitula, but was not flowering yet; in October, the shoots of C. bertolonii were almost faded, whereas some of the hybrids were still blossoming, as was C. acaulon.
The scattered distribution of hybrids C. acaulon × C. bertolonii in Monte Sagro and the exclusive presence of C. bertolonii × C. erisithales in the four studied populations in Alpe Tre Potenze suggest that the hybridization of C. bertolonii with both these species is repeatedly ongoing, which is further confirmed by the above-mentioned herbarium specimens and literary report, as well as by an older herbarium specimen of C. acaulon × C. bertolonii collected probably by A. Bertoloni in Catino dell Sagro on the northern slope of Monte Sagro at the beginning of the nineteenth century and later deposited in the Viennese Natural History Museum (Bertoloni s. a. W 1929(Bertoloni s. a. W /1457Online Resource 15). From a taxonomic point of view, it is somewhat worrying that hybridization with C. acaulon takes place directly on the locus classicus of C. bertolonii-within the Monte Sagro population. However, based on Bertoloni's "Amoenitates," we can conclude that both species have coexisted there for at least 200 years without either going extinct (Bertoloni 1819: 404-405). The neotype from this locality selected by Zanotti and Cristofolini (1985: 35; http:// parla tore. msn. unifi. it/ img72/ FI001 979. jpg) does not show any sign of hybridization.
Although the situation in Monte Sagro illustrates a relatively long coexistence of C. bertolonii with C. acaulon, repeatedly ongoing hybridization, fertility of hybrids C. bertolonii with C. erisithales, and overlap in blossoming and ecological preferences could pose a threat to the genetic integrity of the rare mountain North Apennine endemic C. bertolonii, particularly by genetic erosion from C. erisithales in Alpe Tre Potenze or in the whole Tuscan-Emilian Apennines-as is the case with other rare, diploid, European, narrowly distributed high-mountain thistles, such as C. alsophilum, C. greimleri, C. carniolicum, and C. rufescens Michálková et al. 2018).
Etymology: The epithet is derived from the town Abetone in proximity of Alpe Tre Potenze, where the hybrids between C. bertolonii and C. erisithales were found.