Phylogeography of cliff-dwelling relicts with a highly narrow and disjunct distribution in the western Mediterranean

1538 • AMERICAN JOURNAL OF BOTANY 102 ( 9 ): 1538 – 1551 , 2015 ; http://www.amjbot.org/ © 2015 Botanical Society of America AMERICAN JOURNAL OF BOTANY RESEARCH ARTICLE Mediterranean endemic plant species provide fascinating material for the study of phylogeography ( Th ompson, 2005 ). One of the most interesting regions is the plant biodiversity Baetic-Rifan hotspot ( Médail and Quézel, 1997 ), which includes European and African areas separated by the Mediterranean Sea and the Strait of Gibraltar. Th e high biodiversity of this region is explained by a complex geologic and climatic history that gave rise to highly diverse habitats ( Médail and Quézel, 1997 ; Th ompson, 2005 ). However, this region is characterized from a fl oristic point of view by a high percentage of species common to the European and African ranges (about 75%; Valdés, 1991 ). Th us, the Baetic-Rifan hotspot symbolizes an intricate scenario from a biogeographical perspective ( Molina-Venegas et al., 2013 ). Th e internal zones of the Baetic and Rifan orogens originated from a common terrain located between the Iberian Peninsula and southern France during the Oligocene ( ≈ 30 Ma), which experienced drift and fragmentation into microplates during the Miocene ( Lonergan and White, 1997 ; Rosenbaum et al., 2002 ). When the Baetic-Rifan microplate reached the southwestern Mediterranean (21–18 Ma), it remained separated from the Iberian Peninsula and Africa by the Baetic and Rifan corridors, respectively. Th e positive relief associated with the fold-and-thrust belt at the western frontal zone of the microplate gave rise to a stretch of land that, by the 1 Manuscript received 4 April 2015; revision accepted 20 July 2015. 2 Departamento de Biología Vegetal y Ecología, Universidad de Sevilla, Apartado 1095 41080 Sevilla, Spain; and 3 Departament of Biological Sciences, Sungkyunkwan University, 440–746, Suwon, Korea 4,5 Authors for correspondence: jl[email protected]; sonc[email protected]u, phone +82-31-299-4499, fax +82-31-290-7015 doi:10.3732/ajb.1500152 Phylogeography of cliff -dwelling relicts with a highly narrow and disjunct distribution in the western Mediterranean 1 Jose L. Silva 2,4 , So-Yeon Lim 3 , Seung-Chul Kim 3,5 , and José A. Mejías 2 PREMISE OF THE STUDY: The high biodiversity in the Baetic-Rifan hotspot of Mediterranean region is shaped by complex geological and climatic histories and has been a subject of recent intensive studies. However, very little is known about phylogenetic and biogeographic history of three rare and critically endangered cliff -dwelling species of Sonchus in section Pustulati in this region. METHODS: We investigated the genetic variation and phylogenetic relationships of populations based on nuclear (ITS/ETS) and plastid (3 ′ trnL–ndhJ/ psaI–accD ) DNA sequences, and amplifi ed fragment length polymorphisms (AFLPs). We performed a Bayesian relaxed molecular clock analysis with ITS data to estimate divergence times for major lineages. KEY RESULTS: ITS/ETS and AFLP phylogenies showed high concordance and contrasted with cpDNA data. The divergence between S. masguindalii and S. fragilis was dated at 5.48 Ma, between S. fragilis and S. pustulatus at 3.89 Ma, and between the Baetic and Rifan S. pustulatus at 1.18 Ma. Within each distribution area, AFLP data showed a relatively high genetic structuring and moderate genetic diversity, the latter being impoverished in the Baetic populations. CONCLUSIONS: Our results further confi rm the hybrid origin of S. pustulatus , a critically endangered species. The origin and diversifi cation of lineages appear to have occurred on the temporary land bridge that joined Iberian and North Africa during the Messinian Salinity Crisis (5.96–5.33 Ma) and the subsequent Zanclean fl ood that progressively refi lled the Mediterranean Basin (5.33–3.60 Ma). The only Baetic populations of S. pustulatus most likely originated from the Rifan ones. KEY WORDS Baetic-Rifan region; cliff ecology; disjunct distribution; Gibraltar arc; long-distance dispersal; Messinian-Zanclean ages; relaxed Bayesian molecular clock dating; Sonchus section Pustulati ; vicariance SEPTEMBER 2015 , VOLUME 102 • SILVA ET AL. —PHYLOGEOGRAPHY OF SONCHUS SECTION PUSTULATI • 1539 Middle Miocene (15 Ma), probably joined Africa and Iberia in a position in the middle of the present-day Alboran Sea ( Rosenbaum et al., 2002 ). As the eastward-dipping subduction zone migrated to the west and the Alboran domain underwent progressive back collapse, this land bridge/stretch migrated westward until it reached its current position in the Gibraltar Arc (Late Miocene, 10 Ma). Closure of the Rif and Betic corridors accompanied the process, fi - nally leading to the separation of the Mediterranean and Atlantic seas during the Messinian (6 Ma). Th is promoted the desiccation of the Mediterranean Sea (Messinian Salinity Crisis, 5.96–5.33 Ma; Krijgsman et al., 1999 ; Fauquette et al., 2006 ) and facilitated intercontinental expansion of plant distributions until the opening of the Strait of Gibraltar (5.33 Ma; e.g., Caujapé-Castells and Jansen, 2003 ; Cano-Maqueda et al., 2008 ) and the Zanclean fl ood that progressively refi lled the Mediterranean basin (5.33–3.60 Ma). Th ereaft er, the sea barrier likely made intercontinental migration of plants highly dependent on long-distance seed dispersal and favored vicariance processes in those plants present in both sides before the total refi lling of the Mediterranean, i.e., leading to diff erentiation in pairs of closely related species. Most phylogeographical studies in the Baetic-Rifan region focus upon the role of the Strait of Gibraltar as a bridge or a barrier for species migration and gene fl ow between the northern (European) and southern (African) ranges (see review of Rodríguez-S á nchez et al., 2008 ). Apart from addressing the disjunct populations around the Strait of Gibraltar, phylogenetic/phylogeographic investigation of additional Baetic-Rifan plant taxa that do not occur on both sides of such strait provides us with insights into the role of sea barriers as a limiting factor for range expansion. Incorporating timecalibrated phylogenetic analyses of taxa with disjunct distributions is crucial with regard to fully understanding dispersal, colonization, and isolation processes in the region (e.g., Caujapé-Castells and Jansen, 2003 ; Casimiro-Soriguer et al., 2010 ; Fern á ndez-Mazuecos and Vargas, 2011 ). In addition, analyses of population genetic structure and diversity enable us to determine the degree of connection between them and the regions that potentially acted as the centers of diversifi cation (e.g., Casimiro-Soriguer et al., 2010 ). Furthermore, it is expected that the genetic landscape in species with high habitat specifi city and low dispersal ability, as in many cliff - dwelling species (e.g., García et al., 2002 ; Picó and Riba, 2002 ), is highly structured among their populations. Sonchus pustulatus Willk., S. fragilis Ball, and S. masguindalii Pau and Font Quer constitute a group of rare, long-lived, cliff -dwelling species (Asteraceae), endemic to narrow areas within the Baetic- Rifan hotspot ( Fig. 1a ; Boulos, 1973 ; Silva et al., 2015 ). Th ey belong to the well-supported section Pustulati within the Sonchus subgenus Sonchus ( Boulos, 1973 ; Kim et al., 2007 ), which seems to be an ancient clade within the group ( Kim et al., 2008 ). All three species are distributed in the Rifan range, where they have been considered as very rare ( Fennane and Ibn-Tattou, 1998 ), and S. pustulatus is the only species of the section occurring in the Baetic range, where it is categorized as critically endangered ( Bañares et al., 2004 ; Silva et al., 2015 ). Due to the highly restricted and disjunct distribution, the presumably low dispersal ability, and the ecological particularities of these cliff plants ( Silva, 2014 ), the section Pustulati is a good system to study intercontinental colonization processes and their relationship with the paleogeography of the Baetic-Rifan region. In the current study, we attempted to elucidate the taxonomic position of these species within the section Pustulati and to test whether the origin and diversifi cation of the species are related with the most important past geologic and climatic events occurred in the Baetic-Rifan region. We conducted phylogenetic and phylogeographic studies in Sonchus section Pustulati based on noncoding cpDNA (3 ′ trnL – ndhJ / psaI – accD ) and nrDNA ITS/ETS sequences, and amplifi ed fragment length polymorphism (AFLP). Th e following specifi c questions were addressed: (1) How consistent are phylogenetic analyses with the current taxonomic treatment? (2) Could the origin and diversifi cation of the species be related to the Messinian Salinity Crisis, the subsequent Zanclean fl ood that progressively refi lled the Mediterranean basin and/or the onset of the Mediterranean climate? (3) Could the current disjunct distribution of S. pustulatus have resulted from historical vicariance events rather than recent long-distance dispersal events between both sides of the Baetic-Rifan region? (4) What evidence about recent and current biogeographic trends of these species within each distribution area can be inferred from their current population genetic structuring and known biological/habitat traits? And fi nally, (5) is the narrow distribution range of these rare species associated with low genetic diversity? MATERIAL AND METHODS The species of Sonchus section Pustulati— Th e 19 extant populations of Sonchus section Pustulati are distinguished in fi ve small distribution areas restricted to the Baetic-Rifan Internal Zone ( Silva et al., 2015 ; Fig. 1a ). Th e rare Baetic populations of S. pustulatus are located in the environs of Almeria city ( Fig. 1a , area A, southeastern Spain), whereas the Rifan ones colonize cliff s in the Oued Laou valley and some neighboring valleys ( Fig. 1a , area B, western Rif). Sonchus fragilis is exclusive to the mountains surrounding the city of Tetouan ( Fig. 1a , area C, northwestern Rif). Finally, S. masguindalii is distributed along the coast in the Bokkoya Mountains and some sea cliff s in the environs of Al-Hoceima city ( Fig. 1a , area D, central Rif). Some herbarium specimens of S. pustulatus from northwestern Algeria collected in the 19th century (e.g., herbarium COI), with standard morphological characters of the species, indicate that this species was also present in rocky places on the oceanfront near Ghazaouet (formerly Nemours; Fig. 1a , area E, northwestern Algeria), but its persistence in the area has not been determined. Th ese endemic species are found in a highly fragmented habitat and are specifi cally associated with north-facing cliff s, which are at a low to medium altitude, have an alkaline, mainly limestone substrate, and are located on the coast or less than 20 km from the sea ( Silva, 2014 ; Silva et al., 2015 ). Th ese plants are polycarpic perennials, mostly self-incompatible (except for S. fragilis that has a high incidence of self-compatibility), with generalist entomophilous pollination ( Silva, 2014 ). Achenes have a short-lasting pappus and are heavier than other more widespread Sonchus species ( Silva, 2014 ). Populations sampled and DNA extraction— Leaf tissue samples were collected from 15 to 20 georeferenced individuals per population (using GPS). To obtain representative samples of whole populations and to avoid collecting siblings, we made a linear transect along the cliff s and collected one sample every 10–15 m. Th e leaves were immediately dried and stored in silica gel at room temperature until further processing. Voucher specimens were deposited in the University of Seville herbarium (SEV; Appendix S1, see Supplemental Data with the online version of this article). A total of 316 individuals were used for the molecular analyses. Total genomic 1540 • AMERICAN JOURNAL OF BOTANY DNA was isolated from ca. 1 cm 2 of leaf tissue with the Nucleo Spin Plant II Kit (Macherey-Nagel GmbH and Co. KG, Düren, Germany) according to the manufacturer’s protocol. DNA concentration and quality of each sample were checked on 0.9% agarose gels and quantifi ed in a NanoDrop1000 spectrophotometer (Th ermoScientifi c, Wilmington, USA). DNA was diluted to a fi nal concentration of 30 ng/μL. DNA sequencing— For cpDNA sequencing, 92 individuals were selected (3–5 per population; Table 1 ), plus three individuals of the FIGURE 1 (a) Distribution range of the species of Sonchus section Pustulati , locations of all known populations, and area of occupancy ( Silva et al., 2015 ). Area A (blue), Baetic populations of S. pustulatus ; Rifan populations: area B (red), of S. pustulatus ; area C (orange), of S. fragilis ; area D (green), of S. masguindalii ; and area E, unconfi rmed current presence of S. pustulatus (dark gray). (b) Statistical parsimony network of the six plastid haplotypes found in the section. Circle: haplotype, line: nucleotide substitution, small black dot: inferred intermediate unsampled haplotype. SEPTEMBER 2015 , VOLUME 102 • SILVA ET AL. —PHYLOGEOGRAPHY OF SONCHUS SECTION PUSTULATI • 1541 related species S. palustris , which was used as the outgroup. Th e previous phylogenetic study of Sonchinae ( Kim et al., 2007 ) strongly suggested that S. palustris is the closest sister lineage of section Pustulati , and thus it was used as a sole outgroup. To sequence the ITS/ETS region of nrDNA, we selected three individuals from each population, for a total of 60 individuals including three of the same outgroup species. For cpDNA, of the 34 noncoding regions that were initially surveyed ( Shaw et al., 2007 ), the two regions with the highest potentially informative character value were selected, i.e., the 3 ′ trnL–ndhJ (UAA) [=Tab E] and accD – psaI . PCR amplifi cation of nrDNA ITS and ETS was performed as described by Lee et al. (2005) , except that the ETS region was amplifi ed with the use of the primers ETS1f (5 ′ –CTTTTTGTGCATAATGTATATATAGGGGG–3 ′ ) and 18S –2L (5 ′ –TGACTACTGGCAGGATCAACCAG–3 ′ ) designed by Linder et al. (2000) . Th e same PCR conditions for two cpDNA regions were applied. PCR products were purified with GENEALL EXPIN PCR SV (GeneAll Biotechnology, Seoul, Korea). Sequencing was conducted at the Geno Tech Corp. (Seoul, Korea). Base calling and sequence editing were performed with SEQUENCHER v4.2.2 (Gene Codes, Ann Arbor, Michigan, USA). Sequences were aligned manually with MacClade (v.4.06, Maddison and Maddison, 2003 ). AFLP reaction— For the AFLPs, we analyzed a total of 167 individuals, ranging from 4 to 12 per population ( Table 1 ). We followed the protocol of Vos et al. (1995) with some modifi cations ( Lauterbach et al., 2011 ). A total of 18 AFLP primer combinations were tested for their capacity to detect polymorphisms and their reproducibility of banding patterns in all the Sonchus taxa studied. Th e three most polymorphic primer combinations were used for selective PCR: EcoRI ACA– Tru1l CTG; EcoRI ACA– Tru1l CAG; EcoRI ACC– Tru1l CAA. One sample from each distribution area (3% of the total) was repeated on all diff erent plates to calculate the average reproducibility of the method ( Bonin et al., 2004 ). Th e AFLP fragments were separated on a polyacrylamide gel with an internal size standard (GenomeLab DNA Size Standard Kit 400, Beckman Coulter, Krefeld, Germany) by means of an automated sequencer (CEQ 8000, Beckman Coulter). Data fi les were imported into Genographer software (v1.6.0, J. J. Benham, Montana State University, Bozeman, Montana, USA) and fragments from 75 to 450 bp were identifi ed and scored automatically for their presence and absence. A negative control (no DNA) was included in each plate. Finally, fragments were checked individually, without prior information on their origin, by means of the “thumbnail” function, and a presence/ absence matrix was exported for further analyses. Analyses of cpDNA and nrDNA ITS/ETS sequences— Th ree diff erent phylogenetic analyses were conducted for both cpDNA (two combined noncoding regions) and nrDNA (combined ITS and ETS regions) data sets. We fi rst conducted an equally weighted unordered maximum parsimony (MP) approach ( Fitch, 1971 ) implemented in the program PAUP* v4.0b10 ( Swoff ord, 2002 ). Th e MP analyses included a default heuristic search option with tree-bisection- reconnection (TBR) branch swapping and MULPARS on. Bootstrap support (BS) was calculated by bootstrap analysis from 1000 replicates with the same heuristic options. Secondly, we conducted TABLE 1. Populations of Sonchus section Pustulati and AFLP results. N , population size ( Silva et al., 2015 ); n 1 , n 2 and n 3 , number of analyzed individuals for cpDNA and ITS/ETS nrDNA sequencing and AFLP, respectively. %P, percentage of polymorphic loci; Fr PRI , number of private fragments. DW, rarity index: a , based on all plants analyzed, b , based on populations of S. pustulatus and S. fragilis , and c , based on populations of S. pustulatus . DW values signifi cantly higher ( + ) or lower ( - ) than those to be expected by chance appear in bold. H E , average gene diversity. DW Species and location, m a.s.l. Code Coordinates (N / W) N n 1 n 2 n 3 %P Fr PRI a b c H E ± SE A) Spain, SE Baetic range Sonchus pustulatus 898 15 9 34 0.109 ± 0.056 Bco. San Antonio, Aguadulce, 300 m ANT 36 ° 49.8 ′ / 2 ° 34.4 ′ 484 5 3 12 41.9 1 1.8 2.8 3.2 0.088 ± 0.014 Bco. San Telmo, Almeria, 22 m TEL 36 ° 49.7 ′ / 2 ° 29.0 ′ 112 5 3 11 45.1 1.8 2.6 - 2.9 - 0.101 ± 0.017 Bco. Caballar, Almeria, 100 m CAB 36 ° 50.6 ′ / 2 ° 28.7 ′ 292 5 3 11 47.8 2 1.9 2.6 - 2.9 - 0.093 ± 0.016 B) Morocco, W Rif Sonchus pustulatus 9991 24 15 48 0.167 ± 0.082 Oued Laou valley, 120 m LAUF 35 ° 23.3 ′ / 5 ° 12.6 ′ 1062 5 3 12 50.6 2.3 + 3.6 + 4.2 + 0.160 ± 0.025 Oued Laou valley, 153 m LAUL 35 ° 20.8 ′ / 5 ° 11.0 ′ 540 5 3 11 52.6 2.0 2.9 3.2 0.138 ± 0.023 Oued Laou valley, 150 m LAUO 35 ° 17.4 ′ / 5 ° 13.8 ′ 221 4 3 4 44.3 2.0 3.0 4.3 + 0.145 ± 0.048 Oued Laou valley, 235 m TAL 35 ° 16.0 ′ / 5 ° 13.8 ′ 1980 5 3 10 48.0 — — — 0.146 ± 0.025 Oued Al-Kannar defi le, 244 m KAN 35 ° 13.0 ′ / 5 ° 01.2 ′ 6188 5 3 11 52.2 1.8 2.7 3.2 0.155 ± 0.025 C) Morocco, NW Rif Sonchus fragilis 33 927 18 12 11 0.138 ± 0.074 Montes Ghorghiz, Tetouan, 550 m GHO1 35 ° 32.1 ′ / 5 ° 23.1 ′ 26 985 5 3 6 45.8 1 2.6 + 4.5 + 0.146 ± 0.036 Montes Ghorghiz, Tetouan, 550 m GHO2 35 ° 32.8 ′ / 5 ° 22.6 ′ — 3 3 5 41.9 2.5 + 4.6 + 0.120 ± 0.034 Montes Ghorghiz, Tetouan, 550 m GHO3 35 ° 32.3 ′ / 5 ° 22.3 ′ 661 5 3 — — — — — — Montes Dersa, Tetouan, 275 m DER 35 ° 35.9 ′ / 5 ° 24.6 ′ 6282 5 3 — — — — — — D) Morocco, central Rif Sonchus masguindalii 42 269 35 21 74 0.177 ± 0.086 Torres de Alcal á beach, s.l. TOR 35 ° 09.4 ′ / 4 ° 19.7 ′ 6244 5 3 10 44.3 1 1.8 - 0.146 ± 0.025 Bades beach—Peñón Vélez, s.l. BAD1 35 ° 10.2 ′ / 4 ° 17.9 ′ 6444 5 3 10 50.6 2 2.1 0.139 ± 0.024 Valley to Bades beach, 60 m BAD2 35 ° 09.6 ′ / 4 ° 16.9 ′ 638 5 3 12 43.5 1.8 0.138 ± 0.021 Boumahdi beach, s.l. BOU 35 ° 14.0 ′ / 4 ° 00.7 ′ 19 513 5 3 12 45.1 1.6 0.117 ± 0.018 Cebadilla beach, Al-Hoceima, s.l. ALH 35 ° 14.6 ′ / 3 ° 58.0 ′ 688 5 3 9 42.3 1.5 - 0.096 ± 0.018 Quemado beach, Al-Hoceima, s.l. QUE 35 ° 14.5 ′ / 3 ° 55.5 ′ 6708 5 3 11 52.2 3 2.2 + 0.157 ± 0.025 Sfi ha beach, Al-Hoceima, s.l. SFI 36 ° 12.8 ′ / 3 ° 54.3 ′ 2035 5 3 10 50.2 2.0 0.169 ± 0.029 1542 • AMERICAN JOURNAL OF BOTANY maximum likelihood (ML) analyses, in which optimal models of molecular evolution were chosen with the likelihood ratio test ( Whelan and Goldman, 1999 ) implemented in the program MODELTEST v3.7 ( Posada and Crandall, 1998 ). Model parameters were then imported into PAUP*, and a heuristic search (asis sequence addition, TBR branch swapping, and MULPARS option on) was executed. We also conducted ML bootstrap analyses with 1000 replicates to determine the support for each clade. Th irdly, we performed Bayesian analyses (BI) using the program MRBAYES v3.1.2 ( Huelsenbeck and Ronquist, 2001 ; Ronquist and Huelsenbeck, 2003 ). We calculated the likelihood parameters for the Bayesian analysis using the program MRMODELTEST v2.2 ( Nylander, 2004 ). We ran the Bayesian Markov chain Monte Carlo (MCMC) algorithm for 4 000 000 generations, with four simultaneous chains (three “cold” and one “heated”), starting from random trees and sampling every 100 generations. The program TRACER v1.5 ( Rambaut and Drummond, 2009 ) was used to evaluate the burn-in and to examine log likelihoods, ensuring that the run was in the stationary phase and that adequate eff ective sample sizes (ESS) were attained. A conservative 20% burn-in was removed from the sampled set of trees and a 50% majority-rule consensus tree was generated from the remaining trees. For estimates of Bayesian clade support, we considered strong support for values of posterior probability (PP) ≥ 0.95, moderate support for 0.90 ≤ PP ≥ 0.94, and weak to no support for PP ≤ 0.89. A haplotype network was also constructed from the cpDNA data set. Due to the absence of homoplasy in this data set, based on the consistency index (CI) of the MP analysis (see Results), haplotypes were networked manually using the principle of maximum parsimony. Th e manual network of haplotype was confi rmed by the statistical parsimony implemented in the program TCS version 1.21 ( Clement et al., 2000 ). Gaps were treated as missing data, and the connection limit was set to 95% in accordance with Hart and Sunday (2007) . To calibrate the phylogeny, we performed a Bayesian relaxedclock approach using the uncorrelated-rates model implemented in the program BEAST v1.7.3 ( Drummond et al., 2006 ; Drummond and Rambaut, 2007 ). Our newly generated sequences from Sonchus section Pustulati were combined with the ITS data set of the subtribe Sonchinae available in GenBank, which were generated by Kim et al. (2007) . We identifi ed the best-fi t substitution model for the ITS data (GTR+I+G) with MRMODELTEST v2.2 ( Nylander, 2004 ). Th e input data were compiled in the program BEAUTI v1.7.3 ( Drummond et al., 2006 ; Drummond and Rambaut, 2007 ) with the two prior sets as follows: (1) the age for the most recent common ancestor of genus Dendroseris in the Juan Fernandez Islands; normal prior distribution with a mean of 3.3 Ma and standard deviation of 0.4 (giving a 95% confi dence interval ranging from 2.6–4.0 Ma; Sang et al., 1994 ); (2) the age for the monophyletic clade of Reichardia , early-diverged genus within the Sonchinae; uniform prior distribution with a mean of 8.1 Ma and standard deviation of 2.8 (giving a 95% confi dence interval of 3.6–13.6 Ma; S.-C. Kim, unpublished data). Th e clock model was set to relaxed uncorrelated lognormal, and the Yule process was chosen as the speciation process. Posterior distributions for each parameter were estimated by means of a MCMC run for 40 000 000 generations with parameters logged every 5000 generations. Th e output log fi les were analyzed with TRACER v1.5 ( Rambaut and Drummond, 2009 ) to assess convergence and to confi rm that the eff ective sample sizes (ESS) for all parameters were larger than 200, ensuring that the MCMC had run long enough to produce a valid estimate of the parameters. We discarded 10% of burn-in (i.e., 800 trees) trees based on the TRACER results and produced a maximum credibility tree using TreeAnnotator v1.7.3 ( Drummond et al., 2006 ; Drummond and Rambaut, 2007 ). AFLP analyses— To assess population genetic diversity from the presence/absence matrix, we calculated the percentage of polymorphic fragments (%P) with the program AFLP-SURV v1.0 ( Vekemans, 2002 ) and the number of private fragments (Fr PRI ) with the program FAMD v1.25 ( Schlüter and Harris, 2006 ). Average gene diversity ( H E ) at the population level was computed with the program Arlequin v3.11 ( Excoffi er et al., 2005 ). We determined an additional measure of isolation degree, the Rarity 1 index (equivalent to the frequency of down-weighted marker values; i.e., DW sensu Schönswetter and Tribsch, 2005 ) with R v2.14.1 ( R Development Core Team, 2011 ) and AFLPDAT functions ( Ehrich, 2006 ). Th e values were calculated from 100 000 repetitions by means of the R function Rarity.permut, which implements a permutation approach to assess whether rarity values for populations are higher or lower than what is to be expected by chance, assuming that all markers are distributed in all populations at random. To estimate amongpopulations relationships, we computed Nei’s pairwise genetic distance data and 10 000 bootstrap distance matrices using AFLPSURV v1.0 ( Vekemans, 2002 ). We used these genetic distances to build a neighbor-joining tree (PHYLIP soft ware; Felsenstein, 1993 ). We employed two Bayesian approaches to infer the number of distinct genetic clusters ( K ) in our data using structure v2.3.1 ( Pritchard et al., 2000 ) and the R package Geneland v4.0.2 ( Guillot et al., 2005 ; Guillot and Santos, 2010 ). For the structure analysis, we used the admixture and the correlated allele frequency models. Th is combination most accurately assigns individuals to closely related groups ( Falush et al., 2007 ). Twenty independent runs for each K = 1–18 were performed by means of 75 000 iterations of burn-in followed by 500 000 iterations to ensure convergence of the MCMC. Th e best estimate for the number of clusters ( Δ K , as described by Evanno et al., 2005 ) was computed with STRUCTURE HARVESTER ( Earl and vonHoldt, 2011 ). Th e soft ware CLUMPP, v1.1.2 ( Jakobsson and Rosenberg, 2007 ) was used with the Greedy algorithm and 10 000 random input orders of 20 independent runs to determine the optimal alignment of clusters across individual runs for each K . Th e results from CLUMPP were imported into the program DISTRUCT v1.1 ( Rosenberg, 2004 ) for viewing. For the Geneland analysis, we ran 100 000 MCMC iterations and 20 replications for each value of K , setting the thinning to 100, allowing K to vary from 1 to 18, and with all spatial individual coordinates added. We employed the correlated frequency model, which enables even subtle genetic structuring to be detected ( Guillot, 2008 ). A burn-in period of 200 was computed in the postprocessing. We then calculated the mean PP distribution of the data for each of the 20 runs and selected only the 10 runs with the highest posterior distribution to be considered in the analysis. Th e results of computations were visually checked but were not consistent across runs (10 independent MCMC simulations) in terms of estimated number of population K , so the conclusion was based on the run giving the highest average PP, as suggested by the Geneland soft ware manual. We determined the extent of genetic diff erentiation, measured as F ST (the fi xation index), both for pairs of populations (pairwise F ST ) and for geographically based sets of populations (AMOVA; Arlequin v3.1.1; Excoffi er et al., 2005 ). Th e confi dence intervals of SEPTEMBER 2015 , VOLUME 102 • SILVA ET AL. —PHYLOGEOGRAPHY OF SONCHUS SECTION PUSTULATI • 1543 the F ST values were determined through bootstrapping (20 000 replicates), as implemented in Arlequin. Finally, to detect possible patterns of isolation by distance, we performed Mantel tests with 10 000 permutations using XLSTAT soft ware ( Addinsoft , 2010 ). To this end, we tested the correlation between population pairwise- F ST values and geographic distances for each distribution area. RESULTS cpDNA sequences — The aligned sequences of accD–psaI and 3 ′ trnL (UAA)– ndhJ intergenic spacers were 775 and 723 bp long, respectively. Of 1498 bp of aligned sequences from the two combined intergenic regions, 1475 characters were constant, and 23 characters were parsimony informative. For the accD–psaI intergenic spacer, three sets of indels were found: (1) a 7-bp insertion and a 17-bp direct repeat insertion shared by S. fragilis and S. masguindalii ; (2) a 4-bp direct repeat insertion and an 8-bp direct repeat deletion in all ingroup taxa; and (3) a 7-bp insertion shared by S. pustulatus and S. masguindalii . As for the 3 ′ trnL (UAA)– ndhJ intergenic spacer, 23-bp and 40-bp deletions were found in S. fragilis and S. pustulatus , respectively. Th e MP analysis for cpDNA data provided more than 100 000 equally parsimonious trees with a tree length (TL) of 23, a consistency index (CI) of 1.0 and a retention index (RI) of 1.0. Th e topology of the resulting MP and ML trees was identical and also consistent with the Bayesian analysis results ( Fig. 2 ). Th ese trees showed early divergence of S. pustulatus within the section Pustulati , followed by a clade containing the two other remaining species of S. masguindalii and S. fragilis (MP: 100%BS / ML: 100%BS / BI: 1.00PP). Both S. pustulatus and S. masguindalii turned out to be not monophyletic perhaps due in part to slow mutation rate of two short chloroplast noncoding regions. In the haplotype network ( Fig. 1b ), two haplotypes were found in S. pustulatus , H1 for Rifan populations (LAUF, LAUL, LAUO, TAL, and KAN) and H2 for Baetic populations (ANT, TEL, and CAB). Sonchus fragilis had only one haplotype (H3) and S. masguindalii three: H4 (TOR, BAD1, and BAD2), H5 (BOU, ALH, QUE, SFI, and one individual from BAD2), and H6 (two individuals from SFI). ITS and ETS nuclear ribosomal DNA sequences — Th e aligned ITS and ETS regions were 646 and 393 bp, respectively. Of 1039 aligned sites, 908 were constant, 2 were variable but parsimony uninformative, and 129 were variable and parsimony informative. In the Rifan populations of S. pustulatus , we detected few polymorphic sites in the ITS and ETS sequences, with only four of 393 aligned sites in the ETS (1%) and 12 of 495 aligned sites in the ITS 1 and 2 regions (2.4%). Th e populations of Baetic S. pustulatus , S. fragilis , and S. masguindalii showed no polymorphisms in the ITS and ETS sequences. Th e MP search found more than 100 000 equally parsimonious trees, with a TL of 144, a CI of 0.9583, and an RI of 0.9946 (trees not shown). Th e model test analysis found GTR+G as the best model for nucleotide evolution, and the ML search found a single best tree ( Fig. 3 ). Th e BI tree (not shown) presented the same tree topology as that of the MP and ML. All three diff erent analyses based on FIGURE 2 One maximum likelihood (ML) tree from the combined analysis of plastid accD – psaI /3 ′ trnL (UAA)– ndhJ sequences. Numbers at branches are maximum parsimony bootstrap values/ML bootstrap values/ posterior probability. 1544 • AMERICAN JOURNAL OF BOTANY combined ITS/ETS showed that the section Pustulati is monophyletic (100%BS/100%BS/1.00PP) and is quite divergent from the chosen outgroup species S. palustris . Th e monophyly of section Pustulati was further confi rmed based on the broader phylogenetic analysis of subtribe Sonchinae ( Kim et al., 2007 ; Appendix S2, see online Supplemental Data). Th ey also clearly demonstrated early divergence of S. masguindalii within the section Pustulati (100%BS/100%BS/1.00PP) and a sister relationship between S. pustulatus and S. fragilis (90%BS/93%BS/1.00PP). Within S. masguindalii , populations TOR, BAD1, and BAD2 turned out to be sisters to populations BOU and ALH (86%BS/99%BS/0.99PP), and in turn, this group of fi ve populations was sister to SFI and QUE (100%BS/100%BS/1.00PP). Th ere was no ribotype variation within S. fragilis , whereas three ribotypes were found in each of the S. masguindalii and S. pustulatus species. In S. pustulatus , two and one ribotypes were found in the Rifan and Baetic populations, respectively. Th e most recent common ancestor age of Sonchus section Pustulati was estimated to be approximately 5.48 Ma (3.171–8.766 Ma 95% CI; Fig. 4a ; Appendix S2). Within section Pustulati ( Fig. 4 ), the FIGURE 3 One maximum likelihood (ML) tree from the analysis of nuclear ITS/ETS sequences. Numbers at branches are maximum parsimony bootstrap values /ML bootstrap values/posterior probability. divergence time between S. pustulatus and S. fragilis was approximately 3.89 Ma (2.04–6.48 Ma 95% CI) and between Baetic and Rifan populations of S. pustulatus , approximately 1.18 Ma (0.42–2.41 Ma 95% CI). Th e latter time was similar to the observed divergence between S. masguindalii populations (1.39 Ma; 0.46–2.83 Ma 95% CI). AFLP marker diversity — Th e negative controls showed no amplifi cation. Th e average reproducibility of the AFLP makers, i.e., the average proportion of consistently reproduced bands over all replicates, was 98.4 ± 1.17% (mean ± SE). All individuals (167) exhibited diff erent multilocus genotypes. Th e percentage of polymorphism (%P) varied from 41.9 to 52.6%, and Frag PRI was zero in most populations ( Table 1 ). Th e DW values, either taking into account all studied plants or only those plants within populations of S. pustulatus and S. fragilis , were signifi cantly higher in the S. fragilis populations and, in general, were signifi cantly lower in the Baetic populations ( Table 1 ). Th e genetic diversity values of the Baetic populations of S. pustulatus were low (total H E = 0.109) in comparison with the remaining distribution areas (total H E = 0.138– 0.177), and with Rifan S. pustulatus (total H E = 0.167; Table 1 ). Genetic distances among populations — Two main groups within section Pustulati were recognized in the neighbor-joining phylogram based on Nei’s genetic distance ( Fig. 5a ). One of them comprised the S. pustulatus and S. fragilis species, and the other comprised S. masguindalii . Th e fi rst group was clearly divided into three clusters. Th e populations of S. fragilis were separated from the whole set of S. pustulatus (82% BS), and within this latter set, the Baetic and Rifan populations were also separated (82% BS). Th e Rifan populations of S. pustulatus from the Oued Laou valley were grouped together excluding LAUO, which was clustered apart with KAN (78%), the most isolated population located in the Oued Al-Kannar valley (20 km away, Fig. 1a ). Within the Baetic cluster, the fi rst node separated population ANT from the others (100%). A notable feature within the cluster of S. masguindalii was the presence of an independent clade (100% BS) that included the two easternmost populations, SFI and QUE. STRUCTURE and Geneland results — Th e Bayesian analysis conducted with structure showed that K = 2 was the most likely number of clusters in our data set, followed by K = 7 (online Appendix S3); K = 2 clearly separated all individuals of S. masguindalii from the rest ( Fig. 5b ), and K = 7 separated fi ve groups in which most individuals SEPTEMBER 2015 , VOLUME 102 • SILVA ET AL. —PHYLOGEOGRAPHY OF SONCHUS SECTION PUSTULATI • 1545 FIGURE 4 (a) Chronograms of the evolution of Sonchus section Pustulati based on molecular dating from the nuclear ITS/ETS sequences. Capital letters (A–D) represent the distribution areas of taxa as in Fig. 1a . Node bars represent the 95% highest posterior density intervals for the divergence time estimates of clades. (b–g) Phylogeographic reconstruction of the lineages across the western Mediterranean basin since the Middle Miocene. Capital letters on maps represent the estimated occupied region for the ancestor (Pre-D) and taxa (A–D). Black arrows indicate events of diversifi cation. Hypotheses suggested for the origin of the current disjunction of S. pustulatus : (1) recent long distance seed-dispersal events over the Alboran Sea (gray arrows in g), and (2) short-distance dispersal events across the land bridge that connected Iberia and North Africa (gray arrows in f) followed by historical geographical isolation (g). Paleogeographic information was taken from Krijgsman et al. (1999) , Martín et al. (2001) , Rosenbaum et al. (2002) , Fauquette et al. (2006) , and García-Castellanos et al. (2009) . 1546 • AMERICAN JOURNAL OF BOTANY exhibited very high assignment probabilities ( Fig. 5b ). All individuals of S. fragilis were grouped in the same cluster, and individuals of S. pustulatus were separated into two clusters: a Rifan one and a Baetic one. Th e individuals from the populations QUE and SFI were also separated from the remaining individuals of S. masguindalii , although with a higher level of genetic admixture. Th e two remaining clusters appear in many individuals of diverse origins and may result from noise produced by the sensitivity of the structure algorithm to variation in sample size ( Kalinowski, 2011 ). Th e Bayesian analysis performed with Geneland found K = 6 as the most likely number of clusters in our data set (Appendix S3). Th is analysis also separated individuals according to distribution areas, but giving three (rather than two) clusters for individuals of S. masguindalii : those from the western (populations TOR, BAD1 and BAD2), intermediate (ALH and BOU), and eastern ends (SFI and QUE). Th e splitting between western and intermediate populations could be a consequence of the over substructuring of population nuclei that Geneland generates in cases of isolation by distance ( Cushman et al., 2006 ), although it was not detected by the Mantel test for this species (see below). Partition of genetic variance among individuals and populations — Th e AMOVA results ( Table 2 ) showed a generally very high FIGURE 5 (a) Neighbor-joining analysis of 17 populations of Sonchus section Pustulati based on Nei’s genetic distance data. Numbers by nodes are bootstrap values >50. Colored circles denote main genetic groups inferred from Bayesian clustering for K = 7; and colored squares for K = 2. (b) Genetic structure of populations inferred using the model-based Bayesian algorithm implemented in the software STRUCTURE for K = 7 and K = 2. Each vertical line represents an individual that is divided into K colored segments depicting their membership in each of the K clusters. Populations are separated by vertical black bars. genetic structuring among geographical distribution areas (62% of genetic variation, F ST = 0.70). Nearly half of the genetic variation also occurred among groups both between the Rifan populations of S. fragilis and S. pustulatus (47%, F ST = 0.57) and between the Baetic and Rifan populations of S. pustulatus (40%, F ST = 0.51). All pairwise- F ST values between populations from diff erent distribution areas were high (online Appendix S4). Within each distribution area, the AMOVA results always showed a relatively high genetic structuring, especially in S. masguindalii ( F ST = 0.14– 0.25; Table 2 ). Th e most geographically close populations showed low genetic diff erentiation, e.g., Baetic S. pustulatus TEL-CAB (pairwise- F ST < 0.01) and Rifan S. pustulatus from the Oued Laou valley (0.05–0.09), S. masguindalii SFI-QUE (0.02), and S. masguindalii BAD1-BAD2 (0.09). No significant isolation by distance was detected within the S. masguindalii ( r = 0.34, P = 0.141, N = 7), Rifan S. pustulatus ( r = 0.46, P = 0.183, N = 5) or Baetic S. pustulatus ( r = 0.23, P = 0.685, N = 3) population subsets. DISCUSSION Phylogenetic congruence with current taxonomic treatment in Sonchus section Pustulati — Present phylogenetic results support current taxonomic treatment of the three species that confi rm section Pustulati Boulos of genus Sonchus : S. masguindalii , S. fragilis , and S. pustulatus . Th ese were clearly separated with each of the used molecular markers. Th e high congruence identifi ed between the ITS/ETS and AFLP analyses supports the monophyletic origin of Sonchus section Pustulati ( Kim et al., 2007 , 2008 ). Furthermore, the topological incongruence detected between phylogenies based on maternally inherited chloroplast sequences (cpDNA) and biparentally inherited nuclear sequences (ITS/ETS and AFLP) led to the hypothesis that S. pustulatus originated by interspecifi c hybridization from S. fragilis , as a likely paternal contributor, and a still unknown maternal donor from the root lineage of sections Sonchus / Asperi ( Kim et al., 2008 ). However, alternative hypotheses, such as incomplete lineage sorting of ancient polymorphisms or S. pustulatus simply being the carrier of a captured plastid lineage to explain the incongruence between nuclear and chloroplast data, cannot be completely ruled out. Additional study based on several independent single- or low-copy nuclear loci may sort out those confl icting patterns between chloroplast and nuclear DNA. A pre-Mediterranean clade: Late Miocene-Pliocene origin of the species — In this study, we also found that Sonchus section Pustulati