Pliocene-Pleistocene ecological niche evolution shapes the phylogeography of a Mediterranean plant group

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This is the peer reviewed version of the following article: Benítez Benítez, Carmen, Escudero Lirio, Marcial, Rodríguez Sánchez, Francisco, Martín-Bravo, S.; Jiménez-Mejías, P. (2018). Pliocene-Pleistocene ecological niche evolution shapes the phylogeography of a Mediterranean plant group. .Molecular Ecology. 27 (7), 1696-1713. which has been published in final form at https://doi.org/10.1111/mec.14567. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Use of Self- Archived Versions. This article may not be enhanced, enriched or otherwise transformed into a derivative work, without express permission from Wiley or by statutory rights under applicable legislation. Copyright notices must not be removed, obscured or modified. The article must be linked to Wiley’s version of record on Wiley Online Library and any embedding, framing or otherwise making available the article or pages thereof by third parties from platforms, services and websites other than Wiley Online Library must be prohibited." Page 1 of 48 Molecular Ecology 1 Pliocene-Pleistocene ecological niche evolution shapes the phylogeography of a 2 Mediterranean plant group 3 Benítez-Benítez, C.1*, Escudero, M.2, Rodríguez-Sánchez, F.3, Martín-Bravo, S.1† & 4 Jiménez-Mejías, P.4† 5 1Área de Botánica, Departamento de Biología Molecular e Ingeniería Bioquímica, 6 Universidad Pablo de Olavide, carretera de Utrera km 1 s.n., 41013 Seville, Spain. 7 2Department of Plant Biology and Ecology, University of Seville, Spain. 8 3Departamento de Ecología Integrativa, Estación Biológica de Doñana, Consejo 9 Superior de Investigaciones Científicas, Avda. Américo Vespucio 26, 41092 Sevilla, 10 Spain. 11 4Department of Botany, Smithsonian Institution, Washington, D.C., DC, United States. 12 *Author for correspondence: cbenben[email protected]s 13 †Both authors contributed equally to the senior authorship 14 Running head: Bioclimatic niche evolution in Carex Page 1 of 48 Molecular Ecology 1 15 Abstract 16 Estimating species ability to adapt to environmental changes is crucial to understand 17 their past and future response to climate change. The Mediterranean basin has 18 experienced remarkable climatic changes since the Miocene, which have greatly 19 influenced the evolution of the Mediterranean flora. Here we examine the evolutionary 20 history and biogeographic patterns of two sedge sister species (Carex, Cyperaceae) 21 restricted to the western Mediterranean basin, but with Pliocene fossil record in central 22 Europe. In particular, we estimated the evolution of climatic niches through time and its 23 influence in lineage differentiation. 24 We carried out a dated phylogenetic-phylogeographic study based on seven DNA 25 regions (nDNA and ptDNA) and fingerprinting data (AFLPs), and modelled ecological 26 niches and species distributions for the Pliocene, Pleistocene, and present. Phylogenetic 27 and divergence time analyses revealed that both species form a monophyletic lineage 28 originated in the late Pliocene-early Pleistocene. We detected clear genetic 29 differentiation between both species with distinct genetic clusters in disjunct areas, 30 indicating the predominant role of geographic barriers limiting gene flow. We found a 31 remarkable shift in the climatic requirements between Pliocene and extant populations, 32 although the niche seems to have been relatively conserved since the Pleistocene split of 33 both species. This study highlights how an integrative approach combining different 34 data sources and analyses, including fossils, allows solid and robust inferences about the 35 evolutionary history of a plant group since the Pliocene. 36 Keywords: allopatric speciation, Carex sect. Phacocystis, genetic structure, niche 37 conservatism, niche differentiation, species distribution modelling 38 Molecular Ecology Page 2 of 48 2 39 Introduction 40 The interplay of ecological and evolutionary processes shapes the distribution of 41 species and the dynamics of speciation, dispersal, adaptation and extinction, which 42 results from geological and climatic changes (Webb, Ackerly, McPeek & Donoghue, 43 2002; Lavergne, Mouquet, Thuiller & Ronce, 2010). Therefore, to explain 44 biogeographic and speciation patterns, it is crucial to understand how ecological 45 requirements evolve through time. An increasing number of studies on the evolution of 46 species climatic niches shows that, although negligible niche change may happen since 47 the split of sister lineages (phylogenetic niche conservatism, PNC; Harvey & Pagel, 48 1991; Ackerly, 2003; Wiens & Graham, 2005), there are also examples of ecological 49 innovation (niche shift) since the species’ divergence (Ackerly, 2003). 50 One of the most recognized signatures of niche conservatism is the similarity of 51 ecological preferences among closely related species (Wiens & Graham, 2005; Losos, 52 2008), hypothesis previously identified by Darwin (1859) as a consequence of their 53 relatively recent common ancestry. In addition to shared ancestry, Grafen (1989) 54 proposed that niche similarity is also due to shared environmental constraints and 55 biogeographic history (Wiens & Graham, 2005; Losos, 2008; Crisp & Cook, 2012). The 56 niche conservatism idea has played an important role particularly in relation to climate 57 change impacts (Thuiller, 2003), as climatic tolerance limits the ranges of species 58 (Wiens & Graham, 2005), even at very large geographical and evolutionary scales 59 (Crisp et al., 2009; Mairal, Sanmartín & Pellissier, 2017). Holt (1996) related the 60 absence of evolutionary change with the difficulty for adaptive traits to evolve and 61 stabilize in novel ecological conditions. However, niches are rarely considered 62 completely conserved due to processes of local adaptation (Thompson, 2005, Wiens & 63 Graham, 2005). In particular, the presence of barriers that limit expansion (e.g. Page 3 of 48 Molecular Ecology 3 64 mountain ranges; Taberlet, Fumagalli, Wust-Saucy & Cosson, 1998) may play a key 65 role for niche evolution in allopatric speciation (Wiens, 2004a), since the adaptation to 66 ecological conditions could be associated to the constraints imposed in different 67 microhabitat preferences (Thompson, 2005). Moreover, recent studies have 68 demonstrated that niche differences could not always be related to phylogenetic 69 distances, since in some cases sister species have not conserved the ancestral niche 70 (Revell, Harmon & Collar, 2008; Münkemüller, Boucher, Thuiller & Lavergne, 2015). 71 The assumption of niche conservatism has recently been further challenged based on 72 examples of rapid niche shifts during species diversification (Evans, Smith, Flynn & 73 Donogue, 2009). Only niche shift allows species to colonize and adapt to new habitats 74 ecologically dissimilar from those on their ancestral areas (Spalink et al., 2016), or 75 persist in the same area after an environmental change. There is also evidence of niche 76 shifts derived from biological invasions (Broennimann et al., 2007; Stiels et al., 2015). 77 Large-scale shifts in plant species across biomes have also been revealed (e.g. Koecke 78 et al., 2013; Jara-Arancio et al., 2014). The divergence of ecologically relevant traits 79 may have initially evolved sympatrically in the past, however range shifts resulting in 80 small range overlap can promote high trait divergence (Anacker, 2004) associated to the 81 evolutionary specialization of species to a particular climatic regime (Thompson, 2005). 82 The Mediterranean basin is considered one of the Earth´s hotspots of plant 83 biodiversity regarding species richness and endemism (Médail & Quézel, 1997; Myers 84 et al., 2000). Plant distribution patterns around the Mediterranean basin have been of 85 particular biogeographic interest due to the influence of complex historical geological- 86 climatic processes and extant heterogeneous environmental conditions in this region 87 (Médail & Quézel, 1997; Médail & Diadema, 2009; Thompson, 2005). The origin of the 88 Mediterranean-type climate and its flora is particularly interesting from an evolutionary Molecular Ecology Page 4 of 48 4 89 perspective. The Messinian Salinity Crisis (5.96-5.3 million years ago (Ma)) was a 90 significant period of sub-tropical flora extinction and adaptation to aridity, triggering 91 diversification of Mediterranean lineages (Rodríguez-Sánchez et al., 2008; Fiz-Palacios 92 & Valcárcel, 2013). The progressive global cooling and local aridification already 93 initiated during the late Miocene (Suc et al., 1995; Milne & Abbott, 2002; Thompson, 94 2005) culminated with the onset of the Mediterranean climatic regime during the late 95 Pliocene (3.4-2.8 Ma). This climatic progression modified the composition and structure 96 of the Mediterranean forests, which lost elements from the pre-existing subtropical 97 lauroid flora and began to resemble contemporary vegetation (Thompson, 2005). 98 Remnants of such subtropical vegetation persisted in the form of relict elements mostly 99 at the eastern or western parts of the basin (Médail & Diadema, 2009; Fernández- 100 Mazuecos et al., 2016; Míguez et al., 2017), as well as in Mediterranean islands and 101 Macaronesia (Médail & Quézel, 1997; Médail & Diadema, 2009). The later Quaternary 102 climatic oscillations (starting c. 2.5 Ma), characterised by the alternation of colder 103 (glacial) and warmer (interglacial) periods, strongly affected the distribution and genetic 104 structure of species due to the recurrent range shifts (e.g. Mansion et al., 2008; Gentili et 105 al., 2015). In this case, the Mediterranean basin and Macaronesia served as glacial 106 refugial areas (Vargas, 2007; Médail & Diadema, 2009; Gentili et al., 2015; Mairal, 107 Sanmartín & Pellissier, 2017), where the long-term persistence of isolated sets of 108 populations from ancestral species frequently resulted in the formation of new allopatric 109 lineages (Wiens, 2004b). The role of niche evolution in the biogeographic and 110 diversification patterns of the Mediterranean basin is crucial but still poorly understood. 111 Our study group is a complex of two closely related sister species, Carex reuteriana 112 Boiss. and C. panormitana Guss. (Cyperaceae), allopatrically distributed in the western 113 Mediterranean basin (Fig. 1), where they inhabit creeks and river shores at low and Page 5 of 48 Molecular Ecology 5 114 medium altitudes (20-1900 m; Luceño & Jiménez-Mejías, 2008). Carex reuteriana is 115 endemic to the Iberian Peninsula and northwestern Africa (Luceño & Jiménez-Mejías, 116 2008), and C. panormitana has a restricted distribution in Sardinia and in a few isolated 117 populations in Sicily and Tunisia (Pignatti, 1982; Urbani, Gianguzzi & Ilardi, 1995; 118 Jiménez-Mejías, Martín-Bravo, Amini-Rad & Luceño, 2013). A close phylogenetic 119 relationship between both disjunct taxa has been recently found (Global Carex Group, 120 2016). Currently, two subspecies are accepted for C. reuteriana (Jiménez-Mejías et al., 121 2011): (1) C. reuteriana ssp. reuteriana, which shows a patchy and disjunct distribution 122 in the central and western Iberian Peninsula; and (2) C. reuteriana ssp. mauritanica 123 (Boiss. & Reut.) Jim.-Mejías & Luceño, in the southwestern Iberian Peninsula and 124 northwestern Africa. Geographic isolation has been previously proposed as the main 125 factor driving the diversification of C. reuteriana (Jiménez-Mejías et al., 2011). In 126 contrast, the phylogeography of C. panormitana has not been previously studied. 127 Interestingly, the recent study of fossil evidences has shown that the putative ancestor of 128 C. reuteriana-C. panormitana was widely distributed throughout central Europe and 129 western Asia during the Pliocene (Jiménez-Mejías & Martinetto, 2013; Jiménez-Mejías 130 et al., 2016) indicating the former presence of species of this group in areas where it is 131 now absent. 132 Phylogenetic niche conservatism is the most common pattern in phylogenetic studies 133 of species diversification (Ackerly, 2003; Wiens, 2004; Crisp et al., 2009; Wiens et al., 134 2010). However, several case studies in the Mediterranean basin suggested little niche 135 conservatism (Donoghue, 2008; Rundel et al., 2016). We investigate the evolutionary 136 history of the C. reuteriana-C. panormitana complex using an integrated approach 137 comprising phylogenetic, phylogeographic and dating analyses, combined with species 138 distribution modelling (SDM). In the present study we aim to evaluate the amount of Molecular Ecology Page 6 of 48 6 139 niche conservatism that occurred along the diversification of the C. reuteriana-C. 140 panormitana group in the context of climate changes in the Mediterranean through the 141 Plio-Pleistocene. In addition, we assess the phylogeographic structure to infer the main 142 drivers of differentiation among both species and subspecies and their populations. 143 Materials and Methods 144 Estimation of divergence times 145 DNA amplification and matrix construction 146 Two nuclear (ITS and ETS-1f) and three plastid (psbA-trnH, rpl32-trnLUAG and 147 ycf6-psbM) regions were selected to assess the phylogenetic relationships within the 148 Eurasian group of Carex sect. Phacocystis according to previous results by Jiménez- 149 Mejías (2011). Fifteen samples of 12 species were obtained from silica-dried field- 150 collected material and herbarium samples (Appendix S1). DNA extraction was 151 performed as in Jiménez-Mejías (2011). The primers and amplification procedures for 152 the nuclear and plastid regions followed those in Dragon & Barrington (2009), Shaw et 153 al. (2005) and Shaw, Lickey, Schilling & Small (2007), respectively. Sequence 154 chromatograms were edited using Geneious v. 6.1.8 (Biomatters Ltd., Auckland, New 155 Zealand) and automatically aligned with MUSCLE (Edgar, 2004). Newly obtained 156 edited sequences were included in a matrix with a representative sampling of the genus 157 Carex as compiled by Míguez et al. (2017) and in GCG (2016). To this matrix, we 158 added ITS, ETS1f, psbA-trnH and rpl32-trnLUAG sequences from GenBank (Appendix 159 S1). A final combined matrix consisting of 134 concatenated sequences was constructed 160 with an aligned length of 5211 sites. Two samples of C. reuteriana (one sample per C. 161 reuteriana subspecies) and two samples of C. panormitana (one from Sardinia and 162 another from Sicily) were included. 163 Dating analysis 13 Page 13 of 48 Molecular Ecology 288 SNAPP analyses of AFLP data and the date of origin of the complex that resulting from 289 BEAST results (see Multispecies coalescent AFLP analysis in Results section). The 290 inheritance scalar considered more appropriate was 1 for autosomal nDNA and indels 291 and 0.5 for ptDNA to reflect the expected effective population size of an inherited locus 292 (Qiu et al., 2009; Pettengill & Moeller, 2012). A range of mutation rates was calculated 293 for each locus (excluding indels) based on the 95% highest posterior density interval for 294 the crown node age of the complex. This range was estimated using DnaSP v. 5.1 295 (Librado & Rozas, 2009) to calculate the average number of nucleotide differences 296 between populations per each locus. These ranges were included in the model to 297 estimate parameters related to species. Population model was constrained to allow only 298 migration between current and ancestral sister lineages. Multiple trial runs were 299 performed in IMA2 to establish an adequate set of priors in the final analyses. We 300 performed two different runs with the selected combination of a priori maximum values 301 for the parameters migration (m = 200), population size (q = 50) and time (t = 3) using 302 150 chains and five millions of generations, with a burn-in of one million of 303 generations. In each run 5,000 trees were stored, adding up 10,000 genealogies which 304 were used to calculate the final parameters in demographic scales using the mean 305 mutation rates of the three DNA regions. We used a generation time of two years 306 (which is the minimum for the vast majority of Carex species; Ball & Reznicek, 2002), 307 as well as five and ten years. 308 Species distribution modelling 309 Information on the geographic distribution of Carex reuteriana-C. panormitana and 310 their putative fossil ancestor (C. panormitana-type; Table S1) were used to assess the 311 present ecological niche and its evolution from the Pliocene to present times. The 312 current distribution range of the complex was collected mainly from GBIF 14 Molecular Ecology Page 14 of 48 313 (https://www.gbif.org/) and corrected and complemented with our own field sampling 314 (Appendix S4), bibliographic references (Jiménez-Mejías, Escudero, Chaparro & 315 Luceño, 2007; Gianguzzi, Cusimano, Ilardi & Salvatore, 2013; Urbani, Calvia & 316 Pisanu, 2013) and herbarium specimens (COFC, SEV, and SS herbaria; Thiers, 2015). 317 Unreliable points due to incorrect georeferencing (e.g. populations falling in the sea or 318 outside the known distribution range) or taxonomic issues (e.g. known 319 misidentifications) were removed. The final dataset for extant species, accurately and 320 comprehensively representing their current distribution, was composed of 462 point 321 localities, 433 belonging to C. reuteriana (315 for ssp. reuteriana and 118 for ssp. 322 mauritanica) and 29 belonging to C. panormitana (which included all known 323 populations: six close subpopulations of the only population from Sicily, two 324 populations from Tunisia, and 21 located in Sardinia). The final dataset for the fossil 325 ancestor was composed of 17 Pliocene fossil records identified as C. panormitana-type 326 (Fig. 1; Jiménez-Mejías et al., 2016), which were used to validate the Pliocenic 327 palaeodistribution projection. 328 We retrieved 19 bioclimatic variables (resolution of 2.5 minutes) from WorldClim 329 (http://www.worldclim.org/; Hijmans et al., 2005) for current conditions and projections 330 to the Late Quaternary; and from ecoClimate (resolution of 0.5º; http://ecoclimate.org/; 331 Lima-Ribeiro et al., 2015), for projections to the Pliocene. We also used soil pH (15 cm 332 depth) from ISRIC (http://www.isric.org/) as environmental predictor, since C. 333 reuteriana ssp. reuteriana is acidophilous, thus not indifferent to the substrate. We 334 selected environmental variables based on ecological knowledge of our species and 335 statistical criteria (Mod et al., 2016). Thus, to avoid using highly correlated variables for 336 the modelling, we calculated the absolute value of the correlation matrix converting it 337 into a matrix of distances in the form of a dendrogram. We only selected one variable 15 Page 15 of 48 Molecular Ecology 338 per clade when the branch was below 0.5. Furthermore, we also used the variance 339 inflation factor (VIF), a measure of how much collinearity increases variance in a 340 model, using HH package (Heiberger, 2017) in R (R Core Development Team, 2017). 341 We selected variables with VIF < 5 which had sound biological interpretation for our 342 species. In the end, in addition to soil pH, we used four bioclimatic variables, two 343 representing temperatures (bio2: mean diurnal range, and bio4: temperature seasonality) 344 and two representing precipitation (bio15: precipitation seasonality, and bio16: 345 precipitation of wettest quarter). 346 The maximum entropy algorithm, as implemented in MaxEnt v. 3.3.3 (Elith et al., 347 2011), was used to evaluate the potential distribution of our species complex for: (1) 348 present times; (2) the Late Quaternary, including (2a) the Last Glacial Maximum (LGM, 349 c. 21ka) and (2b) the Last interglacial (LIG, 120-140 ka)); and (3) the middle Pliocene 350 (3 Ma). For model calibration, C. reuteriana was subdivided into two subspecies 351 (Jiménez-Mejías et al., 2011), to account for the retrieved genetic structure (see results). 352 On the contrary, C. panormitana was considered as a single group. For the projection of 353 the suitable range at the middle Pliocene we used a model based on all extant 354 occurrences since in that period the species most likely did not display phylogeographic 355 structure yet, according to estimated divergence times (see results). We used the same 356 background area for all present and past projections (Fig. 4, Fig. S3). Finally, for 357 comparison we also fitted models to the Pliocene fossil records using all 17 occurrences 358 of C. panormitana-type. 359 For each model we ran 100 replicates, using a random partition of 20% of the 360 occurrences data set to test the model. We used the software default parameters, but 361 excluding threshold and hinge features which may lead to overfitting (Merow, Smith & 362 Silander, 2013; Phillips et al., 2017). We evaluated the calibrated models by comparing 16 Molecular Ecology Page 16 of 48 363 response curves, area under the curve (AUC) scores and jackknife tests to select those 364 uncorrelated variables with the highest relative contributions to the model. The retained 365 variables were subsequently used to repeat the modelling analyses for present, Late 366 Pleistocene and Pliocene. 367 In addition, the environmental niche occupied by extant taxa and Pliocene fossil 368 records was visualised through a Principal Components Analysis (PCA) of the retained 369 variables using the “prcomp” function in R. We also used the function “ggbiplot” in 370 ggplot2 package (Wickham, 2009) for mapping individual climatic variables into the 371 multivariate climate space. We compared the fitted response curves for the retained 372 variables between extant taxa and Pliocene fossils. We calculated values of Schoener’s 373 index (D) statistic for niche overlap (Schoener, 1968; Warren, Glor & Turelly, 2008). 374 Niche overlap between genetic groups (see results) was evaluated using the niche 375 similarity test which assesses whether the niches of both groups are more similar than 376 expected by chance (Warren, Glor & Turelly, 2008). We also calculated the niche 377 equivalency that tests whether the niches occupied by two taxa are identical (Warren, 378 Glor & Turelly, 2008; Broennimann et al., 2012). Schoener’s D and niche similarity and 379 equivalency tests were calculated in ecospat package (Di Cola et al., 2017) 380 Results 381 Estimation of divergence times from DNA sequences 382 Diversification of the core Phacocystis was estimated to have occurred during the 383 Late Miocene (mean 7.97 Ma, 95% HPD: 4.93-10.37 Ma; Table 2a; see Fig. S1). The 384 split of the two main lineages (lineages 1 and 2, see Fig. 3a) could have taken place 385 around 6.2 Ma (95% HPD: 4.29-9.02 Ma; Fig. 3a), also coinciding with Late Miocene. 386 The stem node of C. panormitana - C. reuteriana complex was dated at around 4.79 Ma 17 Page 17 of 48 Molecular Ecology 387 (95% HPD: 3.8-6.43 Ma; Fig. 3a, Table 2a) during the late Miocene-early Pliocene, 388 while the crown node fell in the Upper-Middle Pleistocene (mean 0.95 Ma; 95% HPD: 389 0.19-2.48 Ma; Fig. 3a, Table 2a). 390 391 AFLP population genetic analyses 392 The final AFLP matrix was composed of 182 polymorphic loci for the 130 393 sampled individuals from 18 populations (Fig. 1). The PCoA of the whole AFLP dataset 394 (Fig. 2c) revealed three well-defined clusters: 1) C. reuteriana, comprising individuals 395 of subspecies reuteriana and mauritanica segregated at different sides of the cluster 396 (Fig. 2d); 2) Sardinian populations of C. panormitana; and 3) Sicilian and Tunisian 397 populations of C. panormitana, with a slight separation between the samples coming 398 from each of the two areas (Fig. 2e). The three first factors of the PCoA including all 399 samples and taxa accounted for 23.69%, 10.94% and 8.48% of the variation, 400 respectively. The PCoA scatter plot for C. reuteriana ssp. reuteriana did not show 401 geographical structure, whereas the PCoA performed for C. reuteriana ssp. mauritanica 402 revealed a clear subdivision into three subclusters with geographic congruence 403 (Jiménez-Mejías et al., 2011): (1) Iberian populations from Sierra Morena (north of 404 Guadalquivir river valley), (2) Iberian populations from Betic ranges (south of 405 Guadalquivir river valley), and (3) the northern African (Moroccan) populations. No 406 clear geographical pattern was detected among the samples of Sardinian populations. 407 The STRUCTURE analysis of the whole dataset using ∆K method yielded four 408 optimal genetic clusters (K=4) with taxonomic/geographic correspondence (Fig. 2a). 409 The first cluster comprised Sardinian populations of C. panormitana, the second one 410 included Sicilian-Tunisian populations of C. panormitana, the third comprised all 411 populations of the C. reuteriana ssp. mauritanica, and the last cluster incorporated the 18 Molecular Ecology Page 18 of 48 412 populations of the subspecies reuteriana. In this STRUCTURE analysis the admixture 413 found between groups was incidental (Fig. 2a). On the other hand, the results from Ln 414 Pr (X/K) tool displayed six optimal genetic clusters (K=6; Fig. 2b). Subspecies from C. 415 reuteriana were subdivided in two clusters each one: (1) subspecies mauritanica into 416 north and south Guadalquivir valley, (2) and subspecies reuteriana into one cluster 417 constituted by two populations from CW Spain and the other one with the other 418 populations (Fig. 2b). 419 SNAPP analyses yielded a well-resolved topology in which both the two species (C. 420 reuteriana, C. panormitana) and the two subspecies of C. reuteriana were well 421 supported (>0.95 PP; Fig. 3). The two main lineages within C. panormitana (Sardinia, 422 Sicily-Tunisia) received slightly lower support (>0.85 PP; Fig. 3b). The complex C. 423 reuteriana-C. panormitana first differentiated c. 0.95 Ma at the Upper-Middle 424 Pleistocene (95%HPD: 0.663-1.280 Ma; Fig. 3b). The split between the two C. 425 reuteriana subspecies took place around 0.576 Ma ago (95% HPD: 0.428-0.734 Ma; 426 Fig. 3b, Table 2B). Subspecies mauritanica diversified in two lineages (0.436 Ma; 95% 427 HPD: 0.313-0.562 Ma; Fig. 3b, Table 2b): (1) Iberian populations from Sierra Morena, 428 and (2) Iberian populations from the south of Guadalquivir river valley plus the 429 Moroccan populations. The crown node of ssp. reuteriana was dated back to 0.313 Ma 430 (95% HPD: 0.205-0.436 Ma; Fig. 3b, Table 2b). Divergence of the Sicilian-Tunisian 431 lineage from the Sardinian lineage of C. panormitana was estimated at 0.785 Ma (86% 432 HPD: 0.471-1.045 Ma; Fig. 3b, Table 2b).The split between Sicilian and Tunisian 433 populations was dated at 0.529 Ma ago (95% HPD: 0.313-0.755 Ma; Fig. 3b, Table 2B). 434 Finally, Sardinian populations were those that diversified most recently (0.144 Ma; 95% 435 HPD: 0.085-0.212 Ma; Fig. 3b, Table 2b). Effective population sizes from SNAPP 436 analysis were higher for C. reuteriana than for C. panormitana; in turn, Sardinian 19 Page 19 of 48 Molecular Ecology 437 populations belonging to the latter species were estimated to have a smaller effective 438 size than Sicilian-Tunisian populations. However, both subspecies from C. reuteriana 439 showed a similar effective population size. In addition, the common ancestor of the 440 lineages had larger population size than subsequently divergent groups (see Table S3). 441 Since we are aware of the limitations derived from the SNAPP analysis, a detailed 442 section about methodological concerns has been included (see Appendix 7). 443 444 nDNA-ptDNA haplotype analyses 445 Statistical parsimony analysis of plastid data set revealed three haplotypes (H1-H3; 446 Fig. S2a, Table 1). The most common haplotype (H3) was shared between C. reuteriana 447 ssp. mauritanica and C. panormitana (Fig. S2a, Table 1). The other haplotypes were 448 exclusive of C. reuteriana ssp. reuteriana (H2) and C. reuteriana ssp. mauritanica (H1) 449 respectively Fig. S2a, Table 1). The nDNA network displayed ten haplotypes (H1-H10; 450 Fig. S2b, Table 1). No haplotypes were shared between C. panormitana and C. 451 reuteriana. Carex panormitana displayed two haplotypes, the most common (H9) 452 present in all sampled populations, while H10 was found in a few individuals from 453 Sicily (Fig. S2b, Table 1). Carex reuteriana ssp. reuteriana displayed six haplotypes 454 (H1, H3 and H5-H8) and C. reuteriana ssp. mauritanica five haplotypes (H1-H5), the 455 two subspecies sharing three haplotypes. 456 Isolation-with-migration model 457 Mean divergence times estimated for the two ancestors of C. reuteriana-C. 458 panormitana with the IMA2 analysis (Fig. S2c) were more recent than those obtained 459 from the SNAPP analysis based on AFLPs (see below; Fig. 3b). However, taking into 460 account 95% confidence intervals, ages from IMA2 (Fig. S2b) were mostly congruent 461 with those from SNAPP (Fig. 3b, see Appendix S5). Additional results about estimated 20 Molecular Ecology Page 20 of 48 462 values of effective population size and detailed methodological caveats about IMA2 463 have been included in Appendix S6 and S7, respectively. 464 Species distribution modelling 465 The four bioclimatic variables that most contributed to the models, in addition to soil 466 pH, were: (bio2) Mean Diurnal Temperature Range (mean of monthly differences 467 between maximum and minimum temperature), (bio4) Temperature Seasonality 468 (standard deviation of monthly mean temperature*100), (bio15) Precipitation 469 Seasonality (coefficient of variation of monthly precipitation), and (bio16) Precipitation 470 of the Wettest Quarter (see Table S4). AUC values were above 0.9 for all models, which 471 indicate good predictive ability (Swets, 1988). The potential distribution inferred by 472 MaxEnt at present for the different taxa partitions considered within the C. reuteriana- 473 C. panormitana complex (Fig. 4) closely matched the current populations distribution in 474 each case (Fig.1). Projections to the LGM (21ka; Fig. S3a) and LIG (120-140 ka; Fig. 475 S3b) revealed potential distributions similar to the current ones, although with a more 476 restricted distribution for the LIG. Our two independent approaches for inferring middle 477 Pliocene (3 Ma) distributions (i.e., projecting models calibrated with either current 478 occurrences or Pliocene fossils) revealed similar results. In both cases, the potential 479 distribution expanded into higher latitudes (Central Europe) (Fig. 4d-e). 480 Pliocene populations occupied significantly less seasonal climates (both for temperature 481 and precipitation) than extant populations (Fig. 5a). On the other hand, all three extant 482 taxa overlapped strongly in their climatic spaces (Fig. 5b): only C. panormitana differed 483 in occupying climates with lower temperature oscillations and intermediate precipitation 484 seasonality, than C. reuteriana subspecies (Fig. 5b). The PCA revealed that Pliocene C. 485 panormitana-type fossils appeared well differentiated in environmental space from 486 populations of extant species (Fig. 6). With respect to extant taxa, the PCA revealed an Molecular Ecology Page 20 of 48 20 487 overlap between C. reuteriana and C. panormitana (Fig. 5b, Fig. 6). However, all niche 488 similarity tests failed to show niche similarity higher than expected by chance (Table 3). 489 Pairwise statistical comparisons of ecological niches suggested that they are not 490 significantly more similar than expected by chance (Table 3). However, values were far 491 closer to significance (especially for the comparison between ssp. reuteriana and ssp. 492 mauritanica) than for the comparison between fossils and extant taxa. 493 Discussion 494 A vicariant genetic structure shaped by Pleistocene glaciations 495 We have assessed the evolutionary history, ecological niche shifts, and 496 biogeographic relationships in the C. reuteriana-C. panormitana complex constituted 497 by two sister species disjunctly distributed in the western Mediterranean. Both species 498 form a strongly supported monophyletic group (GCG, 2016) which seems to have 499 diverged from its closest relatives during the Pliocene (Fig.3a, Table 2a). The estimated 500 times of diversification within this group were straitened within the Pleistocene (Fig.6, 501 Table 2A), suggesting a role of the Pleistocene glacial-interglacial cycles, as has been 502 inferred for other western Mediterranean plant groups (Terrab et al., 2008; González- 503 Martínez et al., 2010; Tremetsberger et al., 2016; among many others). The genetic 504 structure and estimated divergence ages points to range contraction followed by 505 isolation in different western Mediterranean glacial refugia (Iberian and Italian 506 peninsulas, northwestern Africa, and Tyrrhenian islands; Médail & Diadema, 2009) as 507 the most likely scenario for the speciation and population differentiation within the C. 508 reuteriana-C. panormitana complex (Figs. 2, 3b). The harsher climatic conditions 509 during the glacial periods seem to have resulted in the extinction of the group’s ancestor 510 from Central and Northwestern Europe (North Italy, Germany, Poland, Netherlands), 27 Page 27 of 48 Molecular Ecology 634 3- Ball, P.W., & Reznicek, A.A. (2002). Magnoliophyta: Commelinidae: 635 Cyperaceae. In: Flora of North America, North of Mexico. Carex Linnaeus. 636 Oxford, UK: Oxford University Press, New York, 23, 254–273. 637 4- Bonin, A., Bellemain, A.E., Bronken Eidesen, P., Pompanon, F., Brochmann, 638 C., & Taberlet, P. (2004). How to track and assess genotyping errors in 639 population genetic studies. Molecular Ecology, 13, 3261-3273. 640 5- Boucher-Lalonde, V., Morin, A., & Currie, D.J. (2016). Can the richness- 641 climate relationship be explained by systematic variations in how individual 642 species’ ranges relate to climate? Global Ecology and Biogeography, 25, 527- 643 539. 644 6- Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.H., Xie, D., 645 Suchard, M.A., Rambaut, A., & Drummond, A.J. (2014). BEAST 2: A 646 software platform for Bayesian evolutionary analysis. PloS Computational 647 Biology, 10, e1003537. https://doi.org/10.1371/journal.pcbi.1003537.g002 648 7- Broennimann, O., Fitzpatrick, M.C., Pearman, P.B., Petitpierre, B., Pellissier, 649 L., Yoccoz, N.G., & Guisan, A. (2012). Measuring ecological niche overlap 650 from occurrence and spatial environmental data. Global Ecology and 651 Biogeography, 21, 481–497. 652 8- Broennimann, O., Treier, U.A., Muller-Scharer, H., Thuiller, W., Peterson, 653 A.T., & Guisan, A. (2007). Evidence of climatic niche shift during biological 654 invasion. Ecology Letters, 10, 701-709. 655 9- Brummitt R.K. (2001). World geographical scheme for recording plant 656 distributions, 2nd edition. Carnegie Mellon University, Pittsburgh, 657 Pennsylvania, USA: Hunt Institute for Botanical Documentation 658 10- Bryant, D., Bouckaert, R., Felsenstein, J., Rosenberg, N.A. & Roy Choudhury, 659 A. (2012). Inferring species trees directly from biallelic genetic markers: by 660 passing gene trees in a full coalescent analysis. Molecular Biology and 661 Evolution, 29, 1917-1932. 662 11- Cavender-Bares, J., & Pahlich, A. (2009). Molecular, morphological, and 663 ecological niche differentiation of sympatric sister oak species, Quercus 664 virginiana and Q. geminata (Fagaceae). American Journal of Botany, 96, 665 1690-1702. 666 12- Clement, M., Posada, D., & Crandall, K. (2000). TCS: a computer program to 667 estimate gene genealogies. Molecular Ecology, 9, 1657-1660. 668 13- Crisp, M.D., Arroyo, M.T.K., Cook, L.G., Gandolfo, M.A., Jordan, G.J., 669 McGlone, M.S., Weston, P.H., Westoby, M., Wilf, P., & Linder, H.P. (2009). 670 Phylogenetic biome conservatism on a global scale. Letters, 458, 754-756. 671 14- Crisp, M.D., & Cook, L.G. (2012). Phylogenetic niche conservatism: what are 672 the underlying evolutionary and ecological causes? New Phytologist, 196, 681- 673 694. 674 15- Darriba, D., Taboada, G.L., Doallo, R., & Posada, D. (2012). JModelTest 2: 675 More models, new heuristics and parallel computing. Nature Methods, 9, 772. 676 16- Darwin, C. (1859). On the Origin of Species. Ed. J. Murray. London. 677 17- De Castro, O., Véla, E., Vendramin, G.G., Gargiulo, R., & Caputo, P. (2015). 678 Genetic structure in the Genista ephedroides complex (Fabaceae) and 679 implications for its present distribution. Botanical Journal of the Linnean Society, 680 177, 607-618. 681 18- Di Cola, V., Broennimann, O., Petitpierre, B., Breiner, F.T., D’ Amen, M., 682 Randin, C., Engler, R., Pottier, J., Pio, D., Dubuis et al. (2017). Ecospat: an R 28 Molecular Ecology Page 28 of 48 683 package to support spatial analyses and modeling of species niche and 684 distributions. Ecography, 40, 774-787. 685 19- Donoghue, M.J. (2008). A phylogenetic perspective on the distribution of plant 686 diversity. Proceedings of the National Academy of Sciences, 105, 11549-11555. 687 20- Dragon, J., & Barrington, D.S. (2009). Systematics of the Carex aquatilis and 688 C. lenticularis lineages: geographically and ecologically divergent sister 689 clades of Carex section Phacocystis (Cyperaceae). American Journal of 690 Botany, 96, 1896-1906. 691 21- Drummond, A.J.,& Rambaut, A. (2007). BEAST: Bayesian evolutionary 692 analysis by sampling trees. BMC evolutionary biology, 7, 214. 693 22- Drummond, A.J., Suchard, M.A., Xie, D., & Rambaut, A. (2012). Bayesian 694 Phylogenetics with BEAUti and the BEAST 1.7. Molecular Biology & 695 Evolution, 29, 1969–1973. 696 23- Earl, D.A., & VonHoldt, B.M. (2012). STRUCTURE HARVESTER: a 697 website and program for visualizing STRUCTURE output and implementing 698 the Evanno method. Conservation Genetics Resources, 4, 359-361. 699 24- Edgar, R.C. (2004). MUSCLE: multiple sequence alignment with high accuracy 700 and high throughput. Nucleic Acids Research, 32, 1792-1797. 701 25- Elith, J., Phillips, S.J., Hastie, T., Dudík, M., Chee, Y.E., & Yates, C.J. (2011). 702 A Statistical explanation of MaxEnt for ecologists. Diversity and 703 Distributions, 17, 43-57. 704 26- Evans, M.E.K., Smith, S.A., Flynn, R.S.,& Donoghue, M.J. (2009). Climate, 705 Niche evolution, and Diversification of the “Bird-Cage” Evening Primroses 706 (Oenothera, Sections Anogra and Kleinia). American Society of Naturalists, 707 173, 225-240. 708 27- Fernández-Mazuecos, M., Blanco-Pastor, J.L. & Vargas, P. (2013). A 709 phylogeny of toadflaxes (Linaria Mill.) based on nuclear internal transcribed 710 spacer sequences: systematic and evolutionary consequences. International 711 Journal of Plant Sciences, 174, 234-249. 712 28- Fernández-Mazuecos, M,, Jiménez-Mejías, P., Martín-Bravo, S., Buide, M.L., 713 Álvarez, I., & Vargas, P. (2016). Narrow endemics on coastal plains: Miocene 714 divergence of the critically endangered genus Avellara (Compositae). Plant 715 Biology, 18, 729-738. 716 29- Fiz-Palacios, O., & Valcárcel, V. (2013). From Messinian crisis to 717 Mediterranean climate: A temporal gap of diversification recovered from 718 multiple plant phylogenies. Perspectives in Plant Ecology, Evolution and 719 Systematics, 15, 130-137. 720 30- Gentili, R., Bacchetta, G., Fenu, G., Cogoni, D., Abeli, T., Rossi, G., 721 Salvatore, M.C., Baroni, C., & Citterio, S. (2015). From cold to warm-stage 722 refugia for boreo-alpine plants in southern European and Mediterranean 723 mountains: the last chance to survive or an opportunity for speciation? 724 Biodiversity, 16, 247-261. 725 31- Gianguzzi, L., Cusimano, D., Ilardi, V., & Salvatore R. (2013). Distribution, 726 ecology, vegetation and conservation survey on the relictual population of 727 Carex panormitana Guss. (Cyperaceae) in Sicily (Italy). Webbia, Journal of 728 Plant Taxonomy and Geography, 68, 159-175. 729 32- Global Carex Group. (2016). Megaphylogenetic specimen-level approaches to 730 the Carex (Cyperaceae) phylogeny using regions ITS, ETS, and matK: 731 implications for classification. Systematic Botany, 41, 500-518. 29 Page 29 of 48 Molecular Ecology 732 33- González-Martínez S.C., Dubreuil, M., Riba, M., Vendramin, G.G., 733 Sebastiani, F., & Mayol, M. (2010). Spatial genetic structure of Taxus baccata 734 L. in the western Mediterranean Basin: Past and present limits to gene 735 movement over a broad geographic scale. Molecular Phylogenetics & 736 Evolution, 55, 805-815. 737 34- Grafen, A. (1989). The phylogenetic regression. Philosophical Transactions of 738 the Royal Society of London, Series B: Biological Sciences, 326, 119-157. 739 35- Harvey, P.H., & Pagel, M.D. (1991). The Comparative Method in Evolutionary 740 Biology. Oxford, UK: Oxford University Press. 741 36- Heiberger, R.M. (2017). HH: Statistical Analysis and Data display: Heiberger 742 and Holland. R package v. 3.1-34. https://CRAN.R-project.org/package=H 743 37- Hey, J. (2010b). Isolation with Migration Models for more than two populations. 744 Molecular Biology and Evolution, 27, 905-920. 745 38- Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., & Jarvis, A. (2005). 746 Very high resolution interpolated climate surfaces for global land areas. 747 International Journal of Climatology, 25, 1965-1978. 748 39- Hosseinian Yousefkhani, S.S., Rastegar-Pouyani, E., & Aliabadian, M. 749 (2016). Ecological niche differentiation and taxonomic distinction between 750 Eremias strauchi strauchi and Eremias strauchi kopetdaghica (Squamata: 751 Lacertidae) on the Iranian Plateau based on ecological niche modeling. Italian 752 Journal of Zoology, 83, 408-416. 753 40- Janes, J.K., Miller, J.M., Dupuis, J.R., Malenfant, R.M., Gorrell, J.C., 754 Cullingham, C.I., & Andrew, R.L. (2017). The K=2 conundrum. Molecular 755 Ecology. https:// doi:10.1111/mec.14187 756 41- Jara-Arancio, P., Arroyo, M.T.K., Guerrero, P.C., Hinojosa, L.F., Arancio, & 757 G., Méndez, M.A. (2014). Phylogenetic perspectives on biome shifts in 758 Leucocoryne (Alliaceae) in relation to climatic niche evolution in western 759 South America. Journal of Biogeography, 41, 328-338. 760 42- Jiménez-Mejías, P. (2011). “Taxonomía y Sistemática de las secciones 761 Ceratocystis y Phacocystis del género Carex en Europa y la cuenca 762 mediterránea”. Universidad Pablo de Olavide, Sevilla, Spain: Tesis doctoral. 763 43- Jiménez-Mejías, P., Escudero, M., Chaparro, A.J., & Luceño, M. (2007). 764 Novedades corológicas del genero Carex para la Peninsula Iberica. Acta 765 Botanica Malacitana, 32, 305-309. 766 44- Jiménez-Mejías, P., Escudero, M., Guerra-Cárdenas, S., Lye, K.A., & Luceño, 767 M. (2011). Taxonomic delimitation and drivers of speciation in the Ibero- 768 North African Carex sect. Phacocystis river-shore group (Cyperaceae). 769 American Journal of Botany, 11, 1855-1867. 770 45- Jiménez-Mejías, P., & Luceño, M. (2011). Euro + Med Plantbase. The 771 information resource for Euro-Mediterranean plant diversity. Retrieved from 772 http://ww2.bgbm.org/EuroPlusMed/query.asp 773 46- Jiménez-Mejías, P., Martín-Bravo, S., Amini-Rad, M., & Luceño, M. (2013). 774 Disentangling the taxonomy of Carex acuta s.l. in the Mediterranean Basin 775 and the Middle East: re-evaluation of C. panormitana Guss. and C. kurdica 776 Kük. ex Hand.-Mazz. Plant Biosystem, 148, 64-73. 777 47- Jiménez-Mejías, & P, Martinetto, E. (2013). Toward an accurate taxonomic 778 interpretation of Carex fossil fruits (Cyperaceae): A case study in section 779 Phacocystis in the Western Palearctic. American Journal of Botany, 100, 780 1580-1603. 30 Molecular Ecology Page 30 of 48 781 48- Jiménez-Mejías, P., Fernández-Mazuecos, M., Amat, M.E., & Vargas, P. 782 (2015). Narrow endemics in European mountains: high genetic diversity 783 within the monospecific genus Pseudomisopates (Plantaginaceae) despite 784 isolation since the late Pleistocene. Journal of Biogeography, 42, 1455-1468. 785 49- Jiménez-Mejías, P., Martinetto, E., Momohara, A., Popova, S., Smith, S.Y., & 786 Roalson, E.H. (2016). A commented synopsis of the pre-Pleistocene fossil 787 record of Carex (Cyperaceae). The Botanical Review, 82, 258-345. 788 50- Koecke, A.V., Muellner-Riehl, A.N., Pennington, T.D., Schorr, G., & 789 Schnitzler, J. (2013). Niche evolution through time and across continents: the 790 story of Neotropical Cedrela (Meliaceae). American Journal of Botany, 100, 791 1800-1810. 792 51- Lavergne, S., Mouquet, N., Thuiller, W., & Ronce, O. (2010). Biodiversity 793 and Climate Change: Integrating Evolutionary and Ecological Responses of 794 Species and Communities. Annual Reviews, 41, 321-350. 795 52- Librado, P., & Rozas, J. (2009). DnaSP v5: software for comprehensive 796 analysis of DNA polymorphism data. Bioinformatics, 25, 1451-2. 797 53- Lima-Ribeiro, M.S., Varela, S., González-Hernández, J., Oliveira, G., Diniz- 798 Filho, J.A.F., & Terribile, L.C. (2015). ecoClimate: a database of climate data 799 from multiple models for past, present, and future for Macroecologists and 800 Biogeographers. Biodiversity Informatics, 10, 1-21. 801 54- Losos, J.B. (2008). Phylogenetic niche conservatism, phylogenetic signal and 802 the relationship between phylogenetic relatedness and ecological similarity 803 among species. Ecology Letters, 11, 995-1007. 804 55- Luceño, M., & Jiménez-Mejías, P. (2008). Carex (Cyperaceae). In: 805 Castroviejo S. et al. eds. Flora Iberica. Plantas vasculares de la Península 806 Iberica e Islas Baleares. Cyperaceae-Pontederiaceae, 18, 109-250. 807 56- Maguilla, E., Escudero, M., Waterway, M.J., Hipp, A.L., & Luceño, M. 808 (2015). Phylogeny, Systematics, and trait evolution of Carex section 809 Glaerosae. American Journal of Botany, 102, 1-17. 810 57- Mairal, M., Sanmartín, I., & Pellissier, L. (2017). Lineage-specific climatic 811 niche drives the tempo of vicariance in the Rand Flora. Journal of 812 Biogeography, 44, 911-923. 813 58- Mansion, G., Rosenbaum, G., Schoenenberger, N., Bacchetta, G., Rosselló, 814 J.A., & Conti, E. (2008). Phylogenetic analysis informed by geological history 815 supports multiple, sequential invasions of the Mediterranean Basin by the 816 angiosperm family Araceae. Systematic Biology, 57, 269-285. 817 59- Martín-Bravo, S. & Escudero, M. (2012). Biogeography of flowering plants: a 818 case study in mignonettes (Resedaceae) and sedges (Carex, Cyperaceae). In: 819 Stevens, L. (ed.). Global Advances in Biogeography, 257-290. 820 60- Médail, F., & Quézel, P. (1997). Hot-spots analysis for conservation of plant 821 biodiversity in the Mediterranean Basin. Annals of the Missouri Botanical 822 Garden, 84, 112–127. 823 61- Médail, F., & Diadema, K. (2009). Glacial refugia influence plant diversity 824 patterns in the Mediterranean Basin. Journal of Biogeography, 36, 1333– 825 1345. 826 62- Merow, C., Smith, M.J., & Silander, Jr. J.A. (2013). A practical guide to 827 MaxEnt for modeling species’ distributions: what it does, and why inputs and 828 settings matter. Ecography, 36, 1058-1069. 829 63- Míguez, M., Gehrke, B., Maguilla, E., Jiménez-Mejías, P., & Martín-Bravo, S. 830 (2017). Carex sect. Rhynchocystis (Cyperaceae): a Tertiary subtropical relict Molecular Ecology Page 30 of 48 30 831 in the western Palearctic showing a dispersal derived Rand-Flora Pattern. 832 Journal of Biogeography. https://doi:10.1111/jbi.13027 833 64- Milne, R.I., & Abbott, R.J. (2002). The origin and evolution of tertiary relict 834 floras. Advances in Botanical Research, 38, 281-314. 835 65- Mod, H.K., Scherrer, D., Luoto, M. & Guisan, A. (2016): What we use is not 836 what we know: environmental predictors in plant distribution models. Journal 837 of Vegetation Science, 27, 1308-1322. 838 66- Münkemüller, T., Boucher, F.C., Thuiller, & W., Lavergne, S. (2015). 839 Phylogenetic niche conservatism-common pitfalls and ways forward. 840 Functional Ecology, 29, 627-639. 841 67- Myers, N., Mittermeier, R.A., Mittermeier, C.G., da Fonseca, G.A.B., & Kent, 842 J. (2000). Biodiversity hotspots for conservation priorities. Nature, 403, 853– 843 858. 844 68- Peakall, R., & Smouse, P.E. (2012). GenAlEx 6.5: genetic analysis in Excel. 845 Population genetic software for teaching and research-an update. 846 Bioinformatics, 28, 2537-2539. 847 69- Pettengill, J.B. & Moeller, D.A. (2012). Phylogeography of speciation: 848 allopatric divergence and secondary contact between outcrossing and selfing 849 Clarkia. Molecular Ecology, 21, 4578-4592. 850 70- Phillips, S.J., Anderson, R.P., Dudík, M., Schapire, R.E. & Blair, M.E. (2017). 851 Opening the black box: an open-source release of Maxent. Ecography, 40, 852 887-893. 853 71- Pignatti, S. (1982). Flora d’Italia. Bologna, Italy: Edagricole. 854 72- Pritchard, J.K., Stephens, M., & Donnelly, P. (2000). Inference of population 855 structure using multilocus genotype data. Genetics 155: 945-959. 856 73- Qiu, Y.X., Qi, X.S., Jin, X.F., Tao, X.Y., Fu, C.X., Naiki, A., & Comes, H.P. 857 (2009). Population genetic structure, phylogeography, and demographic 858 history of Platycrater arguta (Hydrangeaeae) endemic to East China and 859 South Japan, inferred from chloroplast DNA sequence variation. Taxon, 58, 860 1226-1241. 861 74- R Core Team. (2017). R: A language and environment for statistical 862 computing. R Foundation for Statistical Computing. Vienna, Austria. 863 Retrieved from https://www.r-project.org/ 864 75- Ramasamy, R.K., Ramasamy, S., Bindroo, B.B., & Naik V.G. (2014). 865 STRUCTURE PLOT: a program for drawing elegant STRUCTURE bar plots 866 in user friendly interface. SpringerPlus. doi: 10.1186/2193-1801-3-431. 867 76- Rambaut, A. (2012). FigTree, tree figure drawing tool v. 1.4.0. Retrieved from 868 http://tree.bio.ed.ac.uk/software/figtree/ 869 77- Rambaut, A., & Drummond, A.J. (2014). Tracer v. 1.6 Retrieved from 870 http://beast.bio.ed.ac.uk/Tracer 871 78- Revell, L.J., Harmon, L.J., & Collar, D.C. (2008). Phylogenetic Signal, 872 Evolutionary Process, and Rate. Systematic Biology, 57, 591-601. 873 79- Rodríguez-Sánchez, F., Pérez-Barrales, R., Ojeda, F., Vargas, P., & Arroyo, J. 874 (2008). The Strait of Gibraltar as a melting pot for plant biodiversity. 875 Quaternary Science Reviews, 27, 2100-2117. 876 80- Rundel, P.W., Arroyo, M.T.K., Cowling, R.M., Keeley, J.E., Lamont, B.B., & 877 Vargas, P. (2016). Mediterranean Biomes: Evolution of their Vegetation, 878 Floras, and Climate. Annual review of Ecology, Evolution and Systematics, 47, 879 383-407. 32 Molecular Ecology Page 32 of 48 880 81- Schoener, T.W. (1968). Anolis lizards of Bimini: resource partitioning in a 881 complex fauna. Ecology, 49, 704-726. 882 82- Shaw, J., Lickey, E.B., Beck, J.T., Farmer, S.B., Liu, W., Miller, J., Siripun, 883 K., Winder, C.T., Schilling, E.E., & Small, R.L. (2005). The tortoise and the 884 Hare II: relative utility of 21 noncoding chloroplast DNA sequences for 885 phylogenetic analysis. American Journal of Botany, 92, 142-166. 886 83- Shaw, J., Lickey, E.B., Schilling, E.E., & Small, R.L. (2007). Comparison of 887 whole chloroplast genome sequences to choose noncoding regions for 888 phylogenetic studies in angiosperms: the Tortoise and the Hare III. American 889 Journal of Botany, 94, 275-288. 890 84- Spalink, D., Drew, B.T., Pace, M.C., Zaborsky, J.G., Li, P., Cameron, K.M., 891 Givnish, T.J., & Sytsma, K.J. (2016). Evolution of geographical place and niche 892 space: Patterns of diversification in the North American sedge (Cyperaceae) 893 flora. Molecular Phylogenetics and Evolution, 95, 183-195. 894 85- Stiels, D., Gaißer, B., Schidelko, K., Engler, J.O., & Rödder, D. (2015). Niche 895 shift in four nonnative estrildid finches and implications for species distribution 896 models. Ibis, 157, 75–90. 897 86- Suc, J.P., Bertini, A., Combourieu-Nebout, N., Diniz, F., Leroy, S., Russo- 898 Ermolli, E., Zheng, Z., Bessais, & E., Ferrier, J. (1995). Structure of West 899 Mediterranean vegetation and climate since 5.3 ma. Acta Zoologica Cracova, 900 38, 3-16. 901 87- Swets, J.A. (1988). Measuring the accuracy of diagnostic systems. Science, 902 240, 1285. 903 88- Swofford, D.L. (2002). PAUP*: Phylogenetic analysis using parsimony (and 904 other methods) v. 4.0b10. Sinauer, Sunderland, Massachusetts, USA. doi: 905 10.1111/j.0014-3820.2002.tb00191. 906 89- Taberlet, P., Fumagalli, L., Wust-Saucy, A.G., & Cosson, J.F. (1998). 907 Comparative phylogeography and postglacial colonisation routes in Europe. 908 Molecular Ecology, 7, 453–464. 909 90- Terrab, A., Schönswetter, P., Talavera, S., Vela, E., & Stuessy, T.F. (2008). 910 Range-wide phylogeography of Juniperus thurifera L., a presumptive 911 keystone species of western Mediterranean vegetation during cold stages of 912 the Pleistocene. Molecular Phylogenetics & Evolution, 48, 94-102. 913 91- Thiers, B. (2015). Index Herbariorum: A global directory of public herbaria 914 and associated staff. New York Botanical Garden's Virtual Herbarium. 915 Retrieved from http://sweetgum.nybg.org/science/ih/ 916 92- Thompson, J.D. (2005). Plant evolution in the Mediterranean. Oxford, UK: 917 Oxford University Press. 918 93- Thuiller, W. (2003). BIOMOD – optimizing predictions of species 919 distributions and projecting potential future shifts under global change. Global 920 Change Biology 9: 1353-1362. 921 94- Tremetsberger, K., Ortiz, M.A., Terrab, A., Balao, F., Casimiro-Soriguer, R., 922 Talavera, M. & Talavera, S. (2016). Phylogeography above the species level 923 for perennial species in a composite genus. AoB Plants, 8, 142- 159. 924 95- Tripp, E.A., & McDade, L.A. (2014). A rich fossil record yields calibrated 925 phylogeny for Acanthaceae (Lamiales) and evidence for marked biases in 926 timing and directionality of intercontinental disjunctions. Systematic Biology, 927 63, 660-684. 33 Page 33 of 48 Molecular Ecology 928 96- Urbani, M., Gianguzzi, L., & Ilardi, V. (1995). Notes on the distribution and 929 ecology of Carex panormitana Guss. (Cyperaceae). Giornale Botanico 930 Italiano, 129, 186. 931 97- Urbani, M., Calvia, G., & Pisanu, S. (2013). Carex panormitana Guss. In: 932 Rossi G, Foggi B, Gennai M, Gargano D, Montagnani C, Orsenigo S, Pedrini, 933 S, eds. Schede per una Lista Rossa della Flora vascolare e crittogamica 934 Italiana, 45, 134–136. 935 98- Vargas, P. (2007). Are Macaronesian Islands refugia of relict plant lineages?: 936 a molecular survey. In: Weiss S., Ferrand N. (eds) Phylogeography of 937 Southern European Refugia. Springer, Dordrecht. 938 99- Warren, D.L., Glor, R.E., & Turelly, M. (2008). Environmental niche 939 equivalency versus conservatism: Quantitative approaches to niche evolution. 940 Evolution, 62, 2868-2883. 941 100- Webb, C.O., Ackerly, D.D., McPeek, M.A., & Donoghue, M.J. (2002). 942 Phylogenies and community ecology. Annual Review of Ecology and 943 Systematics, 33, 475–505. 944 101- Wickham, H. (2009). ggplot2: Elegant Graphics for Data Analysis. Springer- 945 Verlag New York. 946 102- Wiens, J.J. (2004a). Speciation and Ecology revisited: Phylogenetic niche 947 conservatism and the origin of species. Evolution, 58, 193-197. 948 103- Wiens, J.J. (2004b). What is speciation and how should we study it? The 949 American Naturalist, 163, 914-923. 950 104- Wiens, J.J., Ackerly, D.D., Allen, A.P., Anacker, B.L., Buckley, L.B., 951 Cornell, H.V., Damschen, E.I. et al. (2010). Niche conservatism as an 952 emerging principle in ecology and conservation biology. Ecology Letters, 13, 953 1310 – 1324. 954 105- Wiens, J.J., & Graham, C.H. (2005). Niche Conservatism: Integrating 955 Evolution, Ecology, and Conservation Biology. Annual Review of Ecology, 956 Evolution, and Systematics, 36, 519-539. 957 106- Won, Y.J., & Hey, J. (2005). Divergence population genetics of chimpanzees. 958 Molecular Biology and Evolution, 22, 297-307. 959 107- Zitari, A., Tranchida-Lombardo, V., Cafasso, D., Noureddine, H., Scopece, 960 G., & Cozzolino, S. (2011). The disjointed distribution of Anacamptis 961 longicornu in the West-Mediterranean: the role of vicariance versus long- 962 distance seed dispersal. Taxon, 60, 1041-1049. 963 964 Figure legends 965 Figure 1. Populations used in the genetic study. Extant populations are displayed with 966 circles. Carex reuteriana ssp. mauritanica is showed in two different green tones 967 according to the two detected subgroups (see Fig. 2): (A) Sierra Morena populations, 968 and (B) South of Guadalquivir valley plus Moroccan populations. Shaded striped areas 969 depict the approximate known distribution range of the two species according to 970 Jiménez-Mejías et al. (2011). The unshaded populations in Sicily and Tunisia 971 correspond to the only populations known there Occurrences of the Pliocene fossils are 972 displayed using squares. 973 Figure 2. Clustering of individuals by STRUCTURE (see results): each individual is 974 represented by a vertical bar, same color in different individuals indicates that they 975 belong to the same genetic cluster. The four genetic clusters (K=4) retrieved in (a) 976 correspond from left to right to: (1) Sardinian populations and (2) Sicilian-Tunisian 34 Molecular Ecology Page 34 of 48 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 populations belonging to C. panormitana, (3) C. reuteriana ssp. reuteriana and (4) C. reuteriana ssp. mauritanica. In (b) is shown another STRUCTURE analysis results with six genetic clusters (K=6) due to the subspecies mauritanica splitting in two subgroups (South of Guadalquivir valley plus Moroccan populations and Sierra Morena populations) and subspecies reuteriana splitting in other two (two CW Spain populations distinct from the rest). Principal coordinate analysis scatter plots of AFLP data from: (c) the whole dataset, (d) C. reuteriana s.l. and (e) C. panormitana. Populations from C. panormitana are split in three different groups corresponding to each geographic area (see Fig. 1). Figure 3. (a) Maximum-credibility-clade chronogram from the Bayesian divergence time analysis of the combined nuclear-plastid data set showing the study group and related lineages within Carex sect. Phacocystis. Branch lengths correspond to the time scale (Ma). Node ages are given above branches and posterior probability values (>0.80) below branches. Node bars represent the 95% highest posterior density (HPD) intervals for the divergence time estimates of the nodes (see Fig. S1 in supporting information for the full tree). (b) Densitree of the C. reuteriana-C. panormitana complex showing the complete genealogical tree set (thin lines) and consensus tree (thick lines) from SNAPP analysis of AFLP data. Node ages are given above branches (Ma) and the posterior probability below branches. Both values are showed only for nodes with supports equal or greater than 0.8. Figure 4. Climatic suitability and potential distributions (a) at present as predicted by MaxEnt (Phillips et al., 2006) for from left to right: Carex reuteriana ssp. reuteriana, ssp. mauritanica and C. panormitana. (b) The left plot shows the distribution projected for the species complex (Carex reuteriana-C. panormitana) at the middle Pliocene (3 Ma) with a model calibrated with present occurrences, and the right plot displays the potential distribution of the species complex at the middle Pliocene according to the model calibrated with fossil records. Figure 5. Comparison of fitted responses to the four bioclimatic variables of (a) extant species complex (continuous line) vs. Pliocene C.panormitana-type fossils (dashed line), and (b) extant taxa only (dashed line represents C. reuteriana ssp. mauritanica, continuous line ssp. reuteriana, and dotted line C. panormitana). Colours are as in Fig. 1, except for all populations of C. panormitana which are displayed in red. Figure 6. Differences in the environmental niche of extant and fossil populations of Carex reuteriana-C. panormitana species complex. Both plots depict the population values in an environmental space defined by the first two axes from a principal components analysis (PCA) of the four bioclimatic variables (see Fig. 5) and soil pH. (a) Comparison of Carex reuteriana-C. panormitana extant complex (orange) and their Pliocene fossil ancestor (b) Comparison of populations of the extant species complex; coloured dots represent the different taxa/population groups according to the retrieved genetic structure. Page 41 of 48 Molecular Ecology Differences in the environmental niche of extant and fossil populations of Carex reuteriana-C. panormitana species complex. Both plots depict the population values in an environmental space defined by the first two axes from a principal components analysis (PCA) of the four bioclimatic variables (see Fig. 5) and soil pH. (a) Comparison of Carex reuteriana-C. panormitana extant complex (orange) and their Pliocene fossil ancestor (b) Comparison of populations of the extant species complex; coloured dots represent the different taxa/population groups according to the retrieved genetic structure. 190x120mm (300 x 300 DPI) Molecular Ecology Page 42 of 48 Table 1. Geographic location of each sampled population of the C. reuteriana-C. panormitana complex, number of individuals included in each molecular study (AFLPs and DNA sequences), haplotype number (ptDNA and nDNA) indicating the number of individuals included in each one, voucher, and herbarium (acronyms according to Index Herbariorum; Thiers, 2015) where specimens are deposited. Labelling of the populations specifies the taxa (REU = Carex reuteriana ssp. reuteriana; MAU = Carex reuteriana ssp. mauritanica; PAN = Carex panormitana), and the country following TDWG botanical countries nomenclature (Brummitt, 2001) [MOR = Morocco, POR = Portugal, SAR = Sardinia, SIC = Sicily, SPA = Spain, TUN = Tunisia]. One of the Tunisian populations (PAN_TUN-Tub) has only been included in the DNA sequence study. Haplotype numbers are represented in the Figure S2a. Page 43 of 48 Molecular Ecology Taxon/ Population Locality Longitude/ Latitude Ni AFLPs/ Haplotype (ptDNA) Haplotype (nDNA) Voucher/ Herbarium Ni DNA sequences (individuals) (individuals) C. reuteriana s.l. C. reuteriana ssp. reuteriana Portugal, Tras os M. Escudero et al., REU_POR-TM REU_POR-BL Montes, Lamego, Bigorne, Petrarouca Portugal, Beira Litoral, Coimbra, Lousã Spain, Ávila, Sierra de -7.88/41.03 6/5 H2(4) H1(5) H1(4)/ -8.23/40.10 7/5 H2(5) H8(1) H5(2)/ 37ME07 (UPOS- 7374) M. Escudero et al., 60ME07 (UPOS- 7373) J.M. Marín, REU_SPA-Av Gredos, Las Chorreras del Tormes -5.16/40.34 5/5 H2(5) H6(2)/ H7(1) 5504JMM (UPOS- 1004) Molecular Ecology Page 44 of 48 REU_SPA-CcN Spain, Cáceres, Valley of Jerte -5.75/40.22 4/4 H2(3) H1(3)/ H7(1) P. Jiménez-Mejías & I. Pulgar, 57PJM07 (UPOS- 6957) REU_SPA-CcS Spain, Cáceres, Ibor river Spain, H1(3)/ -5.44/39.62 4/5 H2(5) H3(2) H1(3)/ P. Jiménez-Mejías & G.E. Rodríguez , 24PJM13 (UPOS- 5449) P. Jiménez-Mejías REU_SPA-To Toledo, Navalucillos -4.66/39.64 3/5 H2(3) H3(2) & G.E. Rodríguez, 60PJM13 (UPOS- 5479) C. reuteriana ssp. mauritanica Spain, P. Jiménez-Mejías, Sevilla, El 35PJM07 (UPOS-7372) MAU_SPA-Se Ronquillo, -6.17/37.67 7/5 H3(4) H3(5) Rivera de Huelva Spain, Cádiz, MAU_SPA-CaGu El Gastor, Guadalete -5.45/36.88 5/5 H1(5) H4(5) P. Jiménez-Mejías, 34PJM07 (UPOS-7371) river Spain, Cádiz, P. Jiménez-Mejías & I. MAU_SPA-CaAl Alcornocales -5.59/36.55 8/5 H1(5) H3(5) Pulgar, 17PJM07 Natural Park, (UPOS) Page 45 of 48 Molecular Ecology MAU_SPA-J Spain, Jaén, Despeñaperr os -3.06/38.39 8/5 H3(5) H3(5) P. Jiménez-Mejías & L. Reina, 67PJM09 (UPOS) MAU_MOR-Lao MAU_MOR-Lou C. panormitana Tunisia-Sicily Morocco, Tanger, Rif, Oued Laou Morocco, Tanger, Rif, Oued Loukos Tunisia, -5.30/35.14 5/5 H3(5) -5.44/35.03 7/5 H3(5) H1(1)/ H2(3)/ H3(1) H1(1)/ H2(3)/ H5(1) A.J. Chaparro et al., 8AJC05 (UPOS-1637) A.J. Chaparro et al., 3AJC05 (UPOS-1630) P. Jiménez-Mejías & PAN_TUN Jendouba, El Feija National Park 8.31/ 36.49 9/6 H3(5) H9(6) G.E. Rodríguez, 132PJM13 (UPOS-6636) J. Höller PAN_TUN-Tub Tunisia, Tubarka 9.00/37.00 -/1 H3(1) H9(1) (M-0223053) PAN_SIC Italy, Sicily, Fiume Oreto 13.34/ 38.09 14/14 H3(14) H9(10)/ D. Cusimano s.n. (SS) H10(4) Molecular Ecology Page 46 of 48 Sardinia PAN_SAR-Bau Italy, Sardinia, Bau Mela river, Villagrande Italy, Sardinia, 9.42/39.98 9/5 H3(5) H9(5) M. Urbani s.n., 2013 (SS) PAN_SAR-Pira Cantoniera, Pirae´onni, Villagrande Italy, Sardinia, 9.40/ 40.02 10/5 H3(5) H9(5) M. Urbani s.n., 2013 (SS) M. Urbani s.n., 2013 PAN_SAR-Ber Ramacaso river, Berchidda Italy, Sardinia, 9.24/ 40.82 8/5 H3(4) H9(5) (SS) PAN_SAR-Cal Miriacheddu river, Calangianus 9.26/ 40.89 11/5 H3(5) H9(5) M. Urbani s.n., 2013 (SS) Page 47 of 48 Molecular Ecology Table 2. (a) Divergence dates of the most important clades resolved in the divergence time estimation analysis of Carex under an uncorrelated lognormal clock model using the combined matrix of ETS1f, ITS, psbA-trnH, rpl32-trnLUAG and ycf6- psbM regions in BEAST (see Fig. 3a, Fig. S1). Posterior probabilities, mean time to the most common recent ancestor in Ma and 95% highest posterior density (HDP) intervals are shown (see also Fig. S1 in Supporting Information). (b) Divergence dates of the most important clades resolved in the divergence time estimation from AFLP data using a multispecies coalescent approach implemented in SNAPP (see Fig. 3b). Posterior probabilities, mean time to the most common recent ancestor in Ma and 95% HPD intervals are shown. reuteriana (Stem node) Lineage B: C. nigra ssp. transcaucasica, C. trinervis, C. nigra, C. appendiculata, C.cespitosa, C. omskiana (Stem node) 1 3.71 3.51 4.29 b C. reuteriana 1 0.576 0.430 0.732 Iberian lineage (ssp. reuteriana) 1 0.314 0.194 0.435 Afro-Iberian lineage (ssp. mauritanica) 0.97 0.432 0.314 0.558 C. panormitana 0.86 0.745 - - Sardinian lineage 1 0.147 0.009 0.215 Sicilian-Tunisian lineage 0.95 0.525 0.320 0.766 probability (Ma) interval (Ma) Clade Posterior Median 95% HPD Max. Min. a Core Phacocystis 0.62 7.97 4.93 11.28 (2) Lineage A+ Lineage B 0.84 4.79 3.8 6.43 Lineage A: C. panormitana – C. 1 0.95 0.19 2.48 Molecular Ecology Page 48 of 48 Table 3. Pairwise statistical test for comparison of ecological niche overlap between the different genetic and taxonomic groups that compose the current complex, so as with its fossil ancestor. Their statistical significance is represented by p-values (Warren, Glor & Turelly, 2008). Hypotheses: “The niche overlap is more similar/equivalent than at random” Numbers of populations (N) Schoener’s D Niche similarity (1->2) Niche similarity (2->1) Niche equivalency p p p Among groups populations Complex vs. fossils 462/17 0 1 1 - C. reuteriana vs. C. panormitana 433/29 0.193 0.129 0.139 0.525 C. reuteriana ssp. reuteriana vs. ssp. 315/118 mauritanica 0.095 0.04 0.04 1 C. panormitana Sicilian-Tunisian vs. Sardinian 8/21 0.169 0.14 0.158 0.604