Integrative Taxonomy of Armeria arenaria (Plumbaginaceae), with a Special Focus on the Putative Subspecies Endemic to the Apennines

Simple Summary Armeria arenaria is a highly variable Western European species, for which three subspecies are recorded in Italy. Armeria arenaria subsp. arenaria has been reported from Northern Italy, while A. arenaria subsp. marginata and A. arenaria subsp. apennina are considered endemic to the Apennines. The taxonomic value of these two latter taxa is unclear and the actual occurrence of A. arenaria s.str. in Italy has never been addressed. Following an integrated taxonomic approach, in this study we show that all the Italian records of A. arenaria s.str. should be actually referred to A. arenaria subsp. praecox and that only one Northern Apennine endemic taxon can be recognized, namely, A. arenaria subsp. marginata. Abstract Three subspecies of Armeria arenaria are reported from Italy, two of which are considered endemic to the Apennines. The taxonomic value of these two taxa (A. arenaria subsp. marginata and A. arenaria subsp. apennina) is unclear. Moreover, the relationships between A. arenaria subsp. praecox and Northern Italian populations—currently ascribed to A. arenaria subsp. arenaria—have never been addressed. Accordingly, we used an integrated taxonomic approach, including morphometry, seed morpho–colorimetry, karyology, molecular systematics (psbA–trnH, trnQ–rps16, trnF–trnL, trnL–rpl32, and ITS region), and comparative niche analysis. According to our results, French–Northern Italian populations are clearly distinct from Apennine populations. In the first group, there is evidence which allows the recognition of A. arenaria s.str. (not occurring in Italy) and A. arenaria subsp. praecox. In the second group, the two putative taxa endemic to the Northern Apennines cannot be separated, so a single subspecies is here recognized: A. arenaria subsp. marginata.


Introduction
Most of our biological knowledge of plant diversity comes from the foundations laid by alpha taxonomy, which played a crucial role in discovering and documenting plant diversity around the world. Nevertheless, although science is progressing, taxonomists seem to struggle to keep pace with novel methods and approaches. Indeed, hundreds of putative new species are described annually [1], but most of them are still described on qualitative grounds. In such approaches, the information that a taxonomist collects to shape his/her idea about the "species" in question is often obscure [2], so that biases [3] in taxonomist decisions [4] can dramatically affect taxonomic treatment. For instance, the number of species in the genus Armeria varies dramatically under different taxonomic circumscriptions elaborated by different taxonomists [5][6][7]. The subjectiveness of these processes may have contributed to what has been called taxonomic anarchy [8]. Integrated taxonomic approaches aim to address this problem with the consilience principle [2], according to which multiple and complementary approaches (morphology, phylogenetics, cytology, etc.) [9] are used to try falsifying taxonomic hypotheses in a Popperian sense [10]. This represents a step towards an omega taxonomy [11,12] that needs the integration of different skills [13,14].
The genus Armeria Willd. (Plumbaginaceae, Limonioideae) includes up to 95 accepted, mostly Holarctic, perennial species [15]. In Italy, the current knowledge on the taxonomy and systematics of this genus is largely derived from traditional alpha-taxonomic revisions [16,17], which indicate the existence of 23 taxa in 18 species [18]. However, the taxonomic value of some of these taxa is still debated [18], and the picture is further complicated by the fact that species boundaries within Armeria are difficult to establish [6,19] and weak [20][21][22]. In this scenario, it has been demonstrated that homoploid hybrid speciation [21] can play a crucial role in the emergence of new species [23,24], given that all species tested so far are diploid, with 2n = 2x = 18 chromosomes [16,25,26]. The use of nrDNA and (maternally inherited) cpDNA markers helps to elucidate the phylogenetic relationships even under hybridization scenarios [27,28]. Armeria arenaria (Pers.) F. Dietr. complex currently includes 13 subspecies in its whole range [29] and, according to Arrigoni [17], three subspecies occur in Italy: A. arenaria subsp. arenaria, distributed across the Central-Western Alps [18]; A. arenaria subsp. apennina Arrigoni, endemic to the Tuscan-Aemilian Apennines; and A. arenaria subsp. marginata (Levier) Arrigoni, also endemic to the Northern and up to the Central Apennines. Armeria arenaria subsp. praecox (Jord.) Kerguélen ex Greuter, Burdet & G. Long, described from south-eastern France, is reported as doubtfully occurring in Italy. Arrigoni [17] considers A. arenaria subsp. apennina as intermediate between A. arenaria s.str. and A. arenaria subsp. marginata. The same author [17] also claims that there is a series of unclear intermediate forms distinguished by the transition of some putatively diagnostic character states. However, the circumscription of these subspecies is based only on a qualitative morphological approach. All these factors led to the consideration of A. arenaria subsp. apennina and A. arenaria subsp. marginata as two subspecies of uncertain taxonomic value [18].
For these reasons, there is need to use an integrated approach to address the taxonomy [9] of these putative subspecies. To achieve a sound taxonomic circumscription, we performed morphometric analyses, including living populations from type localities, complemented by seed morpho-colorimetry, karyotype asymmetry estimation, molecular systematics, and comparative niche analysis (for similar integrative approaches, see [27,28]). In this study we aim: (1) to test the current taxonomic circumscription; (2) to verify the occurrence in Italy of A. arenaria subsp. praecox; and (3) to clarify the nomenclature of the group.

