Phylogeny, biogeography, and morphological ancestral character reconstruction in the Mediterranean genus Fumana (Cistaceae)

Fumana is a diverse genus of the Cistaceae family, consisting of 21 currently accepted species. In this study, nuclear (ITS) and plastid (matK, trnT‐L) molecular markers were used to reconstruct the phylogeny and to estimate divergence times, including 19 species of Fumana. Phylogenetic analyses (Bayesian Inference, Maximum Parsimony and Maximum Likelihood) confirmed the monophyly of Fumana and did not support the infrageneric divisions previously established. The results support four main clades that group species that differ in vegetative and reproductive characters. Given the impossibility to define morphological characters common to all species within the clades, our proposal is to reject infrageneric divisions. Molecular dating and ancestral area analyses provide evidence for a Miocene diversification of the genus in the north‐western Mediterranean. Ancestral state reconstructions revealed ancestral character states for some traits related to xeric and arid habitats, suggesting a preadaptation to the Mediterranean climate.


Introduction
The Mediterranean Basin is known for its richness in species and high endemism, making it one of the major hotspots for global biodiversity (Médail & Quezel, 1999;Sauquet et al., 2009). The contemporary flora in the Basin has been influenced by tectonic movements and climatic oscillations acting at different spatial and temporal levels (Thompson, 2005), which led to complex patterns of connection-isolation between territories (Rosenbaum et al., 2002;Meulenkamp & Sissingh, 2003;Ree & Sanmartín, 2009). Also, the formation of land bridges between the Tethys and Paratethys seas during the Oligocene and Miocene led to biotic expansions across the Mediterranean (Oosterbroek & Arntzen, 1992;Salvo et al., 2010), and we can still recognize biogeographical links between western and eastern Mediterranean taxa (Médail & Diadema, 2009;Jabbour & Renner, 2011).
The evolutionary path of a plant family can be inferred from the evolution of morphological characters, based on phylogenetic analyses. Certain patterns of character evolution are typical in Mediterranean plant families and may indicate specific adaptations to climatic changes (Ackerly et al., 2002;Ackerly, 2004). Accordingly, small and narrow leaves and a high trichome density have been viewed as adaptations to the increasing dryness and seasonality of the Mediterranean region (Fiz-Palacios et al., 2006;Turini et al., 2010). Recent years have seen a growing interest in both spatial and temporal patterns of diversification and speciation of plant groups in the Mediterranean region. Researchers have also attempted to understand the changes in the morphological characters that have marked the course of evolution of these groups (Guzmán & Vargas, 2005, 2009Salvo et al., 2010). Despite this, there are very few groups of Mediterranean plants that have been studied from both perspectives (Guzmán & Vargas, 2005, 2009; Galbany-Casals diverse and one of the least known genera of the family, and therefore represents an interesting taxon to study the evolutionary processes in the Mediterranean Basin. The morphological differentiation of Fumana is mainly based on the presence of a whorl of sterile stamens and anatropous ovules arrangement (Spach, 1836a(Spach, , 1836b. Studies on morphological characteristics of vegetative and, principally, reproductive traits (inflorescence, stamens, ovules, pollen, and seeds) (Spach, 1836a;Willkomm, 1864;Grosser, 1903;Janchen, 1920;Jean & Pons, 1963;Güemes & Molero, 1993) have led to diverse proposals on the infrageneric organization of the genus that have been accepted to date. Species in the current Fumana genus have been divided in up to three subgenera (subgenera Fumana, Fumanopsis (Pomel) Janch, and Pomelina Maire) by Janchen (1920), Maire (1923), and Güemes & Molero (1993), which were treated in the genus category (genera Fumana, Fumanopsis Pomel, and Pomelina (Maire) Güemes & Raynaud) by Pomel (1860) and Raynaud (1992). To evaluate the coherence of these proposals,  included species of the three taxa in their molecular phylogeny of the Cistaceae. Their results raised doubts about the infrageneric divisions.
To date, there are 21 accepted species of Fumana, most of them having a circum-Mediterranean distribution. The genus is distributed north to south, from the island of Gotland (Sweden) (located in the parallel 57°N) to the Anti-Atlas in southern Morocco and Algeria (along 30°N parallel); and west to east, from Agadir (in the meridian 9°W) to the Urals (60°E meridian) (Grosser, 1903;Janchen, 1920Janchen, , 1925. Unlike other Cistaceae (Cistus and Helianthemum), Fumana is poorly represented in the Mediterranean islands, and missing in the eastern Atlantic oceanic islands (Canary Islands, Azores, Madeira). Only one species (F. procumbens (Dunal) Gren. & Godr.) reaches the Circumboreal region, and three species (F. arabica (L.) Spach, F. laevis (Cav.) Pau, and F. thymifolia (L.) Spach ex Webb) extend to the Saharo-Arabic region (Coode, 1965;Güemes & Molero, 1993). The distribution of the genus covers four biogeographic regions: Mediterranean, Irano-Turanian, Circumboreal, and Saharo-Arabic. The Mediterranean region has been proposed as the main center of diversification, especially the Iberian Peninsula, with 13 species, of which six are endemic to this territory. The Irano-Turanian region, especially the Anatolian peninsula, has been proposed as a secondary center of diversification, with 10 species, of which three are endemic (Janchen, 1920).
Despite previous scientific interest in understanding the evolution of Cistaceae in the Mediterranean region, the evolutionary history of Fumana has never been examined from a phylogenetic viewpoint or within a biogeographic context. The morphological characters of each of the 21 Fumana species are well known, which should allow evaluating their value for infrageneric classification in light of independently conducted molecular phylogenetic analyses. Therefore, firstly, we here perform a phylogenetic study of 19 of the 21 accepted species using nuclear (ITS) and plastid (matK, trnT-L) molecular markers, to test the infrageneric classification previously proposed for Fumana by Janchen (1920), Maire (1923, and Güemes & Molero (1993). Secondly, the diversification of the genus Fumana has been suggested to have started 5.3 myr ago (Aparicio et al., 2017); however, this date was based on estimates derived from investigation of the time of divergence of the main clades of Helianthemum. Therefore, given the reported age of diversification of Fumana (Aparicio et al., 2017), the distribution of current species (Grosser, 1903;Janchen, 1920Janchen, , 1925Coode, 1965;Güemes & Molero, 1993) and the high species diversity in the Iberian Peninsula (Güemes & Molero, 1993), we here also estimate the start of diversification of the genus and its main clades along with a reconstruction of their ancestral areas to test the hypothesis that the genus originated and diversified in the western Mediterranean in the late Miocene to Pliocene. Finally, we reconstruct selected character states in an ancestral state reconstruction to test whether the pattern of morphological evolution in the genus is related with the main paleo-climatic events that have occurred in the Mediterranean Basin. Specifically, we address the following questions: (i) is the phylogenetic reconstruction consistent with the previous systematic subdivisions of the genus (Janchen, 1920;Maire, 1923;Güemes & Molero, 1993)?;(ii) what are the phylogenetic relationships among Fumana taxa?; (iii) when and where did Fumana and its main clades diversify, and which biogeographic processes have affected the distribution patterns of the current species in the Mediterranean region?; (iv) which has been the pattern of evolution of the ancestral characters in the genus? 2 Material and Methods 2.1 Taxonomic sampling A total of 55 Fumana samples, representing 19 of the 21 species currently accepted (Coode, 1965;Heywood, 1968;Greuter et al., 1984;Güemes & Molero, 1993;Güemes, 1999) were used for the study ( Fig. 1; Table 1). Fumana grandiflora Jaub. & Spach and F. oligosperma Boiss. & Kotschy could not be sampled as neither species has been collected since their first description in the 19th century, and DNA extraction from the original herbarium collection was not possible. Species were represented by more than one population, with the exception of F. fontqueri Güemes, F. juniperina (Lax. ex Dunal) Pau, F. lacidulemiensis Güemes, and F. trisperma Hub.-Mor. & Reese, because each has only one known population. The populations were sampled throughout the geographic range of each species (Figs. S1A-S1E) according: to Güemes & Molero (1993, for the western Mediterranean species; to Coode (1965), for the eastern Mediterranean species; and to Heywood (1968), for the species occurring in the north of the Mediterranean. As part of the outgroup, we also newly collected three species of Cistus (one population for each species; Table 1). In addition to these newly collected samples, we included further relatives of Fumana in our phylogenetic analyses, following previous phylogenetic studies Aparicio et al., 2017; for details, see Section 2.3).
2.2 DNA extraction and amplification DNA was extracted from freshly collected leaves, subsequently dried and stored in silica gel, or from leaves taken from herbarium vouchers (Table 1). DNA was extracted with the Speedtools Plant DNA extraction Kit (Biotools, Spain) following the manufacturer's protocol but modifying the lysis step by adding 2-Mercaptoethanol and polyvinylpyrrolidone (PVP) till reaching a final concentration of 0.2% and 3%, respectively. Before DNA extraction, an extra step was added to reduce the excess of polysaccharides, which, based on a preliminary test, were thought to inhibit DNA isolation and amplification. Therefore, 20-30 mg of plant material was grinded and 1 mL of NaCl (5 M) was added. The material was shaken (vibrational frequency 30 Hz, 28.00 agitations per second, 90 s) and then centrifuged at maximum revolution  Güemes & Molero (1993, Coode (1965), andHeywood (1968) is shown in Fig. S1. Symbols represent species and numbers above symbols represent the number of haplotype (see Fig. 4). A, Geographic location of F. arabica and F. fontanesii (species grouped in Clade I). B, Geographic location of F. hispidula, F. juniperina, F. laevipes, F. laevis and F. thymifolia (species grouped in Clade II); C, Geographic location of F. aciphylla and F. bonapartei (species grouped in Clade III); D, Geographic location of F. baetica, F. ericifolia, F. ericoides, F. fontqueri, F. lacidulemiensis, F. paphlagonica, F. paradoxa, F. procumbens, F. scoparia, and F. trisperma (species grouped in Clade IV). Continued speed for 2 min in a standard table-top centrifuge. The NaCl solution was then removed and the steps repeated 2-4 times. DNA amplification of nuclear (ITS) and plastid regions (matK, trnT-L) was carried out by polymerase chain reaction (PCR). The ITS4 and ITS5 standard primers were used to amplify the ITS region (Sun et al., 1994), and the matK intron and the trnT-L intergenic spacer were amplified using the 3914F and 1470R (Johnson & Soltis, 1994) and the a and b primers (Taberlet et al., 1991), respectively. Amplifications were unsuccessful in many samples (15 for ITS, 37 for matK, and 31 for trnT-L), so we designed new internal primers for all regions based upon preliminary results, with two 21/20nucleotid-long internal primers for ITS (ITS-intF: 5′-GTT GCG TGA CGC CCA GGC AG-3′; ITS-intR: 5′-GAG CAC AGC CTC CGT GGC TAG-3′); and two 21/20-nucleotid-long internal primers for the matK region (matK-intF: 5′-GTC AAT TRA ATA AAT GGA TAG-3′; matK-intR: 5′-AGA GGA AGA CTC TTT TAM CC-3′). For the trnT-L region, just one 21-nucleotid-long internal primer was designed (trnTL-intF: 5′-GTA CAT ACG AAT TAC GCA AAC-3′), and then combined with the standard primers a, b, and d from Taberlet et al. (1991) (as reverse primers). DNA was amplified using a FlexCycler (AnalyticJena AG, Jena, Germany) or a 2720 ThermalCycler (AppliedBiosystems, Foster City, USA). After 4 min at 94°C pretreatment, PCR conditions were set as follows: 39 cycles of 1 min at 94°C, 1 min at 45-58°C, and 90 s at 72°C. We added 0.2-0.8 µL of 10 mg/mL BSA (bovine serum albumin) in a total of 20 µL reaction volume in all reactions and 0.2-1 µL DMSO (dimethyl sulfoxide) was only included in reactions for ITS amplification. The PCR products were purified using spin filter columns (QIAquick PCR Purification Kit, California), following the manufacturer's protocol. The cleaned product was then sequenced directly using dye terminators (Big Dye Terminator v. 2.0; Applied Biosystems, LitteleChalfront, UK) following the manufacturer's instructions and run on polyacrylamide electrophoresis gels (7%) using an Applied Biosystems Prism Model 3730 automated sequencer. We could not obtain 7 ITS (F. fontqueri 1, F. laevis 2, and F. paradoxa 2, F. procumbens 4, 5, 7, and 9) and 4 trnT-L sequences (F. procumbens 7, C. albidus, C. clusii, and C. creticus) because the amplification failed. Consequently, the number of newly sequences generated in this study for the ITS, matK, and trnT-L regions was, respectively, 51 (48 of Fumana, 3 of Cistus), 58 (55 of Fumana, 3 of Cistus), and 54 (all of Fumana). Sequences were deposited in GenBank under the accession numbers KJ534083 to KJ534245 (Table 1).

Phylogenetic analyses and plastid haplotype network
To perform phylogenetic analyses, two matrices were constructed: one with 70 sequences of the ITS region, and the other one with 61 sequences of the two concatenated plastid regions (matK, trnT-L). The ITS matrix contained the 48 Fumana (ingroup) and the three Cistus (outgroup) sequences obtained in this study, and also included as outgroups 19 sequences that were obtained from GenBank. Sequences of both nuclear and plastid regions were aligned using MAFFT v.6.822 (Katoh, 2008) hosted on the CIPRES Science Gateway (Miller et al., 2010), inspected and corrected manually on BioEdit v.7.0.9.0 (Hall, 1999) to minimize the number of gaps following the method of Kelchner (2000). Phylogenetic analyses were performed using maximum parsimony (MP), maximum likelihood (ML), and Bayesian inference (BI) separately for nuclear (ITS) and for the concatenated sequences (matK, trnT-L). The robustness for all nodes was estimated with posterior probability (PP) in BI and bootstrap values (BS) in both MP and ML.
For the MP analysis, the dataset was analyzed using equally weighted parsimony in TNT v.1.1 (Goloboff et al., 2008), with a heuristic search and a tree memory of 10 000. Gaps were treated as missing data in all analyses. We chose 1000 replicates of Wagner trees, followed by tree bisection-reconnection (TBR) branch swapping, and saved five trees per replication. A strict consensus tree was then generated. Nodal support was calculated using bootstrap resampling with 1000 replications summarizing the absolute frequency of each group.
To adjust the BI and ML analyses with the most proximate model available, we tested the simplest model of sequence evolution that best fit the sequence data via the bottom-up strategy of hierarchical likelihood ratio test (hLRT) and the Akaike information criterion (AIC, Akaike, 1979) using jModelTest v.0.1.1 (Guindon & Gascuel, 2003;Posada, 2008). The test was run separately for each of the three independent data sets (ITS, matK, trnT-L). The BI analysis was performed using MrBayes v.3.1.2 (Ronquist & Huelsenbeck, 2003), with concatenated plastid sequences treated as partitioned, implementing the best-fitted nucleotide substitution model (GTR + G for ITS; GTR + G for matK; GTR + G for trnT-L). Data matrices were run for 30 × 10 6 generations (with a burning of approximately 10%) on four Markov chain Monte Carlo (MCMC) chains with a temperature of 0.2 with sampling every 1000th generation. Tracer v.1.4  was then used to measure the effective sample size (ESS) of each parameter, which in all cases exceeded 100 and 50% majority-rule consensus tress were constructed. The ML analysis was calculated using RaXML v.7.0.3 (Randomized Axelerated Maximum Likelihood; Stamatakis, 2006) with the GTRGAMMA model and default settings for both nuclear (ITS) and plastid (matK, trnT-L) data sets separately, treating plastid sequence data set as partitioned. Ten ML searches were performed starting from 10 different randomized parsimony trees to obtain the best scoring tree. A standard nonparametric bootstrap with 100 replicates was carried out for internal support using the default estimation algorithm.
Homogeneity between the plastid and nuclear datasets was not tested using the incongruence length difference (ILD) test (Farris et al., 1994) since it is known to be susceptible to both false positives and false negatives (Pirie, 2015). Instead, Fumana accessions grouped in the main clades were inspected and compared between the nuclear and plastid phylogenetic reconstruction. The ITS tree produced a limited resolution but the BI, MP, and ML topology was congruent with those of the plastid tree. Then, we generated one matrix with concatenated ITS, matK, and trnT-L data, which included the sequences obtained in this study of 55 Fumana (ingroup) and of three Cistus (outgroup) accessions plus the sequences obtained from GenBank of the three Dipterocarpaceae accessions (outgroup) used before (see above). The matrix of concatenated ITS, matK, and trnT-L data was used to perform BI, MP, and ML analyses following exactly the same procedure describe above (with concatenated sequences treated as partitioned in BI and implementing the best-fitted nucleotide substitution model).
The number of plastid haplotypes and relationships among them were studied using the concatenated matK and trnT-L sequences of Fumana accesions for the main clades defined in the BI, MP, and ML plastid trees (which resulted in congruence among them). Then, four haplotype networks were constructed using species of the main clades (corresponding to Clade I, II, III, and IV) via the median joining algorithm (Bandelt et al., 1999) in the software PopArt (Population Analaysis with Reticulate Trees; http:// popart.otago.ac.nz).

Divergence time estimates and dispersal-vicariance analyses
Divergence time estimates for the matrix with concatenated ITS, matK, and trnT-L data sets were performed under BI using BEAST 1.6.1 . Xml-files for the BEAST analysis were constructed using BEAUti 1.6.1 (BEAST package) in which the datasets were analyzed under partition-specific models. For the two genetic data sets, we used the GTR + G model as the best fit substitution models. This model was obtained by the model test implemented in jModelTest v.0.1.1 software (Guindon & Gascuel, 2003;Posada, 2008), based on the AIC (Akaike, 1979). The data were analyzed under the uncorrelated lognormal relaxed clock model (UCLD), which is more likely to concede precise estimates than the uncorrelated relaxed clock model that presumes an exponential distribution of the evolutionary rates (Baele et al., 2012). Following the rigorous approach of Vanneste et al. (2014), a pure birth prior (Yule model) was employed in all runs, which assumes a constant speciation rate for each branch of the tree. We employed Yule prior as it is recommended for trees describing the relationships between individuals from different species (BEAST manual, version 1.4). According to Vanneste et al. (2014), we used a uniform prior between 0 and 100 for the Yule birth rate; an exponential prior with mean 0.5 on the rate heterogeneity parameter; an exponential prior with mean 1/3 on the standard deviation of the UCLD clock model; and a diffuse gamma prior with shape 0.001 and scale 1000 on the mean of the UCLD clock model.
The BEAST analysis was calibrated by using the same fossil records as described in . The tree root, consisting of the divergence time of Dipterocarpaceae and Cistaceae, was constrained with a minimum of 23 myr and a maximum of 39 myr, following Wikström et al. (2001). The prior for the age of the root was therefore set to a normal distribution with a mean of 31 myr and a standard deviation of 4.1 myr. We chose a normal distribution as it places higher probability on intermediate dates, providing a more appropriate prior calibration (Ho & Phillips, 2009). The stem of the crown group of the Cistaceae family was constrained using the earliest fossil found of the family of the Cistaceae. This fossil was described as a reproductive structure of Cistinocarpum roemeri Conis (Palibin, 1909) and dated from the Middle Oligocene in Germany (28 myr old). We therefore employed an exponential distribution prior with an offset of 28 myr and a mean of 1.5 myr.
The MCMC post chain was run with for 50 × 10 6 generations (with a burnin of approximately 10%) and sampled every 1000th generation. Tracer v.1.4  was then used to measure the effective sample size (ESS) of each parameter, which in all cases exceeded 100. Trees were then summarized with Tree Annotator v.1.6.1 (Rambaut & Drummond, 2010) as maximum clade credibility, mean node heights and a 0.5 posteriori probability limit. FigTree v.1.3.1 (Rambaut, 2009) was used to visualize the tree.
To reconstruct ancestral areas of distribution, a dispersalvicariance analysis (S-DIVA) was performed using RASP v.2.0 beta (Yu et al., 2010). This method resolves the phylogenic uncertainty of using a collection of trees. DIVA allows for the reconstruction of ancestral distributions without any previous assumptions about the area (Ronquist, 1997), and its use has been recommended under reticulated biogeographical scenarios, such as the Mediterranean Basin (Sanmartín, 2003;Oberprieler, 2005). After discarding 10 000 trees from a BI analysis of the matrix with concatenated ITS, matK, and trnT-L data sets, we employed a subsample of 20 000 trees with the slow ancestral reconstruction option selected to infer ancestral distribution areas. To define the areas, a paleographical criterion was followed (Meulenkamp & Sissingh, 2003) and the selected areas ( Fig. S1F; Table 1) were: A, north-western Mediterranean; B, south-western Mediterranean; C, southeastern Mediterranean; D, north-eastern Mediterranean; E, Eurosiberian. The biogeographical analysis was restricted to a maximum number of five areas, given that this is the maximum number of areas occupied by Fumana procumbens Gren. & Godr., the more widespread species. Outgroups were excluded from the analysis and were coded as "null," according to Yu et al. (2012).

Ancestral state reconstruction analysis
There are 15 morphological characters that have traditionally been considered for circumscription and differentiation in Fumana (Table S1). For the analysis of character evolution we chose six characters that are considered taxonomically important for species circumscription and for species relationships and one character (mucilage secretion in seeds), which has not been considered before, but potentially important to establish differences within Cistaceae (Engelbrecht et al., 2014). Therefore, characters of seed morphology (dispersal unit -seed versus fruit-, number of seeds per fruit, mucilage secretion and seed coat ornamentation), leaf morphology (margin and form) and types of trichomes (presence of glandular trichomes in stems and leaves/presence of non-glandular trichomes in stems and leaves) were analyzed and mapped on a pruned total evidence phylogeny. Character states were determined for each species from fresh and herbarium material. Finally, the complete morphological matrix was performed coding for a total of seven characters. We used the "drop.tip" command of the "ape" software (Paradis et al., 2013) in R v.3.0.1 (R Core Team, 2013) to prune the consensus tree of the concatenated data (ITS, matK, and trnT-L) obtained from the BI analysis, excluding repeated species. Fumana fontqueri was excluded from the analysis as we could not obtain the ITS sequences. This gave the species an imprecise position in the pruned tree and would have misled the outcome of the compete survey. To infer patterns of character evolution, we used the ML function of Mesquite v.2.74 (Maddison & Maddison, 2009) to trace character states on the tree. The "Trace Character History" option was used under the likelihood reconstruction method to display the ancestral state. The maximum likelihood model provides information on genetic branch lengths and uses the Markov k-state one parameter model (Mk1), which assumes a single rate for all transitions between character states (Lewis, 2001). Character states with a significant likelihood for reconstruction were considered the most likely ancestral states (i.e., using the average likelihood decision threshold of 2.0), with a proportional likelihood of 0.88 or higher (Maddison & Maddison, 2009).
ITS sequence data produced limited resolution with polytomies in all MP, BI, and ML phylogenetic analyses (Fig. 2). The BI analysis reached equilibrium after 100 000 generations, using GTR + G as the simplest model. The ITS tree depicted Fumana as a monophyletic group (0.99 PP, 100% BS in MP, and 88% BS in ML) and supported the monophyletic status of two species (F. fontanesii and F. laevipes), since all accessions corresponding to each species were grouped together with moderate to high PP and BS values (F. fontanesii: 1 PP, 95% BS in MP, 96% BS in ML; F. laevipes: 0.99 PP, 74% BS in MP, and 86% BS in ML).
Plastid sequences (matK and trnT-L) gave more resolution in the BI, MP, and ML phylogenetic reconstructions, with Fumana sequences forming a monophyletic group in the three analyses (1 PP, 77% BS in MP, and 70% BS in ML; Fig. S2). BI, MP, and ML analyses of plastid sequences yielded similar topology with BI displaying higher values. Using GTR + G as the simplest model, the BI analysis for the combined matK and trnT-L matrix reached equilibrium after 100 000 generations. Four conspecific accessions (F. aciphylla, F. bonapartei, F. fontanesii, and F. laevipes) formed well supported monophyletic groups in all three analyses (BI, MP, ML). The consensus tree of the three analyses of plastid regions revealed four major clades (named Clade I, II, III, and IV), supported with ≥0.90 PP and ≥70% BS values (with the exception of Clade III: with 0.84 PP and 61% BS in MP). MatK and trnT-L sequences revealed Clades I and II as sister clades, but with low PP (0.78) and no BS support (<50%) in MP, but moderate BS (74%) in ML; Clades III and IV were also sister clades with high PP and BS values (0.99 PP, 77% BS in MP, and 70% BS in ML).
The matrix of concatenated ITS, matK, and trnT-L sequences gave similar BI, MP, and ML results as those of the plastid matrix (Fig. 3), with similar PP and BS values for all clades (except for Clade III in the ML analysis). BI analysis for the combined nuclear and plastid data sets matrix reached equilibrium after 150 000 generations (GTR + G as the simplest model for each data set). Congruent with the BI, MP, and ML plastid tree, the three analyses for the concatenated ITS, matK, and trnT-L sequences resulted in high support for the monophyly of the genus Fumana (1 PP, 99% BS in MP, and 95% BS in ML) and the monophyly of the same four species (F. aciphylla, F. bonapartei, F. fontanesii, and F. laevipes). For all analyses, the tree obtained from the concatenated nuclear and plastid data sets also supported the presence of the same four major clades of species (Clade I, II, III, IV; with the exception of Clade III: with 54% BS in MP and <50% BS in ML), depicting Clades III and IV as sister clades (1 PP, 85% BS in MP, and 85% BS in ML).
For the four major clades of species (Clade I, II, III, IV), the haplotype network analysis revealed a moderate amount of variation at the concatenated plastid matrix (matK and trnT-L) (Figs. 1, 4). The majority of the species (13 species) did not reveal variation of the haplotypes or had haplotypes that differed by one to two mutational steps. However, a high amount of divergence was shown within F. arabica (grouped in Clade I), where each population had a unique haplotype (H1  and H2, Fig. 4) that were separated by 36 mutations. Clades I, II, and III revealed, respectively, 3 (H1 to H3), 6 (H4 to H9), and 4 (H10 to H13) different haplotypes, which were not shared among species within the clade (Fig. 4). Clade IV had 12 different haplotypes (H14 to H25), four of them shared among species within the clade. H14 was shared among accessions of F. fontqueri, F. paphlagonica, and F. procumbens, H21 was shared among accessions of F. ericoides and F. scoparia, H22 and H25 was shared among accessions of F. ericifolia and F. scoparia.
The dispersal-vicariance analysis ( Fig. 5; Table 3) estimated the north-western Mediterranean area as constituting the potential ancestral area, with a probability of 100%. Clades I, II, and IV are supported as having a north-western Mediterranean origin (all three with a probability of 100%), whereas Clade III is estimated to have a north-eastern Mediterranean (100% of probability) origin. For the major clades, our results also support (with ≥0.50% of probability) dispersal events in Clade I, II, and IV, and one vicariance event between Clade III and IV.

Ancestral state reconstruction analysis
Ancestral states of all seven characters were reconstructed for all nodes of the tree (Figs. 6, S3). The character state reconstruction showed that seed number and ornamentation were equivocally reconstructed (Fig. S3). We found a higher likelihood for a nine-seeded state as being ancestral, but it was not strongly supported. Species with six (reticulated) and nine (papillate) seeds, with two exceptions (F. bonapartei with six reticulated seeds and F. fontanesii with nine papillateseeds), were found to be well separated in Clades II and IV, respectively. Threeseeded F. aciphylla and F. trisperma were not located together in the same clade, but were separated in Clades III and IV, respectively. Both the two types of seed ornamentation are present in Clade I as well as in Clade III.  Table 1). Numbers above branches indicate Bayesian posterior probabilities (PP); bootstrap values (BS) from the maximum parsimony (MP) and maximum likelihood (ML) analyses are also indicated above branches (supported clades: ≥0.90 PP; ≥70% BS). A hyphen represents incongruence between BI tree and MP or ML consensus tree.  Table 1). Numbers above branches indicate Bayesian posterior probabilities (PP); bootstrap value (BS) from the maximum parsimony (MP) and maximum likelihood (ML) analyses are also indicated above branches (supported clades ≥0.90 PP; ≥70% BS). A hyphen represents incongruence between BI tree and MP or ML consensus tree.
Ancestral state reconstruction for the leaf margin and leaf form of the species of Fumana (Fig. 6) revealed that nonrevolute leaf margin was treated as the most likely ancestral character state, with a change to revolute leaf margin in three species of Clade II (in F. hispidula, F. thymifolia, and F. laevis). Focusing on the leaf form, a lanceolate leaf form was the most likely ancestral character state, which changed to ericoid in Clade IV, and to ovate in Clade I (in F. arabica) and Clade II (in F. thymifolia). Fumana laevipes in Clade II was the only species shifting to a filiform leaf shape. Glandular trichomes were reconstructed as being the ancestral state and changed twice to non glandular trichomes, once in Clade I (in F. fontanesii) and the other in Clade IV (in F. baetica, F. procumbens, and F. paphlagonica) (Fig. 6).
Seed dispersal has been reconstructed as the most likely state for the ancestor of Fumana while fruit dispersal is shown to be a derived character. The diaspore and dispersal mechanism only changed in three species in Clade IV (in F. baetica, F. procumbens, and F. paphlagonica) from seed to fruit dispersal (Fig. 6). Strong mucilage secretion in seeds has been reconstructed as the most likely ancestral state and is present throughout most of the clades (Fig. 6). There were two changes to a very weak or absent mucilage secretion of seeds in Clade IV (in F. baetica, F. procumbens, and F. paphlagonica) and Clade I (in F. fontanesii).

Phylogenetic analyses and systematic implications
Our phylogenetic reconstruction using ITS sequences of Fumana and several representatives of the remaining seven Cistaceae genera (Fig. 2) provides evidence that this genus forms a monophyletic group within the family. This is coherent with the morphological characters of Fumana, which also clearly differentiate this genus from the rest of the Cistaceae genera (sterile stamens and anatropous ovules arrangement) (Spach, 1836a(Spach, , 1836b. The generic division of Fumana has been accepted by some authors, but not by others (Table S2). Willkomm (1864), Maire (1923, and Güemes & Molero (1993), in their partial revisions of the genus (only with species distributed in the western Mediterranean), and Janchen (1920), in his global review, divided the genus into subgenera, sections, or in the case of Pomel (1860) and Raynaud (1992), in different genera. In contrast, Grosser (1903) did not accept any infrageneric classification. Our results with ITS show a lack of resolution in all phylogenetic analyses (Fig. 2); however, the plastid tree ( Fig. S2) and the combined data tree of nuclear and plastid sequences (Fig. 3) confirm the presence of the four major clades (Clade I, II, III, IV), comprised by the same species. The species grouped in these clades mostly coincide with those included in the sections of Janchen's (1920) infrageneric classification. That is, there is an important coincidence in the composition of species of: Clade I and sect. Platyphyllon Janch.; Clade II and sect. Helianthemoides Willk.; Clade III and sect. Megalosperma Janch.; and Clade IV and sect. Leiosperma Janch. However, there are also some inconsistencies. Specifically, F. fontanesii situated by Janchen (1920) within sect. Leiosperma (section coinciding with our Clade IV), is grouped in the Clade I in our phylogeny. Also F. trisperma, described after Janchen's (1920) taxonomy, according to this author, would be situated in sect. Megalosperma (section coinciding with our Clade III), but in our results it is grouped in Clade IV. These inconsistencies reveal that the four main clades of species found in this work cannot be defined with the characters that distinguish the sections of Janchen (1920). Similarly, the four clades presented here cannot be delimited using the vegetative, embryological, and palynological characters used in the division of the partial revisions of Pomel (1860), Willkomm (1864), or Maire (1923. Consequently, given the absence of morphological characters shared by the species of each clade, our proposal is in agreement with that of Grosser (1903) and is to reject any infrageneric division.
Focusing on the details of the taxonomic analyses, our tree of combined nuclear and plastid data supports the delimitation of four species, which are F. fontanesii, F. laevipes, F. aciphylla, and F. bonapartei, whose morphological characters strongly differ from each other and from the rest of species in the genus. These four species, which have always been considered by taxonomists as independent species, have unique haplotypes in the haplotype network and are also confirmed as monophyletic groups in the Chronogram obtained with Bayesian inference (BI) dating of the combined data (ITS, trnT-L, and matK) and ancestral area reconstruction using S-DIVA analysis. Ancestral area reconstructions and divergence times are only indicated for nodes with BI phylogenetic support (≥0.90 PP). The white squares represent ancestral area reconstructions with a probability ≥0.50%. Letters inside squares correspond to areas (A, north-western Mediterranean; AD, north-western/north-eastern Mediterranean; ABD, north-western/south-western/north-eastern Mediterranean). Contemporary distribution is denoted next to each taxon (A, north-western Mediterranean; B, south-western Mediterranean; C, south-eastern Mediterranean; D, northeastern Mediterranean; E, Eurosiberian). Names of clades or groups are given above squares. Time scale at the bottom in million years ago. phylogenetic reconstruction. However, accessions belonging to 11 species of the 15 species with more than one population sampled did not form monophyletic groups in any phylogenetic reconstruction. Specifically, the absence of monophyly in F. arabica (in Clade I) is consistent with its high haplotypic divergence. The haplotypes of the two populations studied, one continental (Greece) and the other insular (Cyprus), are unique in the species and are separated by more than 35 mutational steps. The continental and insular populations of this species are known to have morphological differences that have led to the delimitation of two varieties: F. arabica var. arabica and F. arabica var. incanescens Hausskn., accepted only in Eastern Mediterranean floras (Coode, 1965;Zohary, 1972;Meikle, 1977). In the light of the results these taxa should be further studied. Fumana laevis, F. hispidula, F. thymifolia (in Clade III), and F. baetica (in Clade IV) did not form monophyletic groups either. However, these species have little or nondivergent haplotypes that are not shared with other species, which would give more support to the delimitation of the species. In contrast, there are three pairs of non-monophyletic species (all in Clade IV) in which the two species of the pair share haplotypes with each other (F. ericifolia-F. paradoxa, F. ericoides-F. scoparia, and F. paphlagonica-F. procumbens). Although the species of these pairs show a clear morphological and ecological differentiation, our data are not valid to determine their recognition as distinct species.

Divergence times and ancestral areas
Based on molecular dating and on the reconstruction of ancestral areas, the diversification of Fumana (crown node) started in the Miocene (20.72 myr ago, 95% HPD: 13.63-28.08) in the Mediterranean (Fig. 5). The origin of the differentiation of the major clades (crown nodes) occurred between 14.25 myr ago (Clade IV, 95% HPD: 7.50-22.75) and 8.98 myr ago (Clade II, 95% HPD: 3.44-16.11). This is consistent with a recent review showing phylogenetic evidence for a Miocene origin of many Mediterranean plant groups with different life forms and biogeographic histories (Vargas et al., 2018) such as Erodium (Fiz-Palacios et al., 2010), Narcissus (Santos-Gally et al., 2012), or Saxifraga sect. Saxifraga (Vargas, 2000;Deng et al., 2015). The evolutionary trajectory of Fumana has probably been impacted by the main paleo-climatic events that have occurred in the Mediterranean Basin since the Miocene. In particular, there is evidence of the existence, perhaps temporally, of a proto-Mediterranean climate by the middle to late Miocene (Rundel et al., 2016), prior to the onset of the Mediterranean climate type (2.8 myr ago; Suc, 1984). The time of differentiation of the four main clades coincides with the period estimated for the proto-Mediterranean climate in the Mediterranean Basin.
The results of the reconstruction of ancestral areas specifically support a north-western Mediterranean ancestor for the genus Fumana ( Fig. 5; Table 3). The number of species and haplotypes as well as the geographic distribution of the current species also suggests that the north-western Mediterranean, probably south-eastern Iberia, was the main center of diversification of the genus and of species differentiation (13 species and 13 haplotypes in southeastern Iberia). Overall, the results allow us to describe two different patterns of spatio-temporal variation. On the one hand, Clade I and Clade II began their diversification, respectively, 11.07 myr ago (95% HPD: 4.64-19.55) and 8.98 myr ago (95% HPD: 3.44-16.11) in the north-western Mediterranean. Later dispersal events of ancestral and/or recent species of these clades could explain why many current species have a wide distribution throughout several regions of the Mediterranean. Alternatively, the sister clades III and IV show a different pattern. The most recent common ancestor to both clades was dispersed from the northwestern to north-eastern Mediterranean more than 17.63 myr ago (95% HPD: 9.99-27.27). In that period, there was a corridor of land in the north of the Mediterranean connecting the Tethys and the Paratethys seas which favored the dispersion of many linages in the Mediterranean region (Sanmartín, 2003;Thompson, 2005). Subsequently, our data supports a process of geographical vicariance between north-eastern and north-western Mediterranean that separated these two sister clades. Clade IV initially differentiated in the north-western Mediterranean, and experienced posterior dispersal events towards other regions of the Mediterranean. In contrast, Clade III began its diversification in the north-eastern Mediterranean, probably without further dispersal events to other areas since the current species are restricted to this area.

Ancestral character states
The onset of diversification of Fumana (20.72 myr ago) coincides with a time of a warm and mainly subtropical climate in the Mediterranean Basin (Utescher et al., 2000;Agustí & Antón, 2002;Böhme, 2003). During this period, desert conditions expanded against tropical forests, making new niches for plants available (Jiménez-Moreno & Suc, 2007). A large development of "warm steppes" in southeast Iberia has been documented in the Miocene (Jimenez-Moreno et al., 2010). There is paleobotanical evidence suggesting that xerophytic Mediterranean vegetation arose in this period (Jiménez-Moreno & Suc, 2007;Barrón et al., 2010;Rundel et al., 2016). These findings are consistent with our character reconstruction results which provide evidence that the most common ancestor of Fumana had leaf traits (narrow leaves with glandular trichomes) that are considered typical xeromorphic adaptations (Rudall, 1980;Moon et al., 2009). Strong mucilage secretion was also an ancestral state in Fumana (Fig. 6). Seeds of many Fumana species produce a sticky and thick mucilaginous layer around the seed coat when they come into contact with water, and this adheres the seeds to the surface they reach upon drying (Grubert, 1974;Engelbrecht et al., 2014). This works as antitelechoric dispersal mechanism adhering the seeds to the ground, which after drying out, remain glued to the soil. It has been proposed that strong mucilage secretion is a selective response to soil erosion and a mechanism for preventing ant predation in open semiarid shrublands (Engelbrecht & García-Fayos, 2012;Engelbrecht et al., 2014). This trait could probably have been a favorable state of character facing erosive conditions in the Mediterranean during the Miocene. The origin of the main four clades in the middle-late Miocene took place without changes in the state of the ancestral characters, with the exception of Clade IV where the shape of the leaf evolved from lanceolate to ericoid. Moreover, a loss of glandular trichomes, a reduction in mucilage secretion, and a change in the dispersal unit evolved within one lineage of Clade IV (Fig. 5, Procumbens group). This clade began to diversify 9.73 myr ago (95%: 4.65-16.09), after the occurrence of an important cooling event (middle Miocene Climate Transition~15 to~13.7 myr) in large parts of the western Mediterranean Basin (Jimenez-Moreno et al., 2010;Hernández-Ballarín & Peláez-Campomanes, 2017). The extant species of this lineage (F. baetica, F. fontqueri, F. paphlagonica, and F. procumbens) currently inhabit moister and fresher environments of high mountains. As these species live in rocky habitats where seedling growth on rocks and stones is impossible, strong mucilage secretion could increase the risk of gluing the seeds onto an unfavorable surface, thus impeding germination and growth. In these cases, dispersal of seeds inside the fruit may increase the possibility of seedling recruitment far from rocks. Regardless of the changes in the traits mentioned above, the evolution of the characters in Fumana cannot be considered as dynamic as that of Cistus, the only genus of Cistaceae previously studied. In Cistus, whose diversification originated after the onset of the Mediterranean climate (1.56 myr ago), a dynamic evolution of characters and an adaptive radiation in the white flowered lineage has been documented . In contrast, in Fumana, whose diversification took place earlier, many of the ancestral states of the characters studied here are still present in most of the current species.
In summary, our results do not resolve phylogenetic relationships at species level in Fumana, which would require the use of a battery of high-resolution molecular markers. However, our work supports Fumana as a genus with its evolutionary history placed far before the onset of the Mediterranean climate. The genus has probably originated in the north-western Mediterranean during the early Miocene in regions that were under extreme dry conditions in an environment of expanding warm steppes. This plant genus gave rise to several species that effectively stayed in the Mediterranean area, in some cases with an extensive distribution and in other cases in very small areas. eastern Mediterranean; D, north-eastern Mediterranean; E, Eurosiberian. Fig. S2. Majority rule consensus tree (50%) from BI analysis based on plastid sequences (trnT-L, matK) of Fumana species. Population numbers are given after species name (see Table 1). Numbers above branches indicate Bayesian posterior probabilities (PP); bootstrap value (BS) from the Maximum Parsimony and Maximum Likelihood analyses are also indicated above branches (supported clades ≥0.90 PP; ≥70% BS). A hyphen represents incongruence between BI tree and MP or ML consensus tree. Fig. S3. Likelihood-based ancestral state reconstruction of seed number and ornamentation. Proportional likelihoods of the most likely state are shown at nodes for all species and clades. Character states are mapped onto the majority rule consensus tree (50%) from BI analysis based on plastid and nuclear sequences (ITS, trnT-L, and matK). Table S1. Morphological characters of the species of the genus Fumana (Dunal) Spach. Table S2. Different systematic treatments of the current species of the genus Fumana (Dunal) Spach (Cistaceae).