Plant evolution in alkaline magnesium-rich soils: A phylogenetic study of the Mediterranean genus Hormathophylla (Cruciferae: Alysseae) based on nuclear and plastid sequences

Habitats with alkaline edaphic substrates are often associated with plant speciation and diversification. The tribe Alysseae, in the family Brassicaceae, epitomizes this evolutionary trend. In this lineage, some genera, like Hormathophylla, can serve as a good case for testing the evolutionary framework. This genus is centered in the western Mediterranean. It grows on different substrates, but mostly on alkaline soils. It has been suggested that diversification in many lineages of the tribe Alysseae and in the genus Hormathophylla is linked to a tolerance for high levels of Mg+2 in xeric environments. In this study, we investigated the controversial phylogenetic placement of Hormathophylla in the tribe, the generic limits and the evolutionary relationships between the species using ribosomal and plastid DNA sequences. We also examined the putative association between the evolution of different ploidy levels, trichome morphology and the type of substrates. Our analyses demonstrated the monophyly of the genus Hormathophylla including all previously described species. Nuclear sequences revealed two lineages that differ in basic chromosome numbers (x = 7 and x = 8 or derived 11, 15) and in their trichome morphology. Contrasting results with plastid genes indicates more complex relationships between these two lineages involving recent hybridization processes. We also found an association between chloroplast haplotypes and substrate, especially in populations growing on dolomites. Finally, our dated phylogeny demonstrates that the origin of the genus took place in the mid-Miocene, during the establishment of temporal land bridges between the Tethys and Paratethys seas, with a later diversification during the upper Pliocene.

The recently revised tribe Alysseae DC. [16][17], although reasonably well delimited, still presents some conflicting evidence from a phylogenetic perspective [1][2][3]16], that only recently has been partially solved [17]. According to Španiel et al. [17], the tribe currently includes 24 genera, with a distribution mainly centered in Eurasia and North Africa, encompassing 277 species of which more than a half belong to the genus Alyssum L. (114 spp.). The following genera by number of species are Odontarrhena C.A.Mey. ex Ledeb. (87 spp.) and Hormathophylla Cullen & T.R.Dudley (10 spp.) [18]. All members of Alysseae are characterized by stellate trichomes, latiseptate or cylindrical (rarely angustiseptate) fruits, with few seeds, winged seeds and a basic chromosome number of x = 8, although smaller or larger derived numbers have also been observed [17] (Alyssum, Odontarrhena, Clypeola, Hormathophylla). Hormathophylla has been the subject of several taxonomic treatments since its original description [19][20][21][22][23] and, although current evidence suggests that Hormathophylla is a monophyletic genus [3,16], these studies are based on a limited sampling of species omitting several taxa of controversial identity.
The genus Hormathophylla contains about 10 species distributed in the western Mediterranean (from northern Algeria to southern France), many of which are narrow endemics from the Baetic ranges in SW Spain (Fig 1). These species can be distinguished morphologically from Alyssum by their white or pink flowers, lateral sepals which are small or not gibbous at all, non-appendiculated filaments and larger fruits between 5-10 mm [3]. Regarding fruit morphology, the genus includes a group of species, such as H. spinosa (L.) P. Küpfer From a cytogenetic perspective, the great variability found in chromosome number (x = 7, 8, 11, 15 [21,24]) in comparison to most of the related genera, which regularly bear a basic chromosome number x = 8, suggests that speciation processes have been accompanied by alterations in chromosome number through polyploidy and dysploidy mechanisms [21]. Küpfer [21], in his thorough study, offered a detailed analysis of possible scenarios of chromosome number evolution in Hormathophylla. He also pointed out the existence of four morphologically different groups in the genus that are in accordance with his hypotheses about the origin of different chromosome numbers: (1)  Apart from the morphological and chromosome variability, a unique trait observed in this genus is the ability to thrive on a wide range of substrates. Rešetnik et al. [16] underscored the importance of the soil in the diversification process of the tribe, supporting the monophyly and sister relationship between the two clades containing most of the Alysseae serpentinophytes: the so-called Bornmuellera and Clypeola clades. Species of Hormathophylla that are part of the Bornmuellera clade can be found growing on a wide range of soil types. They are The number of species in the genus that thrive on dolomites is remarkable, as this is a demanding substrate that shares common features with serpentines [25][26][27]. Whereas serpentinophytes have been thoroughly studied in Alysseae, in which some of the species have evolved to become heavy-metal hyperaccumulators [11,[28][29], little is known about the evolution of species on alkaline soils. The combined presence of high levels of Mg [30-31] and a low proportion of Ca/Mg are the common conditions defining these alkaline habitats, a feature also present in serpentines [25][26]. In addition to soil chemical composition, another factor influencing plant adaptation to this potentially toxic substrate is their parallel occurrence in highly xeric environments, like steppes or rocky habitats [32].
The present study attempts to address several of the questions still unanswered about the evolution of this genus. First, considering that several groups in the tribe Alysseae were previously revealed to be poly-or paraphyletic, it was compelling to establish whether a comprehensive sampling supports Hormathophylla as a monophyletic genus. Secondly, given the previous uncertainty about the generic assignment of H. purpurea (considered by some authors as a separate monotypic genus, Nevadensia Rivas Mart.), it was necessary to elucidate its phylogenetic relationship to other Hormathophylla species. Thirdly, the degree of distinctiveness between pairs of species (e.g. H. cadevalliana/H. longicaulis, H. ligustica/H. saxigena) has not been previously inspected by molecular markers. Finally, to ascertain whether the great morphological variation could be related either to type of substrate or to taxonomic group, we performed a morphometric study of leaf trichomes. Despite the high homoplasy affecting trichome morphology in the tribe Alysseae [10,33], when studied at the generic level, it has been noted that several trichome types are often restricted to well-defined groups of genera within this tribe [11]. The morphology of trichomes was also succesfully used for the distinction of species and subspecies in some intricate species complexes of Alysseae [34][35].
Once the phylogenetic framework was established, a complementary approach was carried out to determine the putative association between substrate and lineages, including the study of species growing on more than one type of soil, the evolution of morphological traits, such as fruit or trichome morphology, and the chromosome evolution. Finally, we used the phylogeny to perform an analysis of biogeographic patterns presenting a spatial-temporal framework by molecular dating to identify patterns compatible with the different paleoclimatic Mediterranean events that could have been involved in diversification episodes over time. To carry out these objectives, we explored a set of regions in two cell genomic compartments. In the chloroplast DNA sequences, we used regions trnL-F and ndhF [16,33,36] and for nuclear we used the ITS region from ribosomal nuclear DNA, previously used in studies of genera related to Hormathophylla [11,16,37] to explore both the phylogenetic relationships among species as well as among different edaphic races or ecotypes [29,38].

Taxon sampling
We sampled all the species that were previously assigned to Hormathophylla with an emphasis on the most taxonomically controversial taxa from the Baetic ranges (Table 1). Each population sample consisted of ten individuals that were collected at a distance of at least ten metrers apart from each other. Plant material was dried in silica gel and stored at room temperature. Mey.) Boiss., respectively) because of their uncertain generic affinities and their possible relationship to Hormathophylla. In this sense, the putative ascription of C. antiatlantica to Hormathophylla was informally suggested by Küpfer (in some annotated herbarium specimens, e.g. MA 121991) as well as by Maire, who highlighted its morphological similarity to H. cochleata and Alyssum sect. Ptilotrichum (C.A.Mey.) Hook.f. [21,39].
To assess monophyly of Hormathophylla and its position in the Alysseae tribe, we used data sets that contained up to 94 accessions belonging to the genus Hormathophylla and 120 more from other species of the Brassicaceae family, the first composed of ribosomal sequences (DS1) and others composed of plastid sequences (DS2) (in those cases where sequences corresponding respectively to ndhF and ITS sequences were available in GenBank). All datasets used are detailed in S1 Table. Analyzed outgroups were chosen according to previous phylogenetic studies [3,10,16,40]. Name, voucher information with herbarium abbreviations according to Thiers [41] and GenBank numbers for all the accessions used are shown in Table 2 (from Hormathophylla genus) and Table 3 (rest of taxa used in the studies). One more data set (DS3) was used that was composed of plastid sequences (trnL-F, trnT-trnL and rpl32-trnl) sequenced in the case of the genus Hormathophylla and species from related genera Fibigia, Brachypus and Alyssoides.

