Nélida Padilla-García

Blanca M Rojas-Andrés

Noemí López-González

Mariana Castro

Sílvia Castro

João Loureiro

Dirk C Albach

Nathalie Machon

M Montserrat Martínez-Ortega

The challenge of species delimitation in the diploid-polyploid complex Veronica subsection Pentasepalae

Keywords:

A reliable taxonomic framework and the identification of evolutionary lineages are essential for effective decisions in conservation biodiversity programs. However, phylogenetic reconstruction becomes extremely diffi-cult when polyploidy and hybridization are involved. Veronica subsection Pentasepalae is a diploid-polyploid complex of ca. 20 species with ploidy levels ranging from 2x to 10x. Here, DNA-ploidy level estimations and AFLP fingerprinting were used to determine the evolutionary history, and species boundaries were reviewed in an integrated approach including also previous data (mainly morphology and sequence-based phylogenetic reconstructions). Molecular analyses were performed for 243 individuals from 95 populations, including for the first time all taxa currently recognized within the subsection. Phylogenetic reconstruction identified four main groups corresponding almost completely to the four clusters identified by genetic structure analyses. Multiple autopolyploidization events have occurred in the tetraploid V. satureiifolia giving rise to octoploid entities in central Europe and north of Spain, whereas hybridization is demonstrated to have occurred in several populations from the Balkan Peninsula. Furthermore, our study has established the taxonomic status of taxa, for the most part recovered as monophyletic. Cryptic taxa within the group have been identified, and a new species, Veronica dalmatica, is fully described. This study highlights the implications of polyploidy in species delimitation, and illustrates the importance to conserve polyploid populations as potential sources of diversification due to evolutionary significance of genome duplications in plant evolution.

Introduction

