Pedicularis rostratospicata subsp. marsica (P. Sect. Rostratae, Orobanchaceae), a New Subspecies from the Central Apennines (Italy)

The new subspecies Pedicularis rostratospicata subsp. marsica is here described based on morphological and molecular analyses. The new taxon is endemic to few localities of the Central Apennines within the Abruzzo, Lazio and Molise National Park (Central Italy). Pedicularis rostratospicata subsp. marsica can be distinguished from the other two currently accepted subspecies (subsp. rostratospicata and subsp. helvetica) by their taller stems, shorter petioles of basal and middle leaves, narrower blades of basal leaves, longer inflorescences with longer internodes and more flowers, and longer calyx lobes. Molecular analysis confirmed the autonomy of the new taxon. Furthermore, the conservation status assessment of the new subspecies according to IUCN categories and criteria is proposed and discussed, and an analytical key to the three subspecies of P. rostratospicata is presented.


Introduction
Pedicularis L. is the largest genus of Orobanchaceae Vent., with more than 600 species worldwide; it occurs mainly in arctic-alpine regions in the northern hemisphere [1,2]. There is no a recent worldwide revision of Pedicularis, or even a comprehensive phylogeny. In Italy, it is represented by 26 taxa (species and subspecies), and among these, only P. elegans Ten. is endemic to the Central and Southern Apennines [3][4][5].
Pedicularis rostratospicata Crantz belongs to a group of species that has usually been called P. sect. Rostratae Benth. (syn. Rhyncholophae), characterized by perennial herbs, with alternate floral and cauline leaves, and flowers with an upper lip that is more or less falcate, terminating in a distinct, usually long beak without teeth [6][7][8]. Currently, two subspecies are accepted: subsp. rostratospicata occurring mainly in the Central and Eastern Alps and in one locality of the Northern Apennines distributed in Germany, Helvetia, Austria, Slovenia and Italy [9] where it is recorded in Trentino-Alto Adige, Friuli-Venezia Giulia and Emilia-Romagna at Mt. Prado [4,10]; subsp. helvetica (Steininger) O.Schwarz occurring mainly in the Central and Western Alps distributed in France, Helvetia, Austria and Italy [9] where it is recorded in Valle d'Aosta, Piemonte, Lombardia, Trentino-Alto Adige and Liguria [4].
In Central Italy, P. rostratospicata is known only for few localities of the Abruzzo, Lazio and Molise National Park in Abruzzo administrative region [11,12]. According to Bruno and Bazzichelli [13], Bazzichelli and Furnari [14], Pignatti [15], and Conti and Bartolucci [16], this taxon deserves further study and actually does not fit among the known subspecies.
A morphological investigation and a molecular analysis have been carried out on living specimens and herbarium material coming from the Central Apennines (Italy) and the Alps,

Morphometric Analyses
A total of 20 morphological characters were selected and scored for 15 dried individuals of Pedicularis rostratospicata s.l. ("MAR") from three localities in the Central Apennines (Abruzzo: Mt. Marsicano, Mt. Petroso-Valle Cupella, and Serra delle Gravare), 19 individuals of P. rostratospicata subsp. helvetica ("HEL") from several localities in the Western and Central Alps (Italy, Austria and Helvetia), and 21 individuals of P. rostratospicata subsp. rostratospicata ("ROS") from the Eastern Alps (Italy, Austria and Slovenia), with a resulting dataset of 55 individuals × 20 variables. Of the morphological characters studied, 18 were quantitative and 2 were qualitative ( Table 1). The quantitative and qualitative data were subjected to a principal coordinate analysis (PCoA) and cluster analysis (CA) using the average linkage method (UPGMA). The similarity matrix was calculated using the Gower coefficient, which is suitable for mixed data [19]. To test the current taxonomic hypothesis, we applied jackknifed linear discriminant analysis (LDA) with the individuals being a priori assigned to the three taxa based on morphological features. The multivariate analyses were carried out with the PAST v.4.11 software package [20].
Each quantitative continuous character was also subjected to univariate analysis (a Kruskal-Wallis test with Bonferroni corrections for multiple comparisons) with SPSS v.25 software [21]. Furthermore, the variability of the analyzed morphological characteristics was described using standard statistical parameters (mean, standard deviation, minimum, maximum, and 25th and 75th percentiles).