DNA isolation, PCR amplification and sequencing
Plant material was desiccated (0.10-0.05 g) and ground to a fine powder with a mixer mill MM400 (Resch, Haan, Germany). Total DNA was extracted using a DNAeasy Plant Mini Kit (Qiagen), and average concentration was estimated with the aid of a spectrofluorometer (Synergy Mx Microplate Reader, Biotek).
The ribosomal region (ITS1-5.8S-ITS2) was amplified in a reaction volume of 20 μL, which contained 10 ng of genomic DNA, 1 μM each primer (C26A and N-NC18S10 [42]), 0.25 μL of Gotaq polymerase (5uds/μl) (Promega, Madison, Wisconsin, USA), 1.5 μL MgCl 2 (25 mM), and 1.6 μL dNTPs (2.5mM). PCR conditions were: 5 min at 95˚C, followed by 35 cycles of 1 min at 95˚C, 1 min at 50˚C and 1 min at 72˚C each, followed by an incubation at 72˚C for 10 min. PCR amplification products were purified using a Genelute PCR cleanup Kit (Sigma-Aldrich, St Louis, MO, USA). Sequencing reactions were performed by means of the BigDye Terminator V3.1 (Applied Biosystems, Foster City, CA, USA) sequencing kit. Sequences were purified following manufacturer recommendations, and using the same primers as in the PCR. Reaction products were run on an ABI 3100 Avant (Applied Biosystems) automatic sequencer.
For cpDNA sequences, three regions were selected, according to the degree of polymorphism shown. These sequences were: trnL-trnF, trnT-trnL and rpl32-trnL. In the phylogenetic analysis of Alysseae, a fragment of the ndhF region was also included [16]. Primers used in the amplification reaction were those proposed by Shaw et al. [43] (in rpl32-trnL fragment), Taberlet et al. [44] (in the case of trnL-trnF and trnT-trnL) and Beilstein et al. [33] (in ndhF fragment). PCR settings were those suggested by Shaw et al. [43,45] with minor modifications.
Purification of amplification products and sequencing reactions were performed as in the ITS region. All sequences were uploaded to Genbank and accession numbers are shown in Table 2.

Sequence analysis and haplotype identification
Forward and reverse sequences were compared, assembled and corrected where necessary using ClustalW [46] and BioEdit V 7.0.5.3 [47], to establish a consensus sequence in each sample. In the case of ribosomal sequence alignment, the use of Muscle [48] was necessary. Following the recommendations of Fuertes-Aguilar et al. [49][50], sites were considered as polymorphic (additive polymorphic sites or APS) when they exhibited more than one nucleotide in the chromatogram, an when the weakest signal was at least 25% of the strongest one. All of the obtained sequences were submitted to NCBI GenBank. The evolution model represented by a phylogenetic tree, does not always adequately describe evolutionary scenarios like hybrid speciation, recombination or chloroplast capture via introgression [51]. Hence, reconstruction of phylogenetic networks may help in the interpretation of relationships [51][52]. For this particular purpose, one analysis was performed with data set DS3, which contained concatenated sequences from the amplified plastid regions (trnL-trnF, trnT-trnL and rpl32-trnL), obtained from the samples of the genus Hormathophylla, and the outgroup from the genera Fibigia and Alyssoides Mill. (in total 96 sequences). Statistical parsimony networks (statistical parsimony algorithm [52]) were constructed using the program TCS V. 1.21 [53]. In this analysis, gaps were treated as missing data.