The delimitation of species boundaries is a classic problem for biologists. Until about seventy years ago, taxonomists have focused on morphological differences between species for their circumscription, rather than on coherence with their evolutionary history. However, since the 1940s, a wider interest in the evolutionary history of organ-isms arose [START_REF] Huxley | The New Systematics[END_REF]. In 1950, Hennig published his theory of phylogenetic systematics giving rise to the origin of cladistics, which revolutionized the field of taxonomy [START_REF] Hennig | Grundzüge einer Theorie der phylogenetischen Systematik[END_REF]. Despite originally considered for the analysis of morphological characters, it is equally suitable for other types of characters that have been used with taxo-nomic purposes during the last decades. Currently, molecular phylogenies, complementing morphological characters, are the key instruments for biologists and biosystematists who try to understand the evolutionary processes that shape the history of species. Nevertheless, evolutionary histories involving radiations or complex processes such as hybridization, introgression and/or polyploidization, complicate phylogenetic reconstruction [START_REF] Naciri | Species delimitation and relationships: the dance of the seven veils[END_REF]. This, together with a lack of morphological differences and uncertainties over reproductive isolation among polyploids and their diploid progenitors, has resulted in taxonomic biases within polyploid complexes [START_REF] Soltis | Autopolyploidy in angiosperms: have we grossly underestimated the number of species?[END_REF][START_REF] Barley | The challenge of species delimitation at the extremes: diversification without morhological change in Philippine sun skinks[END_REF]. Here, the importance of the species concept is fundamental. The biological concept of species proposed by [START_REF] Mayr | Systematics and the Origin of Species, from the Viewpoint of a Zoologist[END_REF] is difficult to apply when working with closely related species in which hybridization and introgression are common. Most plant taxonomists have traditionally relied on morphology to delimit species boundaries (i.e., morphological species concept), whereas others adopted in the last decades a concept based on genetic differ-ences and monophyly (i.e., genetic and phylogenetic species concepts). However, in species groups with frequent hybridization and poly-ploidization, the general lineage concept [START_REF] De Queiroz | Different species problems and their resolution[END_REF][START_REF] De Queiroz | Species concepts and species delimitation[END_REF]) may be more appropriate. According to this concept, species are defined as separately evolving metapopulation lineages, which can be identified by different properties accumulated by that species during the process of speciation (e.g., reproductive isolation, morphological or genetic differences, monophyly, etc.). This general lineage concept has been broadly adopted and promoted the development of an integrative taxonomic approach in which multiple and complementary methods are integrated to delimit species boundaries [START_REF] Dayrat | Towards integrative taxonomy[END_REF]. This ap-proach argues that taxonomic inference should be based on congruence across different types of characters and analyses. When results from different sources of data are incongruent, caution to delimit species is preferable since taxonomic conclusions may have significant implica-tions [START_REF] Carstens | How to fail at species delimitation[END_REF]. For example, regarding conservation and sustainable use of natural resources, in accordance with the Convention on Biological Diversity (https://www.cbd.int/gti/importance.shtml), taxonomy is necessary for effective decision making, because it pro-vides basic understanding about the components of biodiversity. In a world where wild species are increasingly under threat, the conserva-tion status of a taxon can only be correctly evaluated under the light of a clear taxonomic framework [START_REF] Mace | The role of taxonomy in species conservation[END_REF]. The identification and preservation of evolutionary processes is also essential in conservation programs, especially for endangered, rare and endemic species.

In the present study, complementary methodologies are used to address the taxonomic challenges of a study group with a complex evolutionary history, Veronica subsection Pentasepalae Benth. This subsection is a monophyletic lineage within V. subgenus Pentasepalae (Benth.) M.M. Mart. Ort., Albach & M.A. Fisch. (Rojas-Andrés et al., 2015). It has a recent origin (mean crown age 2.8 Mya., [START_REF] Meudt | Is genome downsizing associated with diversification in polyploid lineages of Veronica?[END_REF] and is composed of ca. 20 closely related perennial species distributed in Eurasia and North Africa. Interestingly, the group comprises five species and three subspecies endemic to single countries or sometimes only a small area within one country. Some of them are included in regional, national and/or international Red Lists (Peñas de [START_REF] Peñas De Giles | Veronica tenuifolia subsp. fontqueri (Pau) M.M. Mart. Ort[END_REF][START_REF] Cabezudo | Lista Roja de La Flora Vascular de Andalucía[END_REF]; Alcántara de la [START_REF] De La Fuente | Catálogo de especies amenazadas en Aragón[END_REF][START_REF] Petrova | Red list of Bulgarian vascular plants[END_REF][START_REF] Bilz | Veronica turrilliana. The IUCN Red List of Threatened Species[END_REF] International Union for Conservation of Nature, 2016), although there is a clear lack of information for numerous species that have not yet been carefully evaluated. The most important diversification center of V. subsection Pentasepalae is the Balkan Peninsula. The group is characterized by the presence of a pentapartite calyx with the fifth sepal being significantly smaller, by a capsule usually rounded at the base, and a base chromosome number of x = 8. However, although the group is well defined within Veronica [START_REF] Albach | Chromosome numbers in Veroniceae (Plantaginaceae): review and several new counts[END_REF], the existence of morphologically intermediate forms within the subsection due to overlapping morphological character states makes V. subsection Pentasepalae one of the taxonomically most complicated groups within the genus [START_REF] Albach | Phylogeny of Veronica in the Southern and Northern Hemispheres based on plastid, nuclear ribosomal and nuclear low-copy DNA[END_REF]. Since Bentham described V. subsection Pentasepalae in 1846, numerous taxonomic treatments have been proposed (for a historical review of monographs and Floras, see Rojas-Andrés and Martínez-Ortega, 2016), and several studies based on morphological, karyological or molecular data have tried to elucidate the evolutionary history of the group (e.g., [START_REF] Martínez-Ortega | Palynological study of Veronica sect. Veronica and sect. Veronicastrum (Scrophulariaceae) and its taxonomic significance[END_REF], 2004, 2009;[START_REF] Andrés-Sánchez | Combining traditional morphometrics and molecular markers in cryptic taxa: towards an updated integrative taxonomic treatment for Veronica subgenus Pentasepalae (Plantaginaceae sensu APG II) in the western Mediterranean[END_REF]. In the most recent molecular study, [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF] used nuclear and plastid DNA sequence data to perform a phylogenetic analysis of the subsection. Despite the contributions of that study to the understanding of the evolutionary history of the group, a high degree of incongruence was found between the ITS and plastid DNA trees, probably caused by hybridization and in-complete lineage sorting (ILS). Hence, some questions about the monophyly and the relationships among species remained unresolved. Such questions are unlikely to be answered using a few loci alone, especially considering the prevalence of hybridization and polyploidization in the genus [START_REF] Albach | Incongruence in Veroniceae (Plantaginaceae): evidence from two plastid and a nuclear ribosomal DNA region[END_REF].

The variety of ploidy levels in the subsection, ranging from diploid to decaploid (data previous to 2008 summarized by [START_REF] Albach | Chromosome numbers in Veroniceae (Plantaginaceae): review and several new counts[END_REF][START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF], indicate that polyploidy has been a crucial process in the evolution of the group. Polyploidy or whole-genome duplication (WGD) is a frequent mechanism of evolution and speciation in flowering plants [START_REF] Stebbins | Variation and Evolution in Plants[END_REF][START_REF] Grant | Plant Speciation[END_REF][START_REF] Soltis | Advances in the study of polyploidy since plant speciation[END_REF][START_REF] Soltis | Polyploidy and angiosperm diversification[END_REF][START_REF] Soltis | Polyploidy and genome evolution in plants[END_REF][START_REF] Mayrose | Recently formed polyploid plants diversify at lower rates[END_REF][START_REF] Kellog | Has the connection between polyploidy and diversification actually been tested?[END_REF]. Despite ongoing research regarding the distinction between the types of polyploids [START_REF] Levin | The Role of Chromosomal Change in Plant Evolution[END_REF][START_REF] Soltis | What we still don't know about polyploidy[END_REF][START_REF] Husband | The incidence of polyploidy in natural plant populations: major patterns and evolutionary processes[END_REF][START_REF] Doyle | Double trouble: taxonomy and definitions of polyploidy[END_REF], two main categories are generally recognized based on their origin: (i) autopolyploids that arise within a species, via intraspecific hybridization and duplication of similar genomes (homologous) and (ii) allopolyploids formed by interspecific hybridization and chromosome doubling of genomes that are more or less divergent (homeologous). The prevalence of different types of polyploids in nature has been intensively discussed [START_REF] Müntzing | The evolutionary significance of autopolyploidy[END_REF][START_REF] Stebbins | Types of polyploids: their classification and significance[END_REF][START_REF] Lewis | Polyploidy in species populations[END_REF][START_REF] Parisod | Evolutionary consequences of autopolyploidy[END_REF], and recent studies suggest a parity in the incidence of autopolyploidy and allopolyploidy [START_REF] Barker | On the relative abun-dance of autopolyploids and allopolyploids[END_REF] but see [START_REF] Doyle | Double trouble: taxonomy and definitions of polyploidy[END_REF]. The differentiation between these processes is fundamental to evaluate the importance of polyploidization and hybridization in plant evolution. In this context, the diploid-polyploid complex Veronica subsection Pentasepalae is an excellent model to gain deeper insights into the contribution of these mechanisms to the evolutionary history of angiosperms.

The aim of this study is to clarify the phylogenetic relationships of V. subsection Pentasepalae by analyzing the nuclear genome using Amplified Fragment Length Polymorphism (AFLP). In addition to its use in phylogeographic studies, the AFLP technique is now widely used to infer phylogenetic relationships and to identify hybridization and polyploid events in recently evolved polyploid non-model groups [START_REF] Meudt | Amplified fragment length polymorphism data reveal a history of auto-and allopolyploidy in New Zealand endemic species of Plantago (Plantaginaceae): new perspectives on a taxonomically challenging group[END_REF][START_REF] Reberning | The evolutionary history of the whiterayed species of Melampodium (Asteraceae) involved multiple cycles of hybridization and polyploidization[END_REF][START_REF] Himmelreich | Phylogenetic relationships in the extreme polyploid complex of the New Zealand genus Leptinella (Compositae: Anthemideae) based on AFLP data[END_REF][START_REF] Zozomová-Lihová | Taxonomy and evolutionary history of Alyssum montanum (Brassicaceae) and related taxa in southwestern Europe and Morocco: diversification driven by polyploidy, geographic and ecological isolation[END_REF]. Compared to the previous study by [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF], in addition to using AFLPs, a molecular tool for which markers are distributed throughout the genome, we expand the study to include for the first time all taxa currently recognized within the subsection. Also, we added individuals that are difficult to identify to species and of different ploidy level from miXed-ploidy populations. AFLPs were generated to address the following specific points: (i) The role of auto-and allopolyploidization processes in the evolutionary history of V. subsection Pentasepalae; (ii) The implications of these processes in species delimitation and classification; and (iii) A review of the taxonomic status of the taxa within V. subsection Pentasepalae.

Materials and methods

Plant material

Samples were collected in the field during 2009-2015 except for one population of V. satureiifolia and one population of V. tenuifolia subsp. fontqueri that were collected in 2002. Localities, initial taxonomic assignment, and further information about samples are given in Table S1. Fresh leaf material was collected and stored in silica gel. For V. krylovii, three individuals were included of which two were selected from herbarium material and one from a cultivated specimen in the Botanical Garden of Oldenburg (Germany). Veronica orientalis, which belongs to V. subsection Orientales (Wulff) Stroh of V. subgenus Pentasepalae was chosen as outgroup. The complete data-set comprises 243 individuals from 95 populations (outgroup included) covering the geographic distribution of each taxon (Fig. 1). From each location, 2-3 individuals were included, except for populations with miXed-ploidy levels. In these exceptional cases, two individuals per cytotype were analyzed whenever possible. Initial plant identification was based on the most recent taxonomic treatment (Rojas-Andrés and Martínez-Ortega, 2016), with the exception of some taxa that were not recognized by those authors, but whose names are used here to test their status [i.e., V. crinita f. bosniaca and V. thracica were included under the variation of V. crinita, and V. macrodonta under the variation of V. austriaca subsp. austriaca in that taxonomic treatment]. Material that was difficult to identify was initially catalogued using morphological affinity to other species (e.g., V. affinis linearis). Additionally, V. austriaca subsp. jacquinii/V. orbiculata indicates individuals of intermediate morphology between both species. Vouchers are deposited in the herbaria ALTB, GDA, MGC, OLD, SALA, VANF and WU (herbarium acronyms follow Thiers, 2017).

DNA-ploidy level estimation using flow cytometry

DNA-ploidy levels were estimated by flow cytometry using silica gel dried leaves. Individuals from each sampled population were measured separately. Nuclear suspensions were prepared following the method described by [START_REF] Galbraith | Rapid flow cytometric analysis of the cell cycle in intact plant tissues[END_REF] in which leaf tissue of each in-dividual was chopped together with leaf tissue from an internal stan-dard using a sharp razor blade in a Petri dish containing a buffer so-lution, namely Woody Plant Buffer (WPB; [START_REF] Loureiro | Two new nuclear isolation buffers for plant DNA flow cytometry: a test with 37 species[END_REF]. * * Solanum pseudocapsicum L. (2C = 2.589 pg; [START_REF] Temsch | Genome size in liverworts[END_REF], Zea mays L. 'CE-777' (2C = 5.43 pg; [START_REF] Lysak | Estimation of nuclear DNA content in Sesleria (Poaceae)[END_REF], Pisum sa-tivum L. 'Ctirad' (2C = 9.09 pg; [START_REF] Doležel | Plant genome size estimation by flow cytometry: inter-laboratory comparison[END_REF] and Pisum sativum L. 'Kleine Rheinländerin' (2C = 8.84 pg; [START_REF] Greilhuber | Genome size variation in Pisum sativum[END_REF] were used as internal standards depending on the C-value and standard availability. The suspension of isolated nuclei was filtered through a 48 µm nylon mesh, incubated with RNase to degrade double stranded RNA and stained with a saturating solution of propidium iodide fol-lowing [START_REF] Loureiro | Two new nuclear isolation buffers for plant DNA flow cytometry: a test with 37 species[END_REF] and [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF]. For each individual, one run of 5000 counts was made on a CyFlow SL (Partec GmbH, Münster, Germany; equipped with a 488 nm solid-state laser) or a CyFlow Space (Partec GmbH, Münster, Germany; equipped with a 532 nm solid-state laser). Results were acquired using Partec FloMax software v2.4d (Partec GmbH, Münster, Germany). A proXy of the ho-loploid genome size (2C) was calculated as follows: Veronica 2C nuclear DNA content (pg) = (Veronica G1 peak mean/internal standard G1 peak mean) * genome size of the internal standard. The DNA-ploidy level was estimated for each sample based on the values of the genome size proXy and the available chromosome counts for the studied species [START_REF] Martínez-Ortega | Species boundaries and phylogeographic patterns in cryptic taxa inferred from AFLP markers: Veronica subgen. Pentasepalae (Scrophulariaceae) in the western Mediterranean[END_REF][START_REF] Albach | Chromosome numbers in Veroniceae (Plantaginaceae): review and several new counts[END_REF].

DNA extraction and AFLP genotyping

Total genomic DNA was extracted from silica gel dried material following the CTAB protocol of [START_REF] Doyle | A rapid DNA isolation procedure for small quantities of fresh leaf tissue[END_REF]. The quality of the extracted DNA was checked on 1% TAE-agarose gels and the amount of DNA was estimated using a Nanodrop 2000C Spectrophotometer (Thermo Scientific). All extractions are stored at -80 °C at the "Biobanco de ADN Vegetal" (University of Salamanca, Spain). The AFLP procedure followed the method described by [START_REF] Vos | AFLP: a new technique for DNA fingerprinting[END_REF] with slight modifications. Genomic DNA (ca. 100 ng) was digested with MseI (New England BioLabs) and EcoRI (Fermentas) and ligated to double-stranded adaptors with T4 DNA-Ligase (Thermo Scientific) in a single restriction-ligation reaction for 3 h at 37 °C. Products were di-luted and pre-amplified using primers EcoRI-A (5′ GAC TGC GTA CCA ATT CA -3′) and MseI-C (5′ GAT GAG TCC TGA GTA AC -3′). Taq DNA Polymerase (BIOTOOLS BandM Labs. S.A) was used in the following PCR conditions: 2 min at 72 °C, 29 cycles of 30 s at 94 °C, 30 s at 56 °C and 2 min at 72 °C, and a final extension of 10 min at 72 °C. At this step, the preselective amplified fragments were visualized on 1% TBE-agarose gel. After dilution of pre-selective products, selective amplifi-cations were performed with the following PCR profile: 10 min at 95 °C, 13 cycles of 30 s at 94 °C, 1 min at 65 °C (decreasing 0.7 °C in each cycle) and 2 min at 72 °C, followed by 24 cycles of 30 s at 94 °C, 1 min at 56 °C and 2 min at 72 °C, with a final extension of 10 min at 72 °C. All PCR reactions were performed on an Eppendorf-Mastercycler-Pro thermocycler. Twelve individuals from ten different taxa representing the whole diversity of the final dataset were used to screen a total of eight different combinations of selective primers. Four primer combi-nations were finally selected (Table S2) based on the number and clarity of the peaks, and the polymorphism level observed among in-dividuals, which was checked to be sufficiently variable (i.e., overall genetic similarity among individuals from the same population was higher than that found among individuals from different populations, and much higher than the similarity detected among individuals from different taxa). Final selective PCR products were multiplexed for genotyping using the internal GeneScan 500 LIZ Size Standard (Applied Biosystems) in a multi-capillary sequencer ABI Prism 3730 (Applied Biosystems). Negative controls were run at each step of the process and 4.5% (=11) of the samples were replicated in each independent run of PCR from the same extracted DNA to assess genotyping errors [START_REF] Bonin | How to track and assess genotyping errors in population genetic studies[END_REF][START_REF] Bonin | Statistical analysis of amplified fragment length polymorphism data: a toolboX for molecular ecologists and evolutionists[END_REF][START_REF] Pompanon | Genotyping errors: causes, consequences and solutions[END_REF]. Final error rate was esti-mated after automated scoring according to [START_REF] Bonin | How to track and assess genotyping errors in population genetic studies[END_REF] by comparing the 1/0 matrices obtained for the replicated samples. Dif-ferences detected here could be due either to technical causes and/or to the automated scoring process. The degree of reproducibility of the data set was also investigated by analyzing the placement of replicates in a Neighbor-Joining tree.

Optimization procedure of automated AFLP scoring

Two different protocols were tested for the optimization of scoring parameters: the protocol of [START_REF] Holland | Optimizing automated AFLP scoring parameters to improve phylogenetic resolution[END_REF] and the open-source software optiFLP version 1.54 developed by [START_REF] Arthofer | OptiFLP: software for automated optimization of amplified fragment length polymorphism scoring parameters[END_REF]. The results obtained with these methodologies did not show incongruence or significant differences between them. OptiFLP was chosen for our analyses because of its greater flexibility, faster analysis and the possibility to run the program with its "unsupervised mode", which uses phylogenetic tree's robustness to find settings that maximize the differences between groups of profiles. To use the software designed by [START_REF] Arthofer | OptiFLP: software for automated optimization of amplified fragment length polymorphism scoring parameters[END_REF], electropherograms were first visualized in the software PeakScanner v.1.0 (Applied Biosystems) with all default settings except for a "light smoothing". Samples with poor quality profiles were discarded and AFLP data were exported to the open-source software optiFLP v.1.54 for the optimization of scoring parameters. Subsequently, fragments were automatically scored with tinyFLP v.1.30 [START_REF] Arthofer | TinyFLP and tinyCAT: software for automatic peak selection and scoring of AFLP data tables[END_REF] using the parameters optimized by optiFLP software (Table S3) and the data matrices from the different markers were concatenated using tinyCAT v.1.2 [START_REF] Arthofer | TinyFLP and tinyCAT: software for automatic peak selection and scoring of AFLP data tables[END_REF]. A single scoring procedure was run to create data matrices to be used in subsequent genetic structure analyses.

Phylogenetic analyses

AFLP data were analyzed in a phenetic framework (i.e., distance based clustering), due to the limitations reported for alternative methods commonly used for phylogenetic reconstruction [START_REF] Albach | Amplified fragment length polymorphisms and sequence data in the phylogenetic analysis of polyploids: multiple origins of Veronica cymbalaria (Plantaginaceae)[END_REF][START_REF] Himmelreich | Phylogenetic relationships in the extreme polyploid complex of the New Zealand genus Leptinella (Compositae: Anthemideae) based on AFLP data[END_REF]. With the aim of understanding the phylogenetic relationships between closely related taxa of the complex, and investigating the possible occurrence of hybridization and/or polyploidization processes, a Neighbor-Net was calculated based on Jaccard distances using SplitsTree4 v.4.13.1 [START_REF] Huson | Application of phylogenetic networks in evolutionary studies[END_REF]. Additionaly, Neighbor-Joining (NJ) trees based on Jaccard and Nei-Li distances were also built using Split-sTree4 v.4.13.1 and PAUP * 4.0b10 [START_REF] Swofford | PAUP * : phylogenetic analysis using parsimony ( * and other methods)[END_REF], respectively, to assess the influence of the distance measure on the results. Bootstrap values (1000 replicates) from the NJ tree based on Jaccard distance were transferred to the Neighbor-Net graph.

Genetic structure analyses

The genetic structure was investigated in the entire AFLP dataset, as well as in data subsets using the same conditions. Data subsets were obtained from the partition of the entire dataset according to the four main K = 4 clusters identified during the initial analysis. Since we are not able to corroborate if the populations under study follow the Hardy-Weinberg model, the genetic structure was initially investigated using two different approaches: non-hierarchical K-means clustering [START_REF] Hartigan | Algorithm AS 136: a k-means clustering algorithm[END_REF], which does not assume Hardy-Weinberg equilibrium, and Bayesian clustering analysis based on the MCMC algorithm using Structure v.2.3.4 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. Non-hierarchical K-means clustering was performed using the R script of [START_REF] Arrigo | Origin and expansion of the allotetraploid Aegilops geniculata, a wild relative of wheat[END_REF]. Numbers of K from 1 to 21 were tested and at least twenty independent runs starting from random seeds were performed for each K. To determine the optimal number of genetic clusters, the method of [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] was followed as adapted in [START_REF] Arrigo | Origin and expansion of the allotetraploid Aegilops geniculata, a wild relative of wheat[END_REF]. Bayesian clustering analyses were performed in Structure adopting an admiXture model and assuming correlated allele frequencies among populations [START_REF] Falush | Inference of population structure: extensions to linked loci and correlated allele frequencies[END_REF] according to a methodology for dominant markers [START_REF] Falush | Inference of population structure using multilocus genotype data: dominant markers and null alleles[END_REF]. Twenty replicates were run for each K from 1 to 21 with a burn-in length of 100,000 generations followed by 1,000,000 additional sampled generations. Structure analyses were run on the computer cluster developed by the UMS 2700 OMSI at the MNHN (Muséum National d'Histoire Naturelle, Paris). The optimal K value was determined using Structure Harvester [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF] following the method of [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF]. The output files were exported to CLUMPP v.1.1.2b [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF] to perform an alignment of cluster assignments across the replicate analyses that we visualized afterwards using Distruct v.1.1 [START_REF] Rosenberg | DISTRUCT: a program for the graphical display of population structure[END_REF]. The results obtained with both approaches were independently displayed on a Principal Coordinates Analysis (PCoA) [START_REF] Krzanowski | Principles of Multivariate Analysis[END_REF] based on the Jaccard distance index using NTSYSpc 2.2 (EXeter Software, Setauket, NY; Rohlf, 2005). The percentages of variance explained by the two different clustering methods were also compared by AMOVA analyses [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data[END_REF] performed in Arlequin v3.5 [START_REF] Excoffier | Arlequin (version 3.0): an integrated software package for population genetics data analysis[END_REF][START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under LinuX and Windows[END_REF]. Furthermore, PCO-MC (principal coordinate-modal clustering; [START_REF] Reeves | Accurate inference of subtle population structure (and other genetic discontinuities) using principal coordinates[END_REF] was implemented to test the significance of clusters found in PCoA using PCO-MC software (https://www.ars.usda. gov/plains-area/fort-collins-co/center-for-agricultural-resourcesresearch/plant-germplasm-preservation-research/docs/reeves-pco-mc/). The P-value cutoff was set to 0.9999 and the stability cutoff to 15% to maximize sensitivity to subtle population structure while minimizing type I error (Reeves andRichards, 2009, 2010).

Results

DNA-ploidy level determination

DNA-ploidy level estimations according to flow cytometric analyses are shown in Table S1 and Fig. 1. Ploidy was determined for most samples (94%), but for a few (6%) this was problematic likely due to the age of leaf material. In general, our results were in accordance with previous data with the group harboring diploids (2x), tetraploids (4x), hexaploids (6x) and octoploids (8x). Heterogeneity in DNA-ploidy level within a taxon was found only for V. austriaca subsp. austriaca (4x, 6x), V. austriaca subsp. jacquinii (2x, 6x), V. orbiculata (2x, 4x) and V. rosea (2x, 4x). In most cases, only one DNA-ploidy level was observed per population, with the exception of one population of V. rosea in Algeria (2x, 4x), and two populations of V. orbiculata (2x, 4x) in Bosnia and Herzegovina and Croatia, where two DNA-ploidy levels were observed within the same populations. All populations initially determined as V. austriaca subsp. jacquinii/V. orbiculata were tetraploid except for one population from Montenegro (pop. [START_REF] Holland | Optimizing automated AFLP scoring parameters to improve phylogenetic resolution[END_REF], which was shown to be diploid and was finally ascribed to a new species which is described here (i.e., V. dalmatica sp. nov., see Section 5). Similarly, a population from Bosnia and Herzegovina initially identified as V. affinis jacquinii (finally ascribed to V. dalmatica) was confirmed to be diploid (pop. 20), as well as one population of V. affinis linearis from FYROM (pop. 42, with an a posteriori identification as V. linearis). Most of the individuals initially catalogued as V. affinis kindlii (pop. 32, 33, 35) or V. affinis orsiniana (pop. 51, 52, 53, 54, 55) were diploids, except for two populations labeled as V. affinis kindlii (from Greece and Montenegro, pop. 31 and 34 respectively) that were found to be tetraploids.

Automated scoring of the AFLP data and degree of reproducibility

A total of 1127 polymorphic fragments were scored with the software tinyFLP (Table S2). The error rate per locus obtained for our final Table 1 Overview of the boostrap (BS) values detected by different phylogenetic analyses supporting each taxon included in the study. Taxa are shown according to the four groups identified by the placement of taxa in the Neighbor-Net network. Final taxonomic assignments, DNA-ploidy level and classification of individuals in clusters using different methodologies are indicated. Species recovered as independent clusters in PCO-MC analyses are indicated by a checkmark (✓).

Neighbor-Net Groups

Taxa

Ploidy BS values Clustering

NJ Jaccard NJ Nei-Li K-means K = 2 Structure K = 3 Structure K = 4 PCO-MC GROUP I V. kindlii
A C ─ V. crinita Kit 2x 100.0 100.0 1 A C ─ V. crinita (=V. crinita f. bosniaca) 2x 62.9 61.0 1 A C ─ V. thracica Velen. 2x 100.0 100.0 1 A C ─ V. krylovii Schischk. 2x 100.0 100.0 1 A C ─ V. orbiculata A. Kern 2x, 4x 73.0 78.0 2 A C ✓ V. prostrata L. 2x 100.0 100.0 2 A C ✓ V. turrilliana Stoj. & Stef. 2x 100.0 100.0 2 A C ✓ V. affinis kindlii 2x, 4x
2 C D ─ 2 C D ─ 2 C D ─ 1 A C ─ 1 C D ─ V. sennenii (Pau) M.M. Mart. Ort. & E. Rico 8x V. angustifolia Vahl (Bernh.) 8x V. teucrium L. 8x
data-set optimized with the optiFLP software was on average 2.55%. In the NJ analyses, siX of the eleven replicated samples were placed with their respective original samples, with a bootstrap value > 98%. The other five replicated samples were recovered at least in the same cluster as their respective original samples and others of the same population.

Phylogenetic reconstruction

Phylogenies reconstructed with different distance methods were congruent with one another and supported the monophyly of most of the species previously recognized by Rojas-Andrés and Martínez-Ortega (2016), with high bootstrap values (Table 1). However, since only one population of V. rhodopea was included in our study, monophyly of this particular species could not be confirmed. On the other hand, most of the internal nodes of the NJ trees were not supported by bootstrapping (Fig. S1), and the Neighbor-Net (Fig. 2) showed a high degree of re-ticulation existing in the group. Nevertheless, four main groups (I, II, III and IV) were identified according to the placement of individuals in the network and are presented below.

Group I. This low-supported group comprised five monophyletic diploid species: V. kindlii (BS = 100), V. linearis (BS = 99.9), V. or-siniana (BS = 100), V. rhodopea (BS = 100), and V. teucrioides (BS = 100). Diploid individuals from the south of Italy (ind. 135-140) of uncertain taxonomic identity, which are morphologically similar to V. orsiniana (V. affinis orsiniana), were recovered as monophyletic together with V. kindlii. One tetraploid population initially determined as V. affinis kindlii (ind. 80-82), was recovered as an independent lineage (BS = 100).

Group II. The monophyly of this group was clearly supported (BS = 99). It comprised three species recovered as monophyletic with bootstrap values of 100%. The diploid V. tenuifolia with three sub-species [subsp. tenuifolia, subsp. javalambrensis, subsp. fontqueri] and the tetraploid V. aragonensis are endemic to the Iberian Peninsula. The third species V. rosea, mostly diploid but comprising some tetraploid individuals in a single population, is endemic to North Africa.

Group III. This group included mostly diploid species with highly supported monophyly but there was low support for it as a whole. First, V. krylovii (BS = 100), one of the species representing the subsection in Siberia and Kazakhstan, was recovered as a strongly supported clade. Second, V. prostrata from central Europe was found to be also well resolved (BS = 100), as well as V. turrilliana (BS = 100), an endemic taxon from the border region of Bulgaria and Turkey. Individuals identified as V. crinita were recovered in two distinct lineages (BS = 100). Finally, this group comprised diploid individuals of V. austriaca subsp. jacquinii (BS = 100), and miXed cytotypes (2x, 4x) of V. orbiculata (BS = 71.7), as well as tetraploid individuals of intermediate morphology recorded as V. austriaca subsp. jacquinii/V. orbiculata.

Group IV. Following the initial taxonomic classification, this group included four polyploid taxa: (i) tetraploid and hexaploid cytotypes of V. austriaca comprising three subspecies: subsp. austriaca, subsp. dentata, and subsp. jacquinii; (ii) the tetraploid species V. satureiifolia; (iii) the octoploid V. sennenii, endemic from the north of Spain; and (iv) octoploid individuals of V. teucrium var. teucrium, and var. angustifolia. The monophyly of this group was well supported (BS = 93.9) but not the monophyly of most of the species within it. Phylogenetic analysis (Fig. 2) only supported the monophyly of hexaploid individuals of V. austriaca subsp. jacquinii but with a very low bootstrap value (BS = 59.6). Veronica satureiifolia and V. sennenii were recovered together with octoploid individuals identified as V. teucrium subsp. an-gustifolia (BS = 67.5).

Genetic structure

Following the method implemented in Structure Harvester, Bayesian clustering analysis supported an optimal partition of the subsection in three clusters. On the contrary, non-hierarchical K-means clustering analysis of the same dataset estimated K = 2 as the most likely number of genetic clusters. However, PCoA (Fig. 3) and AMOVA analyses (Table 2) demonstrated that the clustering proposed by Structure explained a higher percentage of the variance among groups than K-means (see Table 2). Accordingly, we here focus on results of Bayesian clustering. High levels of admiXture were found in Bayesian clustering analyses performed at higher values of K (Fig. S2C). It should be pointed out that most taxa included in our analyses (with the ex-ception of polyploids from group IV and of V. teucrioides from group I) were recovered as independent clusters when Structure analyses were performed at K = 20 (Fig. S2C). In addition, an exclusive cluster was found grouping tetraploid populations of V. austriaca subsp. jacquinii/V. orbiculata and V. orbiculata.

The clusters revealed by Structure at K = 4 (Figs. 2 andS2A) generally concurred with the groups identified by the Neighbor-Net with the only exception of V. orsiniana (Table 1), and partially corresponded with geographic regions: cluster A included a group of narrow endemics mostly restricted to the south of the Balkan Peninsula; cluster B com-prised the three well recognized species from the Iberian Peninsula and North Africa; cluster C recovered most diploids from the Balkan Pe-ninsula, V. krylovii from Russia, and V. orsiniana; cluster D included all the polyploid taxa mainly from central Europe and north of Spain.

Additional Bayesian clustering analyses within clusters A, B, C and D estimated an optimal K = 5, K = 3, K = 3 and K = 2, respectively (Fig. S2B). According to the results obtained for cluster A, the four diploid species and the tetraploid population of V. affinis kindlii from Mt. Vermion (pop. 31) were recovered in independent clusters. In cluster B, the three species from the Ibero-North African group (V. aragonensis, V. rosea and V. tenuifolia) were recovered in independent and homogeneous clusters almost without admiXture among them. When analyses were performed within cluster C, one cluster grouped a single species (i.e., V. orsiniana) and the other two clusters divided the taxa from group III in two subgroups. One subgroup included diploids of V. austriaca subsp. jacquinii and V. orbiculata and all intermediate populations displaying high levels of admiXture [V. austriaca subsp. jacquinii/V. orbiculata, V. crinita (=V. crinita f. bosniaca) and V. affinis kindlii]. Another subgroup, also with a certain degree of admiXture, comprised the remaining diploid species from group III [V. crinita, V. crinita (=V. thracica), V. turrilliana, V. prostrata, V. krylovii]. Within cluster D, there was an obvious geographic pattern in which polyploids form two separate clusters, however, with many individuals forming a continuous gradation of proportion scores between the two clusters (Fig. 4). Hexaploid individuals classified as V. austriaca subsp. jacquinii (pop. 11-16) were included in one cluster with a very high proportion score (> 0.99; data not shown). By contrast, tetraploid individuals identified as V. satureiifolia, octoploid populations assigned to V. senneni and affinis individuals, had a proportion score > 0.99 (data not shown) to be defined in a second cluster. Most individuals of V. teucrium var. angustifolia showed a high genetic affinity to this second cluster with low levels of admiXture.

PCO-MC analysis recovered ten species as significant independent clusters: V. aragonensis, V. dalmatica, V. kindlii, V. linearis, V. orbiculata, V. orsiniana, V. prostrata, V. rosea, V. tenuifolia, and V. turrilliana (Table 1).

Discussion

The importance of auto-and allopolyploidization in the subsection and the recurrent formation of polyploids

The diversity of cytotypes and the existence of miXed-ploidy levels within species and populations in the group reveal that polyploidization has occurred likely continuously since the origin of the subsection ca. 2.8 Mya [START_REF] Meudt | Is genome downsizing associated with diversification in polyploid lineages of Veronica?[END_REF]. The pattern of reticulation shown by the Neighbor-Net (Fig. 2) and the high levels of admiXture found by Structure suggest that V. subsection Pentasepalae is composed of species S1).

that are in the initial stages of divergence. Furthermore, based on morphological and (phylo-)genetic intermediacy between potential parental species, hybridization is confirmed in V. subsection Pentasepalae. In addition, ILS cannot be excluded as a cause of the lack of resolution observed for internal nodes of the Neighbor-Net and NJ trees. Nevertheless, phylogenetic analyses demonstrate that polyploid taxa distributed mainly in central Europe (group IV) constitute a very well supported group (Fig. 2). It should be pointed out that the higher number of AFLP fragments present in polyploid individuals may produce a bias in posterior analyses towards the apparent monophyly of most of the polyploids. However, this artificial grouping due to a higher number of AFLP fragments in polyploids can be discarded in our dataset because other polyploid species of the subsection are not recovered within group IV (e.g., V. aragonensis and V. orbiculata). Thus, our study could indicate a common origin of polyploid entities from group IV at the tetraploid level. This hypothesis was previously rejected by [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF] due to the existence of morphological variation among polyploids and to the fact that most polyploid species have a polytopic origin [START_REF] Soltis | Polyploidy: recurrent formation and genome evolution[END_REF]. However, a possible explanation is that hexa-and octoploids have emerged (probably several times independently within each lineage) after a previous differentia-tion of lineages at the tetraploid level. Furthermore, the extinction of diploid or some of the tetraploid ancestors within group IV is also likely, which together with the limitations of the available methodologies hampers the obtention of ancestor-derivative patterns within group IV [START_REF] Stebbins | Chromosomal Evolution in Higher Plants[END_REF][START_REF] Mcdade | Hybrids and phylogenetic systematics II. The impact of hybrids on cladistic analysis[END_REF][START_REF] Buggs | The le-gacy of diploid progenitors in allopolyploid gene expression patterns[END_REF].

Our results suggest that polyploid species in the subsection may have emerged by different processes. Whereas autopolyploidization appears to be the main evolutionary force for some taxa, allopolyploidization also seems to be common. Evidence for autopolyploidization is found, for example, in group II. Tetraploid individuals have been found in a population of V. rosea from Algeria (Table S1, pop. 63), which cluster together with the rest of diploid individuals of the species (Figs. 2 andS2), thus suggesting a recent autopolyploid event occurring within the population.

Analyses further point to an autopolyploid origin of the octoploid V. sennenii from the tetraploid V. satureiifolia. The individuals belonging to V. satureiifolia and V. sennenii are recovered in the same group in the Neighbor-Net without any clear separation between individuals of both species, and they form a homogeneous cluster in the Bayesian clustering analyses (Fig. 4). Flow cytometric analyses (Table S1) have confirmed that both species (and consequently, both ploidy levels) grow in sympatry in the province of Huesca, in the north of Spain (pop. 69,4x;and pop. 73,8x). On the basis of all these results, this is another example of autopolyploid speciation in natural populations [START_REF] Soltis | Autopolyploidy in angiosperms: have we grossly underestimated the number of species?[END_REF]. Likewise, according to the Neighbor-Net, most octoploid individuals determined as V. teucrium var. angustifolia are nested with V. satureiifolia and V. sennenii and have very similar genetic composition in clustering analyses (Fig. 4). Furthermore, one population of V. satureiifolia (pop. 70, 4x) and one of V. teucrium var. angustifolia (pop. 87, 8x) have been found in close proXimity (about 400 m) in the region of Île-de-France. Veronica satureiifolia and V. teucrium var. angustifolia were also shown to share the same cpDNA haplotype [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF]. Consequently, a plausible interpretation is that multiple autopolyploidization events might have occurred in the tetraploid V. satureiifolia giving rise to octoploids that have been identified as V. sennenii in the Iberian Peninsula and V. teucrium var. angustifolia in France. Alternatively, a past continuous distribution area of the octoploid entity and a subsequent fragmentation scenario cannot be discarded, although it seems unlikely considering the distance of 500 km between the French southernmost and the Spanish northernmost populations, and the existence of the Pyrenees in between.

Our results also reveal that at least two episodes of polyploidy have occurred within V. orbiculata (Fig. 2; Group III), although the processes seem to be more complex. Autopolyploid formation has been detected in one population from Bosnia and Herzegovina (Table S1; pop. 43) and might be occurring in other, not surveyed, populations. Tetraploid individuals found in this population (ind. 107, 108) are nested within the diploid individuals (ind. 109-116) with a BS value of 100% (Fig. 2). Moreover, individuals of both cytotypes are recovered together as a significant cluster in PCO-MC and Structure analyses. In contrast, tetraploid individuals (ind. 117-119) from another population in Croatia (Table S1; pop. 45) are recovered well separated from the diploids of the same population (ind. 114-116; BS = 99.5) and are not included together within any significant cluster in PCO-MC analyses. Moreover, an exclusive cluster is found in Structure analyses (Fig. S2C) that groups these tetraploid individuals with tetraploid populations recorded as V. austriaca subsp. jacquinii/V. orbiculata. Thus, these tetraploids are probably the result of an allopolyploidization event. Consequently, V. orbiculata is a further example of a diploid-polyploid species with numerous independent origins of polyploid entities as shown in other species [START_REF] Soltis | Polyploidy: recurrent formation and genome evolution[END_REF][START_REF] Bardy | Disentangling phylogeography, polyploid evolution and taxonomy of a woodland herb (Veronica chamaedrys group, Plantaginaceae s.l.) in southeastern Europe[END_REF][START_REF] Bardy | EXtensive gene flow blurs species boundaries among Veronica barrelieri, V. orchidea and V. spicata (Plantaginaceae) in southeastern Europe[END_REF].

There are other strong arguments of recurrent formation of allopolyploids within group III. Specifically, tetraploid individuals from Bosnia and Herzegovina (Table S1; ind. 42-47) recorded as V. austriaca subsp. jacquinii/V. orbiculata due to their transitional morphology, are recovered by the NJ in an intermediate position between these species (Fig. S1). Thus, a hybrid origin of these populations is suggested with diploid V. austriaca subsp. jacquinii and V. orbiculata as putative parental species. In addition, diploid individuals located in Montenegro labeled as V. affinis kindlii (ind. 83-84), are also recovered in an intermediate position together with a population of V. crinita, putatively belonging to f. bosniaca (ind. 64-66). Furthermore, the position of individual 85 (2x) in the Neighbor-Net suggests that homoploid hybridization could be also an important evolutionary process occurring in the group, which is here demonstrated for the first time and of which V. × gundisalvi may represent an additional example [START_REF] Martínez-Ortega | Species boundaries and phylogeographic patterns in cryptic taxa inferred from AFLP markers: Veronica subgen. Pentasepalae (Scrophulariaceae) in the western Mediterranean[END_REF].

Last, hybridization and/or introgression events may have affected the tetraploid population of V. aff. kindlii located in Mt. Vermion (Greece; Fig. 2, group I, pop. 31). Bayesian clustering analyses showed high levels of admiXture with the polyploid group IV (Fig. S2A). Veronica austriaca subsp. jacquinii is the only species from group IV distributed in this southern area of the Balkan Peninsula. Thus, the position of population 31 in the Neighbor-Net could be influenced by hybridization and/or introgression processes involving V. austriaca subsp. jacquinii and representatives from group I or its ancestors.

The challenging task of delimitating species within a recently diverged diploid-polyploid complex

Species delimitation within recently diverged plant complexes is currently a major challenge for systematists. In general, at this level of lineage separation, phenotypic differences among species may not be evident and a clear phylogenetic signal is not always obtained [START_REF] Federici | DNA barcoding to analyse taxonomically complex groups in plants: the case of Thymus (Lamiaceae)[END_REF]. Consequently, other characteristics (e.g., ploidy levels, differences in habitat, pollinators, phenology, etc.) are important lines of evidence when delimiting species in this recently diverged, phenotypically, and phylogenetically complex groups. This situation requires the adoption of the general lineage concept of species in which different species properties (that have been used as criteria under rival species concepts), serve as lines of evidence to assess lineage delimitation [START_REF] De Queiroz | Species concepts and species delimitation[END_REF].

Identifying biological diversity at the species level is even more challenging when processes such as polyploidy are involved in the evolution of a group. Polyploidy has long been considered a mechanism of direct sympatric speciation [START_REF] Otto | Polyploid incidence and evolution[END_REF][START_REF] Schemske | Understanding the origin of species[END_REF][START_REF] Rieseberg | Plant speciation[END_REF]. However, recent studies suggest that polyploid speciation is not necessarily an instantaneous process [START_REF] Husband | The incidence of polyploidy in natural plant populations: major patterns and evolutionary processes[END_REF]. The formation of unreduced gametes and other biological traits are fundamental in initial stages of polyploid emer-gence and establishment [START_REF] Rieseberg | Plant speciation[END_REF][START_REF] Fowler | Critical factors in the establishment of allopolyplois[END_REF]. Both, the rates of production of unreduced gametes and the successful long-term establishment and spread of new polyploid in-dividuals are affected by genetic and environmental factors [START_REF] Ramsey | Pathways, mechanisms and rates of polyploidy formation in flowering plants[END_REF][START_REF] Comai | The advantages and disadvantages of being polyploid[END_REF][START_REF] Lafon-Placette | Current plant speciation research: unravelling the processes and mechanisms behind the evolution of reproductive isolation barriers[END_REF]. Re-gardless of the timing of the process, it has been estimated that 15% of angiosperm speciation events, and even more in Veronica, are asso-ciated with a ploidy increase [START_REF] Albach | Chromosome numbers in Veroniceae (Plantaginaceae): review and several new counts[END_REF][START_REF] Wood | The frequency of polyploid speciation in vascular plants[END_REF]. Indeed, it has been corroborated the effect of genome duplication on countless features (e.g., reproductive biology, phenotype, physiology, geographical and environmental distributions of cytotypes, genetic, epigenetic and genomic consequences, and so forth; reviewed in [START_REF] Ramsey | Ecological studies of polyploidy in the 100 years following its discovery[END_REF]. Whether caused by ploidy per se, adapta-tion or founder effects and genetic drift, these changes may maintain polyploids as separately evolving metapopulation lineages from their parental taxa, which justify their treatment as separate taxonomic species [START_REF] De Queiroz | Species concepts and species delimitation[END_REF].

Our study has demonstrated that most of the species of V. subsection Pentasepalae are still in the initial stages of divergence. Moreover, autoand allopolyploids have been identified within the group. In this situation, the taxonomic status of diploid and polyploid taxa within V. subsection Pentasepalae is reviewed adopting the general species concept of de [START_REF] De Queiroz | Species concepts and species delimitation[END_REF] and making use of an integrative taxonomic approach. In our case study, no significant differences in habitat preference are observed, experimental data on reproductive biology are not available, and probably many species share pollinators to a great extent. Thus, we have based the decisions of species delimitation on ploidy level, phylogeny and genetic divergence, but also on information from morphology, distribution, ecology, etc., available in Rojas-Andrés and Martínez-Ortega (2016). Furthermore, a conservative approach to taxonomy is preferable when incongruences among different lines of evidence are found [START_REF] Carstens | How to fail at species delimitation[END_REF]. When such situation was encountered, we maintained the last taxonomic treatment of Rojas-Andrés and Martínez-Ortega (2016).

Moreover, we think that populations identified in this study, for which we have not obtained sufficient evidence to be delimited as species (e.g., tetraploid hybrid populations catalogued as V. austriaca subsp. jacquinii/V. orbiculata that cannot be identified as a different species but could potentially evolve as a distinct lineage) would have to be considered as functional units of biological diversity. In these cases, this recognition would help to address future ecological, evolutionary and taxonomic questions [START_REF] Ramsey | Ecological studies of polyploidy in the 100 years following its discovery[END_REF][START_REF] Laport | Out of one, many: the biodiversity considerations of polyploidy[END_REF].

Additionally, we consider that further molecular tools (i.e., mole-cular studies using neutral markers as SSRs; López-González et al., in prep.), morphological data (e.g., traits with potential impact on in-dividual fitness), and deeper ecological and biological information (e.g., environmental distribution analyses among cytotypes, crossing experiments to understand reproductive interactions) are needed to reevaluate whether the species rank is appropriate for some of these taxa (e.g., polyploids from group IV).

Taxonomic considerations

This study provides new insights into the systematics of the polyploid complex V. subsection Pentasepalae. Our analyses support 20 distinct species in the group. The most recent taxonomic treatment available (Rojas-Andrés and Martínez-Ortega, 2016) has been revised and updated (changes summarized in Table S1).

First, the individuals of V. austriaca subsp. jacquinii, included in this study are placed in two separate phylogenetic lineages differentiated by their ploidy levels (diploids vs. hexaploids) (Fig. 2). Monophyly of the diploid individuals, which represent an example of cryptic species in the subsection, is well supported by phylogenetic reconstruction (BS = 99.8), whereas hexaploids are recovered with a low bootstrap value (BS = 59.6). Additionally, PCO-MC and Bayesian clustering analysis recover these populations as significant and independent clusters (Table 1 and Fig. S2C). Indeed, after an exhaustive revision of herbarium specimens, morphological characters corresponding to each of these species have been found (see Section 5). In addition, the dis-tribution area of the diploid cytotypes is restricted to the Adriatic coast of Albania, Bosnia and Herzegovina, Croatia and Montenegro. Based on all these lines of evidence, we consider that diploid individuals of V. austriaca subsp. jacquinii should be recognized at the specific rank as V. dalmatica N.Pad.Gar., Rojas-Andrés, López-González and M.M.Mart.Ort (see Section 5).

Second, the analyses presented here allow the recognition of V. thracica at the species level. Veronica thracica was described by [START_REF] Velenovsky | Dritter Nachtrag zur Flora von Bulgarien[END_REF] to differentiate plants mainly occurring in Bulgaria characterized by white hairy stems and oval-obovate, deeply cordate, almost amplexicaulus leaves. This name was later combined under V. teucrium as subspecies [START_REF] Velenovsky | Flora bulgarica. Descriptio et enumeratio systematica plantarum vascularium in principatu Bulgariae sponte nascentium[END_REF] or variety [START_REF] Maly | Beiträge zur Kenntnis der illyrischen Flora[END_REF] and has been related to V. crinita by other authors (e.g., Watzl, 1910;[START_REF] Peev | Veronica L[END_REF]. These individuals were considered within the variation of V. crinita in the most recent taxonomic treatment due to their morphological similarities (Rojas-Andrés and Martínez-Ortega, 2016). Genetic data now provide evidence that these populations constitute an independent evolutionary lineage differentiated from typical V. crinita described from Hungary (Figs. 2 andS2C) and it represents an additional example of a cryptic species within the subsection [START_REF] Martínez-Ortega | Species boundaries and phylogeographic patterns in cryptic taxa inferred from AFLP markers: Veronica subgen. Pentasepalae (Scrophulariaceae) in the western Mediterranean[END_REF]. Furthermore, after examining the morphological characters of this material, we found that V. thracica has dense tomentose indument on leaves and stems, formed by patent to slightly incurvate hairs that confer a whitish (light green in vivo) color to the plant. In contrast, V. crinita has villous indument on leaves and stems, which is constituted by crooked, generally interwoven hairs that confer a brownish green color to the plant. The leaves are concolor (i.e., upper and lower leaf sides of the same color) and more densely pilose in V. thracica, while they are slightly bicolor (i.e., dark/dive green color of the upper leaf side vs. green color of the underside of the leaf) and comparatively not so densely pilose in V. crinita. We consider that there is enough evidence to recognize these lineages as separate species, and consequently, the recognition at the specific rank of the Bulgarian populations is proposed.

Finally, several taxa have been described under V. teucrium, but only two varieties have been recognized in the last taxonomic treatment of the subsection (Rojas-Andrés and Martínez-Ortega, 2016): V. teucrium var. teucrium and V. teucrium var. angustifolia. These two mostly allo-patric octoploid entities are morphologically distinct. Veronica teucrium var. teucrium is mainly distributed in Germany, Austria, and Bulgaria, while V. teucrium var. angustifolia occurs in France according to our sampling. Moreover, both varieties are recovered in different subclusters in our molecular analyses (Fig. 4). Based on all available data, we consider that both entities should be recognized at the specific level as V. teucrium L. and V. angustifolia (Vahl) Bernh., respectively. Additionally, apart from geographic differentiation, individuals of V. teucrium var. angustifolia are very similar to V. sennenii in size and appearance. Taken together these data would suggest that V. sennenii and V. teucrium var. angustifolia have the same parental origin, although they could have arisen from the same or different autopolyploid events. If they were considered synonyms, the name V. angustifolia (Vahl) Bernh. would prevail at the specific level, according to the principle of priority. But this taxonomic decision should not be firmly adopted until additional exhaustive analyses including more populations of these taxa are performed.

Additionally, the importance of some populations from the south of Italy identified in Flora d'Italia as V. austriaca and their relationship with those from the Balkan Peninsula have been previously highlighted [START_REF] Fischer | Veronica L[END_REF]. However, the identity of these plants (ind. 135-140; labeled in this study as V. affinis orsiniana) has remained unclear for many years. Our analyses confirm the identity of these plants as V. kindlii and show that this species should be considered independent from V. orsiniana or V. austriaca (Table 1, Fig. 2). The name V. kindlii has recently been resurrected to designate those populations from the Balkan Peninsula, which were previously known as V. orsiniana [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF]. Thus, a clear amphi-Adriatic distribution of V. kindlii is now demonstrated. Finally, there is no evidence that the individuals initially identified as V. affinis kindlii belong to V. kindlii. Unfortunately, the taxonomic status of these entities remains unresolved. Additional exhaustive field sampling in order to have a good representation of these unresolved entities and posterior molecular studies could shed some light on the taxonomic identity of these individuals.

Another important outcome is the corroboration of the genetic distinctiveness and monophyly of V. linearis, a diploid endemic species from FYROM that passed unnoticed for many years. This name did not appear in Floras or monographs of Veronica. The plant was initially described as V. kindlii var. linearis by [START_REF] Bornmüller | Zur Flora von Mazedonien[END_REF] and has recently been elevated to the species level based on morphological evidence (Rojas-Andrés et al., 2016;Rojas-Andrés and Martínez-Ortega, 2016). According to our phylogenetic reconstruction (Fig. 2) and PCO-MC analyses (Table 1), V. linearis is recovered as monophyletic within group I and its closest relatives are V. kindlii and V. teucrioides (BS = 74.7 for [V. linearis + V. kindlii + V. teucrioides]). In addition, one population with dubious morphological characters labeled as V. affinis linearis (pop. 42) is recovered within V. linearis. Nevertheless, clustering analyses (Fig. S2B) showed introgression with V. teucrioides, which is in agreement with the dubious determination based on morphological characters.

With regard to the delimitation of varieties and subspecies, the theoretical framework behind their concept is less clear as it is for species. We have attempted to avoid the use of these ranks in the proposed taxonomic changes, but in two cases the data available are not conclusive and the subspecies rank has been retained: (i) Within group II, three subspecies are recognized under V. tenuifolia.

Their different distribution areas and the divergence found in the Neighbor-Net and NJ trees between populations corresponding to each subspecies could indicate reduced gene flow among them (Figs. 2 and S1), as previously showed by studies based on AFLP and morphology [START_REF] Martínez-Ortega | Species boundaries and phylogeographic patterns in cryptic taxa inferred from AFLP markers: Veronica subgen. Pentasepalae (Scrophulariaceae) in the western Mediterranean[END_REF][START_REF] Andrés-Sánchez | Combining traditional morphometrics and molecular markers in cryptic taxa: towards an updated integrative taxonomic treatment for Veronica subgenus Pentasepalae (Plantaginaceae sensu APG II) in the western Mediterranean[END_REF]. However, other analyses do not support the recogni-tion of the subspecies as independent clusters (see PCO-MC and genetic structure results in Table 1, Fig. S2). Likewise, phylogenetic analyses based on nuclear and plastid DNA sequences did not dif-ferentiate among the three subspecies currently recognized [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF]. Due to the incongruences found among dif-ferent sources of data, we suggest to maintain their current formal rank as subspecies. (ii) The subspecific rank has also been retained for some taxonomic entities belonging to group IV (i.e., three subspecies recognized under V. austriaca). The lack of resolution in our AFLP analyses (Table 1, Fig. 4) as well as in nuclear and plastid DNA trees [START_REF] Rojas-Andrés | EXploring the intricate evolutionary history of the diploid-polyploid complex Veronica subsection Pentasepalae (Plantaginaceae)[END_REF] do not support the recognition of current subspecies as different species, nor the unification in a single species. Unfortunately, the phylogenetic relationships among these taxa remain unresolved. These subspecies have been described in the last taxonomic treatment given the morphological and chorological differences among them (Rojas-Andrés and Martínez-Ortega, 2016). Consequently, we have retained the subspecies rank within V. austriaca (subsp. austriaca; subsp. dentata, subsp. jacquinii), at least until future studies clarify the evolutionary history and taxonomy of these polyploid entities.

Conclusions and description of a new species

The exhaustive sampling of V. subsection Pentasepalae, and the use of AFLP fingerprinting together with flow cytometry data provided new insights into the evolutionary history and species delimitation of a taxonomically complex plant group, in which auto-and allopolyploidization appear to be active evolutionary processes, even nowadays. Based on all sources of data currently available, V. subsection Pentasepalae contains at least 20 monophyletic species, five of them narrow endemics. This taxonomic framework is essential to design suitable conservation strategies. Future studies should focus on trying to understand in more detail the role that hybridization has played in the evolution of the subsection, and on ecological factors that make polyploidy so important in plant evolution and speciation. In ac-cordance with the data obtained, a new species is here formally de-scribed: Next, we provide a diagnosis between V. dalmatica and the morphologically closest taxa, as well as a full description of V. dalmatica, which is parallel to the descriptions provided by Rojas-Andrés and Martínez-Ortega (2016). The indument is described according to [START_REF] Beentje | The Kew Plant Glossary: An Illustrated Dictionary of Plant Terms[END_REF]. Two measurements given together refer always to length × width.

Diagnosis: V. dalmatica differs from V. austriaca subsp. jacquinii by its smaller plant size (10-16 vs. 25-50 cm), having shorter stem-hairs (0.3- 0.4 vs. 0.8-1.2 mm), smaller leaves (12-16 × 5-10 vs. 20-30 × 10-20 mm), shorter styles (3-5 vs. 4-7 mm) and tiny capsules with a less deep sinus (up to 0.5 vs. 1.0 mm). Attending to chromosome number, V. dalmatica is diploid (2n = 16) whereas V. austriaca subsp. jacquinii is frequently hexaploid [2n = (32), 48, (64), ( 80)].

Overlapping in ploidy level, V. dalmatica (2n = 16) is morphologically differentiated from V. orbiculata (2n = 16, 32) by having the eglandular hairs of the stem not arranged in two opposite lines along it. Apical shoot leaves are opposite in V. dalmatica, whereas in V. orbiculata they can be opposite, alternate or verticillate by three.

Description: Stems (6) 10-16 [START_REF] Otto | Polyploid incidence and evolution[END_REF] cm long, slightly ascending to erect, covered by eglandular hairs (0.18) 0.30-0.40 (0.81) mm long, incurvate, ± appressed and antrorse, not arranged in 2 opposite lines along the stem; apical shoot bearing (5) 8-11 (16) pairs of leaves. Leaves opposite, (8) 12-16 (20) × (3) 5-10 (16) mm; ovate, obovate, or narrowly to widely trullate; more or less rounded or cuneate at the base; pinnatifid to pinnatisect, with linear-lanceolate to narrowly el-liptic segments, variable in width, entire, revolute to subrevolute, subglabrous or pilose, covered by hairs (0.06) 0.10-0.18 (0.23) mm long, sessile to shortly petiolate. Basal leaves pinnatifid to pinnatisect, segments 0.4-1.0 mm wide; medium leaves (i.e., those situated in the central part of the stem) pinnatipartite to pinnatisect, segments 0.25-1.00 mm wide; uppermost leaves pinnatisect, rarely bipinnatifid, segments 0.25-0.60 mm wide. Leaves of the apical shoot opposite, linear to lanceolate, narrowly elliptic, entire, dentate-serrate to pinnatisect, revolute to subrevolute. Racemes axillary, opposite, exceptionally solitary, bearing (9) 20-40 (48) flowers, loosely to densely arranged; peduncles (2.5) 3.0-7.0 (11) cm long, covered by a non-glandular indument similar to that of the leaves; bracts (1.5) 3.0-5.0 (8.0) mm long, linear, entire, exceptionally pinnatifid to pinnatisect at the base with one or two segments, glabrous or subglabrous, covered by hairs similar to those covering the leaves; pedicels (1.6) 3.0-5.0 (8.5) mm long. Calyx (0.7) 2.0-4.0 (5.0) mm long, with (4) 5 sepals, linear-lanceolate, usually shorter than the capsule, glabrous or subglabrous. Corolla 9-15 mm in diameter, light or dark blue. Capsule (2.0) 3.0-5.0 (6.0) × (2.0) 3.0-4.5 (5.3) mm, glabrous, widely elliptic or widely obovate to very widely depressed ovate-obovate, rounded at the base, slightly emarginated or rounded at apex, sinus up to 0.5 (0.6) mm depth. Style (2.8) 3.0-5.0 (6.0) mm long. Seeds (0.9) 1.3-1.7 × 1.5-1.8 (2.0) mm, ca. 8 per capsule. Notes. The plant is illustrated in Rojas-Andrés and Martínez-Ortega, 2016; Fig. 3(f-i). Apical shoot is not represented. Specimens examined -See AppendiX A.

Fig. 1 .

 1 Fig. 1. Maps of sampling sites. Population codes follow Table S1, symbol shapes represent ploidy level (○ 2x; □ 4x; Δ 6x; ◊ 8x; ☆ miXed-ploidy populations), and colors indicate cluster affiliation. An asterisk (*) indicates missing data for DNA-ploidy level. (A and B) Locations of the 93 populations of Veronica subsect. Pentasepalae analyzed in this study. (C) Detailed distribution map of studied populations from the Balkan Peninsula.

Fig. 2 .

 2 Fig. 2. (A) Neighbor-Net network based on 1127 AFLP scored fragments of 241 individuals of Veronica subsect. Pentasepalae using Jaccard's genetic distances. Individual codes are shown following Table S1. Arcs delimits taxa whose initial taxonomic determination and ploidy are shown. Range of colors of the arcs differentiates the four clusters identified by the Structure analyses. Bootstrap values (BS) > 50% are shown. (B) Bayesian clustering analyses based on the entire AFLP dataset. Four main clusters from a STRUCTURE analyses with K = 4 are represented by different colors. Black lines separate different populations which are indicated below the graph (population codes follow TableS1).

Fig. 3 .

 3 Fig. 3. Principal Coordinate Analysis (PCoA) of the AFLP dataset of 241 individuals based on Jaccard's distances and DCENTER module. AXis 1 and AXis 2 explain 8.36% and 6.20% of the variation, respectively. (A) The two clusters supported by genetic structure analyses using K-means clustering algorithm are represented by colors. (B) Colors indicate the four genetic clusters from K = 4 estimated by Bayesian clustering analyses using Structure.

Fig. 4 .

 4 Fig. 4. Genetic structure analysis based on AFLP data of 58 individuals from 26 populations (corresponding with group IV/cluster D, which comprises most polyploid taxa from central Europe and north of Spain). (A) Bayesian model-based clustering at K = 2 using Structure. Colors represent different clusters (dark blue vs. light blue) and separate individuals of different populations, which are indicated below the graph (codes follow Table S1). Taxonomic names from population are shown above the graph. (B) Part of the Neighbor-Net (presented in Fig. 2) representing the 58 analyzed individuals. (C) Distribution map of populations included in genetic structure analysis of cluster D.

  Veronica dalmatica N.Pad.Gar., Rojas-Andrés, López-González & M.M.Mart.Ort., sp. nov. -Type: Holotype: Bosnia and Herzegovina, Republica Srpska: between Brgat and Trebinje, 42.68289N, 18.28949E, 283 m, 10/VI/2015, clearings in a forest of Carpinus betulus with Paliurus spina-christi. Leg. M. Martínez Ortega, X. Giráldez, N. Padilla and N. López, MO6119 (SALA 157047!).

  Chromosome Number -2n = 16 Habitat. -Dry and stony meadows, steppes, forest glades and shrublands, rocky slopes; usually on calcareous soils; (50) 200-1100 (1400) m above sea level. Distribution. -W Balkan Peninsula; Albania, Bosnia and Herzegovina, Croatia, Montenegro. Etymology -The epithet indicates geographical distribution of the species. Dalmatia is a historical region of the Adriatic Sea ranging from Rab (Croatia) to the Bay of Kotor (Montenegro) including a small area of Bosnia and Herzegovina.

Table 2

 2 Analysis of molecular variance (AMOVA) performed with different grouping approaches. Percentage of variation explained by different methodological groupings (K-means model, Neighbor-Net and Structure algorithm) are shown.

	Clustering approach	K value	Source of variation	Sum of squares	Variance components	Percentage of variation	Statistic
	K-means model	K = 2	Among clusters	569.761	3.18	3.10	Fct = 0.031
			Among populations	16272.464	52.08	50.75	Fsc = 0.524
			Within populations	7198.333	47.36	46.15	Fst = 0.538
	Neighbor-Net	K = 4	Among clusters	3310.346	15.78	15.03	Fct = 0.150
			Among populations	13873.399	42.14	40.12	Fsc = 0.472
			Within populations	6971.500	47.10	44.85	Fst = 0.551
	Structure algorithm	K = 3	Among clusters	2589.107	16.12	14.96	Fct = 0.150
			Among populations	14594.638	44.50	41.31	Fsc = 0.486
			Within populations	6971.500	47.10	43.73	Fst = 0.563
	Structure algorithm	K = 4	Among clusters	3398.245	16.60	15.74	Fct = 0.157
			Among populations	13785.499	41.75	39.59	Fsc = 0.470
			Within populations	6971.500	47.10	44.67	Fst = 0.553
	Structure algorithm	K = 20	Among clusters	9018.235	36.00	34.74	Fct = 0.347
			Among populations	6507.650	21.60	20.85	Fsc = 0.319
			Within populations	5889.167	46. 01	44.41	Fst = 0.556

Acknowledgements

We are grateful to X. Giráldez and all other colleagues for collecting material that was used in this study: A. Abad de Blas, A.C. Torres, D. Pinto-Carrasco, E. Rico, G. Calabrese, J. Fuentes, J. Peñas de Giles, L. Gutiérrez, T. Romero, M. Santos-Vicente, R. Peña-García, S. Andrés-Sánchez, S. Rubio, V. Di Donato and V. Lucía. We thank P. Kosachev, B. Frajman and Julio Fuentes for providing material. We would like to express our gratitude to S. Filoche for providing us information to find V. satureiifolia in Fontainebleau. We thank Teresa Malvar for technical support in the lab, P. Reeves for help running PCO-MC analysis, and to our friends E. Rico and X. Giráldez for their continuous support. We also thank two anonymous reviewers for their helpful comments and suggestions. This work was supported by the Spanish Ministerio de Economía y Competitividad (projects CGL2009-07555 and CGL2012-32574); and the University of Salamanca (Ph.D. grant to NPG co-founded by Banco Santander). FCT with POPH/FSE funds financed SC and MC (starting grant project IF/01267/2013 and SFRH/BD/89910/ 2012).

Appendix B. Supplementary material

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ympev.2017.11.007.