AFLPseq Fingerprinting
Genomic DNA for genetic fingerprinting was extracted in accordance with the CTAB DNA extraction protocol of Doyle and Dickson [22] and Doyle and Doyle [23]. The AFLPseq fingerprinting method was proposed by Dorfner et al. [24] and combines the genome complexity-reducing AFLP approach [25] with the next-generation sequencing (NGS) of resulting AFLP bands using the Nanopore sequencer MinION from Oxford Nanopore Technologies (Oxford, United Kingdom). It provides sequence and single-nucleotide polymorphism (SNP) information for hundreds of anonymous loci from across the whole genome and could be used for population genetic, phylogenetic, and species delimitation studies. It is suited for both well-preserved DNA from silica-gel dried leaf material and degraded DNA from herbarium specimens.
The present AFLPseq study comprised 11 Pedicularis accessions (Table 2), either recently collected, silica gel dried material (one accession of P. rostratospicata subsp. helvetica from Valle d'Aosta administrative region and three accessions of P. rostratospicata subsp. rostratospicata from Trentino-Alto Adige administrative region) or herbarium material housed in APP. The AFLPseq procedure followed the protocol given in Dorfner et al.'s work [24] with the following modifications. In the restriction ligation step, we used a double-digestion procedure with restriction enzymes MseI and EcoRI. After the ligation of MseI and EcoRI adapters (MseI adapter: 5 -GACGATGAGTCCTGAG-3 + 5 -TACTCAGGACTCAT-3 ; EcoRI adapter: 5 -CTCGTAGACTGCGTACC-3 + 5 -AATTGGTACGCAGTCTAC-3 ), we continued with the AFLP genome reduction protocol using primers with 1bp overhangs (MseI-C: 5 -GATGAGTCCTGAGTAAC-3 ; EcoRI-A: 5 -GACTGCGTACCAATTCA) in the pre-selective amplification step and in the selective amplification step with additional 1bp (EcoRI 5 -GACTGCGTACCAATTCAA-3 ) or 2bpoverhangs (MseI 5 -GATGAGTCCTGAGTAACTG-3 ), respectively. The two primers used in the latter amplification step, however, were additionally tailored to include Nanopore barcode adapter sequences at the 5 end of the primers (Mse_CTG_Nanopore_fw: 5 -TTTCTGTTGGTGCTGATATTGCGATGAGTCCTGAGTAACTG-3 ; Eco_AA_Nanopore_rv: 5 -ACTTGCCTGTCGCTCTATCTTCGACTCCGTACCAATTCAA-3 ), as suggested in the 'Ligation sequencing amplicons-PCR barcoding (SQK-LSK109 with EXP-PBC001)' protocol by Oxford Nanopore Technologies, substituting a subsequent ligation of the Nanopore barcode adapter with an additional barcoding PCR. To ensure specific binding with long and tailed primers, a two-step variation of the preselective PCR was conducted (94 • C for 2 min; followed by 30 cycles of 94 • C for 20 s and 72 • C for 2 min; and a final step at 72 • C for 2 min). To every 2 µL of a 1:10 diluted preselective PCR product, 5 µL of Taq DNA Polymerase Master Mix RED, 0.25 µL of each 10 µM tailed selective primer, and 2.5 µL of H 2 O were added. After the selective PCR, the length of the fragments ranged from 200 to 700 bp. All subsequent steps (Nanopore barcode PCR, sample multiplexing, size selection, and preparation of Nanopore sequencing library) followed those of Dorfner et al. [24]. The resulting library was sequenced with MinION using three Flongle flow cells.
For the bioinformatic analysis of the Nanopore read data, we employed a simplified version of the SLANG pipeline [24], developed specifically for assembling and orthologizing loci from Nanopore reads while also facilitating variant calling. In our modified approach, we deviated from the original SLANG pipeline by clustering all reads of all OTUs (Operational Taxonomic Units) simultaneously using a single VSEARCH v.2.15.0 [26] clustering step. The clustering threshold was set at 0.85, indicating a high level of sequence similarity required for clustering, while also accounting for variation and Nanopore sequencing errors. This new approach eliminated the need for the within-sample clustering of individual OTUs, as well as among-sample clustering for locus orthologization, as all OTU reads were already orthologized via read clustering during the initial step of this modified approach. For variant identification, we adopted the reference-based variant-calling approach of SLANG, albeit with the modification that allows all reads from an OTU to be mapped simultaneously to all consensus sequences, deviating from the original methodology, where only reads from the respective cluster were allowed to map to their respective reference sequence. This adjustment was implemented to minimize missing data and enhance data completeness. For mapping reads to a reference sequence, MINIMAP2 v.2.17 [27] was implemented. Subsequent variant extraction was employed in SLANG through BCFTOOLS [28]. Variants were considered valid if their frequency at the corresponding alignment position was equal to or greater than 30%, enabling the detection of heterozygosity in diploid individuals. Additionally, a minimum read depth of 5 for that position was required to ensure robust variant calling. Variants passing the filtering criteria were exported as a SNP alignment from the resulting VCF file through a custom script. Finally, we utilized the Python package PYCKMEANS (accessible at https://github.com/TankredO/pyckmeans; accessed on 6 June 2023) to compute Jukes-Cantor genetic distances and subsequently perform a principal coordinate analysis (PCoA) to visualize the relationships among the Pedicularis accessions.

Morphometric Analyses
The first two axes of the PCoA explain 55% of the total variance. The scatterplot shows on the first two axes a clear separation of the three taxa (MAR, HEL, and ROS) involved in this study, and no overlapping areas among individuals were found (Figure 1). The UPGMA dendogram ( Figure 2) shows two main clusters ("a" and "b"): cluster "a" composed of all the specimens of ROS, and cluster b composed of specimens of HEL and MAR. The latter cluster is further separated into two well-defined subclusters; "1" includes all the specimens referable to MAR and "2" includes all the specimens of HEL. Discriminant analysis resulted in 100% (jackknifed) correct classification of individuals a priori attributed to the three taxa ( Figure 3). No overlap between individuals of MAR, ROS, and HEL was found, showing a clear distinction between these units.  The results of the univariate analysis of continuous characters are summarized in Table 3 and depicted in Figure 4. Based on the Kruskal-Wallis test, we conclude that at least one taxon is different in terms of the following characters (p < 0.01): H, PBLL, BBLL, PMLL, BMLL, BMLW, IL, IBIL, CL, CLL, COL, and NL. The states of three characters (H, CL, and SL) are significantly different between the ROS and the other two taxa (p < 0.01). The states of five characters (PBLL, BBLL, BMLL, CL, and NL) are significantly different between HEL and the other two taxa (p < 0.01). The states of three characters (PMLL, IL, and IBIL) are significantly different between MAR and the other two taxa (p < 0.01). In addition, the state of BMLW is significantly different between MAR and HEL.   12.97 ± 1.23 (10.1-)12.5-13.5(-16) 15 .  Morphological analyses provide evidence that the Apennines populations of P. rostratospicata should be regarded as a new taxon, endemic to Central Italy. In addition, the study of the specimens collected on Mt. Prado (Emilia-Romagna administrative region) housed in FI and SPAL allowed us to review and refer to them as P. rostratospicata subsp. helvetica, as previously and correctly reported by Alessandrini et al. [29]. In subsequent papers [4,30], based on the indications from A. Alessandrini, however, it has been attributed by mistake to P. rostratospicata subsp. rostratospicata.

AFLPseq Fingerprinting
In total, 194,212 reads and 68.66 Mbp were sequenced for the 11 Pedicularis accessions. After read preprocessing, 110,374 reads with lengths between 100 bp and 700 bp passed the Q5 quality filter. With the modified SLANG pipeline at a cluster threshold of 0.85, 320 or-thologous loci were inferred, containing 3139 SNPs. After the calculation of Jukes-Cantor distances, the resulting PCoA plot was received (Figure 5), indicating a clear separation of P. rostratospicata subsp. rostratospicata (A1285-A1287) from other accessions through PCo axis 1 (explaining 73.2% of the total variation), while accessions from P. rostratospicata subsp. helvetica (A1288-A1291) and those from the Central Apennines (A1292-A1295) were separated on PCo axis 2 (explaining only 4.5% of the total variation). These results match the multivariate analysis of morphological characters in the study group (Figures 1 and 2), where P. rostratospicata subsp. rostratospicata is well-separated from the other two taxa, while these show overlapping variation in some of the surveyed characters (Table 2, Figure 4). Following Oberprieler [31] in his argumentation scheme on conceptualizing species rank decisions (the 'Wettstein tesseract'), the present situation of genetic and morphological variation in the P. rostratospicata assembly may be best tackled in a taxonomic treatment either accepting three allopatric subspecies of a single species or dividing it into two species with P. rostratospicata subsp. rostratospicata on the one hand and subsp. helvetica plus the Central Apennines (Abruzzo) populations on the other. Despite the strong genetic and morphological separation of P. rostratospicata subsp. rostratospicata and subsp. helvetica seen in the present analyses, the lack of sympatry (mainly in the Eastern Alps for the former, Western and Central Alps for the latter) and the seemingly similar ecology of the two taxa (both being centered in the Caricion ferrugineae assembly showing the same substrate, precipitation and elevation niche and the same flowering time; [32]) argues for the subspecies rank for the two taxa as long as nothing is known about the strength of reproductive barriers between the two (code '1101' in Figure 2 of Oberprieler 2023 [31]). On the other hand, the morphological distinctness (at least in multivariate analyses) of the allopatrically distributed P. rostratospicata subsp. helvetica and Central Apennines (Abruzzo) populations of the species, along with the only marginal differentiation in genetic respects, is also best captured taxonomically at the subspecific rank (code '0101' or '0111' in Figure 2 of Oberprieler 2023 [31]).

Taxonomic Treatment
Pedicularis rostratospicata subsp. marsica F. Conti & Bartolucci, subsp. nov. (Figure 6). Diagnosis: Pedicularis rostratospicata subsp. marsica is distinguished from the other subspecies by its taller stems, shorter petioles in basal and middle leaves, narrower blades in basal leaves, longer inflorescences with longer internodes and more flowers, and longer calyx lobes (Figure 8).
Habitat: Ledges, screes in alpine belt from 1900 to 2100 m a.s.l., subject to long snow cover. Phenology: Flowering from the end of June to beginning of August, fruiting from the end of July to August.
Distribution: Endemic to the National Park of Abruzzo, Lazio and Molise (SAC IT7110205 "Parco Nazionale d'Abruzzo") in Abruzzo administrative region distributed at Forca Resuni, Mt. Petroso, Valle Cupella, Serra delle Gravare, Mt. Marsicano [11][12][13][14]16,33]. The geographic distribution of P. rostratospicata based on the herbarium specimens examined is shown in Figure 9. Conservation status: Pedicularis rostratospicata subsp. marsica is known to be in only one location (five subpopulations) within the NATURA 2000 network in SAC IT7110205 "Parco Nazionale d'Abruzzo". The extent of occurrence (EOO) is 35,374 km 2 , and the area of occupancy (AOO) is 20 km 2 (cell grid 2 × 2 km), as calculated with the GeoCAT (Geospatial Conservation Assessment Tool) software [34]. The main estimated threats include global warming, which favors the arrival of more thermophilous competitive species, and restricted range, which causes a greater risk of extinction in the event of disease or pest attacks. We have no data to evaluate whether or not the population is in decline and whether or not extreme fluctuations occur. According to IUCN criteria [35], the species is assessed to be near-threatened (NT).

Conclusions
In this work, we conducted a morphometric and molecular study of Pedicularis rostratospicata. We were able to demonstrate that the populations of this species occurring in the Central Apennines should be referred to a subspecies new to science, described here with the name P. rostratospicata subsp. marsica. Finally, we provided the conservation status assessment of the new subspecies according to IUCN categories and criteria (assessed as NT), an identification key for the three subspecies of P. rostratospicata (for dried herbarium specimens) and a distribution map based on the herbarium material examined.

Data Availability Statement:
The data presented in the current study are available within the article. Raw data of the AFLPseq fingerprinting were submitted to GenBank under the BioProject accession number PRJNA984162.