Phylogenetic analysis
Phylogenetic analyses were achieved by using Bayesian Inference (BI) and Maximum Likelihood (ML) inference. A group of data sets were created for different analyses (S1 Table). In the case of ribosomal sequences, we used data sets DS1 and DS4, the first with an extended sampling to cover the whole tribe Alysseae, and the second with the sampling focused on Hormathophylla. In the case of chloroplast sequences, we used two data sets: DS2 for the analysis of the whole tribe and DS3 for the analysis centered on the genus. We also assembled a concatenated matrix with sequences belonging to ribosomal and cpDNA sequences, creating a new data set called DS5.
For the BI phylogenetic analysis, we employed MrBayes 3.1.2 [54][55]. Selection of the bestfit model of nucleotide substitution was carried out according to Bayesian and Akaike Information Criteria (BIC and AIC), as implemented in jModelTest [56]. 5x10 6 generations were produced, with four Markov chains and sampling every 100 th generation. The first 10% of the sampled trees was discarded as burn-in. With the remaining trees, a 50% majority rule consensus tree was generated and posterior probability used as an estimation of clade support.
ML analysis was performed using RAxML 7.0.4 software [57][58]. Phylogenetic analyses were carried out through 1000 fast bootstrap analyses followed by a search for the best resulting ML tree in a single run. Because of the demanding computational requirements, a version of the program located in the web CIPRES-gateway (http://8ball.sdsc.edu:8888/cipres-web/ home) was used. Following recommendations of [59], the General Time Reversible model (GTR) was used with an alpha parameter for the shape of the gamma distribution to account for among-site rate heterogeneity for both datasets. Highly congruent topologies were obtained using the parsimony ratchet [60], following the same method as in [61].

Data analysis for the estimation of divergence times
Two data subsets, obtained from DS1 and DS2 respectively, called DS6 and DS7 were used for the estimation of divergence times. As in the phylogenetic analyses, AIC and BIC criteria were considered to determine the substitution model that best fitted the data sequences. Three independent replicates of the BI analysis of the genus Hormathophylla phylogeny generated trees with similar topology. Only one run was used for the construction of the consensus tree, with 25% of sampled trees deleted as burn-in.
Molecular dating analyses were performed with BEAST v1.8 [62], a program designed to estimate divergence times by means of the Bayesian Markov chain Monte Carlo (MCMC) approach. The tool Beauti [63] was used to edit the input file for BEAST. Thereafter, a relaxed molecular clock model was implemented using the uncorrelated lognormal parametric algorithm [62]. All the analyses were carried out on the assumption of a birth-death speciation model as a tree prior [64], assuming constant rates of speciation-extinction per lineage. Four runs with two chain, were performed for the dating analysis, each with the MCMC chain length set to 75,000,000, so as to extract every 10,000th and sample every 10,000 trees. The results of the BEAST analyses were checked in Tracer 1.5 [65] for model of likelihood and parameters convergence between each run and that each run had reached a stationary state.
Both chains were combined using LogCombiner 1.8, after discarding the first 10% of the sampled trees. Results were considered reliable once effective sample size (ESS) values of all

Scanning electron microcopy study. Trichome morphology
Scanning Electron Microscopy (SEM) was used to examine trichome morphology of selected taxa from the tribe Alysseae. Taxa and voucher specimens are listed in S2 Table. These studies were carried out at the Real Jardín Botánico SEM facility. Small portions of leaves were taken and mounted on a stub using double adhesive tape, and sputter-coated with gold in a Sputter Coater Balzers model SCD 004 with a thickness of about 50-700 μm. The specimens were visualized with a Hitachi S3000N digital electron microscope, operated at 25 kv. Since stellate trichomes are usually dichotomously branched up to three orders and the branch number varied depending on their order, we restricted the term ray to the last order number of ramification. This number has to be counted on the upper side of the leaves, where its presence is usually higher. Indument density and number of indument layers, as well as presence of rugosity in the epidermis, were registered and trichomes classified based on type of branching number of first-order arms, number of succesive divisions per main arm, and number of terminal rays. In Hormathophylla species, a series of measures were taken for the following parameters of stellate trichomes: longest and shortest ray length, and its ratio (roundness index), central disc diameter, trichome ray length between same order branching, branching angle and branch diameter.

Flow cytometry and estimation of DNA ploidy levels
Flow cytometry was used to measure relative nuclear DNA content and infer ploidy levels of the studied populations. Ploidy was first analyzed in plants of Hormathophylla from populations with known chromosome numbers [21,[69][70][71][72] and then used to infer ploidy level in https://doi.org/10.1371/journal.pone.0208307.t002 Evolution on alkaline magnesium-rich soils of the Mediterranean genus Hormathophylla      [74] were used as internal standards. Up to 10 samples per locality were analyzed when possible (Table 1). Fresh and young leaves were dried in silica gel immediately after field collection and were stored at 25˚C [75]. For the isolation of nuclei, desiccated leaves (0.5 cm 2 ) were co-chopped Evolution on alkaline magnesium-rich soils of the Mediterranean genus Hormathophylla with the fresh leaf tissues from a standard individual using a razor blade in a Petri dish with 1 mL of ice-cold Otto I buffer [76]. The resulting suspension was filtered through a 42 μm nylon mesh and incubated for at least 5 min at room temperature. One milliliter of a solution containing Otto II buffer [76] supplemented with 2-mercaptoethanol (2μL/mL) and DAPI (4μg/ ml) was added to the flow-through fraction, and stained for 1-2 min. Intensity of fluorescence of 5000 particles (dyed nuclei) was measured using a Partec Cyflow ML instrument with an HBO-100 mercury arc lamp (Partec GmbH, Munster, Germany). The flow cytometry histograms were evaluated using Partec FLO MAX V. 2.4d (Partec GmbH) software. The reliability of the measurements was assessed by computing the coefficients of variation (CVs) from both the analyzed and the standard samples. All analyses above the CV threshold value of 5% were rejected.

Phylogenetic analysis of the genus Hormathophylla using ribosomal sequences
According to the results of the expanded ITS data set (S1 Fig Table). In the jModeltest analysis, the evolutionary model was selected by both AIC and BIC (S1 Table). The tree topologies obtained in Bayesian inference and Maximum likelihood analyses were similar, although the branch support values and some connections varied slightly. The resulting trees are shown in Fig 2, centered in Hormathophylla species and in S1 Fig with an enlarged sampling of outrelated genera.

Species diversification in Hormathophylla
According to our ITS tree results (Fig 2), Hormathophylla can be considered a monophyletic genus (1.00 PP; 100% BS), which would include H. cochleata and H. purpurea as taxa belonging to this group. The split between two lineages, hereafter called A and B, is clear. Both lineages correspond to the two groups of species previously detected based on basic chromosome numbers and are well supported (lineage A was 0.75 PP; 77% BS and lineage B 0.81 PP; 60% BS).
The internal relationships among the chloroplast haplotypes based on the same haplotype data set (DS3, Table 2) were resolved into one single TCS network (Fig 3). Among the 96 analyzed sequences, we found 50 haplotypes, including those in the outgroup (one for Fibigia suffruticosa and Alyssoides utriculata and two for Fibigia clypeata). We found H. pyrenaica to be the species with haplotypes closest to those present in outgroups, separated by 165 steps from the nearest outgroup. The rest of the species were split into three well-differentiated haplotype groups.  (H27 and H30). In this haplogroup the most frequent haplotype was H20, which was present in three of the populations from H. purpurea, and separated by six steps from the nearest haplotype of H. spinosa. Regarding H. saxigena, its haplotype was separated four steps from the exclusive H. ligustica haplotype and 11 steps from nearest haplotype of H. spinosa.
There is an association between the haplotypes and the type of substrate in which the different populations grow. This association is, in some cases, independent of the taxonomic ascription. Within haplogroup I, we find H. lapeyrouseana, with most of the haplotypic lineages present on gypsum, and the two subspecies of H. cochleata, which are present on calcareous outcrops (limestone and dolomite). Within haplogroup II, we find species that were always present on alkaline soils, with a particular predominance of dolomitic or calcodolomitic substrates (H. reverchonii, H. longicaulis, H. cadevalliana and H. lapeyrouseana), but all the individuals analyzed from populations living on serpentines (H. longicaulis) are in this group. Finally, in haplogroup III, the most diverse, we find species on different substrate affinities, although most of them were predominantly on limestone (H. pyrenaica, H. ligustica, H.  saxigena and H. spinosa). Moreover, H. spinosa and H. purpurea were found on siliceous substrates.

Within Hormathophylla plastid phylogenetic analysis
When the plastid data set (DS3) is used to infer the phylogenetic relationships within Hormathophylla, the resulting tree (Fig 4) exhibits a lower resolution than the topology obtained from ITS ribosomal sequences. However, the topology does not exactly match the two lineages A and B, as identified in the ITS phylogeny (Fig 2). The species included in ITS lineage A grouped into its own clade (0.7 PP; 76% BS) in the plastid phylogeny. The clade is completed with the sequences belonging to the southernmost populations of H. lapeyrouseana. Within   trnL-trnF, trnT-trnL and rpl32-trnL, shown as a majority

Phylogenetic analysis of tribe Alysseae based on the plastid region
The analysis based on data set 3 (DS3, Table 2

Phylogenetic analysis of the combined matrix with nrDNA and cpDNA
The analysis of the concatenated sequences of ribosomal and plastid DNA considering the partitioned data (DS5) generated the consensus tree shown in S3 Fig. First, the monophyly of the Hormathophylla genus was conserved (1.00 PP), within the Bornmuellera clade (1.00 PP). As in the plastid analysis, it was not possible to differentiate species from each of the two ITS groups. Nonetheless, we found that species from lineage A came together in a well supported group. Concerning Cuprella antiatlantica and C. homalocarpa, these species grouped into a lineage nested within the Bornmuellera clade, as a sister group of the genus Bornmuellera. The phylogenetic tree was similar to that obtained from chloroplast sequences, except in the case of the relationships between the Alyssum, Clypeola and Aurinia clades.

Estimation of divergence times
The dated phylogeny, with the estimation of divergence times (Fig 5 and S4 Fig), depending on the genome compartment used in its inference, was based on ribosomal (DS6) or plastid sequences (DS7).
As a result of the parametric analysis of the ribosomal sequences, the split between Hormathophylla and its most related group composed of species from the Fibigia genus indicated The data phylogeny based on the use of plastid sequences varied widely with respect to those obtained with ribosomal sequences. Lineage A did not split well from lineage B, but if we consider the southernmost population of H. lapeyrouseana, this group split from the rest of the species 1.94 Ma (0.89-3.37 Ma 95% HPD) ago. The most recent common ancestor of the Hormathophylla genus was established with a mean age of 4.18 Ma (2.29-6.68 Ma 95% HPD). The split between Hormathophylla and the related group composed of species from the Fibigia genus showed a mean age of 16.14 Ma (9.70-22.88 95% HPD). The split of the Alysseae tribe with respect to the rest of the Brassicaceae was about 35.89 Ma (28.03-44.31Ma 95% HPD).

Leaf trichome morphology
The extraordinary trichome complexity and variability of the leaf blade indumentum required an additional effort of classification. All species from Hormathophylla presented stellate trichomes, peltate (briefly stalked), appresed and with a circular to elliptic outline ranging from 0.3 to 0.5 mm in diameter. Trichomes usually show a more or less protuberant central disc, from which a number of primary and secondary dichotomous branches extend usually in the same plane. Primary rays (four to eight) diverge from a central disc that in most species are reduced to a small central area (e.g. H. purpurea and H. cochleata). The number of last-order branches varies between 13 and 36. The central disc can show tuberculate and irregularly bulged ornamentation, but can also be smooth. All trichomes showed rounded tubercles with a different degree of density and size, that become sparser or absent at the branch tips. We were able to distinguish four types of trichomes among all species, which are present in different species depending on their phylogenetic position (Fig 6, S5 Table).
When we mapped trichome types along the phylogentic trees we found an association between lineages and the morphological types. In ITS lineage A (Fig 2), we detected a single type of trichome (type I). In this type, trichomes were stellate and broadly peltate, with an almost perfectly circular outline (diameter: 0.25-0.35 mm); they were largely overlapping and arranged in more than three layers, so dense that no leaf epidermis could be seen unless the hair cover was entirely removed. Primary rays (7-8) diverge from a rounded and undivided central disc (ca. 2/10 diameter) and are reduced to a brief, stout base bearing a total of 28-32 branch tips, almost regularly dichotomous. Rays were sparsely covered by tuberculate Evolution on alkaline magnesium-rich soils of the Mediterranean genus Hormathophylla protrusions, that become sparser and disappear at the branch tips. These type I trichomes have not been observed in any other member of the Alysseae tribe.
With respect to ITS lineage B (Fig 2), on the contrary than in lineage A, we observed trichome differences between the species nested there. Trichomes were very dense, stellate and broadly peltated, usually with a rounded outline, and variable diameter from 0.3 to 0.4 mm. We found that primary branches were longer than those in lineage A, and tubercles were larger, more prominent and thicker, with secondary branching occurring farther from the center. Within this lineage, we could distinguish three trichome types classified depending on diameter, number of rays, presence of membranaceous structures among branches, and central disc diameter. Type II is present in H. lapeyrouseana, and the two subspecies of H. cochleata. The diameter and central disc were larger than the rest of the species (0.41 to 0.44 mm and 0.7 to 0.9 mm, respectively), with branches making up an almost circular outline, and the presence of membranaceous structures among them. Type III appears in H. purpurea and H. pyrenaica. This type of trichomes is similar to the types found in other Alysseae genera, such as Alyssum or Clypeola. Trichomes showed a low number of tubercles, lower number of branch tips (12 to 14), and an elongated central disc, with branches limited by an elliptic outline. Type IV is detected in H. pyrenaica, H. spinosa, H. saxigena and H. ligustica. In this type, trichomes showed an almost perfectly circular outline (diameter: 0.30-0.35 mm); they were largely overlapping and arranged in more than three layers, and so dense that no leaf epidermis could be seen unless the hair cover was entirely removed. Primary rays (4-7) diverge from a rounded and undivided central disc (ca. 1.5-2/10 diameter) and are reduced to a short, stout base bearing a total of 30-36 branch tips, almost regularly dichotomous. Rays were densely covered by rounded tubercles, even at the branch tips.

Chromosome variation and ploidy level in Hormathophylla
For most of the species we detected the same ploidy level (tetraploid) in most of the studied specimens (up to 10 individuals per population; Table 1). We found that all analyzed individuals from H. purpurea, H. spinosa, H. cochleata, H. reverchonii and H. saxigena exhibited a constant level of ploidy in all studied populations. We would like to note that in the particular case of samples belonging to the Morocco populations of the species H. spinosa, we did not find 2x individuals, even though they had been reported previously [77]. The occurrence of two different levels of ploidy in the same species has been recorded, in the cases of H. cadevalliana and H. longicaulis, not associated with any geographic pattern. Octoploid individuals (8x) were found in populations 1, 3 and 4 from H. longicaulis and 13 from H. cadevalliana. For these species, tetraploid individuals (4x) were found in populations 3, 5, 6, 7, and 8 in H. longicaulis and 9, 10, 11 and 12 in H. cadevalliana, which means that in some populations the two levels of ploidy co-occur. For example, in the H. longicaulis population 3, we found nine octoploid and one tetraploid individual. In both species, this is the first time that 8x level is reported for H. cadevalliana, and the first time 4x level is reported for H. longicaulis.

Origin and biogeography of Hormathophylla
Our results indicate that the Hormathophylla clade split from its sister group, composed by species from the genus Fibigia, in the Late Miocene (10.5 Ma) (Fig 5). During this period, the existence of a northern Mediterranean land corridor connecting the seas of Tethys and Paratethys [77][78][79][80] favored dispersal of a number of lineages from east to west throughout the Mediterranean region. These lineages underwent progressive isolation developing eastwest disjunctions, because of the interruption of the corridor by marine transgressions in Late Miocene and the structuration of the Alps (10 Ma) [81]. Such an east-west disjunct distribution within the Bornmuellera clade is observed between lineages present in western-Mediterranean and eastern-Mediterranean, sometimes extended in central Asia (S5 Fig). This is a biogeographic pattern described in insects as the so called Kiermack disjunctions [82] but also in Mediterranean plants like in the tribe Anthemidae [83] and the genus Jasione L. [84].
Therefore, the Hormathophylla ancestor established in western Mediterranean diversified into two different lineages in the upper Pliocene (about 3.8-4.18 Ma ago depending on whether we consider ITS or plastid sequences). However, it is not until the transition between the Late Pliocene and Gelasian in the early Pleistocene (2.2 and 2.9 Ma for lineages A and B, respectively), during the onset of the Mediterranean climate (3.4-2.9 Ma), when the rapid diversification of Hormathophylla species took place, fundamentally in the Baetic ranges. This general pattern of rapid Baetic diversification in plant taxa is well documented by both paleobotanical [85][86][87][88][89] and molecular studies [90]. The onset of the mediterranean climate was a significant environmental change caused by the establishment of rainy temperate seasons and drought in summer months. Based on dated phylogenies of several plant genera, e.g. Dianthus L. [91] or Cistus L. [92], there was an increase in rates of diversification with the establishment of the Mediterranean climate with dry summers.
The small area of distribution of most of the species from the genus, (H. reverchonii, H. cochleata subsp. baetica, H. purpurea or H. pyrenaica), all of them with a Pleistocene crown age, would additionally suggest that fragmenting patterns have played an important role in the biogeographical diversification of Hormathophylla. The influence of the Pleistocene stadialinterstadial dynamics in the southern Europe ranges has already been described in other genera like Armeria Willd. and Anthyllis L. [93][94][95][96].
The consequences of speciation by geographical isolation affected mostly in Baetic ranges (as well as in northern Morocco and southern France), followed by secondary contacts, with expansion and contraction of species distribution reflecting climatic fluctuations of the Quaternary period [97][98]. During the glacial period, broad areas at middle altitudes were covered by cold steppe and tundra biomes, providing suitable habitats for mountain plants [99]. It is possible that interglacial warming caused populations to retract to the top of the mountains [94,100]. The glacial/interglacial cycles in southern Iberian ranges helped to shape the present distribution of genetic lineages, which are geographically contiguous with each species adapting to a different substrate [101][102][103].

The genus Hormathophylla is monophyletic
Our analyses confirm the monophyly of Hormathophylla, resolving two conflicting questions regarding the limits of the genus: the separation of H. purpurea as a monotypic genus and the generic ascription of Alyssum antiatlanticum to Hormathophylla. H. purpurea had been proposed as a separate monotypic genus, Nevadensia (Nevadensia purpurea (Lag. and Rodr.) Rivas Mart.) [22] based on the morphological and ecological differences existing between H. purpurea and the remaining species of Hormathophylla. This view has been followed in some recent publications [104]. However, Küpfer [21] who initially considered that Ptilotrichum purpureum (Lag. & Rodr.) Boiss. should be treated as Alyssum purpureum Lag. & Rodr., later transfered this species to Hormathophylla due to the color of flowers (pink), unusual in Alyssum [18]. Our plastid phylogenetic analyses not only confirm Küpfer's view, but also reveal a close relationship of this species with H. spinosa, the only other species with pink petals (Fig 4). The ITS phylogeny reveals the placement of this species in lineage B, and is totally congruent with its basic chromosome number x = 8 detected by Contandriopoulos [21].
The second potentially disruptive element in the monophyly of Hormathophylla was Alyssum antiatlanticum. A. antiatlanticum is an isolated species in the genus that Maire [39] assigned to Alyssum sect. Psilonema (C.A.Mey.) Hook. f. because of its elongated cylindrical nectaries. The putative assignment of A. antiatlanticum to Hormathophylla was informally proposed by Küpfer based on some annotated specimens (e.g., MA 121991), as well as by Maire who highlighted its morphological proximity to H. cochleata (Coss. & Durieu) P.Küpfer and Alyssum sect. Ptilotrichum (C.A.Mey.) Hook.f. [21,39]. Our phylogenetic analyses constitute the genetic basis for the recently published genus Cuprella, which includes C. antiatlantica and C. homalocarpa and reject its phylogenetic placement within Hormathophylla [17].
The obtained phylogenies reveal the early split of the two lineages according to its basic chromosome number, distinguishing those belonging to lineage A, derived from x = 7, from those belonging to lineage B, probably derived from x = 7 and x = 8 via hybridization and aneuploidy [21]. This differentiation is also supported by a separation in trichome morphology. Whereas species belonging in lineage A (Type I) were very similar, species in lineage B encompass three types (II, III and IV) (S5 Table). Regarding the basic chromosome number x = 8, which is by far the most common in the tribe outside of Hormathophylla (171 vs 7 species), we postulate that species grouped in lineage A with 7 as a basic chromosome number are derived from a common ancestor with x = 8.

Relationships between H. cadevalliana and H. longicaulis
Lineage A was composed of H. reverchonii, H. cadevalliana and H. longicaulis. The low phylogenetic divergence, reflected in the low resolution between species belonging to lineage A, is supported by the almost uniform trichome type I (Fig 5) morphology observed, where no distinction in trichome types is found between H. cadevalliana and H. longicaulis, and the indumentum of both species is very similar to that of H. reverchonii (S5 Table).
An explanation for the lack of differentiation in ITS topology might suggest the occurrence of introgression among some of the taxa. In the case of H. cadevalliana and H. longicaulis, ITS sequences obtained appear intermixed in the tree (Figs 2 and 4) with the presence of additivity patterns (S3 Table). All the evidence seems to indicate that there is either a past or an ongoing hybridization process between the two species in part of their overlapping areas of distribution. Such an admixture pattern is more evident in some of the southern ranges (Baza and Filabres) where populations from both species are even sympatric. The diversity observed in ploidy level remarkably only appears in those localities where sympatric populations from two species occur, e.g. H. longicaulis population 3 groups individuals with 4x-and 8x-levels of ploidy. Likewise, other H. cadevalliana populations (14) exhibit a similar additive pattern with H. reverchonii (S3 Table). In addition, plastid haplotypes from all three species in this lineage are all clustered in the haplotype network (Fig 3). We could not find any correlation between morphological variation and level of ploidy within the same species, which would preclude the possible taxonomic recognition of infraspecific taxa based on that criterion.
If introgression between the two species is actually occurring, the lack of resolution observed in the ITS phylogenetic tree could be explained by the inclusion of sequences in a process of incomplete concerted evolution. The occurrence of individuals with ITS sequences rich in additive polymorphic sites [50] can also be favored by the presence of early generation hybrid individuals [49], in which concerted evolution of multi-copy genes is retarded by a lack of homologous loci recombination. Remarkably, in the genus Hormathophylla, only in H. longicaulis and H. cadevalliana has the presence of octoploid individuals been demonstrated by chromosome counts and ploidy estimations via genome size [21] ( Table 1, this study). A plausible scenario in which hybridization could arise was along the succesive cycles of expansion and contraction due to cycles of glaciation that could favor the hybridization processes as previously demonstrated for other plant groups in the same soutwestern Iberian ranges [102][103].

Fruit morphology. Genetic differentiation between northern and southern populations in H. lapeyrouseana and in relation to adaptation to dolomitic substrates
Not many morphological traits support the species ITS lineages identified in Hormathophylla. All species with cochlear fruit were present within lineage B, although they could have originated separately twice: in H. spinosa and H. cochleata subsp. cochleata, H. cochleata subsp. baetica and H. lapeyrouseana clade, respectively (Fig 2). This result is in line with the considerations made by Al-Shehbaz et al. [2] in relation to the use of the fruit morphology as a taxonomical character in the family Brassicaceae. This author suggests that important changes in fruit morphology can occur faster and independently of other characters, subject to a considerable convergence and sometimes taxonomically irreproducible. There are many studies that suggest that differences of only a few genes can cause substantial variations in shape, size and dehiscence of the fruit in the family Brassicaceae, e.g. in Arabidopsis, genes involved in formation of fruit (FRUITFUL, MADS-box, SHATTERPROOF) can modify shape (e. g., ratio length/width) and dehiscence type [105][106][107][108].
H. cochleata subsp. cochleata, H. cochleata subsp. baetica and H. lapeyrouseana constitute a well supported group both with nuclear and plastid markers (with the exception of the southernmost localities of H. lapeyrouseana). Differentiation between H. cochleata subsp. cochleata and H. cochleata subsp. baetica has been subjected to taxonomic studies, being considered as the same [21] or vicariant species [109]. Our results support a synthetic vision, as we will consider later. Furthermore, molecular data establish a close relationship between these species and H. lapeyrouseana. This is based on the similarity of plastid haplotypes that could be caused by a process of recent hybridization, or to an incomplete lineage sorting, given the great differences that exist between their basic chromosome numbers.
On the other hand, the plastid lineage of the southernmost populations was more related to species belonging to lineage A than lineage B (within haplogroup II). This could be due to an incomplete lineage sorting in plastid sequences, as there is no clear evidence of ancestral hybridization (between species from lineages A and B, proposed by Küpfer [21]), as there are no ribosomal additivities shared. Nonetheless, this should not be ruled out, nor should the existence of horizontal chloroplast transference (haplogroup II, with a high number of plants present on dolomite), which could be associated somewhat to adaptation to special substrates.
Other existing relationships, according to the phylogeny obtained through plastid DNA, are between H. spinosa and H. purpurea, H. ligustica and H. saxigena, all with the same chromosome number. This finding could indicate the existence of a common ancestor, or in the case of H. spinosa, incomplete lineage sorting due to the strong relationship among their haplotypes. H. spinosa, has a larger distribution area, where lineages could refuge and re-establish gene flow during unfavorable periods [110][111].
The lineage of H. saxigena and H. ligustica is consistently recovered both in plastid and ribosomal trees (Figs 2 and 4), which is interesting as they are species with a stenochoric and partially sympatric distribution. Furthermore, the presence of different common additivities and symplesiomorphies indicate that they could have been differentiated as species recently. Their plastid lineage is related to the H. spinosa group. We also find the case of H. purpurea included among haplotypic lineages of H. spinosa. It is not clear if its origin may be due to horizontal transfer (equal chromosome number) or an incomplete lineage sorting process. This could be favored by the partially sympatric distribution of these species in Sierra Nevada, with H. purpurea remaining at the top of the mountains, surrounded by H. spinosa, allowing the existence of gene flow.
Finally, H. pyrenaica, a narrow endemic, with uncertain affinities seems to be related to H. cochleata subsp. cochleata and H. cochleata subsp. baetica based on the ribosomal evidence, but not at the haplotypic level, where it seems to constitute a separate lineage, and is the nearest haplotype to the sequence of outgroups (Fig 3). This may indicate that it is a case of retention of ancestral haplotypes due to incomplete lineage sorting.

The relevance of alkaline substrates in the tribe Alysseae and Hormathophylla
The diversification in the process of speciation resulting from the establishment of the Mediterranean climate (Fig 5) is associated in time and space (in the Baetic ranges) with the colonization of different types of substrate and with the parallel geographical isolation [112][113]. This is clear in the case of species diversification in ITS lineage A (2.17 Ma), all of which are dolomiticolous and restricted to Baetic ranges [25]. Edaphic adaptation has been noted as an important factor in the speciation of different plant groups [95,[114][115][116][117]. This has already been suggested in other genera belonging to the tribe Alysseae, e.g. Phyllolepidum or Odontarrhena [11,37]. Adaptation to different substrate types in Hormathophylla involved a strong selective pressure that could promote morphological differentiation and the development of reproductive barriers, preventing genetic flow between disjunct divergent populations [118][119]. One of the most remarkable ecological traits associated with the tribe Alysseae is its ability to grow on alkaline, serpentine, dolomite and limestone substrates (S6 Table). This adaptability is present in most of the recently studied lineages [16]. Particularly in the group of genera sister to Hormathophylla (1.00 PP; 100 BS) composed of Lutzia, Irania, Clastopus, Physoptychis, Pterygostemon, Degenia, Acuston, Alyssoides, Resetnikia, Brachypus and Fibigia, many are found on serpentine outcrops [120][121]. This group is sister to another one encompassing Phyllolepidum Trinajstic and Bornmuellera, which is also prevalent in alkakine substrate (S1 Fig). In total, if we register soil preferences, species growing on alkaline soils constitute a large majority in the tribe, with at least 63 growing on ultramafic, serpentine and volcanic soils and more than 70 on limestone, gypsum and dolomites. Moreover, among all, these genera have species that grow on alkaline substrates with high levels of Mg (both for serpentine and dolomite).
One of the most studied cases of adaptation to extreme edaphic conditions in the Alysseae are nickel hyperaccumulators [11,[28][29], which concentrate in two genera (Odontarrhena and Bornmuellera [122]). Mapping this ecophysiological trait on the ribosomal phylogeny, these species are grouped into two clades along with a high number of taxa able to grow on ultramafic substrates and species tolerant to xeric substrates with high levels of Mg, belonging to 7 out of the 23 genera of Alysseae (i.e., Bornmuellera, Clypeola, Fibigia, Hormathophylla Meniocus, Odontarrhena and Physoptychis). On the other hand, the species exclusively present on dolomite belong to the genera Hormathophylla, and Phyllolepidum. This relationship between dolomites and serpentines has common factors. The ability to grow on serpentines reveals the adaptation to cope with relatively high concentrations of Ni and other heavy metals. However, in addition to this limiting condition in plant adaptation to ultramafic substrates, there are associated factors developed in the same ecophysiological mechanisms, such as the ability to grow in highly xeric environments [32] in the presence of high levels of Mg [30-31] and a low Ca/Mg ratio, which in many cases are more restrictive than Ni levels [123]. These two characteristics also define soils developed on dolomites [25,124]. The high number of species that grow on particular substrates present in certain lineages seems to adjust to the hypothesis by which some taxa have the availability to thrive on Mg-rich substrates. This seems to be the key preadaptive character in serpentine specialization as demonstrated in other angiosperm families [25,[125][126]. The evolutionary relevance between substrate properties and diversification processes have been demonstrated by large-scale studies in the Alps where substrate properties (siliceous or calcareous) are fundamental to explaining modeling of genetic patterns with a large number of species [115].
The adaptation of plant lineages to alkaline or ultramafic substrates has been observed in other Mediterranean regions of the world, in the case of the genera Lasthenia Cass. [127] or Ceanothus L. [128]. Likewise, this plasticity can be found in other families diversified in the Mediterranean Basin, such as Cistus L., Cistaceae [129], or Onosma L., Boraginaceae [11], where species associated with a wide range of alkaline substrates show that obligate serpentinophytes share evolutionary ancestors with non-serpentinicolous taxa [11,125], as far as they grow on substrates with high levels of Mg.
There is a strong support in the plastid phylogeny for a derivation of serpentinicolous populations from dolomiticolous in H. longicaulis, in this case not associated with speciation. However, the transitions between gypsicolous and dolomiticolous habitat is associated with speciation, as in H. cochleata to H. lapeyrouseana (Fig 3). Therefore, it seems relevant that the comparative study of plastid genomes through phylogenetic methods could shed light on the evolution of edaphic adaptations. Two species in particular, H. longicaulis and H. lapeyrouseana, are especially relevant due to their biedaphic behavior. Such biedaphic behavior has been recorded in other dolomiticolous species such as Convolvulus compactus Boiss. or Jurinea pinnata DC., [124,130], which are present in serpentine and gypsum, respectively. It seems possible that parallel adaptative mechanisms (in addition to tolerance to environmental stress) work on these three types of substrate [26,131]. As in Cecchi et al. [37], we believe that the Bornmuellera clade, in particular the genus Hormathophylla, is an ideal system to perform experimental comparative research on their ability to thrive in metal accumulation conditions and the tolerance to magnesium in plants.

The evolution of trichome morphology in Hormathophylla: Taxonomy or ecology?
The observed differences detected among examined taxa, in terms of shape and density of trichomes, raise the question whether these micromorphological fruit and epidermal structures relate to species/lineages or are associated with the ecology and substrate on which the species grows, in all cases xeric habitats (dolomites, serpentines, gypsum, rocky places, etc.). The results of trichome morphological analysis in the Alysseae tribe by Beilstein et al. [10,33] including species of Hormathophylla revealed that, considering Brassicaceae, stellate dendritic trichomes exhibit a homoplastic behavior with multiple origins across different lineages in the family. Our results confirm the lack of exclusive trichome types in Hormathophylla, and lead to the general conclusion that only closely related species retain similar trichome types.
The presence of the lepidote, stellate and fasciculate trichomes have been linked to functions like atmospheric water absorption, protection from sunlight radiation, and cold insulation in high-mountain habitats [132][133]. In Hormathophylla, no strict association was found between the trichome type and substrate beyond the taxonomic boundaries. In contrast within Hormathophylla, some types are strongly associated with groups of species, for example, only type I trichomes from lineage A species (H. reverchonii, H. longicaulis and H. cadevalliana, with a basic chromosome number of 7), in which trichome morphology can even be used to distinguish them from the rest of the tribe. The morphological features (i.e., thickness terminal branch, the ratio length/thickness, and the number of primary branching (seven or eight)) can be used as diagnostic characters for these species.
Among the other species belonging to the genus, an indication that trichome morphology agrees with species boundaries is that the type II trichome characterizes the two species H. cochleata and H. lapeyrouseana. Both species belong to an allopolyploid lineage derived from a cross of species with a basic chromosome number of x = 7 and x = 8. Interestingly, the trichome size parameters (diameter, thickness of primary rays) of this allopolyploid lineage revealed a significantly more robust structure not found in any other species in the genus, which suggests that this is the result of hybrid vigour traits (S5 Table).