Sampling
In total, we selected 12 populations (Table 1)  (1) to include all the type localities of the four taxa putatively occurring in Italy (FB, LA, LL, and MB-acronyms as in Table 1); (2) to include other populations explicitly cited in [19]: AA, BO, BR, MC, MP, and TV; and (3) to also include a lowland (GA) and the easternmost (PS) populations in Italy. Table 1. Taxa and populations of Armeria arenaria sampled in this study, according to the current taxonomic hypothesis [17]. "Code" corresponds to population acronyms used elsewhere in the manuscript. * = type locality. "Voucher" refers to the specimens stored at Herbarium Horti Botanici Pisani (PI) and freely available for consultation at http://erbario.unipi.it/, accessed on 10 July 2022. See also Figure 6 for the geographical localisation of the sampled populations.

Code
Current Taxonomic  For each population, about 20 flowering individuals were sampled. The number of flowering scapes was counted in the field, whereas pictures were taken to assess the colour of the flowers and involucres of each plant. In total, 229 specimens were collected, and herbarium vouchers were prepared. All vouchers are stored at Herbarium Horti Botanici Pisani (PI), and high-resolution images are freely available for consultation at http://erbario.unipi.it/, accessed on 10 July 2022 (codes in Table 1). Concerning molecular systematics, dried leaves were picked from a subset of three individuals for each population and put in a paper bag with silica gel. Ripe fruits were also collected from the same populations. Seeds were dried at room temperature for two months and cleaned in the Germplasm Bank, Department of Biology of the University of Pisa, using sieves and Agriculex CB-1 Column Seed Cleaner complemented by manual cleaning.

Morphometric Analysis
In total, 49 qualitative and quantitative morphological characters ( Table 2, see also Supplementary Materials for details concerning the calyx) were studied, with a resulting dataset of 223 individuals × 49 variables. Macroscopic measures were taken with a digital calliper (error ± 0.1 mm), whilst microscopic and calyx measurements [30] (Table 2 and Figure S1) were taken through bar-scaled pictures with a Fiji 2.1.0 [31]. To provide a more objective means of counting the number of leaf veins, free-hand transversal sections of Biology 2022, 11, 1060 4 of 21 leaves were prepared. We considered as a "vein" each fascicule composed of the xylem and phloem surrounded by sclerenchyma. The anatomy of summer leaves was surveyed under a Leitz Diaplan light microscope at 40× ( Figure S2). We considered as "involucral bracts" those from the capitulum involucre, as "spikelet bracts" those subtending each spikelet, and as "bracteoles" those under each flower. To take into account the internal variability of the capitulum, we measured a spikelet collected from the middle of the capitulum ("inner spikelet") and a spikelet in contact with the inner involucral bract ("outer spikelet").  All statistical analyses were conducted in R Studio (version 3.6.2) [32]. To test the suitability of the data for factor analysis, the Kaiser-Meyer-Olkin test (MSA = 0.86, psych package [33] and Bartlett sphericity test (p < 0.001, REdaS package [34]) were performed successfully on the correlation matrix. Since there were mixed variables, Gower distance in the FD package [35] with Podani correction [36] was used, whilst Cailliez correction [37] was applied due to the violation of the triangle inequality (i.e., the matrix was not Euclidean). On such a dissimilarity matrix, Principal Coordinate Analysis (PCoA) in the ape package [38] was used to explore the dataset. Graphs were plotted with the ggplot2 package [39]. Oneway ANOSIM in PAST (version 4.09) [40] was used to test the null hypothesis of no difference between groups in the Gower dissimilarity matrix. To test the current taxonomic hypothesis and other alternative groupings based on our results, we applied jackknifed Linear Discriminant Analysis (LDA) in the MASS (for plotting) and Predpsych (to obtain the confusion matrix [41]) packages. Qualitative variables were converted into numbers with integer encoding. Using the PredPsych package, Cohen's Kappa coefficient was estimated for each grouping hypothesis. K coefficient is a measure of how the classification results compare with the values assigned and is generally thought to be a more robust measure than simple percentage agreement calculation, since it considers the possibility of agreement occurring by chance [42]. It ranges from 0 to 1 and K values greater than 0.75 may be taken to represent excellent agreement beyond chance [43].
Each character was statistically tested. For all the quantitative characters, normality was tested with the Shapiro-Wilk test. Normal and non-normal data were checked for homoscedasticity with the Bartlett test and Levene-Brown-Forsythe test, respectively. After checking statistical assumptions, normal and log-normal quantitative characters were compared among groups with one-way ANOVA and the post hoc Tukey-Kramer test (homoscedastic data) or Welch's ANOVA and the post hoc Games-Howell test (heteroscedastic data). On the contrary, non-normal characters were tested with the Kruskal-Wallis and multiple comparisons test of Wilcoxon-Mann-Whitney (homoscedastic data) or a permutation test implemented using the pairwisePermutationTest function of the rcompanion package [44] (heteroscedastic data). To control the family-wise error rate of multiple comparisons, Holm's correction was applied to all the tests. Qualitative nominal and ordinal characters were tested with the R function pairwiseNominalIndependence and pairwise-OrdinalIndependence based on Fisher's exact test and implemented in the rcompanion package [44]. All the tests were considered significant with α< 0.01. The number of statistically significant differences for a variable among population pairs was counted and a pairwise triangular matrix was built. Descriptive statistics for each group were calculated using the describeBy function in the psych package [33].

Seed Morpho-Colorimetric Analysis
For a sample of 100 seeds per accession (cleaned from the fruiting calyx and the membranous pericarp), digital images were acquired using a flatbed scanner (Epson Perfection V550) with a digital resolution of 1200 dpi. When an accession had fewer than 100 seeds, the analysis was carried out on the whole batch available. The system worked with 2D images: seeds were randomly disposed on the scanner tray, so that they did not touch one another, and covered using a box with white paper followed by a box with black paper to avoid interference from environmental light. The images were processed using the software package ImageJ (version 1.52b) (Available online: http: //rsb.info.nih.gov/ij. (accessed on 11 March 2022), and the descriptors of seed-size, shape, and colour features were measured and analysed. A plugin, Particles8 [45] (Available online: https://blog.bham.ac.uk/intellimic/g-landini-software/. (accessed on 11 March 2022), was used to measure 20 colorimetric and 26 morphometric features. This plugin was further enhanced by adding algorithms that can compute the Elliptic Fourier Descriptors (EFDs) for each analysed seed, thus increasing the number of independent variables. Following Terral et al. [46] and Sarigu et al. [47], to minimize measurement errors and optimize the efficiency of shape reconstruction, 20 harmonics were used to define the seed boundaries, obtaining 78 additional variables that were useful to discriminate between the studied seeds. In total, 124 morphometric and colorimetric characters were measured for each seed [48]. Statistical analyses were performed with the software SPSS release 16 (SPSS 16.0 for Windows; SPSS Inc., Chicago, IL, USA) by applying stepwise Linear Discriminant Analysis (LDA).

Karyological Analysis
Seeds were germinated in Petri dishes with 1% agar at 25 • C in an alternating 12/12 h dark/light photoperiod. After about 4 days, radicles emerged, and seedlings were removed from the seed incubator and kept at 4 • C for 24 h in a fridge, then we followed the Feulgen staining protocol. Root tips were pre-treated with 0.4% colchicine for 3 h and then fixed in Carnoy fixative solution for 1 h. After hydrolysis in HCl 1 N at 60 • C for 8 min, the root tips were stained in leuco-basic fuchsine for 2 h; root tips were squashed in a solution of aceto-orcein on a microscope slide.
Chromosomes were observed with a Leitz Diaplan microscope at 100× and pictures were taken with a Leica MC-170HD camera using Leica LAS-EZ 3.0 imaging software. At least four good metaphase plates were measured for each population. Lastly, chromosome numbers and karyological variables, such as THL (Total Haploid Length), M CA (Mean Centromeric Asymmetry), CV CL (Coefficient of Variation of Chromosome Length), and CV CI (Coefficient of Variation of Centromeric Index) were obtained from each plate with MATO 1.1 (version 20210101) [49]. Since all the karyological variables were normal and homoscedastic, they were statistically tested with One-way ANOVA and the post hoc Tukey-Kramer test for more than 3 groups or with two sample t-tests when comparing 2 groups.
Sequences were visually inspected and aligned using the ClustalW algorithm [53] implemented in BioEdit (version 7.2.5) [54] with the default values. An incongruence length difference (ILD) test was carried out in Nona (version 2.0) [55], as a daughter process of Winclada (version 1.00.08) [56], to test the putative incongruence of nuclear and chloroplast partitions prior to combination; default values were used for the analysis. A nucleotide evolution model was calculated for each of the five sequenced regions using jModelTest (version 2.1.10) [57], and the best fitting model was chosen over the others using the Bayesian Information Criterion (BIC) [58]. A Bayesian phylogenetic tree was inferred in MrBayes (version 3.2.6) [59] in two simultaneous, independent runs with the following settings: 2,000,000 generations of MCMC sampling every 2000 generations, and four runs (three cold and one hot). Convergence and mixing were evaluated in Tracer (version 1.7.2) [60]. The consensus Bayesian tree was visualized in FigTree (version 1.4.2) [61]. The best evolution models were K80 [62] for ITS and F81 [63] for chloroplast markers.

Comparative Niche Analysis
Occurrence data for the studied taxa were retrieved directly in the field, from SI-LENE (French National Mediterranean Botanical Conservatory of Porquerolles) (Available online: http://flore.silene.eu/index.php?cont=accueil. accessed on 17 May 2022), GBIF (Available online: https://www.gbif.org. accessed on 17 May 2022), and Wikiplantbase #Italia [64], with a total of 496 points. To test for differentiation in environmental space, we represented and quantified niche overlap using the PCA-based method developed by Broennimann et al. [65]. The Schoener's D index, which ranges from 0 (no overlap) to 1 (full overlap), was used to measure niche overlap [66]. We used niche similarity tests [67] to assess whether the ecological niches of the taxa were more similar than expected at random from their geographical ranges. Niche similarity tests compare the environmental conditions occupied by taxa, taking into account the environmental conditions that are available in the geographic area occupied by each taxon. Briefly, the observed climatic niche overlap between two taxa was compared, with the overlap measured between the niche of one taxon and the randomized niche of the other taxon. This randomized niche was obtained by randomly sampling occurrence points in buffer areas of 10 km around occurrences (the 'background area').

Nomenclature and Distribution
Currently accepted names, basionyms, and homotypic synonyms within Armeria arenaria and its subspecies studied here were taken from the Med-Checklist [68]. Information about the herbaria in which the original material could be stored was derived from [69]. Accordingly, we digitally examined the following herbaria: B, FI, L, LY, M, MPU, P, and SLA (herbarium acronyms follow Thiers [70]). Once we had elaborated the new taxonomic scheme, we used our identification key to assess the geographical distribution of the recognised taxa by checking the herbarium materials stored at ANC, APP, B, CAME, FI, HLUC, MJG, MW, P, and RO. This material was then georeferenced and used to build the map in Figure 6.

Morphometry
The first two axes of the PCoA explain 49% of the total variance. Along the first axis, there is a clear separation of four Apennine populations (AA, LA, MB, and MC, on the right side of Figure 1). Hereafter, we will refer to this group of four populations as "marginatoid". Another group uniting populations from northern Italy (BO, BR, GA, MP, PS, and TV) and France (FB and LL) emerged. We will refer to this group hereafter as "arenarioid" (on the left side of Figure 1). The MP population, initially attributed to A. arenaria subsp. apennina, clearly falls among arenarioid plants.
Biology 2022, 11, x FOR PEER REVIEW 9 of 23 "arenarioid" (on the left side of Figure 1). The MP population, initially attributed to A. arenaria subsp. apennina, clearly falls among arenarioid plants. Along the second axis, the topotypical population of A. arenaria subsp. arenaria (FB) shows a slight separation from the other arenarioid populations (Figure 2). One-way ANOSIM showed that there was, indeed, a significant difference between FB and the rest of the arenarioid populations (BO, BR, GA, MP, PS, TV, and LL) (R = 0.6573, p = 0.001), confirming the separation shown along the second axis of PCoA. LDA performed on the current taxonomic hypothesis (Table 3) obtained an 87% correct classification and K = 0.8. The lowest value of sensibility was scored by A. arenaria subsp. marginata (77.7%), followed by A. arenaria subsp. apennina (85.4%). The percentage of correct classifications and K increased to 99% and 0.9696, respectively, when comparing arenarioid with marginatoid plants. Along the second axis, the topotypical population of A. arenaria subsp. arenaria (FB) shows a slight separation from the other arenarioid populations (Figure 2). One-way ANOSIM showed that there was, indeed, a significant difference between FB and the rest of the arenarioid populations (BO, BR, GA, MP, PS, TV, and LL) (R = 0.6573, p = 0.001), confirming the separation shown along the second axis of PCoA. LDA performed on the current taxonomic hypothesis (Table 3) obtained an 87% correct classification and K = 0.8. The lowest value of sensibility was scored by A. arenaria subsp. marginata (77.7%), followed by A. arenaria subsp. apennina (85.4%). The percentage of correct classifications and K increased to 99% and 0.9696, respectively, when comparing arenarioid with marginatoid plants.  Table 1 and Figure 6. To further investigate the morphological variation within arenarioid plants, we carried out a pairwise comparison using univariate statistical analyses on single characters Figure 2 shows that the highest number of pairwise differences (94) was found between FB and all the other arenarioid populations. The population that shows the second most number of differences is PS (72).
Accordingly, we set up two new alternative grouping hypotheses in both of which marginatoid plants (AA, LA, MB, and MC) were combined in a single group. In the first grouping hypothesis (I), we tested FB, together with all the Northern Italian populations as belonging to the same taxon (as in the current taxonomic hypothesis) against the single population LL, which is the topotypical population of A. arenaria subsp. praecox. In the second grouping hypothesis (II), we tested LL as belonging to the same taxon as all the other Northern Italian arenarioid populations against the single population of FB (which corresponds to A. arenaria s.str.). The performance of LDA was 96% (K = 0.925) under  Table 1 and Figure 6. Table 3. Confusion matrix of the LDA based on the 49 morphometric characters, assuming the current taxonomic hypothesis of Armeria arenaria subspecies as a priori groups, as proposed by Arrigoni [17]. Rows show the membership of each a priori established group, whereas columns show the membership predicted by the classification model. To further investigate the morphological variation within arenarioid plants, we carried out a pairwise comparison using univariate statistical analyses on single characters. Figure 2 shows that the highest number of pairwise differences (94) was found between FB and all the other arenarioid populations. The population that shows the second most number of differences is PS (72).
Accordingly, we set up two new alternative grouping hypotheses in both of which marginatoid plants (AA, LA, MB, and MC) were combined in a single group. In the first grouping hypothesis (I), we tested FB, together with all the Northern Italian populations, as belonging to the same taxon (as in the current taxonomic hypothesis) against the single population LL, which is the topotypical population of A. arenaria subsp. praecox. In the second grouping hypothesis (II), we tested LL as belonging to the same taxon as all the other Northern Italian arenarioid populations against the single population of FB (which corresponds to A. arenaria s.str.). The performance of LDA was 96% (K = 0.925) under grouping hypothesis I and 98% (K = 0.968) under grouping hypothesis II. The two most important qualitative characters are provided in Tables S1-S3, whereas mean (± standard deviation) values of the quantitative morphological characters for each population are provided in Tables S4 and S5.

Seed Morpho-Colorimetry
The LDA performed on the current taxonomic hypothesis of a priori groups gave an overall cross-validated classification performance of 51.7% (Table 4). Armeria arenaria subsp. marginata showed the highest percentage of discrimination, with values of 71.5%, while the lowest (36.3 %) was detected in A. arenaria subsp. apennina (Table 4). Table 4. Confusion matrix of the LDA based on the seed morpho-colorimetric dataset (percentages of correct classification), assuming the current taxonomic hypothesis of Armeria arenaria subspecies as a priori groups, as proposed by Arrigoni [17]. Rows show the membership of each a priori established group, whereas columns show the membership predicted by the classification model. The second LDA, contrasting arenarioid and marginatoid plants, provides an overall percentage of 86% correct classification, with high discrimination performance for the two groups (Table 5). Table 5. Confusion matrix of the LDA based on the seed morpho-colorimetric dataset (percentages of correct classification), assuming "arenarioid" and "marginatoid" groups for Armeria arenaria populations. Rows show the membership of each a priori established group, whereas columns show the membership predicted by the classification model. According to two alternative grouping hypotheses derived from the morphometric analysis, FB was tested, together with all Northern Italian populations, as belonging to the same group, against the single population LL (hypothesis I), and LL was tested as belonging to the same group as all other Northern Italian arenarioid populations against the single population FB (hypothesis II). The discriminant analysis provided an overall percentage of classification of 77.3% and 84.4% for hypotheses I and II, respectively (Table 6, Figure 3). In hypothesis I, high discrimination performance was obtained for LL (78.8%) and marginatoid plants (81.3%). Concerning hypothesis II, higher performances, ranging from 91.0% (in FB) to 81.5% (in marginatoid plants), were detected.
Under grouping hypothesis I, keeping all the Italian arenarioid populations as A. arenaria subsp. arenaria and contrasting them with LL and with marginatoid plants, one-way ANOVA revealed that THL (F = 8.056; p < 0.001) and MCA (F = 10.52, p < 0.001) were significantly different, but not CVCL and CVCI. A post hoc Tukey-Kramer test showed that MCA and THL values differed significantly (p < 0.001) between marginatoid and Italian arenarioid plants (+FB). In contrast, MCA and THL were not significantly different between LL and the marginatoid group, or between the two arenarioid groups.
Under grouping hypothesis II, grouping all the Italian arenarioid populations with LL, contrasting them with FB and with marginatoid plants, One-way ANOVA revealed that THL (F = 8.158; p < 0.001), MCA (F = 11.03; p < 0.001), and CVCI (F = 5.221; p < 0.01) were significantly different among the three groups, while no difference was found in CVCL values.
A post hoc Tukey test showed that MCA and CVCI values differed significantly between marginatoid plants and FB at p < 0.01, whereas THL differs significantly at p < 0.001 between marginatoid plants and all Italian arenarioid plants (+ LL) (Figure 4). In contrast,

Karyotype Structure and Asymmetry
All the studied populations were diploid, with 2n = 2x = 18 chromosomes. They showed medium-sized (4.68 ± 0.64 µm), mostly metacentric (48.6%) or submetacentric (50.7%), chromosomes (see also Figures S3 and S4). One-way ANOVA revealed that all four karyological indices showed no significant differences among the four subspecies as circumscribed according to the current taxonomic hypothesis. However, the arenarioid plants (n = 45) showed significantly lower M CA (t = −4.52, df = 59, p < 0.001) and THL (t = −4, df = 59, p < 0.001) values when compared to marginatoid plants (n = 16). Lower, but not significantly different, values were also observed in CV CL and CV CI . Mean (± standard deviation) values for the karyological indices at population level are provided in Table S6.
Under grouping hypothesis I, keeping all the Italian arenarioid populations as A. arenaria subsp. arenaria and contrasting them with LL and with marginatoid plants, one-way ANOVA revealed that THL (F = 8.056; p < 0.001) and M CA (F = 10.52, p < 0.001) were significantly different, but not CV CL and CV CI . A post hoc Tukey-Kramer test showed that M CA and THL values differed significantly (p < 0.001) between marginatoid and Italian arenarioid plants (+FB). In contrast, M CA and THL were not significantly different between LL and the marginatoid group, or between the two arenarioid groups.
Under grouping hypothesis II, grouping all the Italian arenarioid populations with LL, contrasting them with FB and with marginatoid plants, One-way ANOVA revealed that THL (F = 8.158; p < 0.001), M CA (F = 11.03; p < 0.001), and CV CI (F = 5.221; p < 0.01) were significantly different among the three groups, while no difference was found in CV CL values.
A post hoc Tukey test showed that M CA and CV CI values differed significantly between marginatoid plants and FB at p < 0.01, whereas THL differs significantly at p < 0.001 between marginatoid plants and all Italian arenarioid plants (+ LL) (Figure 4). In contrast, there was no significant difference between Italian arenarioid plants (+ LL) and FB in any of the studied karyological indices. there was no significant difference between Italian arenarioid plants (+ LL) and FB in any of the studied karyological indices.

Molecular Systematics
The number of phylogenetically informative characters obtained from the amplification of the five markers was 36, corresponding to approximately 1.5% of the entire alignment. The markers that showed the highest number of informative characters were the intergenic spacers trnL-rpl32 and ITS, with 13 and 11 phylogenetically informative characters (Table S9). The results of the ILD test showed that all plastid markers were congruent (p > 0.05, Table S10). On the contrary, increasing the number of replicates (up to 100), all the pairwise combinations were congruent except for ITS and trnL-rpl32, which turned incongruent at p = 0.01 (Table S10). Indeed, removing the trnL-rpl32 marker, the ITS and the resulting concatenated plastid matrix become congruent (p = 0.09). However, since the topology of the concatenated trees with and without trnL-rpl32 were not in conflict, we decided to retain the full matrix, which was 2337 bp long.
The Bayesian concatenated consensus unrooted tree is shown in Figure 5. Arenarioid populations are split into two main clades but are collectively well distinct from marginatoid (+ arenarioid PS) populations. The former main clade is more variable and encompasses accessions from the French populations FB and LL (forming a clade), accessions from MP, a clade with the accessions from GA, TV, BO (the latter in a separate clade), as well as those from BR, which do not form a monophyletic group. The second main clade contains two clades and the accessions from PS with an unresolved position. Separate ITS and plastid phylogenies are provided in Figures S5 and S6.

Molecular Systematics
The number of phylogenetically informative characters obtained from the amplification of the five markers was 36, corresponding to approximately 1.5% of the entire alignment. The markers that showed the highest number of informative characters were the intergenic spacers trnL-rpl32 and ITS, with 13 and 11 phylogenetically informative characters (Table S9). The results of the ILD test showed that all plastid markers were congruent (p > 0.05, Table S10). On the contrary, increasing the number of replicates (up to 100), all the pairwise combinations were congruent except for ITS and trnL-rpl32, which turned incongruent at p = 0.01 (Table S10). Indeed, removing the trnL-rpl32 marker, the ITS and the resulting concatenated plastid matrix become congruent (p = 0.09). However, since the topology of the concatenated trees with and without trnL-rpl32 were not in conflict, we decided to retain the full matrix, which was 2337 bp long.
The Bayesian concatenated consensus unrooted tree is shown in Figure 5. Arenarioid populations are split into two main clades but are collectively well distinct from marginatoid (+ arenarioid PS) populations. The former main clade is more variable and encompasses accessions from the French populations FB and LL (forming a clade), accessions from MP, a clade with the accessions from GA, TV, BO (the latter in a separate clade), as well as those from BR, which do not form a monophyletic group. The second main clade contains two clades and the accessions from PS with an unresolved position. Separate ITS and plastid phylogenies are provided in Figures S5 and S6.

Comparative Niche Analysis
Schoener's D values were generally low, ranging from 0 to 0.208. In particular, A. arenaria subsp. praecox was the subspecies showing a niche that overlapped less with those of the other subspecies, according to the current taxonomic hypothesis ( Table 7). The lack of significance in the similarity test indicated that the low niche overlap values were due to habitat availability in the background areas rather than an effect of habitat selection. Taken together, these results suggest differences in optimal niche positions without niche shift.

Comparative Niche Analysis
Schoener's D values were generally low, ranging from 0 to 0.208. In particular, A. arenaria subsp. praecox was the subspecies showing a niche that overlapped less with those of the other subspecies, according to the current taxonomic hypothesis ( Table 7). The lack of significance in the similarity test indicated that the low niche overlap values were due to habitat availability in the background areas rather than an effect of habitat selection. Taken together, these results suggest differences in optimal niche positions without niche shift. Table 7. Results of niche similarity tests in environmental spaces among the different taxa and circumscription hypotheses of Armeria arenaria. Backgrounds were defined by applying 10 km buffer zones around the occurrence points. Current taxonomic hypothesis as stated by Arrigoni [17], (I) = first alternative grouping hypothesis, (II) = second alternative grouping hypothesis. ns = not significant.

Discussion
All our results concur in highlighting that the current taxonomic hypothesis available for Armeria arenaria is no longer supported. Starting from the marginatoid plants, there is no morphometric support at all for distinguishing the two taxa as proposed by Arrigoni [17]. Moreover, the idea that A. arenaria subsp. apennina represents a taxon somehow intermediate between A. arenaria subsp. marginata and A. arenaria subsp. arenaria [17] is only supported by their climatic requirements. Nevertheless, it also should be noticed that the two putative marginatoid taxa show the highest values of niche overlap detected. There is no karyological difference between the two marginatoid taxa, but together they show higher M CA and THL values with respect to the arenarioid plants. Phylogenetically, all marginatoid plants form a highly supported clade, in which the accessions of the two putative subspecies are intermingled. A single alpine arenarioid population (PS) is placed phylogenetically close to the marginatoid plants, suggesting that the genetic differentiation between arenarioid and marginatoid plants occurred only recently and may be derived from incomplete lineage sorting or gene flow. Despite this, morphological evidence fully places PS among arenarioid plants. Accordingly, we deem that the maintenance of the subspecific rank for marginatoid with respect to arenarioid plants is appropriate. Concerning arenarioid plants, they share a set of morphological and karyological features. Altogether, our data also support the maintenance of the subspecific rank for Armeria arenaria subsp. praecox with respect to A. arenaria subsp. arenaria, albeit with different circumscriptions, since all the Italian arenarioid populations agree much better with A. arenaria subsp. praecox than with A. arenaria subsp. arenaria. Indeed, from a morphometric point of view, the FB population (A. arenaria s.str.) shows the highest number of pairwise differences among all the other arenarioid populations, and it also shows the smallest seeds.
As a consequence, we exclude Armeria arenaria subsp. arenaria from the Italian flora, in favor of A. arenaria subsp. praecox, so that the range of the former subspecies is now reduced to Portugal, Spain, and France [68]. We cannot rule out that the range of that subspecies could be further narrowed in the future, given that this taxon is "conceived as a mixed bag that includes the variability of the rest of the populations" [25,71]. Armeria arenaria subsp. praecox has been only doubtfully recorded for Italy so far [17]. However, we clearly show that Italian arenarioid plants have a morphology highly overlapping that of the typical A. arenaria subsp. praecox (Figure 1). Indeed, the highest values of correct classification and K obtained by the discriminant analyses conducted for morphology and seed morpho-colorimetry were found when all the Italian arenarioid populations were grouped with LL and not with FB (which corresponds to the typical A. arenaria s.str.). Italian arenarioid populations are also phylogenetically more closely related to LL than to FB in the plastid tree ( Figure S6). The possible occurrence in Italy of other subspecies occurring in Southern France can be excluded based on the comparison of our data with those published by Baumel et al. [72], Tison et al. [73], and Tison et De Foucault [74] (data not shown).
Geographically and climatically, the marginatoid plants from the Northern Apennines are replaced in the central-western alpine and perialpine areas by A. arenaria subsp. praecox, which is in turn replaced by A. arenaria subsp arenaria in Central-Northern France ( Figure 6). plastid tree ( Figure S6). The possible occurrence in Italy of other subspecies occurring Southern France can be excluded based on the comparison of our data with those pu lished by Baumel et al. [72], Tison et al. [73], and Tison et De Foucault [74] (data n shown).

Conclusions
In this work, we conducted an integrative taxonomic study of Armeria arenaria in Italy. On the basis of nomenclatural, morphometric, seed morpho-colorimetric, karyological, molecular, and comparative niche evidence we were able to demonstrate that the current taxonomic setting for this species is no longer supported. Specifically, we proved that A. arenaria subsp. apennina is a heterotypic synonym of A. arenaria subsp. marginata and that all the previous records of A. arenaria subsp. arenaria for Italy should be attributed to A. arenaria subsp. praecox. Finally, we also provided an identification key for dried herbarium specimens to facilitate the identification of these taxa. The same key was used to reconstruct the distribution of the three subspecies based on 81 herbarium specimens.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology11071060/s1, Table S1, Shape of the leaf apex, according to the second grouping hypothesis; Table S2, Presence of teeth on summer leaf margins, according to the second grouping hypothesis; Table S3, Number of veins on summer leaves based on a cross section, according to the second grouping hypothesis; Table S4, Mean and standard deviation values of quantitative characters for the twelve populations studied; Table S5, Median and inter-quartile range of quantitative discrete characters for the twelve populations studied; Table S6, Means and standard deviations for the karyological indices in Armeria arenaria from the twelve studied populations; Table S7, List of molecular markers and their primers; Table S8, PCR settings used; Table S9, Number of phylogenetically informative characters per marker and for the concatenated matrix; Table S10, IDL test results; Figure S1, Schematic drawing of a calyx and related morphometric characters; Figure S2, Schematic drawings of selected summer leaf cross sections from the type localities of the four putative Armeria arenaria subspecies; Figure S3, Haploid idiograms (x = 9) of the twelve populations studied; Figure S4, Selected metaphasic plates of the twelve populations studied; Figure S5. Bayesian unrooted consensus phylogenetic tree of ITS matrix; Figure S6. Bayesian unrooted consensus phylogenetic tree of plastid matrix.