Environment and biodiversity: a multidisciplinary approach to dynamic patterns in the Iberian Peninsula

. Environment and biodiversity A multidisciplinary approach to dynamic patterns in the Iberian Peninsula José Pedro Rodrigues Tarroso Gomes Supervisors: Paulo Célio Alves José Carlos Brito Rachid Cheddadi Porto 2013 . Environment and biodiversity A multidisciplinary approach to dynamic patterns in the Iberian Peninsula José Pedro Rodrigues Tarroso Gomes Tese apresentada à Faculdade de Ciências da Universidade do Porto para a obtenção do grau de Doutor em Biodiversidade, Genética e Evolução Porto 2013 Este trabalho foi apoiado pela Fundação para a Ciência e Tecnologia (FCT) através da atribuição da bolsa de doutoramento (SFRH/BD/42480/2007) Nota Prévia Na elaboração desta tese, e nos termos do número 2 do Artigo 4odo Regulamento Geral dos Terceiros Ciclos de Estudos da Universidade do Porto e do Artigo 31odo D.L. 74/2006, de 24 de Março, com a nova redação introduzida pelo D.L. 230/2009, de 14 de Setembro, foi efetuado o aproveitamento total de um conjunto coerente de trabalhos de investigação já publicados ou submetidos para publicação em revistas internacionais indexadas e com arbitragem científica, os quais integram alguns dos capítulos da presente tese. Tendo em conta que os referidos trabalhos foram realizados com a colaboração de outros autores, o candidato esclarece que, em todos eles, participou ativamente na sua conceção, na obtenção, análise e discussão de resultados, bem como na elaboração da sua forma publicada. To the memory of my mother. To Anna and Inês. Acknowledgements I fell back asleep some time later on And I dreamed the perfect song It held all the answers, like hands laid on I woke halfway and scribbled it down And in the morning what I wrote I read It was hard to read at first but here’s what it said Eid ma clack shaw Zupoven del ba Mertepy ven seinur Cofally ragdah — BILL CALLAHAN,Eid ma clack shaw I am able to say now that I understand the pathological psychosis inherent to the process of writing a thesis. References are a very important part of the scientific method, providing to the reader the source of information and the path to find it in the massive network of scientific works. Looking for references is, thus, a game to play most likely during all the scientific career. However, in the process of writing this thesis, the quest for the reference has expanded from the scientific sphere to the “real” world. In fact, I felt all the climate events during the process of starting, building and writing a PhD as psychological oscillations. This thesis had a last glacial maximum (∼4 yr before thesis, hereafter BT) when I didn’t knew exactly how to start. A trigger event was required, soon followed by a warming trend, sun shining on a wide extent of the brain, ice gone and, thus, less albedo and more grey matter. Dating of this event was not easy, but around 3.5 yr BT I went to Montpellier to start working with Rachid Cheddadi. I am deeply grateful to him for showing me the joy of doing science, even if it was with fossil pollen! He is enthusiastic about science and that is contagious. I was able to deal with the complex pollen database and climate reconstruction methods until the onset of the Younger Dryas. It was not easy. Things can go wrong but persistence in a cryptic refugium can xvi FCUP analysis revealed that biodiversity patterns were largely affected by these processes, with areas of high velocities of species composition change, but also areas with low velocities related to long species persistence. By the end of the current century, major changes on species richness patterns are predicted that will pose many conservation challenges. The application of Simapse to the study of the hybrid zone at micro-scale resulted in an innovative approach to the study of ecological divergence in the range limits of the three vipers. The results allowed a full description of the contact zone, both at genetic and environmental levels, with a characterization of the population structure. Different ecological requirements from each species were found and a transitional ecotone was suggested where the hybrids were identified. The integrative approach followed here provided exhaustive examples describing macroand micro-scale processes occurring in Iberian Peninsula, confirming the potential of merging results from different research fields. The results presented here are expected to support the conservation effort on Iberian herpetofauna and to have impact on future studies of hybrid zones within an ecological context. However, the potential of the developed tools is not limited to the analyses carried here, and it is also expected the expansion of the applicability to other domains related to biodiversity. Résumé L’approche multidisciplinaire pour l’étude des modèles de biodiversité offre des analyses détaillées qui sont souvent difficiles à atteindre dans le cadre d’études monodisciplinaires. Ceci est démontré ici avec une approche intégrative combinant des données issues de différents domaines de recherche incluant l’écologie du paysage, l’évolution moléculaire, la paléoécologie et la modélisation du climat. La péninsule ibérique offre d’excellentes conditions pour développer ces études car il y a un fort taux d’endémisme, une variabilité du climat exceptionnel et une abondance de données disponibles issues de différents domaines. L’approche intégrative présentée ici est étendue avec différentes échelles temporelles et spatiales, couvrant une large période de la fin du Quaternaire, (15.000 à 3.000 ans avant le présent) ainsi que le siècle actuel avec différentes questions abordées à des échelles différents (macro et micro). Dans le cadre de ce sujet de recherche, deux nouveaux logiciels ont été développés (chapitre 3): 1) E-Clic (section 3.1) est un convertisseur de données de prédiction climatiques futures à des formats usuels utilisés dans l’écologie du paysage et 2) Simapse (section 3.2) est un outil de modélisation de la niche écologique qui met en œuvre une méthode statistique basée sur les réseaux de neurones artificiels pour modéliser la distribution des espèces. La méthode de reconstruction de variables climatiques à partir de données polliniques fossiles a également été adaptée et améliorée pour la quantification de treize variables représentant trois périodes de temps différentes de la fin du Quaternaire (chapitre 4, section 4.1). L’intégration de toutes les données et les méthodes a été réalisée à deux échelles différentes (chapitre 5). Dans un contexte macro-échelle, j’ai analysé la dynamique de la composition spécifique de l’herpétofaune ibérique dans le passé et j’ai effectué une simulation prédictive pour le siècle actuel sous la tendance d’un réchauffement climatique (section 5.1). Une analyse micro-échelle a été réalisée avec l’étude de la divergence écologique dans une zone de contact entre les espèces de trois vipères (Vipera latastei,V. aspis et V. seoanei) en utilisant la modélisation de la niche xviii FCUP écologique et l’analyse multi-traits (nucléaire, mitochondrial et morphologique; section 5.2). Le climat ibérique à la fin du Quaternaire a été caractérisé par une tendance générale au réchauffement avec des transitions brusques et des conséquences importantes sur l’organisation spatiale des variables climatiques. L’analyse de l’évolution du climat a permis de classerla péninsule ibérique en quatre zones distinctes qui partagent une même dynamique climatique (température et précipitations). Les résultats de l’analyse macro-échelle a révélé que les patrons de biodiversité ont été largement affectés par ces processus climatiques, avec des zones ayant des vitesses élevées de changement de la composition des espèces, mais aussi des zones avec des vitesses faibles liées à la longue persistance des espèces. À la fin de ce siècle, des changements importants dans les patrons de la richesse spécifique sont prédits et ils posent de nombreux problèmes de conservation. L’utilisation de Simapse dans l’étude micro-échelle de la zone de contact entre les trois espèces de vipères a permis une nouvelle approche pour aborder la divergence écologique dans les limites de leur distribution géographique. Les résultats ont permis la description génétique et écologique de la zone de contact et la caractérisation de la structure de la population. Des exigences écologiques différentes de chaque espèce ont été mises en évidence et un écotone où les hybrides ont été identifiés a été suggéré. L’approche intégrative suivie ici a fourni des exemples décrivant les processus macro-et micro-échelle qui se produisent dans la péninsule ibérique, ce qui confirme le potentiel de la combinaison des résultats issus de différents domaines de la recherche. Les résultats présentés ici ont pour but de soutenir l’effort de conservation de l’herpétofaune dans la péninsule ibérique et d’avoir un impact sur les futures études de zones hybrides dans un contexte écologique. Le potentiel des outils développés dans cette thèse ne se limite pas aux analyses effectuées ici, et une extension est également prévue de son applicabilité à d’autres domaines liés à la biodiversité. Contents 1 Introduction 1 1.1 Climateevolution............................... 1 1.1.1 General trends in climate oscillations . . . . . . . . . . . . . . . . 1 1.1.2 Europeanclimate........................... 3 1.2 Species sensitivity to climate change . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Consequences of past climate change . . . . . . . . . . . . . . . 5 1.2.2 Predicted impacts of future changes . . . . . . . . . . . . . . . . 8 1.3 Evolutionary patterns in the Iberian Peninsula . . . . . . . . . . . . . . . 9 1.3.1 Characteristics of Iberia . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.2 Climate change and landscape dynamics . . . . . . . . . . . . . 10 1.3.3 Cryptic refugia in a glacial refugium . . . . . . . . . . . . . . . . 12 1.4 Reconstructing the past, predicting the future . . . . . . . . . . . . . . . 13 1.4.1 Pollen as a proxy of the past . . . . . . . . . . . . . . . . . . . . 13 1.4.2 Quantitative reconstructions of climate . . . . . . . . . . . . . . . 15 1.4.3 Predicting future climate . . . . . . . . . . . . . . . . . . . . . . . 16 1.5 Ecological niche modelling as a tool to study niche dynamics . . . . . . 20 1.5.1 Niche concept and traits . . . . . . . . . . . . . . . . . . . . . . . 20 1.5.2 Ecological niche-based modelling . . . . . . . . . . . . . . . . . 22 1.5.3 Artificial neural networks . . . . . . . . . . . . . . . . . . . . . . . 24 1.5.4 Ensemble modelling and model evaluation . . . . . . . . . . . . 27 1.6 References .................................. 29 2 Objectives 47 2.1 Generalobjectives.............................. 47 2.2 Detailed objectives and structure of the thesis . . . . . . . . . . . . . . . 47 3 Analyzing species’ distributions with ecological niche modelling 51 3.1 E-Clic - Easy climate data converter . . . . . . . . . . . . . . . . . . . . 51 xx FCUP Contents 3.1.1 Abstract................................ 51 3.1.2 Easy climate data converter . . . . . . . . . . . . . . . . . . . . . 52 3.1.3 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.1.4 References.............................. 57 3.2 Simapse - Simulation Maps for Ecological Niche Modelling . . . . . . . 59 3.2.1 Abstract................................ 59 3.2.2 Introduction.............................. 59 3.2.3 Simapse................................ 60 3.2.4 Input data and general options . . . . . . . . . . . . . . . . . . . 60 3.2.5 Outputresults............................. 63 3.2.6 Example................................ 64 3.2.7 Discussion .............................. 65 3.2.8 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.2.9 References.............................. 66 3.3 Evaluating ecological niche models with virtual and real species . . . . 68 3.3.1 Abstract................................ 68 3.3.2 Introduction.............................. 69 3.3.3 Methods................................ 71 3.3.4 Results ................................ 75 3.3.5 Discussion .............................. 79 3.3.6 Conclusions.............................. 84 3.3.7 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.3.8 References.............................. 84 4 Reconstruction of past Iberian climate 89 4.1 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. . 89 4.1.1 Abstract................................ 89 4.1.2 Introduction.............................. 90 4.1.3 Methods................................ 91 4.1.4 Results ................................ 96 4.1.5 Discussion .............................. 99 4.1.6 Conclusions..............................102 4.1.7 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.1.8 References..............................103 5 Spatial dynamic patterns in the Iberian Peninsula at two scales 109 5.1 Velocity of biodiversity change: a case study in a global hotspot . . . . . 109 FCUP xxi Contents 5.1.1 Abstract................................109 5.1.2 Introduction..............................110 5.1.3 Methods................................112 5.1.4 Results ................................116 5.1.5 Discussion ..............................119 5.1.6 Finalremarks.............................126 5.1.7 References..............................127 5.2 Hybridization at an ecotone: ecological and genetic barriers between threeIberianvipers..............................134 5.2.1 Abstract................................134 5.2.2 Introduction..............................134 5.2.3 Material and Methods . . . . . . . . . . . . . . . . . . . . . . . . 137 5.2.4 Results ................................143 5.2.5 Discussion ..............................148 5.2.6 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . 154 5.2.7 References..............................154 6 General Discussion 161 6.1 Applications for ecological niche modelling . . . . . . . . . . . . . . . . . 161 6.2 Macro-scale: Integrating space and time in biodiversity studies . . . . . 163 6.3 Micro-scale: environmental divergence in an hybrid zone . . . . . . . . 165 6.4 Futureprospects...............................166 6.5 References ..................................168 7 Conclusions 173 8 Appendices 177 Appendix A 179 Appendix B 187 Appendix C 195 Appendix D 199 Appendix E 233 List of Tables 1.1 Predicted temperatures for 2100 from different emission scenarios. . . . 19 3.1 Overview of Simapse’s general options . . . . . . . . . . . . . . . . . . 62 3.2 Ecological niche variables data and usage by virtual species . . . . . . 71 4.1 Origin and description of the data sources of fossil pollen . . . . . . . . 93 5.1 Environmental factors used in the ecological niche modelling . . . . . . 141 5.2 Assignment of analysed individuals to each cluster after classification ofCMPs....................................143 5.3 Number of individuals assigned to each species . . . . . . . . . . . . . 148 D.1 Detailed model results for the current species distributions . . . . . . . . 200 E.1 Sampled individuals with morfological identifications, cluster membership probabilities and mtDNA results . . . . . . . . . . . . . . . . . . . . 233 List of Figures 1.1 Major orbital perturbations defining long term climate cycles . . . . . . . 2 1.2 Extension of the ice sheet during the last glacial maximum . . . . . . . . 4 1.3 The Iberian Peninsula and Balearic island with major topographic features 10 1.4 Relationship between pollen percentages and taxa density . . . . . . . 17 1.5 IPCC family of emissions scenarios . . . . . . . . . . . . . . . . . . . . 18 1.6 Biotic-Abiotic-Movement diagram representing driving forces of species’ distributions.................................. 21 1.7 Typical structure of a simple artificial neural network with backpropagation algorithm and one hidden layer . . . . . . . . . . . . . . . . . . . . . 27 3.1 E-Clic built-in graphical user interface . . . . . . . . . . . . . . . . . . . 53 3.2 E-Clic interface as a toolbox . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3 Average daily temperature forecasts for the years 2025, 2050, 2075 and 2100...................................... 55 3.4 Graphical user interface from Simapse . . . . . . . . . . . . . . . . . . . 61 3.5 Consensus model of the virtual species presence . . . . . . . . . . . . . 64 3.6 Framework of the study and general procedures to use artificial neural networks ................................... 72 3.7 Effects of the number of neurons in the single hidden layer on network structure.................................... 75 3.8 AUC values resulting from different network structures . . . . . . . . . . 76 3.9 Proportion of presences and absences per predictive class for the RVS 77 3.10 Sample size effects on model prediction accuracy . . . . . . . . . . . . 78 3.11 Sample size effects on AUC estimates for the GVS and SVS . . . . . . 79 3.12 Proportion of presences and non-presences for the real species . . . . 80 4.1 Study area with sample points . . . . . . . . . . . . . . . . . . . . . . . 92 4.2 Example of the influence of pollen abundance on the PDF . . . . . . . . 95 2 FCUP Introduction period of 400 ka and 100 ka cycle; 2) obliquity or axial tilt is the angle between the Earth’s rotational axis and the plane formed by its orbit and changes with a periodicity of 41 ka; and 3) axial precession, with a period of 23 ka, is related to the orientation of the axis. The interaction between these orbital components determines the insolation on Earth which is the amount of energy that the planet receives. Climate evidences from the geological record are in close agreement with the predicted pace of climate (Hays et al., 1976; Zachos et al., 2001; Cheddadi et al., 2005). Other responses of climate have been documented and attributed to Earth’s geographical and topological features and phenomena, as the concentration of atmospheric greenhouse gases (Pielke et al., 1998; Zachos et al., 2001; Clark et al., 2012). As seen, these large-scale factors dominate the periodic nature of glacial and interglacial cycles. Instead of a smooth transition from cold to warm conditions, the period after a glaciation event also inherits this oscillating nature with abrupt climate transitions (Alley et al., 2003; Cheddadi et al., 2005; Clark et al., 2012). This has been particularly evident since the last post-glacial process until the present, when a chain Fig. 1.1 – Major orbital perturbations defining long term climate cycles by changing the amount of radiation reaching the planet, as predicted by the Milankovitch theory. Orbital eccentricity refers to the shape of Earth’s obit around the Sun, responsible for long cycles of 400 and 100 ka. The axial obliquity is the angle formed by Earth axis and it orbit and has a medium period of 41 ka. The axial precession refers to the orientation of the axis and controls shorter cycles of 23 and 19 ka. Combinations of these periods are responsible for major insolation changes on Earth. Adapted from Zachos et al. (2001). FCUP 3 Climate evolution of cold and warm events took place (Bond et al., 1993, 1997; Dansgaard et al., 1993; von Grafenstein et al., 1999; Renssen & Isarin, 2001; Clark et al., 2012). During the glacial period, the Dansgaard-Oeschger events (D-O) were responsible for warming trends with a rhythm of 1,470 years that were followed by a longer colder period (Bond et al., 1993; Dansgaard et al., 1993; von Grafenstein et al., 1999). These events are generally followed by oceanic Heinrich events, resulting in fast cold-to-warm transitions with massive discharges of icebergs in the North Atlantic (Bond et al., 1993). The Heinrich event 1 (H1), simultaneous to the cold period of the Oldest Dryas (OD; 18.0 - 14.6 ka), ended with the Bølling-Allerød warm interstadial. The last D-O event includes a period of cold climate between ∼12.9 to 11.6 ka known as the Younger Dryas (YD), and considered as a Heinrich event by some authors (Bond et al., 1993), marking the beginning of the Holocene. The Holocene, although more stable, had also a ∼1470 year cycle responsible for the major shifts to cold during this epoch, including the smaller 8.2 ka event (Bond et al., 1997; von Grafenstein et al., 1999; Heiri et al., 2004). However, this period of time is clearly marked by the Holocene thermal optimum (HTO; ∼11 to 5 ka) with increasing summer temperatures resulting by the combination of orbital forcing and feedback coming from the melting of the Laurentide ice sheet (Seppä & Birks, 2001; Renssen et al., 2009). 1.1.2 European climate Several reconstructions depicted most climate oscillations in Europe, exposing a divergent spatial pattern of climate-related variables (Birks & Ammann, 2000; Lotter et al., 2000; Renssen & Isarin, 2001; Seppä & Birks, 2001; Davis et al., 2003; Cheddadi & Bar-Hen, 2009). At local and regional scales, reconstructions of temperature and precipitation values are usually based on fossil pollen data (e.g. Davis et al., 2003; Cheddadi & Bar-Hen, 2009) or fossil remains of invertebrates (e.g. Tóth et al., 2012). Using fossil pollen, Davis et al. (2003) found different responses in Europe for mean annual temperature: north and central western Europe have a ∼4oC anomaly between 12 ka and the present, while southern western and eastern Europe were generally warmer, with anomalies peaking at -3 and 3oC, respectively. Major differences in Europe were found in winter temperature, being the summer temperatures very stable. Similar results were found by Cheddadi & Bar-Hen (2009), describing a gradient between colder north-western and warmer south-eastern Europe. The succession of events after the LGM had vast repercussions in Europe. The retraction of the ice after the LGM exposed a great area of land (Fig. 1.2), counter- 4 FCUP Introduction balanced by the increase of sea level (Peltier, 1994). The cold events that took place in this period had a wide impact on the summer temperatures. During the OD there was an expansion of the glacial vegetation indicating lower temperatures in southern Europe (Naughton et al., 2007) and colder summers in the east (Tóth et al., 2012). During the YD there was a fast decline to ∼9oC (Birks & Ammann, 2000; Lotter et al., 2000; Renssen & Isarin, 2001) and also precipitation suffered a drastic decrease (Dormoy et al., 2009). The HTO has a strong latitudinal variation. Evidences of warmer and more humid summers were found in northern Europe (Seppä & Birks, 2001; Davis et al., 2003), but the reconstructions for southwestern Europe point toward cooler summers. At 8.2 ka, during a brief temporal period (∼600 years), the temperature suddenly decreased about ∼1oC (Heiri et al., 2004). During the past decades, the scientific community converged efforts on the asFig. 1.2 – Extension of the ice sheet during the last glacial maximum (∼21 ka). The grey area corresponds to the ice cover during the LGM, when the Weichselian ice sheet covered a great extent of northern Europe, and the Laurentide ice sheet covered most of northern North America. Light green depicts the land area exposed during this period due to the drop of sea level. Based on data from Peltier (2004). FCUP 5 Species sensitivity to climate change sessment of other factors shaping global climate, mainly those related to the impact of human activities. Climate changes are driven by changes in the insolation, which, by its turn, is controlled by orbital parameters and Earth features, such as topography, plate tectonics, oceanic circulation, surface albedo and concentration of greenhouse gases (Pielke et al., 1998; Zachos et al., 2001; Clark et al., 2012). The interactions between atmosphere and terrestrial ecosystems are also an important factor modelling Earth’s climate. This relation is largely controlled by feedbacks processes triggered by water evaporation from plant transpiration, carbon cycle and human activities shaping the landscape (Pielke et al., 1998). Larger-scale anthropogenic influence on the European landscape has its roots in the development of agriculture, especially during the late-Quaternary (Berglund, 2003; Jalut et al., 2009). More recently, the escalating concentration of greenhouse gases in the atmosphere has already been reflected in an increase of temperatures with an anthropogenic origin (IPCC, 2007b). The human climate forcing has a predicted increase of severity for the next century with wide physical and biological consequences (IPCC, 2007a; Rosenzweig et al., 2008). In summary, the global climate since the LGM was a succession of cold and warm events. We have seen that some were very drastic and intense, whereas others persisted for longer periods. What is the spatial impact of the velocity of climate change? This question arises frequently due to the wide impact on biodiversity, and answers point towards a simple explanation: complex topographies are more prone to refrain climate change while flat areas have higher velocities (Loarie et al., 2009; Sandel et al., 2011). Which consequences may interactions between climate and topography have on the general patterns of biodiversity? 1.2 Species sensitivity to climate change 1.2.1 Consequences of past climate change The major biological patterns are largely controlled by the climate (Parmesan & Yohe, 2003; Araújo & Rahbek, 2006; Araújo et al., 2008; Willis & MacDonald, 2011) and the velocity of climate change is a major force behind the global patterns of biodiversity. It has already acted in the past, shaping the current distribution of endemic species (Sandel et al., 2011) and it is predicted to be an extremely severe shaping force during the next century, given the scenarios of climate change (Loarie et al., 2009; Burrows et al., 2011). The velocity of climate change since the LGM is noticeably higher at northern latitudes, reaching a peak of ∼170 m/yr on Canada (Sandel et al., 2011). This corresponds to the area that was mostly covered by the Laurentide ice sheet 6 FCUP Introduction during the last glacial period (see Fig. 1.2). Despite the relatively lower velocities, Europe has a south-to-north gradient of increasing velocities (Sandel et al., 2011), highlighting the prominent role of biodiversity refugia of the Mediterranean basin. The milder climate of southern European latitudes allowed several species to persist during the harsh conditions of the last glaciation (Weiss & Ferrand, 2007). Sandel et al. (2011) related high levels of endemism (in the author’s view, small-ranged species of amphibians, mammals and birds) with lower velocities of climate change, which creates an interesting connection to the world hotspots of biodiversity, like the above case of the Mediterranean basin (Myers et al., 2000). The velocity of climate change is known to have a strong impact on biodiversity, however, the mechanisms behind Quaternary extinctions are not fully understood. The interactions between Quaternary climate change and the sprout of human populations are, nevertheless, related to most extinctions during the last glaciations (Burney & Flannery, 2005; Koch & Barnosky, 2006; Wroe et al., 2006; Brook et al., 2008; Nogués-Bravo et al., 2008, 2010; Barnosky et al., 2010). Although there are some discrepancies between authors on which factor is more influential in this interaction (see, for instance, Burney & Flannery, 2005; Wroe et al., 2006), there is a consensus that the combined action of climate warming and human activities probably led many species to extinction. On one side, temperature increase since the LGM restricted the area where some cold-adapted species were able to thrive (Nogués-Bravo et al., 2008, 2010). Additionally, humans, either directly by active hunting or indirectly by inducing major landscape changes and introducing predators, may have caused drastic reductions of population sizes (Burney & Flannery, 2005; Koch & Barnosky, 2006). For instance, through the analysis of potential niche of the woolly mammoth, Mammuthus primigenius, in different time-frames, Nogués-Bravo et al. (2008) suggested that the survival of this species would have been possible in very small patches and, thus, attributed the extinction to the synergistic effects of climate change and human activities. In Eurasia, extinction processes in the Quaternary are responsible for the loss of 36% of the mega-fauna genera (Barnosky et al., 2010) and are divided in two stages: between 46 ka and 22 ka, with the loss of warm-adapted species; and between 12 ka and 8 ka, during the Pleistocene-Holocene transition, with the extinction of cold-adapted species (Barnosky et al., 2010). The impact of climate change in species extinctions is an important factor modulating biodiversity during the Quaternary. Nevertheless, the migration response or resilience of extant species to climate change defines the current biodiversity patterns (Svenning & Skov, 2007; Araújo et al., 2008; Sandel et al., 2011). During the climate FCUP 7 Species sensitivity to climate change oscillations of the Quaternary, species experienced several contractions, expansions and extirpation of their ranges to track suitable conditions (Taberlet et al., 1998; Hewitt, 2000; Davis & Shaw, 2001; Hewitt, 2004; Taberlet & Cheddadi, 2002). In general, migration paths in Europe follow a south-to-north direction, corresponding to the spread of species from the Mediterranean area, which broadly served as refugia, despite the recent ramifications of this term (Ashcroft, 2010). The analogous climate areas between present and the LGM in Europe revealed interesting relations with the broad pattern of refugia areas (Ohlemüller et al., 2012). These authors divided Europe in source and sink areas depending on whether they had currently wide areas of analogous or non-analogous climate. The southern Europe, mostly below 45oof latitude, revealed an extremely high source potential, along with some northern areas in France and Germany. Sink potential was found to the east of 10o, covering most of Eastern Europe. Despite the importance of the latitudinal movement of the species when tracking the changing climate, the altitudinal gradient was equally important (Davis & Shaw, 2001). Evidences of plants migrations during the late Quaternary are abundant from the paleorecord. In Europe, plant assemblages during the LGM exhibited a gradient from southern steppe to northern tundra that evolved mainly to forest-like biomes during the middle Holocene (Elenga et al., 2000; Prentice et al., 2000; Cheddadi & Bar-Hen, 2009; Clark et al., 2012). Herbaceous taxa dominated the pollen record during the LGM and higher densities of tree pollen appear in southern Europe at 13 ka and increased pronouncedly all over Europe until 5 ka (Cheddadi & Bar-Hen, 2009). Migrations routes drawn from fossil records and phylogenetic analysis of species from the genus Pinus and for Fagus sylvatica located several potential refugia in southern Europe (Cheddadi et al., 2006; Magri et al., 2006). Along with the southern European peninsulas, cryptic refugia was also detected for tree species in central and northern Europe (Cheddadi et al., 2006; Magri et al., 2006; Magri, 2010; Svenning et al., 2008). The same pattern of species’ persistence and migration routes is found on animals (Taberlet et al., 1998; Hewitt, 2000, 2004). The paradigmatic cases of the grasshopper (Chorthippus parallelusi), the hedgehog (Erinaceus europeus) and the bear (Ursus arctos) reveal expansion from southern peninsulas, with differences on the extent of expansion from the source in each case (Hewitt, 2000). Others, like small mammals and vipers, revealed cryptic refugia north to the classic Mediterranean area (Ursenbacher et al., 2006; Provan & Bennett, 2008; Fløjgaard et al., 2009; Vega et al., 2010). The oscillations in species’ ranges have consequences for interspecific genetic 8 FCUP Introduction structure (Hewitt, 1999, 2000, 2004; Petit et al., 2003). The harsh conditions during the glacial period that forced temperate species to reduce their ranges to the southern Europe, also promoted vicariance effects, reducing gene flow and increasing the differentiation between lineages (Taberlet et al., 1998). The horizontal colonization pattern is expressed through wider areas than the vertical, which is restricted to mountain ranges and have less available area for expansion processes. The fast expansion of a population creates a gradient of high to low genetic diversity from source areas to the edge of the distribution (Hewitt, 1999, 2000, 2004). Although species respond differently to climate change, most populations suffered strong bottlenecks and founder effects, recovering from a small number of individuals with a rapid expansion (Taberlet et al., 1998; Hewitt, 1999, 2000, 2004). This resulted in a loss of genetic diversity at the leading edge of the range expansion (Hewitt, 1999). Moreover, the establishment of new populations may hamper the introduction of new individuals and, thus, prevent the increase of diversity (Hewitt, 1999). In plants, due to the dispersal mechanisms by seeds and pollen, adaptation was found important both in the edges of expansion and in the core of the population distribution (Davis & Shaw, 2001). As a result of the expansions that took place after the LGM, several suture zones are found in Europe, which correspond to hybrid zones where previously isolated lineages met in secondary contact (Taberlet et al., 1998; Hewitt, 1999, 2000, 2011). Important suture zones in Europe are coincident with major topographic features like the Alps and the Pyrenees. Populations dispersing from the Iberian and the Balkan Peninsulas, met in a generally flat area corresponding to France and Germany. It is interesting to notice that climate change velocity is higher in flat areas (Loarie et al., 2009; Sandel et al., 2011), probably assisting longer dispersal events. A last common hybrid zone in Europe resides in the Scandinavian Peninsula, where species from several refugia met. The exact location of suture zones depends on the interactions between intervening organisms, and tend to float more on flat areas (Taberlet et al., 1998). 1.2.2 Predicted impacts of future changes During the late Quaternary, species responded rapidly to climate oscillations. Will biodiversity be able to cope with future climate change velocity? The predicted climate change for the next century raises several conservation issues (Skov & Svenning, 2004; Thomas et al., 2004; Araújo & Rahbek, 2006; Araújo et al., 2006; IPCC, 2007a; Petit et al., 2008; Loarie et al., 2009; Willis & Bhagwat, 2009; Carvalho et al., 2010a,b; Dawson et al., 2011; Hof et al., 2011; Schloss et al., 2012). The impact of recent climate change (over the last century) is already noticeable on several species (Root FCUP 9 Evolutionary patterns in the Iberian Peninsula et al., 2003; Parmesan, 2006; Moritz et al., 2008; Tingley et al., 2009) and is predicted to result in multiple extinctions (Parmesan, 2006) and, thus, raise the need for conservation measures. There are multiple operators involved on global extinctions besides climate change (Brook et al., 2008; Willis et al., 2010), and lessons from the past, through the study of fossil evidence and climatic reconstructions, show high resilience of biodiversity to past changes, giving a new breath to biodiversity conservation (Willis et al., 2010). However, the authors advert that future changes have a very different nature, given the higher velocity predicted. In fact, Sandel et al. (2011) predicts a difference of more than two orders of magnitude between past (since LGM) and future climate change velocities. Thus, dispersal ability will be one of the most influential processes defining species’ success in the next century (Schloss et al., 2012), as it was in the past. Biodiversity conservation under such scenarios of rapid climate change is difficult but several avenues of actions have been proposed for conservation at a regional level (Carvalho et al., 2010a), either by bearing in mind the balance between cost and effectiveness (Carvalho et al., 2010b) or by incorporating evolutionary processes into conservation design (Klein et al., 2009). 1.3 Evolutionary patterns in the Iberian Peninsula 1.3.1 Characteristics of Iberia The climate dynamics patterns observed in Europe are also seen at lower spatial scales in the Iberian Peninsula. Although it is considered part of a southern European refugium with a high potential of being a source area for species spreading (Ohlemüller et al., 2012), climate was not stable in Iberia. High amplitude climatic oscillations were also frequent in Iberia during the late Quaternary, following the general warming trend from harsh conditions and promoting shifts of the species’ ranges within the Iberian Peninsula. All major climatic events of the Quaternary are discernible in the Iberian margin fossil record with wide impact on land biodiversity (Sánchez-Goñi et al., 2000; Roucoux et al., 2005; Naughton et al., 2007; Fletcher & Sánchez-Goñi, 2008; Fletcher et al., 2010). The Iberian Peninsula has peculiarities that make it a unique area in Europe. Located in the south-western part of Europe, the peninsula is partially isolated from the rest of Europe with a strong geographical barrier of the Pyrenees mountains (Fig. 1.3). Similarly, to the south, it is bathed by the Mediterranean Sea with the relatively small Strait of Gibraltar (∼14km wide) separating it from the African continent (Fig. 1.3). Such isolation features hampers the dispersal of several species, particularly non- 10 FCUP Introduction flying animals or short-dispersal plants, though the strait and the Pyrenees mountains are known to have been permeable to species dispersal several times in the past (e. g. Petit et al., 2003; Carranza et al., 2006; Arroyo et al., 2007). The high level of endemism in Iberia is comparable to those exhibited by the other peninsulas of the Mediterranean Basin and was recognized as part of a global biodiversity hotspot (Myers et al., 2000). The Iberian topographic features include two plateaus that dominate the landscape north and south of the Central Mountain System (Fig. 1.3). Other mountain systems include the northern Cantabric mountains, the southern Baetic system, the Sierra Morena mountains, and the eastern Iberian system (Fig. 1.3). The peninsula exhibits two major bioclimatic zones: the Mediterranean that covers most of the peninsula, except northern areas, and the mountain ranges where Atlantic bioclimate dominates. 1.3.2 Climate change and landscape dynamics During the LGM, the Iberian Peninsula was cold and dry (Naughton et al., 2007; Fletcher et al., 2010), however, these conditions were extreme nearly after the LGM, at the OD (Roucoux et al., 2005; Naughton et al., 2007; Fletcher et al., 2010). During this period, Iberian landscapes were dominated by steppe and tundra vegetation (Sánchez-Goñi et al., 2000; Carrión et al., 2010a), with constrained forests in the northwest and south (Carrión, 2002; Muñoz-Sobrino et al., 2006; Carrión et al., Fig. 1.3 – The Iberian Peninsula and Balearic island with major topographic features. FCUP 11 Evolutionary patterns in the Iberian Peninsula 2010a), though the LGM had enough humidity to support the development of shrub vegetation (Fletcher & Sánchez-Goñi, 2008). The late OD witnessed the expansion of pioneer trees like the case of Betula in the north (Muñoz-Sobrino et al., 2006). The BA is marked by the increasing temperature and precipitation, during which warm and moist conditions were propitious for the advance of forests over steppe landscapes (Naughton et al., 2007; Fletcher & Sánchez-Goñi, 2008; Carrión et al., 2010a). The expansion of pine and oak forests is very intense during this period (Fletcher et al., 2007; Naughton et al., 2007; Carrión et al., 2010a), at least until the onset of the YD, when abrupt drops of temperature and precipitation took place. Temperature reconstructions show that southwest Europe had a severe transition from the YD to Holocene, similarly to northern and central Europe, but much colder than the southeast (Davis et al., 2003). This pattern was confirmed by different temperature reconstruction methods from fossil pollen (Cheddadi & Bar-Hen, 2009). The severe decrease of temperatures on the YD is not noticeable in all Iberian pollen sites, though some reveal extreme sensitivity to the climatic event (Carrión et al., 2010a). Marine pollen records show that the Iberian Peninsula experienced a strong retraction of temperate forests during this period, particularly oak, with the return to steppe landscape dominated by herbaceous plants (Naughton et al., 2007; Fletcher & Sánchez-Goñi, 2008; Carrión et al., 2010a). Nevertheless, Betula trees expanded slightly during this stage (Muñoz-Sobrino et al., 2006; Naughton et al., 2007). The transition to the Holocene has similarities to the BA period. Cold-to-warm changes marked the end of the YD and the beginning of the new period. A striking pattern observed in temperature reconstructions is the absence of the HTO in Iberia, while it was evident in other European areas (Davis et al., 2003; Cheddadi & BarHen, 2009). However, the expansion of Quercus,Corylus and Alnus in northern Iberia (Naughton et al., 2007) and of other thermophile vegetation, including species from the genera Quercus and Olea, in central and southern Iberia, lead the general increase of forest cover (Dorado-Valiño et al., 2002; Fletcher et al., 2007), indicating an rise of temperature and moisture availability. The substitution of well established pine by oak forests is often regarded as a result of climate change and/or a modification of the fire regime (Carrión et al., 2010a). After 5 ka there was a general expansion of green shrublands, with a decrease of Quercus and Pinus forests (Dorado-Valiño et al., 2002; Fletcher et al., 2007; Naughton et al., 2007; Fletcher & Sánchez-Goñi, 2008; Carrión et al., 2010a). The cold events of low amplitude at 10.1, 9.3, 8.2 and 7.4 ka are reflected in forest declines (Fletcher et al., 2010). Despite the impact of climate dynamics in the Iberian landscape, the anthro- 18 FCUP Introduction quantify the impact on climate (e.g. clouds, land heterogeneity, greenhouse gases, ice cover, aerosols, vegetation and others), have been sequentially added to the modelling equations, contributing to the increase of their complexity, more accurate results, and more processing time (IPCC, 2007a). To face the uncertainties of global economic trends and the effects of human activities on greenhouse gases concentrations, the IPCC has built emission scenarios (SRES; IPCC, 2000) with four categories. Such scenarios rely on driving forces of greenhouse emissions, such as demographic evolution and socio-economic features (Fig. 1.5; Table 1.1). The underling storyline behind the A1 scenario family describes a homogeneous world with reduction of major per capita income differences between regions, and intensive economic growth, and an increasing population until mid-century, followed by a decrease. This scenario is further divided in three groups accordingly to major energy source trends: 1) A1T, intensive use of fossil fuels; 2) A2T, use of Fig. 1.5 – Families of emissions scenarios accordingly to IPCC (adapted from IPCC, 2000). Vertical axis represents trends from environmental-enemy to environmental-friendly economies. Horizontal axis denotes gradients from globalised to regionalised worlds. FCUP 19 Reconstructing the past, predicting the future Table 1.1 – Predicted temperatures for 2100 from different emission scenarios based on the IPCC 4th assessment report (IPCC, 2007b). Values given in Celsius degrees (oC ) represent the multi-model global averages of surface warming in relation to average historical climate (1980-1999) and the ±1 standard deviation range of individual model annual averages (between squared brackets). Economic based Environmental based Increasing globalization A1B A1T A1FI B1 2.8 oC 2.2oC 4.0 oC 1.9 oC [1.7, 4.4] [1.4, 3.8] [2.4, 6.4] [1.1, 2.9] Increasing regionalization A2 B2 3.4 oC 2.2 oC [2.0, 5.4] [1.4, 3.8.] non-fossil energy; 3) A1B, a balanced energy source. The B1 scenario family is also described by a storyline with high emphasis on a homogeneous world with population trends as described previously, but including an option for clean energies, less materialistic societies and economies focused on services and information. The A2 and B2 families of scenarios describe an extremely heterogeneous world with continuously growing population, with lower rates in B2. Both scenarios draw regionally controlled economies, but A2 has a greater fragmentation of per capita income and technology change than other storyline. On the other hand, B2 scenario is more environmentally friendly, with local solutions to sustainability, resulting in a more diversified technology than the A2 scenario. The spatial resolution of GCMs has been increasing, following also the same trend of computational power due to more demanding calculations (IPCC, 2007a). The resolution has increased almost 5 times since the first IPCC’s assessment report models during 1990 (∼500 km resolution). Better resolutions allow depicting regional trends, but biodiversity patterns operate at smaller scales, demanding even higher resolutions (e.g. Araújo et al., 2005; Peterson et al., 2007; Waltari et al., 2007; Brito et al., 2009; Carvalho et al., 2010a,b; McKelvey et al., 2011). The increase of the spatial resolution, or downscaling, is a common process for current weather data, with several interpolating algorithms available. The thin plate smoothing splines algorithm with latitude, longitude and altitude as independent variables is known to produce good results with climate data (Hijmans et al., 2005; Kriticos et al., 2011). However, past and future climate estimates have less resolution to depict relationships between climate and independent variables, especially elevation. To overcome this problem, a simple method known as statistical downscaling has been used to produce higher resolution maps from past and future climate layers (Waltari et al., 2007; Tabor & Williams, 2010; McKelvey et al., 2011). Since a linear relationship with elevation is 20 FCUP Introduction difficult to achieve, anomalies to the present are calculated at coarse resolutions and interpolated to a finer spatial scale, usually with splines, and reprocessed to full climate values by adding the present data at the same finer scale (Tabor & Williams, 2010). The anomalies have wider relations with space and, thus, are easier to interpolate. Moreover, the relationship with elevation is usually well established with the available data from weather stations. This method allows building spatial layers with enough resolution to establish spatial relationships with biodiversity distribution, although with an increase of uncertainty (Tabor & Williams, 2010). 1.5 Ecological niche modelling as a tool to study niche dynamics 1.5.1 Niche concept and traits The niche concept is transversal to many subjects discussed above. For instance, building past environmental layers of climate using plant taxa is based on the concept of climate niche and its stability over time. Also, species’ distributions are largely maintained by niche preferences and their expansion and contraction are also niche driven. Additionally, niche can play an important role in evolutionary processes like speciation (Holt, 2003; Wiens & Graham, 2005). Many species’ related processes are based on the concept of niche and such relation raises broad questions about its definition and the number of niche types a definition may hold. The term niche was first applied by Grinnell (1917) to define the set of environmental conditions where a species is able to thrive and reproduce. Later, Elton (1927), in his description of animal communities, defined the niche as the interacting factors that act as forces maintaining the distribution. This is a more functional view of the niche, where species relations with competitors and resources play a more prominent role. These two type of niches, currently known as Grinnellian and Eltonial niches, respectively (Wiens & Graham, 2005; Soberón, 2007; Wiens et al., 2009), operate on different spatial scales due to the nature of the variables acting in each concept (Soberón, 2007). The variables used to construct the Grinnellian niche, the scenopoetic variables (Hutchinson, 1978; Soberón, 2007; Wiens et al., 2009), are measured at large scales (climate data, for instance); whereas the bionomic variables are related to the Eltonian niche (Hutchinson, 1978; Soberón, 2007; Wiens et al., 2009), and due to their biological nature (e.g. competition for resources), operate mostly at local scales (Soberón, 2007). Three decades after the Eltonian niche definition, Hutchinson (1957) summarized the niche concept in the currently most commonly used form: a hypervolume FCUP 21 Ecological niche modelling as a tool to study niche dynamics of n-dimensions where a species is able to persist indefinitely. The dimensions are the niche-related variables meaningful for the species. The Hutchinson niche concept is further divided in the fundamental and realized niches. The first is, in fact, the ndimensional hypervolume that describes the ecological requirements of the species, whereas the second is a section of that hypervolume where the species exists de facto due to other biological pressures such as the presence of competitors and resources availability (Hutchinson, 1957; Guisan & Zimmermann, 2000; Guisan & Thuiller, 2005). These niche concepts are similar but they rarely completely overlap (Fig. 1.6; Soberón & Peterson, 2005). The complete picture of the niche should, however, include a time component. This was later introduced by incorporating the source-sink theory to the niche concept (Pulliam, 2000), and also the range of dispersal ability (Holt, 2003). These concepts add a new interacting sphere to the niche concept representing the area where species may disseminate given their dispersal ability (Fig. 1.6). Fig. 1.6 – Biotic-Abiotic-Movement diagram representing driving forces of species’ distributions. The study area G represents the scenopoetic and bionomic space containing the species’ niche and accessible areas. Aand JF represent the fundamental niche of the species with positive growth rates and composed by Grinnellian variables; B is the space within Eltonian variables where species can compete and coexist with others; and Mdenotes the area accessible by the species. JRis the sum of the area where species is found (JO) plus the area where species has potential to be found ( ˜ JO). The area where species is found, irrespective of growth rate is JSS, and assumed to be similar to M. Adapted from Soberón (2007), using author’s niche notation for simplicity. 22 FCUP Introduction The tendency to maintain the fundamental and/or realized niche over time is called niche conservatism (NC; Wiens & Graham, 2005; Pearman et al., 2008; Peterson, 2011). The conservation of the niche is central to the analysis of biodiversity patterns due to its broad impact on evolution and, consequently, on species’ distributions. For instance, NC is a major force acting in allopatric speciation (Wiens, 2004, 2008; Wiens & Graham, 2005; Pearman et al., 2008; Glor & Warren, 2011). Vicariance events are related to NC over time, by forcing population splits under conditions that are outside the original niche of the population (Wiens, 2004; Wiens & Graham, 2005; Wiens et al., 2009). In its basic definition, a barrier to species dispersal is a patch of unsuitable or suboptimal habitat dividing a population (Wiens & Graham, 2005). Thus, species dispersal is also controlled by NC positively, by allowing individuals to migrate within their niche, and negatively by creating a resistance to the movement in areas outside the species’ niche range. Time is an influential factor of niche stability (Pearman et al., 2008) and evidences of NC come from very different time-scales (Peterson, 2011). On one hand, short to medium term (100to 106years) phenomena including species invasions and geographical shifts orchestrated during the Quaternary, show a trend to preserve the niche of the species (Peterson, 2011). On the other hand, evolutionary scale events between sister species or distant species occur at higher time-scales and involve phylogenetic hypotheses to be tested (Peterson, 2011). Several studies have demonstrated the existence of labile niches while others have found strong evidence of NC between species (for a review on NC see Pearman et al., 2008). Additionally, cryptic NC was found between evolutionary lineages of lizards with stable fundamental niches hidden within different realized niches (Schulte et al., 2012). Detecting NC may be also hampered by the choice of niche predictors used in each study due to different degrees of stability for different variables (Rödder & Lötters, 2009). Nevertheless, NC is an assumption for ecological niche modelling and should be addressed especially with model transference to other temporal and spatial scales (Wiens & Graham, 2005; Wiens et al., 2009; Petitpierre et al., 2012; Warren, 2012). 1.5.2 Ecological niche-based modelling Ecological niche-based modelling (ENM) is the exercise of capturing a species’ niche based on the locations of presence or presence and absence data and a set of environmental variables. It is based on correlative approaches to predict the distribution of a species and it is often compared to a process-based approach (Guisan & Zimmermann, 2000; Kearney, 2006; Soberón & Peterson, 2005), where the relation between FCUP 23 Ecological niche modelling as a tool to study niche dynamics the species and environmental factors is determined mechanistically to find the range where fitness is maximised. Kearney (2006) argued that only mechanistic models can help defining the niche of the species, whereas ENMs should use the term ’distribution model’ preferably to ’niche model’. However, the purpose of ENMs is in fact the assessment of the niche (Warren, 2012). Nevertheless, the presence data of a species is restricted to the realized niche, and there is a lack of consensus in the scientific community about which niche, fundamental or realized, is the result of an ENM (Araújo & Guisan, 2006). A simple answer to this problematic may reside in the complexity of the models being built (Sillero, 2011), that is related to the chosen algorithm (from climate envelope based to machine learning) and to the amount of ENVs used. Simple models, particularly those with few predictors, tend to capture the part of the niche that is conserved (Peterson, 2011) which corresponds to the fundamental niche, because conservatism applies to the Grinnelian space (Soberón, 2007). On the other hand, very complex models tend to over-fit the realized niche. However, there is no estimate available on the correct number of variables to be used in ENM (Peterson, 2011). Ecological niche modelling is a modern tool and its usage is transverse to research fields of landscape ecology and genetics. It is a major tool in conservation planning (Brito et al., 2009; Torres et al., 2010), especially under scenarios of climate change (Araújo et al., 2006; Carvalho et al., 2010a,b), but also fundamental in studies of niche evolution (Peterson et al., 1999; Rödder & Lötters, 2009; Tingley et al., 2009; McCormack et al., 2010; Petitpierre et al., 2012; Schulte et al., 2012). Finding and quantifying relationships between species presence data and environmental data requires an algorithm capable of discovering patterns. Several algorithms are available to this purpose, ranging from very simplistic envelope models to regressionbased analyses, and the more complex machine-learning techniques (Guisan & Zimmermann, 2000; Elith et al., 2006). Algorithms for modelling are grouped based on whether they use only species presence data or they require data on localities of both species presence and absence. Due to the ambiguous nature of absences (Lobo et al., 2010), these are often substituted by artificially generated data (pseudoabsences; Zaniewski et al., 2002; Pearce & Boyce, 2005). Envelope models require presence-only data to build the hypervolume niche based on the locations where the species is present (Carpenter et al., 1993; Elith et al., 2006; Elith & Leathwick, 2009). Niche estimations with presence-only data can be also achieved using a modified principal component analysis, a technique implemented in the Ecological-Niche Factor Analysis (ENFA; Hirzel et al., 2002). Regression-based methods are extensively used to model presence and ab- 24 FCUP Introduction sence data. The aim of regression models is to find relationships between one or more predictors and the presence and absence data. An important advantage over the envelope methods is the quantification of predictor importance. These methods include linear regression, which tries to find linear solutions with maximum likelihood explaining the species data, but may generate, in turn, unreal values (negative probabilities or higher than 1.0; Guisan & Zimmermann, 2000). More commonly used is the extension of linear models to other possible statistical distributions found in the generalized linear models by means of a link function (Guisan & Zimmermann, 2000). The logit function is widely applied link function in the literature to model data that approximate a binomial distribution, as is the case of binary species presence/absence data. Generalized additive models are extensions of these models to fit non-parametric data and better describe non-linear relationships between the environmental predictors and species data (Guisan & Zimmermann, 2000; Guisan et al., 2006; Elith & Leathwick, 2009). With the increasing availability of larger datasets of both species occurrence and environmental data, the relationships between species and environment become more complex to model and demand efficient and powerful techniques. Machine learning methods are very efficient at pattern-detection, but also, highly demanding in terms of computer processing. Several algorithms are described for the purpose of finding species distribution’s patterns, including genetic algorithms, classification and regression trees, maximum entropy, and artificial neural networks. The latter has been proved efficient in detecting very complex patterns that are usually found in ecological relations. It offers a multitude of advantages, however, it application to ecological modelling is still lacking dedicated tools as other commonly used algorithms. 1.5.3 Artificial neural networks Artificial neural networks (ANN) have several advantages over other approaches. First, it is a simple method with broad literature available. The extensive application in ecological research, perhaps inspired by the biological motivation behind the artificial neuron, generated a diverse and broad literature explaining the ANNs in great detail, resulting in an easy to understand methodology. Secondly, ANNs allow different settings of the learning parameters and reshaping of the network (explained in detail below), creating a vast array of learning options that the user may explore within an ensemble forecasting framework (Araújo & New, 2007; Olden et al., 2008). Thirdly, the train and test error assessment allow reducing model overfitting and thus control the generalization ability of the trained network (Dimopoulos et al., 1999; Özesmi et al., FCUP 25 Ecological niche modelling as a tool to study niche dynamics 2006b). Fourthly, ANNs accept continuous or discrete data, both for dependent and independent variables, thus, models may be made using binary or continuous presence data (e.g. abundance) and benefiting from the panoply of possible ecological niche variables available. Other features usually attributed to ANNs include the ability to process complex, non-linear relations and non-parametric data, found regularly in ecological datasets, with a high prediction success even in the presence of noise (Lek et al., 1996; Dimopoulos et al., 1999; Lek & Guégan, 1999; Spitz & Lek, 1999; Pearson et al., 2002; Özesmi et al., 2006a,b; Olden et al., 2008). The development of ANNs was motivated by the efficient pattern interpretation capabilities of the brain and the introduction of the mathematical model of a simplified neuron by McCulloch & Pitts (1943). During the last decade, a relevant number of scientific papers reviewed the application of ANNs in ecology and the analysis of ANNs’ internal behaviour (Olden & Jackson, 2002; Gevrey et al., 2003, 2006a,b; Özesmi et al., 2006a; Park & Chon, 2007), revealing a growing interest on ANNs, despite the restricted use by a small circle of ecologists with a background in informatics (Olden et al., 2008). The interest of the biological scientific community in ANNs, with pioneering applications using species’ data during the 1990s (Lek et al., 1996; Tan & Smeins, 1996; Mastrorillo et al., 1997), triggered new developments of the method in an ecological context, especially with sensitivity analyses of the model output (Lek & Guégan, 1999; Olden & Jackson, 2002; Gevrey et al., 2006b; Özesmi et al., 2006a). The ecological application of ANNs is still widening, encompassing multiple research fields such as predictive modelling (Özesmi & Özesmi, 1999; Özesmi et al., 2006b), plant ecology (Hilbert & Muyzenberg, 1999), impact assessment (Spitz & Lek, 1999), and potential effects of climate warming (Pearson et al., 2002; Araújo et al., 2006; Xavier et al., 2010), supporting the extreme potential of this method. With the increase of computational power, the demand for spatially explicit predictions has also increased. The ANNs spectrum is constituted of various types of networks, varying in structure and learning methods. One common type of ANN is the feed-forward neural network with back-propagation learning (BPN; Lek & Guégan, 1999). This network has a layered structure of neurons connecting the inputs to an output through one or several hidden layers (Fig. 1.7). The BPN has been used to model ecological systems due to its efficient learning ability with complex relationships between dependent and independent variables and also due to its simple nature, which makes it easy to understand (Lek & Guégan, 1999; Özesmi et al., 2006a). The artificial neuron constitutes 26 FCUP Introduction the unit of the BPN. The output of a neuron with n inputs is defined as: (1) Output =f(net) = f( n X i=1 wixi) where f(net)represents the activation function and wirepresents the weight connecting to input xi. The activation function, usually a linear or a sigmoid function (2), squashes the sum of the products of all connecting weights and the respective neuron output in the previous layer to a value that is passed to the next layer of neurons by the connecting weights. (2) f(x) = 1 1 + e−x The learning of the network occurs during the training stage where vectors of variables for each target (presence or absence in the ecological data) are presented to the network, in a process called supervised learning. The values of the connecting weights are initialized to random values and the input vectors of targets and respective variables vectors are propagated through the network producing an output. The error of the output in relation to the target value of the input (i.e. the value that corresponds to the set of inputs of the training data) is assessed and back-propagated in the network. Usually, the error processed by the network is half the sum of squared errors given by the formula: (3) E=1 2X d∈D (td−od)2 where D is the training set, tdthe target value and odthe output of the network for the d sample in the training data. The learning in BPN occurs by gradient descent when the weight’s change is updated accordingly to a set of learning rules. In a BPN, the learning rules are the back-propagation algorithm: (4) 4wij =−ησE σwij +α4wt−1 ij where wij is the weight connecting a neuron in layer i to neuron in layer j, ηrepresents the learning rate and αthe momentum. The last two parameters are commonly used to control the learning ability. The learning rate affects the convergence of the network by defining the allowed quantity of weight’s change (Özesmi & Özesmi, 1999; Olden & Jackson, 2002). While small learning rates will slowly converge to a global minimum, high values will oscillate in the error surface, and may fail in finding the global minimum. This parameter highly interacts with momentum which defines how much of the last weight change is incorporated in the new weight change, thus giving a direction of progress on the error surface (Lek et al., 1996; Spitz & Lek, 1999; Olden & Jackson, 2002). The process of feeding the network with training data is repeated to minimise FCUP 27 Ecological niche modelling as a tool to study niche dynamics the error until the network is fully trained (Lek & Guégan, 1999). Once a network is trained with the combination of input data, independent variables and targets, it can be used to predict to other locations where the same set of variables is available. 1.5.4 Ensemble modelling and model evaluation The use of different models allows a better assessment of uncertainty due to the variability in the outputs (Pearson et al., 2002; Araújo & New, 2007). Stacking several models allows defining a consensus prediction by averaging and to compute a standard deviation as a measure of prediction uncertainty. The ensembling can be done at several levels: across replicates of the same algorithm with resampling methods or across algorithms. This method reduces the over fitting trend of very complex models and increases accuracy (Breiman, 1996). Assessment of model fitting is done by calculating the rates of correct classification. Most common metrics of algorithm performance (e.g. the receiver operating characteristic curve and its area) are based on the confusion matrix where real values are Fig. 1.7 – Typical structure of a simple artificial neural network with backpropagation algorithm and one hidden layer. The artificial neural network is a structure of layers with neurons (dark grey circles) connected by weights (black arrows). The most basic structure has three layers: the input layer where values from different variables are fed to the network; the hidden layer where neurons sum the information from the previous layer and apply an activation function, sending the result to the neurons on the next layer; and the output layer where the results of the network are returned to the user. The first stage of training the network is to propagate the input data to the output neuron (black arrows) and calculate the error at the output. The second stage is to backpropagate the error to each neuron in the previous hidden layers (red arrow). The third stage consists in adapt the weights to new values based on the calculated error with a learning algorithm (green arrow). After this process, an artificial neural network is considered trained and is able to predict to different combination of values of the same set of variables used in input. 34 FCUP Introduction J., Williams, S., Wisz, M.S. & Zimmermann, N.E. (2006) Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29, 129–151. Elith, J. & Leathwick, J.R. (2009) Species Distribution Models: Ecological Explanation and Prediction Across Space and Time. Annual Review of Ecology, Evolution, and Systematics 40, 677–697. Elton, C.S. (1927) Animal ecology. The Macmillan Company, London. Fletcher, W.J., Boski, T. & Moura, D. (2007) Palynological evidence for environmental and climatic change in the lower Guadiana valley, Portugal, during the last 13 000 years. The Holocene 17, 481–494. Fletcher, W.J., Sanchez-Goñi, M.F., Peyron, O. & Dormoy, I. (2010) Abrupt climate changes of the last deglaciation detected in a Western Mediterranean forest record. Climate of the Past 6, 245–264. Fletcher, W.J. & Sánchez-Goñi, M.F. (2008) Orbitaland sub-orbital-scale climate impacts on vegetation of the western Mediterranean basin over the last 48,000 yr. Quaternary Research 70, 451–464. Fløjgaard, C., Normand, S., Skov, F. & Svenning, J.C. (2009) Ice age distributions of European small mammals: insights from species distribution modelling. Journal of Biogeography 36, 1152–1163. Garrigues, T., Dauga, C., Ferquel, E., Choumet, V. & Failloux, A.B. (2005) Molecular phylogeny of Vipera Laurenti, 1768 and the related genera Macrovipera (Reuss, 1927) and Daboia (Gray, 1842), with comments about neurotoxic Vipera aspis aspis populations. Molecular Phylogenetics and Evolution 35, 35–47. Geraldes, A., Carneiro, M., Delibes-Mateos, M., Villafuerte, R., Nachman, M.W. & Ferrand, N. (2008) Reduced introgression of the Y chromosome between subspecies of the European rabbit (Oryctolagus cuniculus) in the Iberian Peninsula. Molecular Ecology 17, 4489–4499. Gevrey, M., Dimopoulos, I. & Lek, S. (2003) Review and comparison of methods to study the contribution of variables in artificial neural network models. Ecological Modelling 160, 249–264. Gevrey, M., Dimopoulos, I. & Lek, S. (2006a) Two-way interaction of input variables in the sensitivity analysis of neural network models. Ecological Modelling 195, 43–50. FCUP 35 References Gevrey, M., Lek, S. & Oberdorff, T. (2006b) Utility of sensitivity analysis by artificial neural network models to study patterns of endemic fish species. Ecological Informatics: Scope, techniques and applications (ed. F. Recknagel), chap. 14, pp. 293–306, Springer, Berlin, 2nd edn. Glor, R.E. & Warren, D. (2011) Testing ecological explanations for biogeographic boundaries. Evolution 65, 673–683. Godinho, R., Mendonça, B., Crespo, E.G. & Ferrand, N. (2006) Genealogy of the nuclear beta-fibrinogen locus in a highly structured lizard species: comparison with mtDNA and evidence for intragenic recombination in the hybrid zone. Heredity 96, 454–463. Gómez, A. & Lunt, D. (2007) Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography of southern European refugia (eds. S. Weiss & N. Ferrand), pp. 155–188, Springer Netherlands. Grinnell, J. (1917) The niche-relationships of the California Thrasher. The Auk 34, 427–433. Guisan, A., Lehmann, A., Ferrier, S., Austin, M., Overton, J.M.C., Aspinall, R. & Hastie, T. (2006) Making better biogeographical predictions of species’ distributions. Journal of Applied Ecology 43, 386–392. Guisan, A. & Thuiller, W. (2005) Predicting species distribution: offering more than simple habitat models. Ecology Letters 8, 993–1009. Guisan, A. & Zimmermann, N.E. (2000) Predictive habitat distribution models in ecology. Ecological Modelling 135, 147 – 186. Hays, J., Imbrie, J. & Shackleton, N. (1976) Variations in the Earth’s orbit: pacemaker of the ice ages. Science 194, 1121–1132. Heiri, O., Tinner, W. & Lotter, A.F. (2004) Evidence for cooler European summers during periods of changing meltwater flux to the North Atlantic. Proceedings of the National Academy of Sciences of the United States of America 101, 15285–15288. Hewitt, G. (2000) The genetic legacy of the Quaternary ice ages. Nature 405, 907–13. Hewitt, G.M. (1999) Post-glacial re-colonization of European biota. Biological Journal of the Linnean Society 68, 87–112. 36 FCUP Introduction Hewitt, G.M. (2004) Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society of London. Series B, Biological sciences 359, 183–195. Hewitt, G.M. (2011) Quaternary phylogeography: the roots of hybrid zones. Genetica 139, 617–638. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G. & Jarvis, A. (2005) Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25, 1965–1978. Hilbert, D.W. & Muyzenberg, J.V.D. (1999) Using an artificial neural network to characterize the relative suitability of environments for forest types in a complex tropical vegetation mosaic. Diversity and Distributions 5, 263–274. Hirzel, A., Hausser, J., Chessel, D. & Perrin, N. (2002) Ecological-niche factor analysis: how to compute habitat-suitability maps without absence data? Ecology 83, 2027– 2036. Hof, C., Levinsky, I., Araújo, M.B. & Rahbek, C. (2011) Rethinking species’ ability to cope with rapid climate change. Global Change Biology 17, 2987–2990. Holt, R. (2003) On the evolutionary ecology of species’ ranges. Evolutionary Ecology Research 5, 159–178. Hutchinson, G. (1957) Concluding remarks. Cold Spring Harbor Symposia on Quantitative Biology, vol. 42, pp. 415–427. Hutchinson, G. (1978) An Introduction to Population Ecology. Yale University Press, New Haven, 1st edn. IPCC (2000) Special report on emissions scenarios: a special report of Working Group III of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK. IPCC (2007a) Climate Change 2007: Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK. IPCC (2007b) Climate Change 2007: The physical science basis. Contribution of working group I to the fourth assessment report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK. FCUP 37 References Jaarola, M. & Searle, J.B. (2004) A highly divergent mitochondrial DNA lineage of Microtus agrestis in southern Europe. Heredity 92, 228–234. Jalut, G., Dedoubat, J.J., Fontugne, M. & Otto, T. (2009) Holocene circumMediterranean vegetation changes: Climate forcing and human impact. Quaternary International 200, 4–18. Kearney, M. (2006) Habitat, environment and niche: what are we modelling? Oikos 115, 186–191. Klein, C., Wilson, K., Watts, M., Stein, J., Berry, S., Carwardine, J., Smith, M.S., Mackey, B. & Possingham, H. (2009) Incorporating ecological and evolutionary processes into continental-scale conservation planning. Ecological applications : a publication of the Ecological Society of America 19, 206–217. Koch, P.L. & Barnosky, A.D. (2006) Late Quaternary Extinctions : State of the Debate. Annual Review of Ecology, Evolution, and Systematics 37, 215–252. Kriticos, D.J., Webber, B.L., Leriche, A., Ota, N., Macadam, I., Bathols, J. & Scott, J.K. (2011) CliMond: global high-resolution historical and future scenario climate surfaces for bioclimatic modelling. Methods in Ecology and Evolution 3, 53–64. Kühl, N., Gebhardt, C., Litt, T. & Hense, A. (2002) Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction. Quaternary Research 58, 381–392. Lek, S., Delacoste, M., Baran, P., Dimopoulos, I., Lauga, J. & Aulagnier, S. (1996) Application of neural networks to modelling nonlinear relationships in ecology. Ecological Modelling 90, 39–52. Lek, S. & Guégan, J.F. (1999) Artificial neural networks as a tool in ecological modelling, an introduction. Ecological Modelling 120, 65–73. Loarie, S.R., Duffy, P.B., Hamilton, H., Asner, G.P., Field, C.B. & Ackerly, D.D. (2009) The velocity of climate change. Nature 462, 1052–1055. Lobo, J.M., Jiménez-Valverde, A. & Hortal, J. (2010) The uncertain nature of absences and their importance in species distribution modelling. Ecography 33, 103–114. Lobo, J.M., Jiménez-Valverde, A. & Real, R. (2007) AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography 17, 145–151. 38 FCUP Introduction López de Heredia, U., Carrión, J.S., Jiménez, P., Collada, C. & Gil, L. (2007) Molecular and palaeoecological evidence for multiple glacial refugia for evergreen oaks on the Iberian Peninsula. Journal of Biogeography 34, 1505–1517. Lotter, A., Birks, H., Eicher, U., Hofmann, W., Schwander, J. & Wick, L. (2000) Younger Dryas and Allerød summer temperatures at Gerzensee (Switzerland) inferred from fossil pollen and cladoceran assemblages. Palaeogeography, Palaeoclimatology, Palaeoecology 159, 349–361. Macarthur, D. (2012) Face up to false positives. Nature 487, 427–428. Magri, D. (2010) Persistence of tree taxa in Europe and Quaternary climate changes. Quaternary International 219, 145–151. Magri, D., Vendramin, G.G., Comps, B., Dupanloup, I., Geburek, T., Gömöry, D., Latałowa, M., Litt, T., Paule, L., Roure, J.M., Tantau, I., van der Knaap, W.O., Petit, R.J. & de Beaulieu, J.L. (2006) A new scenario for the quaternary history of European beech populations: palaeobotanical evidence and genetic consequences. The New Phytologist 171, 199–221. Manabe, S. & Bryan, K. (1969) Climate calculations with a combined oceanatmosphere model. Journal of the Atmospheric Sciences 26, 786–789. Martínez-Freiría, F., Santos, X., Pleguezuelos, J.M., Lizana, M. & Brito, J.C. (2009) Geographical patterns of morphological variation and environmental correlates in contact zones: a multi-scale approach using two Mediterranean vipers (Serpentes). Journal of Zoological Systematics and Evolutionary Research 47, 357–367. Martínez-Freiría, F., Sillero, N., Lizana, M. & Brito, J.C. (2008) GIS-based niche models identify environmental correlates sustaining a contact zone between three species of European vipers. Diversity and Distributions 14, 452–461. Martínez-Freiría, F., Brito, J.C. & Avia, M.L. (2006) Intermediate forms and syntopy among vipers (Vipera aspis and V. latastei) in Northern Iberian Peninsula. Herpetological Bulletin 97, 14–18. Mastrorillo, S., Lek, S., Dauba, F. & Belaud, A. (1997) The use of artificial neural networks to predict the presence of small-bodied fish in a river. Freshwater Biology 38, 237–246. FCUP 39 References McCormack, J.E., Zellmer, A.J. & Knowles, L.L. (2010) Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation? Insights from tests with niche models. Evolution 64, 1231–1244. McCulloch, W. & Pitts, W. (1943) A logical calculus of the ideas immanent in nervous activity. Bulletin of Mathematical Biology 5, 115–133. McKelvey, K.S., Copeland, J.P., Schwartz, M.K., Littell, J.S., Aubry, K.B., Squires, J.R., Parks, S.A., Elsner, M.M. & Mauger, G.S. (2011) Climate change predicted to shift wolverine distributions, connectivity, and dispersal corridors. Ecological Applications 21, 2882–2897. Médail, F. & Diadema, K. (2009) Glacial refugia influence plant diversity patterns in the Mediterranean Basin. Journal of Biogeography 36, 1333–1345. Melo-Ferreira, J., Boursot, P., Randi, E., A. Kryukov, Suchentrunk, F., Ferrand, N. & Alves, P.C. (2007) The rise and fall of the mountain hare. Molecular Ecology 16, 605– 618. Miraldo, A., Hewitt, G.M., Paulo, O.S. & Emerson, B.C. (2011) Phylogeography and demographic history of Lacerta lepida in the Iberian Peninsula: multiple refugia, range expansions and secondary contact zones. BMC Evolutionary Biology 11, 170. Moritz, C., Patton, J.L., Conroy, C.J., Parra, J.L., White, G.C. & Beissinger, S.R. (2008) Impact of a century of climate change on small-mammal communities in Yosemite National Park, USA. Science 322, 261–264. Muñoz, J., Gómez, A., Green, A.J., Figuerola, J., Amat, F. & Rico, C. (2008) Phylogeography and local endemism of the native Mediterranean brine shrimp Artemia salina (Branchiopoda: Anostraca). Molecular ecology 17, 3160–77. Muñoz-Sobrino, C., Ramil-Rego, P. & Gómez-Orellana, L. (2006) Late Würm and early Holocene in the mountains of northwest Iberia: biostratigraphy, chronology and tree colonization. Vegetation History and Archaeobotany 16, 223–240. Myers, N., Mittermeier, R., Mittermeier, C., Da Fonseca, G. & Kent, J. (2000) Biodiversity hotspots for conservation priorities. Nature 403, 853–858. Naughton, F., Sanchez-Goñi, M., Desprat, S., Turon, J.L., Duprat, J., Malaizé, B., Joli, C., Cortijo, E., Drago, T. & Freitas, M. (2007) Present-day and past (last 25000 years) marine pollen signal off western Iberia. Marine Micropaleontology 62, 91–114. 40 FCUP Introduction Nogués-Bravo, D., Rodríguez, J., Hortal, J., Hortal, J., Batra, P. & Araújo, M.B. (2008) Climate change, humans, and the extinction of the woolly mammoth. PLoS Biology 6. Nogués-Bravo, D., Ohlemüller, R., Batra, P. & Araújo, M.B. (2010) Climate predictors of late quaternary extinctions. Evolution 64, 2442–2449. Ohlemüller, R., Huntley, B., Normand, S. & Svenning, J.C. (2012) Potential source and sink locations for climate-driven species range shifts in Europe since the Last Glacial Maximum. Global Ecology and Biogeography 21, 152–163. Olalde, M., Herrán, A., Espinel, S. & Goicoechea, P.G. (2002) White oaks phylogeography in the Iberian Peninsula. Forest Ecology and Management 156, 89–102. Olden, J.D. & Jackson, D.A. (2002) Illuminating the “black box”: a randomization approach for understanding variable contributions in artificial neural networks. Ecological Modelling 154, 135–150. Olden, J.D., Lawler, J.J. & Poff, N.L. (2008) Machine learning methods without tears: a primer for ecologists. The Quarterly review of biology 83, 171–193. Özesmi, S. & Özesmi, U. (1999) An artificial neural network approach to spatial habitat modelling with interspecific interaction. Ecological Modelling 116, 15–31. Özesmi, S., Tan, C. & Özesmi, U. (2006a) Methodological issues in building, training, and testing artificial neural networks in ecological applications. Ecological Modelling 195, 83–93. Özesmi, U., Tan, C., Özesmi, S. & Robertson, R. (2006b) Generalizability of artificial neural network models in ecological applications: Predicting nest occurrence and breeding success of the red-winged blackbird Agelaius phoeniceus.Ecological Modelling 195, 94–104. Park, Y. & Chon, T. (2007) Biologically-inspired machine learning implemented to ecological informatics. Ecological Modelling 203, 1–7. Parmesan, C. & Yohe, G. (2003) A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42. Parmesan, C. (2006) Ecological and Evolutionary Responses to Recent Climate Change. Annual Review of Ecology, Evolution, and Systematics 37, 637–669. FCUP 41 References Pearce, J. & Boyce, M. (2005) Modelling distribution and abundance with presenceonly data. Journal of Applied Ecology 43, 405–412. Pearman, P.B., Guisan, A., Broennimann, O. & Randin, C.F. (2008) Niche dynamics in space and time. Trends in Ecology & Evolution 23, 149–158. Pearson, R.G., Dawson, T.P., Berry, P.M. & Harrison, P.A. (2002) SPECIES : A Spatial Evaluation of Climate Impact on the Envelope of Species. Ecological Modelling 154, 289– 300. Peltier, W. (1994) Ice age paleotopography. Science 265, 195–201. Peterson, A.T., Soberón, J. & Sánchez-Cordero, V. (1999) Conservatism of Ecological Niches in Evolutionary Time. Science 285, 1265–1267. Peterson, A.T. (2011) Ecological niche conservatism: a time-structured review of evidence. Journal of Biogeography 38, 817–827. Peterson, A.T., Pape¸s, M. & Eaton, M. (2007) Transferability and model evaluation in ecological niche modeling: a comparison of GARP and Maxent. Ecography 30, 550–560. Petit, R.J., Aguinagalde, I., de Beaulieu, J.L., Bittkau, C., Brewer, S., Cheddadi, R., Ennos, R., Fineschi, S., Grivet, D., Lascoux, M., Mohanty, A., Müller-Starck, G., Demesure-Musch, B., Palmé, A., Martín, J.P., Rendell, S. & Vendramin, G.G. (2003) Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300, 1563–1565. Petit, R.J., Hu, F.S. & Dick, C.W. (2008) Forests of the past: a window to future changes. Science 320, 1450–2. Petitpierre, B., Kueffer, C., Broennimann, O., Randin, C., Daehler, C. & Guisan, a. (2012) Climatic Niche Shifts Are Rare Among Terrestrial Plant Invaders. Science 335, 1344–1348. Pielke, R.A., Avissar, R., Raupach, M., Dolman, A.J., Zeng, X. & Denning, A.S. (1998) Interactions between the atmosphere and terrestrial ecosystems: influence on weather and climate. Global Change Biology 4, 461–475. Prentice, C., Guiot, J., Huntley, B., Jolly, D. & Cheddadi, R. (1996) Reconstructing biomes from palaeoecological data: a general method and its application to European pollen data at 0 and 6 ka. Climate Dynamics 12, 185–194. 42 FCUP Introduction Prentice, I., Jolly, D. & BIOME 6000 Participants (2000) Mid-Holocene and glacialmaximum vegetation geography of the northern continents and Africa. Journal of Biogeography 27, 507–519. Provan, J. & Bennett, K.D. (2008) Phylogeographic insights into cryptic glacial refugia. Trends in Ecology & Evolution 23, 564–71. Pulliam, H. (2000) On the relationship between niche and distribution. Ecology Letters 3, 349–361. Rebelo, H., Froufe, E., Brito, J.C., Russo, D., Cistrone, L., Ferrand, N. & Jones, G. (2012) Postglacial colonization of Europe by the barbastelle bat: agreement between molecular data and past predictive modelling. Molecular ecology 21, 2761– 74. Reimer, P.J., Baillie, M.G.L., Bard, E., Bayliss, A., Beck, J.W., Blackwell, P.G., Ramsey, C.B., Buck, C.E., Burr, G.S., Edwards, R.L., Friedrich, M., Grootes, P.M., Guilderson, T.P., Hajdas, I., Heaton, T.J., Hogg, A.G., Hughen, K.A., Kaiser, K.F., Kromer, B., McCormac, F.G., Manning, S.W., Reimer, R.W., Richards, D.A., Southon, J.R., Talamo, Turney, C.S.M., van der Plicht, J. & Weyhenmeyer, C.E. (2009) IntCal09 and Marine09 radiocarbon age calibration curves, 0–50,000 years cal BP. Radiocarbon 51, 1111–1150. Renssen, H. & Isarin, R.F.B. (2001) The two major warming phases of the last deglaciation at 14.7 and 11.5 ka cal BP in Europe: climate reconstructions and AGCM experiments. Global and Planetary Change 30, 117–153. Renssen, H., Seppä, H., Heiri, O., Roche, D.M., Goosse, H. & Fichefet, T. (2009) The spatial and temporal complexity of the Holocene thermal maximum. Nature Geoscience 2, 411–414. Rödder, D. & Lötters, S. (2009) Niche shift versus niche conservatism? Climatic characteristics of the native and invasive ranges of the Mediterranean house gecko (Hemidactylus turcicus). Global Ecology and Biogeography 18, 674–687. Root, T.L., Price, J.T., Hall, K.R., Schneider, S.H., Rosenzweig, C. & Pounds, J.A. (2003) Fingerprints of global warming on wild animals and plants. Nature 421, 57– 60. Rosenzweig, C., Karoly, D., Vicarelli, M., Neofotis, P., Wu, Q., Casassa, G., Menzel, A., Root, T.L., Estrella, N., Seguin, B., Tryjanowski, P., Liu, C., Rawlins, S. & Imeson, A. FCUP 43 References (2008) Attributing physical and biological impacts to anthropogenic climate change. Nature 453, 353–357. Roucoux, K., Abreu, L.D., Shackleton, N.J. & Tzedakis, P.C. (2005) The response of NW Iberian vegetation to North Atlantic climate oscillations during the last 65kyr. Quaternary Science Reviews 24, 1637–1653. Sánchez-Goñi, M.F., Turon, J.L., Eynaud, F. & Gendreau, S. (2000) European Climatic Response to Millennial-Scale Changes in the Atmosphere–Ocean System during the Last Glacial Period. Quaternary Research 54, 394–403. Sandel, B., Arge, L., Dalsgaard, B., Davies, R.G., Gaston, K.J., Sutherland, W.J. & Svenning, J.C. (2011) The influence of Late Quaternary climate-change velocity on species endemism. Science 334, 660–664. Schloss, C.A., Nunez, T.A. & Lawler, J.J. (2012) Dispersal will limit ability of mammals to track climate change in the Western Hemisphere. Proceedings of the National Academy of Sciences of the United States of America 109, 8606–8611. Schulte, U., Hochkirch, A., Lötters, S., Rödder, D., Schweiger, S., Weimann, T. & Veith, M. (2012) Cryptic niche conservatism among evolutionary lineages of an invasive lizard. Global Ecology and Biogeography 21, 198–211. Seppä, H. & Birks, H. (2001) July mean temperature and annual precipitation trends during the Holocene in the Fennoscandian tree-line area: pollen-based climate reconstructions. The Holocene 11, 527–539. Sequeira, F., Alexandrino, J., Rocha, S., Arntzen, J..W. & Ferrand, N. (2005) Genetic exchange across a hybrid zone within the Iberian endemic golden-striped salamander, Chioglossa lusitanica.Molecular ecology 14, 245–254. Sillero, N. (2011) What does ecological modelling model? A proposed classification of ecological niche models based on their underlying methods. Ecological Modelling 222, 1343–1346. Skov, F. & Svenning, J. (2004) Potential impact of climatic change on the distribution of forest herbs in Europe. Ecography 27, 366–380. Soberón, J. & Peterson, A.T. (2005) Interpretation of models of fundamental ecological niches and species’ distributional areas. Biodiversity Informatics 2, 1–10. Soberón, J. (2007) Grinnellian and Eltonian niches and geographic distributions of species. Ecology letters 10, 1115–1123. 50 FCUP Objectives measuring species compositional change. The comparative analysis of past and future changes in species richness provides new insights to the spatial location of glacial refugia for amphibians and reptiles, and identifies areas facing potential dramatic impacts of future climate change. In the second manuscript of this chapter we analyse evolutionary processes occurring at a micro-scale in the Iberian Peninsula. A contact zone between three viper species is analysed using Simapse (chapter 3) to assess the role of environmental factors in shaping the dynamics and genetic structure of the populations. We address three main questions: What is the genetic structure of the focal taxa?; Is such genetic structure associated to ecological segregation between taxa?; and Is the ecological transition associated with barriers at the several species traits? To answer these questions we develop new methods in the scope of landscape genetics that provide an exhaustive picture on how environmental factors influence the maintenance of the contact zone. In chapter 6 I provide a general discussion on the subjects addressed in the previous chapters, emphasizing the general achievements and suggesting future research lines. Finally, the chapter 7 summarises the major conclusions of this dissertation. Chapter 3 Analyzing species’ distributions with ecological niche modelling I call our world Flatland, not because we call it so, but to make its nature clearer to you, my happy readers, who are privileged to live in Space. — A SQUARE [EDWIN A. ABBOT], Flatland 3.1 E-Clic - Easy climate data converter1 3.1.1 Abstract There are an increasing number of studies that are now focusing on the influence of climate change on species’ distributions. However, access to predictive climatic datasets for future scenarios is difficult due to their specific formats and/or the need to be geographically downscaled. The TYN dataset is freely available to users and provides a synthetic format with several climatic models and IPCC future climate scenarios. Moreover, the CRU historical dataset (1901 – 2000) is also available which allows users to create baseline models for current climatic variables. E-Clic is a free, user-friendly software package that offers three different ways to convert these two datasets into a spatially explicit raster format which is compatible with the most common geographic information systems and usable on different platforms. 1This section refers to the published article: Tarroso, P. & Rebelo, H. (2010) E-Clic - Easy climate data converter. Ecography,33, 617–620. 52 FCUP Analyzing species’ distributions with ecological niche modelling 3.1.2 Easy climate data converter Climatic variables provide key input parameters when modelling species’ distributions (Guisan & Zimmermann, 2000). Despite the range of different variables that are frequently used in ecological modelling, climate data is commonly used due to its continuous influence in shaping and limiting species’ distributions, and also because of its correlation with other biologically meaningful variables (e.g. land cover) (Huntley & Webb III, 1989; Thomas et al., 2004; Thuiller et al., 2004). The changing dynamics of climatic forces are highly correlated with the shift of species’ distributions, as shown in the past by the repeated moving patterns of populations during the climate oscillations of the glacial and interglacial periods (Huntley & Webb III, 1989). The current forecast of rapid climate change will force a shift in species’ distributions and, in some cases, will increase their risk of extinction (Parmesan & Yohe, 2003; Thomas et al., 2004). Predicting changes in species’ distributions and the likelihood of their extinction is an important measure in understanding the possible impacts of climate change on biodiversity (Hannah et al., 2002). Recently, several studies have been carried out on a number of species, including plants, birds, marine animals, amphibians and reptiles (Hawkes et al., 2007; Lemoine et al., 2007; Thuiller et al., 2005; Araújo et al., 2006). Consequently, achieving accurate distribution models using present day climatic variables is the first step in projecting these models for future predicted climate scenarios. This will allow for more informed decisions regarding the possible loss of biodiversity. As one of the greatest political, social and scientific concerns of our time, research on climate change has generated a multidisciplinary interest which has produced a wide range of models that combine different social-economic scenarios to predict several climatic variables in the future (Nakicenovic & Swart, 2000). Climatic models are available at different geographic scales but usually have a coarse resolution and can be stored in a difficult format that is not easily integrated into the most commonly used geographic information systems (GIS) (e.g. ArcMap) and ecological modelling software (e.g. MaxEnt). The TYN dataset (Mitchell et al., 2004) provides several downscaled models that are frequently used to study the impact of climate change. Specifically, it consists of monthly data broken down into five variables (precipitation, daily mean temperature, diurnal temperature range, vapour pressure, cloud cover) for five global circulation models (CGCM2, CSIRO mk 2, DOE PCM, HadCM3, ECHam4) over 100 year intervals (from 2001 to 2100). It also includes four special reports on emission scenarios (A1FI, A2, B1 and B2) resulting from the Intergovernmental Panel on Climate Change (IPCC) meeting (Nakicenovic & Swart, 2000). These high-resolution datasets are freely available from the Climate Research Unit web page FCUP 53 E-Clic - Easy climate data converter ( http://www.cru.uea.ac.uk/cru/data/hrg.htm ) at two different scales: the TYN SC 1.0 with European data downscaled to 10’ resolution and the TYN SC 2.0 with global coverage at a 0.5oresolution. Although new future climate models and datasets occasionally emerge, the TYN dataset remains as a stable and complete source of future climate predictions for ecological studies. In addition to the TYN dataset, the CRU TS 1.0 dataset (Mitchell et al., 2004) is also freely available to researchers and comprises of historical climate data ranging from 1901 to 2000 for the same variables, coverage and resolution as TYN SC 1.0. This data is highly useful because it provides researchers the opportunity to widen the time frame of their studies and to test the accuracy of their models if a species range has been well documented over the last century. The equivalent historical dataset for the TYN SC 2.0, the CRU TS 2.10 (Mitchell & Jones, 2005) is also available, including four new variables (monthly average daily maximum temperature, monthly average daily minimum temperature, wet day frequency, frost day frequency). Here we present a new software program to convert public climatic datasets (TYN SC and CRU TS) into formats that are more commonly used and therefore can be directly utilized in GIS. ‘Easy climate data converter’ (E-Clic) ( http://purl.oclc. org/eclic/ ) is a free, open-source application which is written in Python and can be used in different operating systems such as Windows, Mac OS or Linux (Fig. 3.1). It can read both TYN SC and CRU TS datasets and is able to convert them to ASCII raster maps - a format easily readable in most GIS and ecological modelling software (e.g. Arcgis, Idrisi, Maxent). Fig. 3.1 – E-Clic built-in graphical user interface with all the available options to convert from TYN or CRU data to a raster format, running in (from front to back) Linux, Mac OS and Windows. 54 FCUP Analyzing species’ distributions with ecological niche modelling Prior to data conversion, users should download and decompress the TYN and/or CRU datasets to a chosen folder. The CRU dataset is available in two different forms; one single file that includes all the CRU TS range data, and as separate files containing data for each decade. The single file that includes all data must be used. Unpacking of the data from the files done by E-Clic follows the procedure described online with the dataset (see dataset Internet site for more details). In short, it consists of searching the data for the defined period of time, converting it into real units given by the model/scenario, and then testing it with the allowed minimum and maximum limits for the chosen variable. The user can choose one model, scenario and variable from all of the TYN supporting data. The output is the average model for the chosen time period and, optionally, a raster for each month within the same period. The output created by E-Clic is easily integrated into most GIS packages and modelling software. The ability to query a spatial and temporal dataset and create raster data from them avoids, in most cases, the need to post-process the raster datasets. Although there are other software packages available to convert TYN SC and CRU TS data into different formats (see the dataset web page for more software and Solymosi et al. 2008), they lack some of the features presented in E-Clic such as a user-friendly interface, ability to be used on different platforms and/or the option to export data into a raster format. As E-Clic is written in Python, it benefits from the multi-platform availability of this programming language (Fig. 3.1). Moreover, it is not dependent upon external python modules and can be run directly after python installation ( http://www.python.org ), which is already available as default on some platforms (e.g. most of Linux distributions and Mac OS). In addition, this software also offers three different interfaces to convert data, depending on the user’s preferences. Fig. 3.2 – E-Clic interface using ArcGIS software as a toolbox. All options are available and the extent may be automatically defined. FCUP 55 E-Clic - Easy climate data converter E-Clic may be used at the command line with data inputted in a strict sequence format (see the online instructions in E-Clic website for more details on how to use E-Clic at the command line). The command line feature may be used to build quick batches to produce large quantities of data. The built-in graphical user interface (Fig. 3.1) is displayed when the user runs E-Clic without any additional command, the interface is launched displaying all the possible options making the process of data conversion Fig. 3.3 – Average daily temperature forecasts for the years 2025, 2050, 2075 and 2100 for four available scenarios (A1FI, A2, B1 and B2). The plots represent the frequency of pixels with the same value of temperature. The year 2025 is plotted with a solid line (—) 2050 with a dashed line (- - -), 2075 with a dash and dot line (- . -) and 2100 with a dotted line (.. . ). 56 FCUP Analyzing species’ distributions with ecological niche modelling between formats relatively simple. Finally, it may also be integrated into ArcGIS (ESRI, Redlands, CA, USA) as a toolbox (Fig. 3.2) where the graphical interface provided by the GIS package can be used to access all of E-Clic’s functions. When it is used within the GIS or when it detects a valid license, it will also output the raster in a geoTIF format. All result files are saved to the chosen output directory and the names are explicit to the data they contain. We predict that the main application of the output raster maps created by E-Clic will be in species distribution modelling for future climate scenarios. We present here a simple example of the output maps for an area covering the Iberian Peninsula and part of North Africa (Fig. 3.3). With TYN_SC_1.06 data, we have built annual mean values for four daily mean temperature scenarios (ranging from the more extreme A1FI, A2, B2 to the less severe B1) for the years 2025, 2050, 2075 and 2100. We synthesized the data using histograms where the frequency of map cells with the same value is plotted against the temperature. The results show that all the scenarios predict an increasing number of cells with higher temperatures with time and, as expected, the A1FI scenario has the greatest increase in temperature. The increase of the maximum values of the predicted mean temperatures for this area between 2025 and 2100 reaches the highest value of 6.7oC with the A1FI scenario, whereas the lowest is 2.7oC for the B1 scenario. Even though there are a limited number of climatic data conversion software packages available, the majority of them require some knowledge of computer programming. Since the majority of users will not be familiar with this type of programming, especially those who come from different scientific fields such as ecology, accessing climatic data can be difficult. Here we present a user-friendly application that allows anyone with the most basic GIS knowledge to easily access climatic data for future scenarios. Given that climate change studies are becoming more important, we hope that this software will help more researchers to have access to the core datasets for their studies. 3.1.3 Acknowledgments The authors are grateful to Jon Flanders for his extensive review of the manuscript that greatly improved its quality. The authors are also grateful to James Harris, Anna Perera and two anonymous referees for their comments and suggestions. PT and HR were funded by the ’Fundação para a Ciência e Tecnologia’ doctoral grants SFRH/BD/42480/2007 and SFRH/BD/17755/2004, respectively. FCUP 57 E-Clic - Easy climate data converter 3.1.4 References Araújo, M.B., Thuiller, W. & Pearson, R.G. (2006) Climate warming and the decline of amphibians and reptiles in Europe. Journal of Biogeography 33, 1712–1728. Guisan, A. & Zimmermann, N.E. (2000) Predictive habitat distribution models in ecology. Ecological Modelling 135, 147 – 186. Hannah, L., Midgley, G. & Millar, D. (2002) Climate change-integrated conservation strategies. Global Ecology and Biogeography 11, 485–495. Hawkes, L.a., Broderick, a.C., Godfrey, M.H. & Godley, B.J. (2007) Investigating the potential impacts of climate change on a marine turtle population. Global Change Biology 13, 923–932. Huntley, B. & Webb III, T. (1989) Migration: species’ response to climatic variations caused by changes in the earth’s orbit. Journal of Biogeography 16, 5–19. Lemoine, N., Schaefer, H.C. & Böhning-Gaese, K. (2007) Species richness of migratory birds is influenced by global climate change. Global Ecology and Biogeography 16, 55–64. Mitchell, T.D., Carter, T.R., Jones, P.D., Hulme, M. & New, M. (2004) A comprehensive set of high-resolution grids of monthly climate for Europe and the globe: the observed record (1901-2000) and 16 scenarios (2001-2100). Mitchell, T.D. & Jones, P.D. (2005) An improved method of constructing a database of monthly climate observations and associated high-resolution grids. International Journal of Climatology 25, 693–712. Nakicenovic, N. & Swart, R. (2000) Special report on emissions scenarios: a special report of Working Group III of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK. Parmesan, C. & Yohe, G. (2003) A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42. Solymosi, N., Kern, A., Maróti-Agóts, A., Horváth, L. & Erdélyi, K. (2008) TETYN: An easy to use tool for extracting climatic parameters from Tyndall data sets. Environmental Modelling & Software 23, 948–949. 58 FCUP Analyzing species’ distributions with ecological niche modelling Thomas, C.D., Cameron, A., Green, R.E., Bakkenes, M., Beaumont, L.J., Collingham, Y.C., Erasmus, B.F.N., de Siqueira, M.F., Grainer, A., Hannah, L., Hughes, L., Huntley, B., van Jaarsveld, A.S., Midgley, G.F., Miles, L., Ortega-Huerta, M.A., Peterson, A.T., Phillips, l.L. & Williams, S.E. (2004) Extinction risk from climate change. Nature 427, 145–148. Thuiller, W., Araújo, M. & Lavorel, S. (2004) Do we need land-cover data to model species distributions in Europe? Journal of Biogeography 31, 353–361. Thuiller, W., Lavorel, S., Araújo, M. & Sykes, M. (2005) Climate change threats to plant diversity in Europe. Proceedings of the National Academy of Sciences 102, 8245–8250. FCUP 59 Simapse - Simulation Maps for Ecological Niche Modelling 3.2 Simapse - Simulation Maps for Ecological Niche Modelling2 3.2.1 Abstract 1. Artificial neural networks (ANN) are known for their powerful predictive power in the analysis of both linear and non-linear relationships. They have been successfully applied to several fields including ecological modelling and predictive species’ distributions. 2. Here we present Simapse – Simulation Maps for Ecological Niche Modelling, a free and open-source application written in Python and available to the most common platforms. It uses ANNs with back-propagation to build spatially explicit distribution models from species data (presence/absence, presence-only, and abundance). 3. The main features include the automatic production of replicates with different subsampling methods and total control of ANN structure and learning parameters. 4. Simapse uses common text formats as main input and output and provides assessment of variable importance and behaviour and measurement of model fitness. Keywords: artificial neural networks, ecological niche modelling, species distribution, Python 3.2.2 Introduction Artificial Neural Networks (ANNs) have been used in scientific fields where pattern recognition is a primary need and its application to biological systems had a considerable growth during the past decades (Lek et al., 1996; Lek & Guégan, 1999; Özesmi et al., 2006a). As a machine learning algorithm, it is with no surprise that ANNs are increasingly being used in ecological studies. These algorithms are usually seen as more powerful to deal with complex ecological datasets than other methods (Brosse & Lek, 2000; Olden et al., 2008; Özesmi et al., 2006b; Pearson et al., 2002). A common type of ANN is the feed-forward neural network with back-propagation learning (BPN; Lek & Guégan, 1999). This network has a layered structure of neurons connecting the inputs to an output through one or several hidden layers (see supplementary material A.1 for more details). The BPN has been used to model ecological systems due to its efficient learning ability and to its simple nature that makes it easy to understand (Lek & Guégan, 1999; Özesmi et al., 2006a). The unit of the BPN is an artificial neuron with an activation function, usually linear or sigmoid. It squashes the 2This section refers to the published article: Tarroso, P., Carvalho, S. & Brito, J.C. (2012) Simapse - Simulation Maps for Ecological Niche Modelling. Methods in Ecology and Evolution,3, 787-791 66 FCUP Analyzing species’ distributions with ecological niche modelling spectively). We thank A. Townsend Peterson and an anonymous reviewer for the helpful comments to a preliminary version of the manuscript. 3.2.9 References Araújo, M.B. & New, M. (2007) Ensemble forecasting of species distributions. Trends in Ecology & Evolution 22, 42–7. Brosse, S. & Lek, S. (2000) Modelling roach (Rutilus rutilus) microhabitat using linear and nonlinear techniques. Freshwater Biology 44, 441–452. Davis, J. & Goadrich, M. (2006) The relationship between Precision-Recall and ROC curves. Proceedings of the 23rd International Conference on Machine Learning pp. 233–240. Dimopoulos, I., Chronopoulos, J., Chronopoulou-Sereli, A. & Lek, S. (1999) Neural network models to study relationships between lead concentration in grasses and permanent urban descriptors in Athens city (Greece). Ecological Modelling 120, 157–165. Dimopoulos, Y., Bourret, P. & Lek, S. (1995) Use of some sensitivity criteria for choosing networks with good generalization ability. Neural Processing Letters 2, 1–4. Fu, L. & Chen, T. (1993) Sensitivity analysis for input vector in multilayer feed forward neural networks. IEEE International Conference on Neural Networks, pp. 215–218. Gevrey, M., Dimopoulos, I. & Lek, S. (2003) Review and comparison of methods to study the contribution of variables in artificial neural network models. Ecological Modelling 160, 249–264. Gevrey, M., Dimopoulos, I. & Lek, S. (2006a) Two-way interaction of input variables in the sensitivity analysis of neural network models. Ecological Modelling 195, 43–50. Gevrey, M., Lek, S. & Oberdorff, T. (2006b) Utility of sensitivity analysis by artificial neural network models to study patterns of endemic fish species. Ecological Informatics: Scope, techniques and applications (ed. F. Recknagel), chap. 14, pp. 293–306, Springer, Berlin, 2nd edn. Lek, S., Delacoste, M., Baran, P., Dimopoulos, I., Lauga, J. & Aulagnier, S. (1996) Application of neural networks to modelling nonlinear relationships in ecology. Ecological Modelling 90, 39–52. FCUP 67 Simapse - Simulation Maps for Ecological Niche Modelling Lek, S. & Guégan, J.F. (1999) Artificial neural networks as a tool in ecological modelling, an introduction. Ecological Modelling 120, 65–73. Lobo, J.M., Jiménez-Valverde, A. & Real, R. (2007) AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography 17, 145–151. Olden, J.D. & Jackson, D.A. (2002) Illuminating the “black box”: a randomization approach for understanding variable contributions in artificial neural networks. Ecological Modelling 154, 135–150. Olden, J.D., Lawler, J.J. & Poff, N.L. (2008) Machine learning methods without tears: a primer for ecologists. The Quarterly Review of Biology 83, 171–193. Özesmi, S. & Özesmi, U. (1999) An artificial neural network approach to spatial habitat modelling with interspecific interaction. Ecological Modelling 116, 15–31. Özesmi, S., Tan, C. & Özesmi, U. (2006a) Methodological issues in building, training, and testing artificial neural networks in ecological applications. Ecological Modelling 195, 83–93. Özesmi, U., Tan, C., Özesmi, S. & Robertson, R. (2006b) Generalizability of artificial neural network models in ecological applications: Predicting nest occurrence and breeding success of the red-winged blackbird Agelaius phoeniceus.Ecological Modelling 195, 94–104. Pearson, R.G., Dawson, T.P., Berry, P.M. & Harrison, P.A. (2002) SPECIES: A Spatial Evaluation of Climate Impact on the Envelope of Species. Ecological Modelling 154, 289– 300. Peterson, A.T., Papes, M. & Sober, J. (2007) Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecology 3, 63–72. 68 FCUP Analyzing species’ distributions with ecological niche modelling 3.3 Evaluating ecological niche models with virtual and real species3 3.3.1 Abstract Ecological niche models are widely used to study species’ distributions based on presence data collected in the field or gathered from diverse collections. Several algorithms are available and testing an ENM is usually done comparatively between different methods. Since different algorithms tend to answer different questions, this comparison is not trivial to execute. Here we present a framework to make exhaustive tests of an ENM and to assess its predictive ability under controlled conditions. For that purpose we use a set of virtual species representing different ecological demands and a set of real species based on distributional atlases. Both data sets offer an accurate description of the complete distribution of each species, however, the virtual species allows a full control by providing both presence and absence areas. To perform the models we used an artificial neural network (ANN) with nine commonly used environmental variables. We test the performance of the ANN with varying network structures and with different sample sizes. Models for the virtual species with pseudo-absences were compared to the control model with the accurate virtual presences and absences. The performance of models for the real species were assessed with the real distribution and compared with the results of the virtual species, under the same conditions. The main finding is that ANNs are accurate in predicting the distribution of both virtual and real species. The structure of the network has little implications for the final models, but simple networks have generally good results and are preferred. Low sample size produces accurate models and the spatial uncertainty of the resulting models tends to increase with the decreasing number of available samples. The models of the real species were found to correctly predict presence when compared with the full distribution information. The framework that we have used here can be extended to other ecological modelling algorithms to fully assess their performance. It allows to test the implications of the set of parameters typical of each algorithm, and to assess the performance under different scenarios of data availability. 3This study is presently submitted to an international journal in the science-citation index: Tarroso, P. & Brito, J.C. Evaluating ecological niche models with virtual and real species. FCUP 69 Evaluating ecological niche models with virtual and real species 3.3.2 Introduction The spatial context of a species distribution is self-evident and its analysis, either in current time or in other time periods, is a central aspect of biogeography. A research line in this field is the study of the relation between species’ locations and the environment. Ecological niche modelling (ENM) has been widely used to study relationships between species’ distribution and ecological niche variables (ENV) and to predict occurrence in areas outside the sampling locations (Guisan & Zimmermann, 2000; Ferrier, 2002; Guisan & Thuiller, 2005; Elith & Leathwick, 2009). The wide range of applications of these tools reveals the importance of ENMs in current research. It comprises diverse areas as biodiversity conservation (Brito et al., 2009, 2011; Carvalho et al., 2011a), evolutionary biology (Martínez-Freiría et al., 2009; Tomovi´ cet al., 2010), human archeology (Banks et al., 2008), and predictions of climate change effects on biodiversity distribution (Carvalho et al., 2010, 2011b; Rebelo et al., 2010). To predict distributions, ENMs rely on algorithms that relate species occurrence with the ENVs. Comparisons of different algorithms have suggested that machine learning methods, such as artificial neural networks or maximum entropy, are more robust than other traditional methods (Segurado & Araujo, 2004; Elith et al., 2006). Artificial neural networks (ANN) were defined as a promising method for habitat distribution modelling (Guisan & Zimmermann, 2000), but until recently their use was restricted to a small circle of ecological modellers (Olden et al., 2008). Nonetheless, research with ANN is active (Özesmi et al., 2006; Park & Chon, 2007; Olden et al., 2008; Tarroso et al., 2012) and the method has been applied during the last two decades to study several ecological and environmental processes (Lek & Guégan, 1999; Lek et al., 1996; Araújo et al., 2006; Özesmi et al., 2006; Xavier et al., 2010). An advantage of ANN over other methods is the ability to fit non-linear relationships resulting in an high predictive power which can virtually approximate any continuous function (Gevrey et al., 2003; Olden et al., 2008), and thus rendering it an adequate method to model complex ecological data (Barry & Elith, 2006). ANN have been shown to outperform other methods (Segurado & Araujo, 2004). The unit of ANN is the biologically inspired artificial neuron, which processes an input signal and yields an output accordingly to its internal activation function (usually a sigmoid function). These neurons are fully connected by weights on a network structured in layers. The most simple form is a three layered network: the first layer corresponds to the input neurons (with the same number as ENVs available); the second layer is known as hidden layer with an arbitrary number of neurons; and the third is the output layer with one neuron that returns the final result. As a supervised machine learning method, the error of the network 70 FCUP Analyzing species’ distributions with ecological niche modelling is assessed with training data and is propagated backwards through the network to adjust the value of the connecting weights (back-propagation learning). The repetition of this process for several iterations produces a trained network with minimal error for the given input that was used to predict the presence of a species in a specific location with the same set of ENVs. The evaluation of the performance of an ENM is based on the comparative analysis of multiple algorithms with single or multiple species (e.g. Segurado & Araujo, 2004; Elith et al., 2006). This task is not trivial due to the particular design or operation of an algorithm affecting the possible range of commission errors (Peterson et al., 2007), or to the impact of different species strategies (i.e. generalist or specialist), resulting in different occurrence extents that has consequences on the general model performance (Araújo & Williams, 2000; Segurado & Araujo, 2004). Thus, the knowledge of the full distribution of the species would be a valuable asset for a complete analysis of model performance. Although an exhaustive record of distribution data is nearly unattainable with a real species, a virtual species provides an efficient method to overcome this difficulty (Hirzel et al., 2001; Guisan et al., 2006; Zurell et al., 2010). Virtual species enlighten model performance under different circumstances, as the researcher is aware of the complete virtual reality. However, caution is needed when analysing virtual species results. Despite the immense data it potentially provides, the success of ENMs in the virtual world does not guarantee the same success in the real world (Zurell et al., 2010). As such, modelling species distribution with distribution data from atlases constitutes a complementary alternative to evaluate ENM performance. These presence data compendiums usually provide details of sampling effort and are based on a uniform strategy over a defined geographical area and time period constituting an important resource for biodiversity conservation (Robertson et al., 2010) and extensive species data for modelling purposes. A framework to assess the performance of a given algorithm for ENM combining results from virtual species and atlases data are lacking in the scientific literature, especially in the case of ANN. Species’ distribution datasets are usually composed of presence-only data. Absences are associated to high uncertainty (Lobo et al., 2010), therefore most studies are limited to the presence locations. Simulating absences, generally referred as pseudo-absences, is a common method to overcome this limitation. Nevertheless, indulgent usage of artificially generated absence data may have negative impacts on the final model (Chefaoui & Lobo, 2008; Vanderwal et al., 2009; Lobo et al., 2010). Much has been done to evaluate ENMs strategies; however, a framework that provides positive and negative controls to analyse the performance of an ENM using FCUP 71 Evaluating ecological niche models with virtual and real species Table 3.2 – Ecological niche variables (ENV) general information and usage by the generalist virtual species (GVS), specialist virtual species (SVS), and the random virtual species (RVS). GVS (n=4612) SVS (n=778) RVS (n=3764) ENV description Code Units Range Function Mean±sd Function Mean±sd Mean±sd Elevation elv m 0–2989 logistic 590.7±363.6 – 300±280.2 645±416.3 Annual precipitation ap mm 227–1755 logistic 586.2±199.9 – 1075.2±216.5 653.7±256.7 Precipitation of the driest month pdm mm 0–110 normal 13.1±9.9 – 32.6±14.3 18.1±15.9 Precipitation seasonality ps 11–76 logistic 42±13.4 – 38.5±11.9 39.6±14.1 Precipitation of the wettest month pwm mm 32–270 logistic 78.8±31 – 145.9±34.3 86.1±35.9 Slope slp o0–68.2 logistic 4.7±5.1 – 6.8±7.1 6.1±7.1 Temperature annual range tar oC 13.5–34.1 logistic 27.5±3.8 logistic 18.6±2.2 26.7±4.1 Minimum temperature of the warmest month mtwm oC 11.4–36.5 logistic 29.6±3.1 normal 23.2±1.3 28.5±3.9 Minimum temperature of the coldest month mtcm oC -10.3–8.5 logistic 2.1±2.9 normal 4.6±1.9 1.8±3.0 species with different niche requisites and to test the impact of pseudo-absence data it is still lacking. The combination of virtual species with real data provide an opportunity to test ENM thoroughly with a design that allows comparable models to be built with effective controls of their performance. The objective of this study is to develop a testing framework for evaluation of an ENM. The general objective is divided in three key aspects: 1) comparison between presence-only and presence and absence data; 2) assessment of the impact of algorithm parameters; and 3) effect of sample size. We use presence datasets of virtual species, to test both algorithm parameters and sample size influence, and datasets of real species distribution from atlases to confirm them. 3.3.3 Methods We developed a framework to extensively test the ANN predicting ability with virtual and real species (Fig. 3.6). The study area included the Iberian Peninsula covering an extent between the coordinates 43.8oN, 36.1oN, 9.5oW and 3.3oE with a mainland area of 580 000km2. A set of ENVs including climatic and topographic data (Table 3.2) were used for modelling both real and virtual species. The climatic ENVs were obtained from Worldclim ( www.worldclim.org ; Hijmans et al., 2005) with the exception of slope which was derived from altitude data. All ENVs were resampled to 0.09oto match the resolution of the presence data for real species, resulting in a total of 7686 cells. 72 FCUP Analyzing species’ distributions with ecological niche modelling Species data We created three virtual species with different ranges of presence response simulating different niche requisites: 1) a generalist species (GVS) with broad distribution in the Iberian Peninsula; 2) a specialist species (SVS) with restricted and temperaturedependent range; and 3) a random virtual species (RVS) with distribution independent of the used ENVs (appendix B.1). To create both GVS and SVS we used the averaged partial niche coefficient for the chosen variables (Hirzel et al., 2001). Each partial niche coefficient, defining the virtual species response, was extracted from the ENV using either a normal or a logistic function with a range between 0 and 1 (table 3.2). All responses were weighted equally to obtain a distribution with continuous Fig. 3.6 – Framework of the study and general procedures to use artificial neural networks. The input data fed to Simapse was divided in virtual species (top) and real species (bottom). Virtual species provides a context to extensively test the effects of structure of the neural network and sample size. This is compared to presence and absence data which has a dual role (dotted arrows): positive control - it provides full virtual reality to the modelling process; and negative control – it provides a random reality which has no relation with ENVs. The access to full virtual reality allows the calculation of a real AUC (rAUC, with correct virtual absences plus the input presences) as opposed to the final AUC (fAUC, with the pseudo-absences) of the consensus model. Train and test AUC are processed (trAUC and teAUC, respectively). Three real species with presence-only data are used. FCUP 73 Evaluating ecological niche models with virtual and real species value. A threshold was then used to convert the presence probability to a binary presence/absence distribution from which datasets of presence-only and presenceabsence locations were randomly extracted. All virtual species were created using R (R Development Core Team, 2012) with the package rgdal (Keitt et al., 2010). The GVS was created to make a general use of all ENVs range and to have a broad spatial extent in the study area. It had a normal response to the precipitation of the driest month (PDM) and logistic responses to all other ENVs (table 3.2, appendix B.1). An arbitrary threshold of 0.75 was chosen to create a binary distribution for the GVS. We generated a list of 400 random locations equally split between presence and absence datasets. We also generated presence-only datasets with 25, 50, 75, 100, 125, 150, 175 and 200 locations. The SVS was built to be more restricted in its distribution and to have specific thermal requirements. It has a logistic response to the temperature annual range (TAR) and normal response to the maximum temperature of the warmest month (MTWM) and minimum temperature of the coldest month (MTCM; table 3.2, appendix B.1). The binary distribution was created with a threshold of 0.5. The generated datasets include a set of 100 presence and 100 absence locations, plus presence-only datasets with 20, 40, 60, 80 and 100 locations. The RVS was generated by assigning a random number between 0 and 1 to each cell of the study area. The binary distribution of the RVS was produced with a threshold of 0.5 and 400 locations were extracted, equally divided between presences and absences (appendix B.1). The real species used in this study included three lacertid lizards (Lacerta schreiberi,Lacerta agilis, and Psammodromus algirus) which exhibit different distributions and ecological requirements in Iberia. Both L. schreiberi and L. agilis are specialists and occupy different spatial extents, the first ranging in larger areas of the humid north-western Iberia (n=1152; 15% of study area) and the latter restricted to a few high altitude cells in the Pyrenees (n=11; 0.14%). P. algirus covers wider habitat types, occurring throughout most of the study area (n=4596; 60%). The distribution data for each species were obtained by digitalisation of the atlases of Spain and Portugal (Pleguezuelos et al., 2004; Loureiro et al., 2008). Modelling procedure All models were built using the ANN implemented in Simapse (Tarroso et al., 2012). This software uses an ANN with back-propagation learning algorithm to produce potential species’ distributions from binary data. All models were ran with random subsampling with 30 repetitions, from which a final averaged consensus model was built 74 FCUP Analyzing species’ distributions with ecological niche modelling using the standard deviation between models as a measure of uncertainty (Fig. 3.6). Each repetition was trained with 25 burn-in and 2500 training iterations (500 reports after a sequence of 5 internal iterations), and tested with 25% of the data chosen randomly from the input dataset to assure proper generalization by reducing overfitting (Dimopoulos et al., 1995). At each reported iteration, the value of the area under curve (AUC) of the receiver operating characteristic curve (ROC) was saved together with the sum of squared error of the network. A repetition was only accepted in the final consensus model if, at some point of the training process, a network was reached with both train and test AUC≥0.7 with real species, GVS and SVS. For RVS, the network was accepted with AUC≥0.5. We tested the influence of network structure with the GVS and SVS by running models with different structures, ranging from 1 to 18 neurons (double the available ENVs) in a single layer. These models were fed with the virtual presence-only data: 200 presence locations for the GVS and 100 locations for the SVS. Pseudo-absences were generated in the same number as presences. The learning rate (LR) was optimized heuristically for each species. In the case of GVS, the LR was set to 0.01, for SVS it was set to 0.1. Momentum was kept constant at 0.1 in all models. For each model we calculated the average train and test AUCs, (trAUC and teAUC, respectively) and the final AUC (fAUC) value of the consensus model with the presences and pseudo-absences as given by Simapse. The real AUC (rAUC) of the same model was also processed after the substitution of the pseudo-absences by the correct virtual absences (see Fig. 3.6). To evaluate the influence of the number of available presences on the final model, we produced models for the GVS and SVS presence-only datasets with a single layer of nine neurons (same number as ENVs available). The LR value for the GVS was set to 0.1 for 25 and 50 presence locations and to 0.01 for 75, 100, 125, 150, 175 and 200 locations. The chosen LR value for the SVS was 0.6 for the dataset of 20 presence locations, 0.2 for the 40 presence and 0.1 for 60, 80 and 100 presence locations. The trAUC, teAUC, fAUC and rAUC values were calculated. As a positive control for the models, we produced models with presence-absence datasets from all the virtual species (Fig. 3.6), which allowed evaluating the effect of the pseudoabsences on final model quality. The RVS, with no relation with the environment, acts as a negative control, to assess the ability to fit only noisy data. The models of the real species followed the standard modelling procedures. The locations were input to Simapse with the ENVs to produce an averaged consensus prediction. We have chosen a structure with a single hidden layer with 9 neurons (same number as ENVs). The chosen LR was different for each species to maximize the gradient descent of the FCUP 75 Evaluating ecological niche models with virtual and real species network error. Therefore, for the L. agilis,L. schreiberi and P. algirus, the LR value was set to 0.7, 0.005, and 0.0005, respectively. 3.3.4 Results Virtual species predictions As expected, the three virtual species have different niche occupancies (table 3.2) resulting in different spatial extents. The GVS is widespread in southern Iberia (4612 cells), while the SVS occupies the mountainous northern areas (778 cells). The RVS is present throughout all study area (3764 cells) without any visible clustering (appendix B.1). Fig. 3.7 – Effects of the number of neurons in the single hidden layer on network structure. The generalist virtual species (GVS) is represented by solid regression line and circles and the specialist virtual species (SVS) by dashed regression line and triangles. Ratios of cells correctly classified as presence or absence (a), false presences as presences (b), and false absences as absences (c) are given. 82 FCUP Analyzing species’ distributions with ecological niche modelling may increase the probability of falling in a non-sampled presence, thus creating difficulties to the learning algorithm and consequently affecting the AUC values (Lobo et al., 2007; Chefaoui & Lobo, 2008). On the other side, modelling with real presences and absences showed a great accuracy with a steep ANN learning path, as is common with presence/absence methods (Zaniewski et al., 2002; Barry & Elith, 2006). Models for the SVS (10% of the total area; see table 3.2) exhibited high accuracies. The high AUC values when compared to rAUC demonstrate an opposite picture for the GVS. The probability of a pseudo-absence to fall in a location with a real absence was higher, leading to an effective learning and, thus, higher AUC values. In fact, the differences between modelling with real absences or pseudo-absences are unnoticeable for the SVS when analysing the AUC values achieved with the different strategies. The strategy of comparing models based on the same presence dataset using accurate virtual absences or pseudo-absences allowed to compare AUC values. The reduced uncertainty of the virtual absences was reflected in the stability of the rAUC, but randomly chosen absences resulted in fAUC oscillation. This exposes the relation of AUC values based on pseudo-absences with other factors like the relative area of occurrence of the species and the area where are pseudo-absences are generated. This topic has been subject of intense debate in recent literature (e.g. Engler et al., 2004; Lobo et al., 2007; Peterson et al., 2007; Chefaoui & Lobo, 2008; Vanderwal et al., 2009; Lobo et al., 2010) but further discussion is beyond the objectives of the present research. ANN parameters The effect of network structure on the models was analysed by changing the number of neurons available in the single hidden layer. The classification success of the predictive models was very high with all tested network architectures. Interestingly, simple networks, i.e. lower number of neurons, achieve a classification success comparable to more complex networks with the tested datasets. Although the process of selecting network architecture is usually heuristic, modellers should keep the simplest network as it is less prone to overfitting and, thus, is able to predict better in new regions (Dimopoulos et al., 1995; Özesmi & Özesmi, 1999; Özesmi et al., 2006). The problem of overfitting in very complex networks is related to the trend of noise fitting and leading to a lack of predictability (Özesmi & Özesmi, 1999). The slight decrease of predictive success with increasing complexity of the network observed in the GVS reflects this problem: its intricate pattern of presences and pseudo-absences is fitted FCUP 83 Evaluating ecological niche models with virtual and real species along with noise when using a higher number of neurons, inflating the number of false absences. On the other hand, the restricted range and specificity to a few ENVs of the SVS make it less prone to noise presence. Modelling in this case, is very effective with a small increase of predictability accompanying the increasing complexity of the network structure. Sample size Modelling the GVS and SVS with different number of presences resulted in high classification rates, even with the lowest number of presences. Even so, additional presences are extremely informative, as shown by the decrease of false presences and false absences in the SVS and GVS, respectively (Fig. 3.10) but also by the decrease of spatial uncertainty (appendix B.2). The increase of the available samples is usually directly proportional to the accuracy of the predictions (Barry & Elith, 2006; Hernandez et al., 2006; Wisz et al., 2008) and the extra information contained in additional presences allows to better depict the relations between species’ locations and the predictors. The increasing complexity of these relations are data hungry (Barry & Elith, 2006) as seen by the unstable results of the GVS with low number of samples or, on the other side, the rapid achievement of very high correct classification values with small sample size by the SVS. The increasing number of presences in this study has a parallel increase in number of pseudo-absences and, as seen before, the presence uncertainty related to these extra locations in GVS results in lower AUC values when measured with the generated pseudo-absences (Fig. 3.11), even in the presence of accurate models. Real species The final step in the algorithm testing was to submit it to the scrutiny of real data. The predictive models for the three species were very accurate, with well defined presences and absences for L. agilis and L. schreiberi (Fig. 3.12a, 3.12b). Although the latter has a broader extent, both species have specific niche requirements and the results are similar to those of the SVS. In this case of a very restricted distribution in a larger region of study, generalization of the niche to other areas occurs as a side effect of using pseudo-absences (Engler et al., 2004; Vanderwal et al., 2009; Lobo et al., 2010). This results in very accurate prediction of the presences (all presences were correctly predicted) but areas where the species was not detected also have high values of prediction (i.e. high commission errors). Regarding P. algirus, with a large number of presences widespread throughout 84 FCUP Analyzing species’ distributions with ecological niche modelling the study area, we achieved less modelling success with an overlap of presence and non-presence area in the predictive space (Fig. 3.12). Most of the non-presence areas resulting from an atlas with intensive and systematic field work can be accepted as a highly probable absence. However, real species distributions are complex, as opposed to the more simple relationships of the virtual world, and despite the importance of the environment as a shaping factor, there are a multitude of other effects, like human and other sources of natural disturbance, limiting or extending the occurrence of species (Barry & Elith, 2006). In this study we do not use a human related predictor or other source of potential disturbance, therefore some of these absences fall in the high predictive space, hampering the detection of those gaps. 3.3.6 Conclusions The use of virtual and real data allowed testing Simapse with knowledge of the full, though simplistic, virtual reality and the more complex relationships of the real world. The predictive performance of the ANN when modelling both datasets revealed to be extremely high. We demonstrated the efficacy of the algorithm with simple architectures and with an extensive range of sample sizes. When absences are available with high confidence, presence and absence modelling strategies should be prefererred over presence-only data. The framework we developed here extends previous methods for assessing ENM performance and we expect that it will be used in future algorithm comparisons. 3.3.7 Acknowledgments PT and JCB are funded by Fundação para a Ciência e Tecnologia: SFRH/BD/42480/2007 and Programa Ciência 2007, respectively. 3.3.8 References Araújo, M.B., Thuiller, W. & Pearson, R.G. (2006) Climate warming and the decline of amphibians and reptiles in Europe. Journal of Biogeography 33, 1712–1728. Araújo, M. & Williams, P. (2000) Selecting areas for species persistence using occurrence data. Biological Conservation 96, 331–345. Banks, W., Derrico, F., Peterson, a., Vanhaeren, M., Kageyama, M., Sepulchre, P., Ramstein, G., Jost, a. & Lunt, D. (2008) Human ecological niches and ranges during the LGM in Europe derived from an application of eco-cultural niche modeling. Journal of Archaeological Science 35, 481–491. FCUP 85 Evaluating ecological niche models with virtual and real species Barry, S. & Elith, J. (2006) Error and uncertainty in habitat models. Journal of Applied Ecology 43, 413–423. Brito, J., Fahd, S., Geniez, P., Martínez-Freiría, F., Pleguezuelos, J. & Trape, J.F. (2011) Biogeography and conservation of viperids from North-West Africa: An application of ecological niche-based models and GIS. Journal of Arid Environments 75, 1029– 1037. Brito, J.C., Acosta, A.L., Álvares, F. & Cuzin, F. (2009) Biogeography and conservation of taxa from remote regions: An application of ecological-niche based models and GIS to North-African canids. Biological Conservation 142, 3020–3029. Carvalho, S.B., Brito, J.C., Crespo, E.G., Watts, M.E. & Possingham, H.P. (2011a) Conservation planning under climate change: Toward accounting for uncertainty in predicted species distributions to increase confidence in conservation investments in space and time. Biological Conservation 144, 2020–2030. Carvalho, S.B., Brito, J.C., Crespo, E.J. & Possingham, H.P. (2010) From climate change predictions to actions - conserving vulnerable animal groups in hotspots at a regional scale. Global Change Biology 16, 3257–3270. Carvalho, S.B., Brito, J.C., Crespo, E.J. & Possingham, H.P. (2011b) Incorporating evolutionary processes into conservation planning using species distribution data: a case study with the western Mediterranean herpetofauna. Diversity and Distributions 17, 408–421. Chefaoui, R.M. & Lobo, J.M. (2008) Assessing the effects of pseudo-absences on predictive distribution model performance. Ecological Modelling 210, 478–486. Dimopoulos, Y., Bourret, P. & Lek, S. (1995) Use of some sensitivity criteria for choosing networks with good generalization ability. Neural Processing Letters 2, 1–4. Elith, J., Graham, C.H., Anderson, R.P., Dudík, M., Ferrier, S., Guisan, A., Hijmans, R.J., Huettmann, F., Leathwick, J.R., Lehmann, A., Li, J., Lohmann, L.G., Loiselle, B.A., Manion, G., Moritz, C., Nakamura, M., Nakazawa, Y., Overton, J.M., Peterson, A.T., Phillips, S.J., Richardson, K., Scachetti-Pereira, R., Schapire, R.E., Soberón, J., Williams, S., Wisz, M.S. & Zimmermann, N.E. (2006) Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29, 129–151. Elith, J. & Leathwick, J.R. (2009) Species Distribution Models: Ecological Explanation and Prediction Across Space and Time. Annual Review of Ecology, Evolution, and Systematics 40, 677–697. 86 FCUP Analyzing species’ distributions with ecological niche modelling Engler, R., Guisan, A. & Rechsteiner, L. (2004) An improved approach for predicting the distribution of rare and endangered species from occurrence and pseudoabsence data. Journal of Applied Ecology 41, 263–274. Ferrier, S. (2002) Mapping spatial pattern in biodiversity for regional conservation planning: where to from here? Systematic Biology 51, 331–63. Gevrey, M., Dimopoulos, I. & Lek, S. (2003) Review and comparison of methods to study the contribution of variables in artificial neural network models. Ecological Modelling 160, 249–264. Guisan, A., Lehmann, A., Ferrier, S., Austin, M., Overton, J.M.C., Aspinall, R. & Hastie, T. (2006) Making better biogeographical predictions of species’ distributions. Journal of Applied Ecology 43, 386–392. Guisan, A. & Thuiller, W. (2005) Predicting species distribution: offering more than simple habitat models. Ecology Letters 8, 993–1009. Guisan, A. & Zimmermann, N.E. (2000) Predictive habitat distribution models in ecology. Ecological Modelling 135, 147 – 186. Hernandez, P.A., Graham, C.H., Master, L.L. & Albert, D.L. (2006) The effect of sample size and species characteristics on performance of different species distribution modeling methods. Ecography 29, 773–785. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G. & Jarvis, A. (2005) Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25, 1965–1978. Hirzel, A., Helfer, V. & Metral, F. (2001) Assessing habitat-suitability models with a virtual species. Ecological Modelling 145, 111–121. Keitt, T.H., Bivand, R., Pebesma, E. & Rowlingson, B. (2010) rgdal: Bindings for the Geospatial Data Abstraction Library. Lek, S., Delacoste, M., Baran, P., Dimopoulos, I., Lauga, J. & Aulagnier, S. (1996) Application of neural networks to modelling nonlinear relationships in ecology. Ecological Modelling 90, 39–52. Lek, S. & Guégan, J.F. (1999) Artificial neural networks as a tool in ecological modelling, an introduction. Ecological Modelling 120, 65–73. FCUP 87 Evaluating ecological niche models with virtual and real species Lobo, J.M., Jiménez-Valverde, A. & Hortal, J. (2010) The uncertain nature of absences and their importance in species distribution modelling. Ecography 33, 103–114. Lobo, J.M., Jiménez-Valverde, A. & Real, R. (2007) AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography 17, 145–151. Loureiro, A., Ferrand, N., Carretero, M.A. & Paulo, O. (2008) Atlas dos Anfíbios e Répteis de Portugal. Instituto de Conservação da Natureza e Biodiversidade, Lisboa. Mackenzie, D.I. & Royle, J.A. (2005) Designing occupancy studies: general advice and allocating survey effort. Journal of Applied Ecology 42, 1105–1114. Martínez-Freiría, F., Santos, X., Pleguezuelos, J.M., Lizana, M. & Brito, J.C. (2009) Geographical patterns of morphological variation and environmental correlates in contact zones: a multi-scale approach using two Mediterranean vipers (Serpentes). Journal of Zoological Systematics and Evolutionary Research 47, 357–367. Olden, J.D., Lawler, J.J. & Poff, N.L. (2008) Machine learning methods without tears: a primer for ecologists. The Quarterly Review of Biology 83, 171–193. Özesmi, S. & Özesmi, U. (1999) An artificial neural network approach to spatial habitat modelling with interspecific interaction. Ecological Modelling 116, 15–31. Özesmi, S., Tan, C. & Özesmi, U. (2006) Methodological issues in building, training, and testing artificial neural networks in ecological applications. Ecological Modelling 195, 83–93. Park, Y. & Chon, T. (2007) Biologically-inspired machine learning implemented to ecological informatics. Ecological Modelling 203, 1–7. Peterson, A.T., Pape¸s, M. & Eaton, M. (2007) Transferability and model evaluation in ecological niche modeling: a comparison of GARP and Maxent. Ecography 30, 550–560. Pleguezuelos, J.M., Márquez, R. & Lizana, M. (2004) Atlas y libro rojo de los anfibios y reptiles de España. Organismo Autónomo de Parques Nacionales. R Development Core Team (2012) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. 88 FCUP Analyzing species’ distributions with ecological niche modelling Rebelo, H., Tarroso, P. & Jones, G. (2010) Predicted impact of climate change on European bats in relation to their biogeographic patterns. Global Change Biology 16, 561–576. Robertson, M.P., Cumming, G.S. & Erasmus, B.F.N. (2010) Getting the most out of atlas data. Diversity and Distributions 16, 363–375. Segurado, P. & Araujo, M.B. (2004) An evaluation of methods for modelling species distributions. Journal of Biogeography 31, 1555–1568. Tarroso, P., Carvalho, S.B.S. & Brito, J.C. (2012) Simapse - Simulation Maps for Ecological Niche Modelling. Methods in Ecology and Evolution 3, 787–791. Tomovi´ c, L., Crnobrnja-Isailovi´ c, J. & Brito, J.C. (2010) The use of geostatistics and GIS for evolutionary history studies: the case of the nose-horned viper (Vipera ammodytes) in the Balkan Peninsula. Biological Journal of the Linnean Society 101, 651–666. Vanderwal, J., Shoo, L., Graham, C. & Williams, S. (2009) Selecting pseudo-absence data for presence-only distribution modeling: How far should you stray from what you know. Ecological Modelling 220, 589–594. Wisz, M.S., Hijmans, R.J., Li, J., Peterson, a.T., Graham, C.H. & Guisan, a. (2008) Effects of sample size on the performance of species distribution models. Diversity and Distributions 14, 763–773. Xavier, R., Lima, F.P. & Santos, A.M. (2010) Forecasting the poleward range expansion of an intertidal species driven by climate alterations. Scientia Marina 74, 669–676. Zaniewski, A.E., Lehmann, A. & Overton, J.M.C. (2002) Predicting species spatial distributions using presence-only data: a case study of native New Zealand ferns. Ecological Modelling 157, 261–280. Zurell, D., Berger, U., Cabral, J.S., Jeltsch, F., Meynard, C.N., Münkemüller, T., Nehrbass, N., Pagel, J., Reineking, B., Schröder, B. & Grimm, V. (2010) The virtual ecologist approach: simulating data and observers. Oikos 119, 622–635. Chapter 4 Reconstruction of past Iberian climate The world, my friend Govinda, is not imperfect or confined at a point somewhere along a gradual pathway toward perfection. No, it is perfect at every moment. — HERMANN HESSE,Siddhartha 4.1 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. 1 4.1.1 Abstract The climate in the Iberian Peninsula has evolved since the last glacial maximum, promoting distributional shifts of most species. This relation between glacial and postglacial climate and species composition is used here to quantify the climate change in the Iberian Peninsula using several fossil pollen records widespread throughout the study area. We have reconstructed spatial layers (1ka interval) of January minimum temperature, July maximum temperature and minimum annual precipitation using a method based on probability density functions and covering the time period between 15ka and 3ka. In order to evaluate the spatial evolution of climate, we used a functional principal component analysis. Using a clustering method we have identified areas that share similar climate evolutions during the studied time period. The spatial reconstructions show a highly dynamic pattern in accordance with climatic events. The four 1This study is presently submitted to an international journal listed in the science-citation index: Tarroso, P., Carrión, J., Dorado-Valiño M., Queiroz, P., Santos, L., Valdeolmillos-Rodríguez, A., Alves, P. C., Brito, J. C., & Cheddadi, R. Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. 90 FCUP Reconstruction of past Iberian climate cluster areas we found exhibit different climate evolution over the studied period. The clustering scheme is compatible with multiple refugial areas. 4.1.2 Introduction The pattern of present biodiversity is the result of a dynamic process driven at a broad temporal scale by geological events and climatic oscillations. Currently, several threats to biodiversity are identified and the largest implications are distributional shifts (Parmesan & Yohe, 2003) and, more dramatically, species extinction (Hewitt, 2000). The biodiversity hotspots retain high levels of endemism and are considered as the best candidates for preserving species diversity in the future (Myers et al., 2000). The Mediterranean basin was identified as a biodiversity hotspot and palaeoenvironmental studies showed that some areas within this region have played the role of refugia to diverse ecosystems over several hundreds of millenia (Wijmstra, 1969; Wijmstra & Smith, 1976; Van der Wiel & Wijmstra, 1987a,b; Tzedakis et al., 2002). Often, those areas where species have persisted during glacial periods are referred to as glacial refugia (Bennett & Provan, 2008; Carrión et al., 2010a; Hewitt, 2000; Hu et al., 2009; MacDonald et al., 2008; Willis et al., 2010) and the high levels of diversity found at species level is corroborated at molecular level (Hewitt, 2000; Petit et al., 2003). Understanding the past processes that affected biodiversity is fundamental for species conservation decisions facing the future expected global climate changes (Anderson et al., 2006; Willis et al., 2010). Refugia have been generally defined based on species survival with an evident relationship to climatic components (Hewitt, 2000; Bennett & Provan, 2008; Cheddadi & Bar-Hen, 2009; Médail & Diadema, 2009), nevertheless, the term has been used recently with multiple definitions (Bennett & Provan, 2008; Ashcroft, 2010). The biological definition of refugia is related to the physiological limits of a species that under an increasingly stressing environment are forced to shift to suitable areas. Palaeoenvironmental data and molecular analysis have proven useful to locate species ancestral distributions and migration routes (Petit et al., 2003; Cheddadi et al., 2006). However, the locations and extension range of putative refugia still lack spatial consensus and a quantification of its dynamic nature. Fossil pollen is an appropriate proxy for quantitative reconstruction of past climate variables (Webb et al., 1993; Cheddadi et al., 1997; Guiot, 1997; Davis et al., 2003; Cheddadi & Bar-Hen, 2009). Using proxy data to derive a definition of refugia in terms of climate, which has an obvious spatial context, may provide a suitable answer to this problem. Climate oscillations in Europe during the last 15,000 years exhibited latitudinal FCUP 91 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. and longitudinal variations (Cheddadi et al., 1997; Davis et al., 2003; Roucoux et al., 2005; Cheddadi & Bar-Hen, 2009; Carrión et al., 2010a). During the last glacial maximum (LGM), several species found refugia in the southern peninsulas (Hewitt, 2000; Tzedakis et al., 2002; Petit et al., 2003; Weiss & Ferrand, 2007; Bennett & Provan, 2008; Hu et al., 2009; Médail & Diadema, 2009; Ohlemüller et al., 2012). The Iberian Peninsula, with a milder climate than the northern European latitudes (Renssen & Isarin, 2001; Carrión et al., 2010a; Pérez-Obiol et al., 2011) served as a refugium to several species that persisted in this area during the LGM. The current patterns of high diversity in the Iberian Peninsula derive partially from this role during harsh glacial conditions and it is the main reason why this peninsula is considered as an important part of the Mediterranean biodiversity hotspot (Médail & Quézel, 1999; Cox et al., 2006). Although it may be seen globally as static refugia, the vegetation and climate dynamics in Iberia reveal a quite complex picture. The abrupt climate changes that occurred after the LGM (Heinrich event 1, Bølling-Allerød interstadial and the Younger Dryas) did not impact uniformly across the peninsula but rather with distinct velocities and patterns as shown by palynologic data (Roucoux et al., 2005; Naughton et al., 2007; Pérez-Obiol et al., 2011) with consequences on biodiversity patterns as described by molecular studies (Branco et al., 2002; Godinho et al., 2006; Weiss & Ferrand, 2007; Miraldo et al., 2011). The Iberian Peninsula is, therefore, a unique area to study the climate dynamics during the late-Quaternary, with immense variation and of paramount importance in terms of biodiversity conservation. Our main objective in this study is to define areas within the Iberian Peninsula (Balearic Islands included) that share similar climate evolution and susceptible to serve as a potential refugium area to temperate trees. We reconstructed spatially three climate variables and quantified their changes between 15 ka and 3 ka BP, with a 1,000 years interval. Then, using robust statistical methods, we defined geographical areas that have undergone similar climate changes and analysed their spatial dynamics throughout the late Quaternary. 4.1.3 Methods The area for the spatial reconstruction extends throughout the land area of the Iberian Peninsula and the Balearic islands (Fig. 4.1). The method used to produce past climate grids is based on probability density functions (PDF) and requires both fossil pollen records and full distribution of modern plant taxa. The raw fossil pollen data were gathered from author’s contribution and from the European Pollen Database ( www.europeanpollendatabase.net ). We checked each site to fit a quality criteria 98 FCUP Reconstruction of past Iberian climate from -10.6oC to -4.3oC, Tjul ranging between 20.3oC and 21.2oC and the wettest with Pmin between 32mm and 41mm. The cluster C2 (26% of the total area) occupies most of the southern Iberian Plateau. It holds the warmest values for Tjan varying between -4.8oC and 1.0oC and Tjul between 26.1oC and 27.7oC. The dissimilarities between clusters C3 and C4 (each 26% of the total area) occur in Tjan, having the former higher temperatures (-4.9oC to -0.1oC) than the latter (-6.6oC to -1.5oC). Pmin is similar within the C2 and C4 (19 to 25mm and 19 to 26mm, respectively) while C3 had higher values (23 to 28mm). Fig. 4.4 – Hierarchical cluster analysis of the functional PCA components of Tjan, Tjul and Pmin in the last 15ka found in the study area. The top dendrogram represents the size of the clusters of similar climate evolution and the relations between them. Numbers correspond to the cluster. FCUP 99 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. 4.1.5 Discussion The climate of the last 15 ka in the Iberian Peninsula was dynamic, with oscillations of temperature and precipitation that had a major impact on the location, extent and evolution of the glacial refugia during the post-glacial period. Nonetheless, the reconstructed overall trend is a noticeable warming after the 15 ka that results from the increase of the summer insolation in the northern hemisphere (Berger, 1978). An evident pattern that strikes from the results presented here is the homogeneous Tjan throughout the Iberian Peninsula at 15ka, with 2oC difference between clusters. This pattern evolved towards a widening of the temperature range during the Holocene. Despite the stability within each cluster, Tjul and Pmin markedly divide the peninsula since 15ka (Fig. 4.5, appendix C.1). Fossil pollen data provide a record of vegetation changes which constitutes a valuable proxy for reconstructing past climate changes. The major constrain we found Fig. 4.5 – Minimum and maximum temperatures of January and July, respectively, and minimum annual precipitation during the last 15 ka. The solid line represents the average climate in the study area. The remaining lines are the average of each cluster found: C1: short dash line; C2: dotted line; C3: dash-dot line and C4: long dashed line. 100 FCUP Reconstruction of past Iberian climate in the Iberian Peninsula was the low number of sequences available according to our quality criteria for a robust spatial climate reconstruction, both in terms of sampling resolution and number of C14 dates. Nevertheless, the resulting climate reconstruction is consistent with other studies (Davis et al., 2003; Cheddadi & Bar-Hen, 2009), placing the Iberian Peninsula in a particular case in Europe in terms of sensitivity to climatic stadials and interstadials like the Oldest Dryas (OD), Bølling-Allerød interstadial (BA) and Younger Dryas (YD). This study shows that Tjan and Tjul had different patterns over the last 15ka. We identify areas where climate had different patterns and also depict regions of similar climate conditions throughout the last 15 ka. Although we have reconstructed climate with a 1 ka period, the climate stadials and interstadials have a spatial imprint that is visible in the spatial time-series. The OD (∼18 to 14.7 ka) is characterized in Iberia by a vegetation changes compatible with cold and humid conditions followed by a warming trend (Naughton et al., 2007). The OD is followed by the warmer BA (∼14.7 to 12.9 ka). Our results show a similar pattern, with colder condition between 15 ka and 14 ka, followed by a warming trend until 13 ka, with different clusters exhibiting different sensibility to these fluctuations (Fig. 4.5), and showing extreme contrasts for Tjan (Figs. 4.3, 4.5). Precipitation values are low between 15 ka and 14 ka except in mountainous areas, comprised mostly in the first cluster, which remain relatively more humid than the other areas (Figs. 4.3, 4.5). As described earlier in Europe (Renssen & Isarin, 2001; Heiri et al., 2004), Tjan shows more abrupt changes than Tjul. The cold to warm transitions that occurred at ∼14.7 and 11.5 ka (Renssen & Isarin, 2001; von Grafenstein et al., 1999) in Europe had a spatial impact that is well depicted in Tjan. The increasing variability of Tjan after ∼14 ka is related to the expansion of trees from glacial refugia which have modified the albedo (Cheddadi & Bar-Hen, 2009). The BA is followed by the cold YD (∼12.9 to 11.6 ka) during which we observe a reduction of the warmer areas at 12 and 11 ka and an expansion of the more arid ones (Fig. 4.5). When Tjan exhibited a fast increase at 9 ka, Tjul remained stable. This is in agreement with Davis et al. (2003), who showed a similar Summer and Winter temperature evolution for South-western Europe. Nevertheless, a closer look at the scale of Iberia allows to observe that some areas had different vegetation compositions (Carrión et al., 2010a) and vegetation is known to respond to both long term climate trends as well as to abrupt changes (Roucoux et al., 2005). The Holocene warm period (approximately between 8.2 and 5.6 ka, depending on where in Europe) is characterized by increasing summer temperatures (Seppä & FCUP 101 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. Birks, 2001), being more evident in Northern Europe and the Alps and simultaneous with a cooling at lower latitudes (Davis et al., 2003). Our results show this dichotomy within the Iberian Peninsula where two areas (clusters 2 and 3) exhibit a warming between 9 and 5 ka while the other areas record a cooling pattern. This cooling at 8 and 7 ka is likely to be related to the impact of the ∼8.2 ka cold event (Heiri et al., 2004; von Grafenstein et al., 1999). Concerning the precipitation, there is evidence of a wetter climate between 9 ka and 6 ka which confirms what was previously known for the southern European lowlands (Cheddadi et al., 1997). The behaviour of the reconstructed variables after 5ka is likely to be influenced by non-natural ecosystem changes due to human activities such as the forest degradation that begun in lowlands and later in mountainous areas (Carrión et al., 2010b), increasing the difficulty of its interpretation. These human impacts add confounding effects in the fossil pollen record and may lead to reconstructed lower temperatures at 5 ka. On the other hand, human impact at larger scales, capable of leaving noticeable imprints on landscape were likely to happen later (Carrión et al., 2010b) and, furthermore, there are evidences of a cooling and drier stage after 5 ka, marking the end of the Holocene warm period in Europe (Seppä & Birks, 2001), and particularly in the Iberian Peninsula (Dorado-Valiño et al., 2002). The northern areas of the Iberian Peninsula had high values of precipitation during the last 15 ka. An interesting pattern arises at 9 ka, in the mid and lower latitudes, where we observe a shift of precipitation intensity. Mid-latitudes registered the highest differences in precipitation, with a beginning of Holocene increasingly humid. On the other hand, the southernmost latitudes became increasingly arid, with the exception of the 3 ka time slice. This could be due to the increasing human activity with a greater impact on the ecosystems (Carrión et al., 2010b) and their feedback on climate (Pielke et al., 1998). The variability of the climate change during the last 15ka in the Iberian Peninsula had an obvious impact on the presence of putative refugia, migrating pathways of species from these refugia and on the overall recolonisation processes during the postglacial period within the Iberian Peninsula. During this period, climate favoured migrations and expansion processes that culminated in secondary contacts for several lineages previously isolated in patches of suitable habitat (Branco et al., 2002; Godinho et al., 2006; Weiss & Ferrand, 2007; Miraldo et al., 2011). The clustering scheme (Fig. 4.4) is consistent with the molecular evidence of a network of putative refugia within Iberia (Weiss & Ferrand, 2007). Refugia have been associated with climate and habitat stability, with both playing complementary roles (Ashcroft, 2010). 102 FCUP Reconstruction of past Iberian climate However, as shown by large scale landscape analysis (Carrión et al., 2010a,b) and climate reconstructions (Davis et al., 2003; Cheddadi & Bar-Hen, 2009), both have a strong dynamic nature in the Iberian Peninsula, and promoted the formation of patches of suitable habitat during harsh conditions. The highly structured population that many species exhibit in the Iberian Peninsula have contributed decisively to the idea of refugial diversity (Hewitt, 2000; Weiss & Ferrand, 2007). Overall, the information included in the multidimensional climate data allowed us to define putative refugia as areas that shared a similar climate evolution during the late-Quaternary and where temperature and precipitations are suitable to support the survival of temperate trees. The range of climate changes over the past 15,000 years has largely affected the recolonization process of species. The defined clusters can be associated with potential isolation or dispersal events of species throughout the studied time span. Particularly, the third cluster (Fig. 4.4) includes areas that have already been described as glacial refugia for several animal and plant species (Weiss & Ferrand, 2007, see chapter 5 for a review of refugia in Iberian Peninsula). In the area represented by this cluster, the climate evolved with a lower amplitude (∼4oC) and was less sensitive to extreme fluctuations than the other clusters, which is compatible with the persistence of species in these areas. The southern plateau, mostly comprised in the second cluster (Fig. 4.4), recorded also mild conditions which are often associated with southern refugia but a rapid feedback to late-Quaternary events may have prevented persistence or recolonisation processes. 4.1.6 Conclusions The reconstruction of past climates using biological data is an invaluable resource for the study of dynamics of glacial refugial areas. Although there is a limited number of available sites and time range coverage, the spatial combination of fossil pollen data provides a continuous record with a climate signal that can be translated into spatially explicit analysis of climate dynamics. The reconstructed climate variables for the post-glacial period show different patterns of evolution but clearly marked by the lasting impact of episodic climatic events. The Iberian Peninsula had areas that shared similar climate evolution during the lateQuaternary. Some areas that we have identified as potential glacial refugia are consistent with those areas where genetic diversity is highest and which are often considered as refugial areas for several animal and plant species. The analysis of these areas and the related climate provides new insights about the dynamics of refugia through time and space which helps a better understanding of FCUP 103 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. the evolution of biodiversity hotspots both at the species and the intraspecific levels. Thus, such study on the Iberian Peninsula has an obvious interest for conservation issues, especially under the expected future climate change. 4.1.7 Acknowledgments PT was funded with a PhD grant (SFRH/BD/42480/2007) and JCB has a contract (Programme Ciência 2007), both from Fundação para a Ciência e Tecnologia. JC contribution was funded by the project Paleoflora y Paleovegetación ibérica, Plan Nacional de I+D+i, Ref. CGL-2009-06988/BOS. LS acknowledges the contribution of M. C. Freitas and C. Andrade (University of Lisbon) who provide the cores. The authors would like to acknowledge all contributors of the European Pollen Database and the Global Biodiversity Information Facility for making their datasets publicly available to the scientific community. We are very grateful to Basil Davis, for his kind support and comments. We also thank William Fletcher and Maria Sanchez-Goñi for data contribution and comments, and also Penélope González-Sampériz contributions. Their contributions greatly improved the quality of the manuscript. This is an ISEM-contribution noxx-xxxx 4.1.8 References Anderson, N.J., Bugmann, H., Dearing, J.A. & Gaillard, M.J. (2006) Linking palaeoenvironmental data and models to understand the past and to predict the future. Trends in Ecology & Evolution 21, 696–704. Ashcroft, M.B. (2010) Identifying refugia from climate change. Journal of Biogeography 37, 1407–1413. Bennett, K. & Provan, J. (2008) What do we mean by ‘refugia’? Quaternary Science Reviews 27, 2449–2455. Berger, A. (1978) Long-term variations of caloric insolation resulting from the Earth’s orbital elements. Quaternary Research 9, 139–167. Branco, M., Monnerot, M., Ferrand, N. & Templeton, A.R. (2002) Postglacial dispersal of the European rabbit (Oryctolagus cuniculus) on the Iberian Peninsula reconstructed from nested clade and mismatch analyses of mitochondrial DNA genetic variation. Evolution 56, 792–803. Carrión, J.S., Fernández, S., González-Sampériz, P., Gil-Romera, G., Badal, E., Carrión-Marco, Y., López-Merino, L., López-Sáez, J.A., Fierro, E. & Burjachs, F. 104 FCUP Reconstruction of past Iberian climate (2010a) Expected trends and surprises in the Lateglacial and Holocene vegetation history of the Iberian Peninsula and Balearic Islands. Review of Palaeobotany and Palynology 162, 458–475. Carrión, J., Fernández, S., Jiménez-Moreno, G., Fauquette, S., Gil-Romera, G., González-Sampériz, P. & Finlayson, C. (2010b) The historical origins of aridity and vegetation degradation in southeastern Spain. Journal of Arid Environments 74, 731–736. Cheddadi, R., Yu, G., Guiot, J., Harrison, S. & Prentice, I.C. (1997) The climate of Europe 6000 years ago. Climate Dynamics 13, 1–9. Cheddadi, R. & Bar-Hen, A. (2009) Spatial gradient of temperature and potential vegetation feedback across Europe during the late Quaternary. Climate Dynamics 32, 371–379. Cheddadi, R., Vendramin, G.G., Litt, T., François, L., Kageyama, M., Lorentz, S., Laurent, J.M., de Beaulieu, J.L., Sadori, L., Jost, A. & Lunt, D. (2006) Imprints of glacial refugia in the modern genetic diversity of Pinus sylvestris.Global Ecology and Biogeography 15, 271–282. Cox, N., Chanson, J. & Stuart, S. (2006) The status and distribution of reptiles and amphibians of the Mediterranean Basin. IUCN, Gland, Switzerland and Cambridge, U.K. Davis, B.A.S., Brewer, S., Stevenson, A.C., Guiot, J. & Data Contributors (2003) The temperature of Europe during the Holocene reconstructed from pollen data. Quaternary Science Reviews 22, 1701–1716. Dorado-Valiño, M., Rodríguez, A.V., Zapata, M.B.R., García, M.J.G. & Gutiérrez, I.D.B. (2002) Climatic changes since the Late-glacial/Holocene transition in La Mancha Plain (South-central Iberian Peninsula, Spain) and their incidence on Las Tablas de Daimiel marshlands. Quaternary International 93-94, 73–84. Furrer, R., Nychka, D. & Sain, S. (2012) fields: Tools for spatial data. URL http: //CRAN.R-project.org/package=fields . Godinho, R., Mendonça, B., Crespo, E.G. & Ferrand, N. (2006) Genealogy of the nuclear beta-fibrinogen locus in a highly structured lizard species: comparison with mtDNA and evidence for intragenic recombination in the hybrid zone. Heredity 96, 454–463. FCUP 105 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. Guiot, J. (1997) Palaeoclimatology: Back at the last interglacial. Nature 388, 25–27. Heiri, O., Tinner, W. & Lotter, A.F. (2004) Evidence for cooler European summers during periods of changing meltwater flux to the North Atlantic. Proceedings of the National Academy of Sciences of the United States of America 101, 15285–15288. Hewitt, G. (2000) The genetic legacy of the Quaternary ice ages. Nature 405, 907– 913. Hicks, S. (2006) When no pollen does not mean no trees. Vegetation History and Archaeobotany 15, 253–261. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G. & Jarvis, A. (2005) Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25, 1965–1978. Hu, F.S., Hampe, A. & Petit, R.J. (2009) Paleoecology meets genetics: deciphering past vegetational dynamics. Frontiers in Ecology and the Environment 7, 371–379. Keitt, T.H., Bivand, R., Pebesma, E. & Rowlingson, B. (2012) rgdal: Bindings for the Geospatial Data Abstraction Library. URL http://CRAN.R-project.org/package= rgdal . Kühl, N., Gebhardt, C., Litt, T. & Hense, A. (2002) Probability Density Functions as Botanical-Climatological Transfer Functions for Climate Reconstruction. Quaternary Research 58, 381–392. Laurent, J.M., Bar-Hen, A., François, L., Ghislain, M. & Cheddadi, R. (2004) Refining vegetation simulation models: from plant functional types to bioclimatic affinity groups of plants. Journal of Vegetation Science 15, 739–746. MacDonald, G., Bennett, K., Jackson, S., Parducci, L., Smith, F., Smol, J. & Willis, K. (2008) Impacts of climate change on species, populations and communities: palaeobiogeographical insights and frontiers. Progress in Physical Geography 32, 139–172. Médail, F. & Diadema, K. (2009) Glacial refugia influence plant diversity patterns in the Mediterranean Basin. Journal of Biogeography 36, 1333–1345. Médail, F. & Quézel, P. (1999) Biodiversity hotspots in the Mediterranean Basin: setting global conservation priorities. Conservation Biology 13, 1510–1513. 106 FCUP Reconstruction of past Iberian climate Miraldo, A., Hewitt, G.M., Paulo, O.S. & Emerson, B.C. (2011) Phylogeography and demographic history of Lacerta lepida in the Iberian Peninsula: multiple refugia, range expansions and secondary contact zones. BMC Evolutionary Biology 11, 170. Myers, N., Mittermeier, R., Mittermeier, C., Da Fonseca, G. & Kent, J. (2000) Biodiversity hotspots for conservation priorities. Nature 403, 853–858. Naughton, F., Sanchez Goñi, M., Desprat, S., Turon, J.L., Duprat, J., Malaizé, B., Joli, C., Cortijo, E., Drago, T. & Freitas, M. (2007) Present-day and past (last 25000 years) marine pollen signal off western Iberia. Marine Micropaleontology 62, 91–114. Ohlemüller, R., Huntley, B., Normand, S. & Svenning, J.C. (2012) Potential source and sink locations for climate-driven species range shifts in Europe since the Last Glacial Maximum. Global Ecology and Biogeography 21, 152–163. Parmesan, C. & Yohe, G. (2003) A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42. Pebesma, E.J. (2004) Multivariable geostatistics in S: the gstat package. Computers & Geosciences 30, 683–691. Pérez-Obiol, R., Jalut, G., Julia, R., Pelachs, A., Iriarte, M.J., Otto, T. & HernandezBeloqui, B. (2011) Mid-Holocene vegetation and climatic history of the Iberian Peninsula. The Holocene 21, 75–93. Petit, R.J., Aguinagalde, I., de Beaulieu, J.L., Bittkau, C., Brewer, S., Cheddadi, R., Ennos, R., Fineschi, S., Grivet, D., Lascoux, M., Mohanty, A., Müller-Starck, G., Demesure-Musch, B., Palmé, A., Martín, J.P., Rendell, S. & Vendramin, G.G. (2003) Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300, 1563–1565. Pielke, R.A., Avissar, R., Raupach, M., Dolman, A.J., Zeng, X. & Denning, A.S. (1998) Interactions between the atmosphere and terrestrial ecosystems: influence on weather and climate. Global Change Biology 4, 461–475. R Development Core Team (2012) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Ramsay, J.O., Wickham, H., Graves, S. & Hooker, G. (2012) fda: Functional Data Analysis. Ramsay, J. & Silverman, B. (2005) Functional data analysis. Springer, New York, 2nd edn. FCUP 107 Spatial climate dynamics in the Iberian Peninsula since 15 000 Yr BP. Renssen, H. & Isarin, R.F.B. (2001) The two major warming phases of the last deglaciation at 14.7 and 11.5 ka cal BP in Europe: climate reconstructions and AGCM experiments. Global and Planetary Change 30, 117–153. Roucoux, K., Abreu, L.D., Shackleton, N.J. & Tzedakis, P.C. (2005) The response of NW Iberian vegetation to North Atlantic climate oscillations during the last 65kyr. Quaternary Science Reviews 24, 1637–1653. Seppä, H. & Birks, H. (2001) July mean temperature and annual precipitation trends during the Holocene in the Fennoscandian tree-line area: pollen-based climate reconstructions. The Holocene 11, 527–539. Tzedakis, P.C., Lawson, I.T., Frogley, M.R., Hewitt, G.M. & Preece, R.C. (2002) Buffered tree population changes in a quaternary refugium: evolutionary implications. Science 297, 2044–2047. Van der Wiel, A.M. & Wijmstra, T.A. (1987a) Palynology of the 112.8-197.8 m interval of the core Tenaghi Philippon III, middle Pleistocene of Macedonia. Review of Palaeobotany and Palynology 52, 89–108. Van der Wiel, A.M. & Wijmstra, T.A. (1987b) Palynology of the lower part (78-120 m) of the core Tenaghi Philippon II, Middle Pleistocene of Macedonia, Greece. Review of Palaeobotany and Palynology 52, 73–88. von Grafenstein, U., Erlenkeuser, H., Brauer, A., Jouzel, J., Johnsen, S.J. & von Grafenstein, U. (1999) A Mid-European Decadal Isotope-Climate Record from 15,500 to 5000 Years B.P. Science 284, 1654–1657. Webb, R.S., Anderson, K.H., Webb, T. & Others (1993) Pollen response-surface estimates of late-Quaternary changes in the moisture balance of the northeastern United States. Quaternary Research 40, 213–227. Weiss, S. & Ferrand, N. (eds.) (2007) Phylogeography of Southern European Refugia. Springer, Dordrecht. Wijmstra, T.A. & Smith, A. (1976) Palynology of the middle part (30–78 metres) of the 120 m deep section in northern Greece (Macedonia). Acta Botanica Neerlandica 25, 297–312. Wijmstra, T. (1969) Palynology of the first 30 m. of a 120m. deep section in northern Greece. Acta Botanica Neerlandica 18, 511–527. 114 FCUP Spatial dynamic patterns in the Iberian Peninsula at two scales and intercepted was forced at zero. This procedure allows to adequate climate projections to current climate when using different sources, maintaining a similar degree of variance. The obtained calibration values were used to adjust all successive decades, generating more reliable estimates of the variables. Ecological niche modelling We used species and climate data as input variables to assess the niche of reptiles and amphibians in the Iberian Peninsula by means of ecological niche modelling. The chosen algorithm was the Artificial Neural Networks (ANN) implemented in the software Simapse v1.1 (Tarroso et al., 2012). The ANN are a machine-learning algorithm known to be very efficient in detecting patterns maintained by non-linear relationships between variables, even in the presence of some noise, as is often the case of ecological models (Olden et al., 2008; Tarroso et al., 2012). ANN attractive features are reflected in the extensive application in ecology, ranging from simple ecological models (Lek et al., 1996), to algorithm ensemble approaches in conservation studies (Carvalho et al., 2010), hybrid zone analysis (see Section 5.2) and climate reconstruction (Cheddadi et al., 1997). It consists on a layered structure of interconnected processing nodes (i.e. artificial neurons) that is able to learn patterns during the training stage, by adjusting the connections between nodes with the purpose of minimizing the output error after each iteration. When fed with species data and ENVs, it will be able to predict presence probability to other locations with the same set of ENVs. We used the current climate ENVs to assess the niche of the 30 selected species. Parameters were optimized heuristically and maintained transversally to all studied taxa. These include 5,000 iterations (with cycles of 1 internal iteration) after 15 burn-in iterations and a structure of the hidden layer with 6 nodes. Learning rates were adapted to presence data complexity (Table D.1): difficulty to model increases with data clustering and smaller learning rates are used (Wilson & Martinez, 2001). The number of pseudo-absences was calculated based on 20% of the available locations without presence data. We built a consensus model after running Simapse for 10 replicates using random sub-sampling. All other values were maintained at Simapse’s default. The final probability model was converted to binary using the 10th percentile value of the model predictions for the presences. Although this threshold classifies 10% of presences as false negatives, it provides accurate potential distributions (Brito et al., 2009). Models for each species were evaluated with receiver operating characteristic curve (ROC) and the respective area under the curve (AUC). The niche models FCUP 115 Velocity of biodiversity change: a case study in a global hotspot using current climate data were projected to the 13 layers of past climate and to 10 decades of predicted future climate, based on 8 combinations of models and scenarios. The projected models of the 30 species, for each projection layer were added to obtain projections of species richness per time layer. Velocity calculation Current climate models were used to assess the species niches and velocities were calculated based on species richness of projected models. Two different velocities measures (Fig. 5.1) were processed for past and future projections independently: 1) the velocity of pixel value maintenance (Vm) sums the minimum distance that the richness value of each pixel has to travel from one layer to the next layer in time and is divided by the total time; and 2) the velocity of richness shifts (Vs) sums the minimum distance a value of richness has to travel in each layer, arranged by chronological order, and divided by the total time. Both velocities have different purposes: while the Vm measures how far a value of richness has to travel after reaching a location, the Fig. 5.1 – Illustration of the velocities metrics. The velocity of pixel maintenance (Vm) corresponds to the cell value maintenance throughout time (a). It is measured the shortest distance between one cell at time t and the nearest cell at the next time layer (t+1) holding an equivalent value of species richness. The velocity of richness shift (Vs), corresponds to the shift of pixel value in space (b). The distance here is the tracking of the species richness value at time t in all subsequent time layers (t+n), irrespective of the original cell location. Distances measure in km are summed and divided by years to calculate the velocities in km.yr-1 units. Species values are tracked with a tolerance value of 1 species. 116 FCUP Spatial dynamic patterns in the Iberian Peninsula at two scales Vs describes how far did each value of richness travelled from its original location. The tracking of species richness has a tolerance of one species, meaning that the distance is found by locating the value with the richness plus or minus one. This procedure allows the possibility of positive infinite velocities: if a value being tracked is not found in the next time layer, an infinite distance is attributed. Geographical coordinates were transformed using a gnomonic projection, adjusted to the centroid of the Iberian Peninsula area, for calculation of velocities in standard units (km.yr−1). Since different future climate scenarios have different parametrizations of physical climate properties and impacts of society behaviours, multi-model averaging may result in deviations and physically impossible scenarios (Knutti et al., 2010). Thus, we projected each combination of model and scenario individually and summarized the results by calculating the minimum and the maximum predicted by the four scenarios, creating a best and worst case scenario for each climate model. All analysis were performed and automated with R language (R Development Core Team, 2012) using the packages ’rgdal’ (Keitt et al., 2012) with the Geospatial Data Abstraction Library (GDAL Development Team, 2011), ’ROCR’ (Sing et al., 2009) and ’RSQLite’ (James, 2011) to interact with an SQLite ( www.sqlite.org ) database storing all the data. Additional processing and automation were done using python ( www.python.org ) and shell scripts. Means are reported with one standard deviation. 5.1.4 Results The models produced for the present climatic conditions had a high average AUC (0.850 ±0.114, appendix Table D.1 and Fig. D.1), indicating high accuracy. However, some species, particularly those with wide distributions, reported low AUC values. There is a general increasing trend between achieved AUC and the proportion of clustered absences (i.e. non-presence locations; appendix D.2), relating low values of AUC with higher scattering of absences’ distribution. The threshold chosen for each model (appendix D.1) classifies some presences located at the edges of the predicted distributions as false negative (appendix D.3). The calibration of the future climate datasets with the present data resulted in less drastic changes of climate (appendix D.4). For 2100, calibrated predictions points to a Tjul ∼0.7oC lower than original models, in average. For Tjan is shown an increase of ∼0.5oC, while for Pmin, calibrated models indicate more ∼5.6mm of precipitation, for average values. Both calculated velocities generated similar spatial patterns (Fig. 5.2). The Vm have slightly higher average values than Vs (Vm = 0.0098 ±0.0090 km.yr−1and Vs FCUP 117 Velocity of biodiversity change: a case study in a global hotspot = 0.0093 ±0.0086 km.yr−1), indicating that pixel values travel longer distances every temporal layer than the original value has to shift during the 12,000 years. Major differences between velocities are located at the western part of the Iberian Peninsula where the values of each pixel are preserved through time, but original species richness (at 15,000 years BP) travelled long distances. Central area of the Peninsula reveal high velocities of change, with pockets of low velocities within higher altitudes. Higher velocities and average richness are generaly located between 500 and 1,000 m altitude (Fig. 5.3), however, with high variance in species richness through the studied time range. The mountain areas higher than 1,000m shown a trend to decrease to near zero velocities and with low changes in species richness. Areas with higher altitudes (>1,000m) had a low to average richness during the studied period, with stable values of species number (appendix D.5). Most diversity was found below 1250 m and, particularly, below 500 m, with high variance in the number of species. The average richness had an absolute decrease from 15,000 (17.3 ±3.5) to 3,000 years BP (14.7 ±3.9). However, this process had an oscillating nature and the minimum was registered at 9,000 years BP (13.7 ±3.0; Fig. 5.4). Locations with highest variance in species number through the studied time period were those that registered highest number of species at the oldest age (Fig. 5.4). The general result of the velocities based on predicted future climate layers is an increase of two orders of magnitude in relation to those of the past (Fig. 5.5). The average Vs is higher than the Vm, however HadCM3 model (minima: Vm = 0.6891 Fig. 5.2 – Velocities of species composition change: Velocity of pixel maintenance (Vm) and velocity of richness shift (Vs) estimated for the period between 15,000 and 3,000 years BP. 118 FCUP Spatial dynamic patterns in the Iberian Peninsula at two scales ±0.9566 km.yr−1< Vs = 1.371 ±1.4446 km.yr−1; maxima: Vm = 1.5479 ±1.6452 km.yr−1< Vs 2.0168 ±2.4781 km.yr−1) exhibit higher differences than CSIRO2 (minima: Vm = 0.5084 ±0.6249 km.yr−1< Vs = 0.5632 ±0.6408 km.yr−1; maxima: Vm = 1.1754 ±0.9742 km.yr−1< Vs = 1.392 ±1.2041 km.yr−1). The climate model HadCM3 also shows a higher number of cells with infinite values (Fig. 5.5), and, minimum values between scenarios also display infinite velocities, predicting that values of richness in these locations will be lost in all available scenarios. The value of richness is predicted to be most unstable at elevations between 750 and 1000 m and near sea level for the next 100 years, and also values at these locations will have to move longer distances (Fig. 5.6). Nevertheless, the two climate models provide different effects of both velocities calculated per altitude. The impact of HadCM3 model is seen more generally at altitudes bellow 1750 m, with many locations reaching high velocities. The highest values of average species richness for the next century are restricted to middle altitudes, and with high levels of variance, indicating changes in number of species at these locations (appendix D.6). The average values of species richness were predicted to decrease in all combinations of models and scenarios (Fig. 5.7). The scenario A1FI showed the highest differences (CSIRO2: 12.7 ±4.0 km.yr−1, 9.4 ±2.6 km.yr−1; HadCM3: 12.8 ±4.0 km.yr−1, 6.2 ±1.4 km.yr−1, for 2010 and 2100 values, respectively for each model), followed by the A2 (CSIRO2: 13.0 ± 4.2 km.yr−1, 10.0 ±3.0 km.yr−1; HadCM3: 12.8 ±4.0 km.yr−1, 7.4 ±1.6 km.yr−1). The scenarios B1 (CSIRO: 12.9 ±4.2 km.yr−1, 11.6 ±3.2 km.yr−1; HadCM3: 12.8 Fig. 5.3 – Relation between velocities of pixel maintenance (Vm) and richness shift (Vs) with altitude in the Iberian Peninsula. Circle size corresponds to the variance of species richness (n. of species/cell) registered for each cell in all time layers. FCUP 119 Velocity of biodiversity change: a case study in a global hotspot ±4.0 km.yr−1, 9.6 ±2.7km.yr−1) and B2 (CSIRO: 12.8 ±4.2 km.yr−1, 11.0 ±3.0 km.yr−1; HadCM3: 12.8 ±4.1 km.yr−1, 9.1 ±2.5 km.yr−1) were very similar, but with differences in severity. The locations predicted to have the highest variance in number of species for the next century are those with a high number of species currently. 5.1.5 Discussion The issue of model transferability has been extensively discussed on scientific literature (Randin et al., 2006; Peterson et al., 2007b; Phillips, 2008; Wenger & Olden, 2012; Heikkinen et al., 2012). Although the question raised by Thuiller et al. (2008) about the relation of model complexity and accuracy is still unanswered, there are arguments pointing to both directions (Stockwell, 2006; Peterson, 2007). Simple models, however, are more prone to capture a well conserved part of the niche, despite being less effective capturing the realized niche, due to the high generalization resulting from a low dimensional niche assessment. The overfitted models were shown to have a negative effect on their transferability (Randin et al., 2006). The framework used here seems appropriate for the study of niche extrapolation to other time frames: 1) the Artificial Neural Networks (ANN) were shown before as an appropriate algorithm to study niche transference (Wenger & Olden, 2012; Heikkinen et al., 2012); 2) although the over parametrization of ANN may require time to fine tune models, it allowed a good control of model overfitting; 3) using few niche limiting variables allowed to build a more relaxed model, supporting the conclusions by Rödder & Lötters (2009); and Fig. 5.4 – Evolution of species richness for the period between 15,000 and 3,000 years BP. Solid line represents the average species richness in the study area and blue band is one standard deviation. Dotted and dashed lines corresponds to the average values of locations in the first and forth quartiles of variance. 120 FCUP Spatial dynamic patterns in the Iberian Peninsula at two scales 4) using 30 species allowed to balance the species-dependent transferability (Randin et al., 2006). By using an efficient machine-learning algorithm with a low number of predictors we were able to build accurate models of the species’ distributions, with low number of either false positives (FP) and negatives (FN). The number of FN is obviously dependent on the threshold chosen for each model. The exclusion of some presences of the model is based on two assumptions: 1) species niche is also dispersal dependent (Soberón, 2007) and individuals may be found outside their fundamental niche; and, at less extent, 2) translocation of animals, either directly or indirectly, has occurred as result of human activities in Iberia (Pleguezuelos et al., 2004). These individuals are, therefore, treated as outliers for the transferable niche and are mainly located at the borders of the distribution (appendix D.3). On the other hand, the number of FP, although also dependent on the threshold, is related to the predictive ability of the model, finding potential presences outside the original locations, i.e., generalization of the model. As observed by the predictive distribution area, some species have higher tendency to be generalists, with predicted ranges extending the area where species were found (appendix D.3). Distributional restricting forces like biotic factors (e.g. presence of competitors; Soberón, 2007), may be among the underlying causes of the high difference between predicted areas and realized niche, as given by presFig. 5.5 – Velocities for predicted scenarios of future climates. The velocity of pixel maintenance (Vm) and velocity of richness shift (Vs) for climate models CSIRO2 and HadCM3 are summarized with minimum and maximum values achieved for four emissions scenarios (A1FI, A2, B1 and B2). White area represents location where infinite velocities were found. FCUP 121 Velocity of biodiversity change: a case study in a global hotspot ence’s locations. However, due to the simplicity of the predictors used in this study, we cannot assess if species with high generalization would be limited by an abiotic factor other than those used. On other cases, especially on species with wider distributions (e.g. Podarcis hispanica,Rhinechis scalaris and others), FP is inflated by scattered absences. Several factors may be creating these isolated absences, including human disturbance in the habitat or evasive species which are difficult to sample. The existence of many isolated absences was shown to be an influential factor lowering the AUC values of the models (appendix Table D.1 and Fig. D.2). This metric can provide misleading results with presence-only modelling (Lobo et al., 2007; Peterson et al., 2007a) but it is reliable when analysed with care, being a widely used metric in ecological studies (Randin et al., 2006; Peterson et al., 2007a; Svenning et al., 2008; Brito et al., 2009; Carvalho et al., 2010; Heikkinen et al., 2012). However, the extensive fieldwork supporting the atlases, which are the basis of the present work, assures a higher degree of confidence on absences data, especially when clustered, and, thus, AUC is expected to show a higher accuracy. The calculation of two velocities measuring richness maintenance and shift, have provided insightful results of the past dynamics and predicted future change in species assemblages. Other studies analysing the velocity of climate change, have based their calculations on linear temporal gradients of temperature and precipitation Fig. 5.6 – Relation between velocities and altitude for future climate predictions. Emission scenarios for CSIRO and HadCM3 climate models are summarized with minimum and maximum values achieved for velocity of pixel maintenance (Vm) and velocity of richness shift (Vs). Red dots at the right of each plot correspond to estimated infinite velocities. 122 FCUP Spatial dynamic patterns in the Iberian Peninsula at two scales changes (Loarie et al., 2009; Sandel et al., 2011). Predictions of temperature and precipitation for the next 100 years have a linear increase (Loarie et al., 2009) and, due to the short time scale, no major fluctuations are expected. Past changes, on the contrary, are usually analysed for longer time periods, capturing intense oscillations within longer trends (Zachos et al., 2001). This is the case of the warming trend after the last glacial maximum until the present days, that registered major cold and warm periods (see Section 4.1; Willis & MacDonald, 2011). Excluding these events from the analyses forces lower predictions of the climate change velocity, related to a spurious linear warming trend since the end of the glacial period (Sandel et al., 2011). The oscillating nature of the past climate warming is, therefore, eliminated by linear estimations of the temporal component of velocities. Our dataset is not possible to analyse with such velocity calculation: species richness oscillates with great variance during the 13,000 years BP studied for the past climate (Fig. 5.3, appendix D.5). Fitting a line to these data would cancel the oscillations in species richness and provide erroneous estimates of velocity. In fact, in extreme cases, predicted expansions and contractions of species ranges would generate such variation of species richness through time that a linear fitting results in a flat slope, indicating no change in species composition. Fig. 5.7 – Predicted evolution of species richness for the current century under four different emissions scenarios. Solid line corresponds to the average species richness in the Iberian Peninsula with one standard deviation (blue band). Dotted and dashed lines represent the average value of species richness for locations in the first and forth quartiles of variance. FCUP 123 Velocity of biodiversity change: a case study in a global hotspot Velocity of changes in the past The past changes on species numbers in the Iberian Peninsula were processed at different velocities. The period between 15,000 and 3,000 years BP has experienced dramatic climate oscillations, as part of the process of glacial to interglacial transition (see Section 4.1; Willis & MacDonald, 2011). The Iberian Peninsula underwent through these changes with a pronounced spatial dynamics (see Section 4.1) that are reflected in the different velocities at which species compositions operate. Most of Central Iberia underwent dramatic changes in this context, and areas of low velocity are scattered throughout the study area. Sandel et al. (2011) points to a maximum climate change velocity of 0.0148 km.yr−1in the Iberian Peninsula which is faster than the average velocity of species composition change estimated in this study. This may indicate that species assemblages and individual species, since each one was modelled independently, had some degree of resilience to climate change in the past. The climate in the Quaternary conducted species distributional shifts, finally shaping the current diversity pattern of reptiles and amphibians (Araújo et al., 2008). The combination of two velocities provides insightful results: areas that preserve quite effectively the species richness are those retaining the original value (measured by Vs) nearby its location (measured by Vm). Areas that share these conditions are found mostly scattered around the coastal line, with few exceptions in the Central System mountain range. Ohlemüller et al. (2012), based on analogous climate locations between the current climate and at the LGM, classified most of southern Iberian Peninsula as high potential source of species for European colonizations. However, as also seen by the different velocities, the climate is a general proxy for species, indicating areas with suitable conditions that may hold species. Here, we follow a nearly opposite reasoning, by using species data and by assuming that areas with slower changes in species number are indicative of climate conditions that favoured species persistence and, thus, preserved overall diversity. Our results support not only the idea of southern Iberian Peninsula as a source area (i.e. with low velocities), but we also add the northern area, with lower velocities, as potential species source for expansion and recolonization processes. The interest of the ecological niche modelling community to derive past distributions of species has gradually increased in the last few years (Waltari et al., 2007; Svenning et al., 2008; Graham et al., 2010). These tools, especially when integrated with molecular analysis, have potential to provide spatial insights of putative refugia for species (Kozak et al., 2008; Hickerson et al., 2010) under a static (e.g. Waltari et al., 2007; Svenning et al., 2008) or dynamic perspective (e.g. Graham et al., 2010). The 226 FCUP Appendix D Continued from previous page. FCUP 227 Appendix D Continued from previous page. 228 FCUP Appendix D Fig. D.5 – Relation between average species richness and altitude for past model results. Circle size is proportional to the variance of species richness throughout time (see fig. 5.3). FCUP 229 Appendix D Fig. D.6 – Velocities of change in species composition for all combinations models and emission scenarios. Velocity of pixel maintenance (Vm) and velocity of richness shift (Vs) are shown. White areas in the map correspond to infinite velocity values. Relations between altitude and velocity and also between average species richness and altitude are shown for each model/emission scenario. Circle size is proportional to the variance by location of species richness through time. 230 FCUP Appendix D Continued from previous page. FCUP 231 Appendix D Continued from previous page. 232 FCUP Appendix D Continued from previous page. Appendix E Supplementary material for section 5.2 Table E.1 – Sampled individuals with morfological identifications, cluster membership probabilities and mtDNA results. Rows marked with an asterisk are reference individuals used in Hardy-Weinberg equilibrium test. Number morphology pop1 pop2 pop3 NlaIII XspI HW 1Vipera aspis 0.03 0.962 0.008 VA VA 2Vipera aspis 0.004 0.991 0.006 VA VA 3Vipera aspis 0.005 0.986 0.009 VA VA 4Vipera aspis 0.005 0.991 0.005 VA VA 5Vipera aspis 0.007 0.983 0.01 VL VL 6Vipera aspis 0.007 0.986 0.008 VA VA 7Vipera aspis 0.015 0.979 0.006 VA VA 8Vipera aspis 0.007 0.989 0.005 VA VA 9Vipera aspis 0.212 0.78 0.007 VL VL 10 Vipera aspis 0.004 0.992 0.004 VA VA 11 Vipera aspis 0.004 0.993 0.003 VL VL 12 Vipera aspis 0.085 0.901 0.014 VL VL 13 Vipera aspis 0.015 0.939 0.046 VA VA 14 Vipera aspis 0.004 0.978 0.017 VA VA 15 Vipera aspis 0.006 0.982 0.012 VA VA 16 Vipera aspis 0.031 0.961 0.008 VA VA 17 Vipera aspis 0.058 0.936 0.007 VA VA 18 Vipera aspis 0.006 0.985 0.008 VA VA 19 0.009 0.985 0.006 VA VA 20 0.009 0.983 0.008 VA VA 21 Vipera latastei 0.984 0.008 0.008 VL VL 22 Vipera latastei 0.986 0.009 0.005 VL VL 23 Vipera latastei 0.482 0.339 0.179 VL VL 24 Vipera latastei 0.085 0.899 0.016 VL VL 25 Vipera latastei 0.991 0.004 0.005 VL VL 26 Vipera latastei 0.988 0.005 0.007 VL VL 27 Vipera latastei 0.992 0.004 0.004 VL VL 28 Vipera latastei 0.991 0.005 0.004 VL VL 29 Vipera latastei 0.969 0.02 0.011 VL VL 30 Vipera latastei 0.977 0.007 0.017 VL VL 31 Vipera latastei 0.982 0.013 0.005 VL VL 32 Vipera latastei 0.976 0.016 0.008 VL VL Continued on next page 234 FCUP Appendix E Table E.1 – Continued from previous page Number morphology pop1 pop2 pop3 NlaIII XspI HW 33 Vipera latastei 0.99 0.005 0.005 VL VL * 34 Vipera latastei 0.979 0.015 0.006 VL VL 35 Vipera latastei 0.992 0.004 0.004 VL VL 36 Vipera latastei 0.987 0.006 0.007 VL VL 37 Vipera latastei 0.973 0.021 0.006 VL VL 38 Vipera latastei 0.991 0.005 0.004 VL VL * 39 Vipera latastei VL VL 40 Vipera latastei 0.989 0.007 0.004 41 Vipera latastei 0.989 0.005 0.005 VL VL 42 Vipera latastei 0.988 0.006 0.006 VL VL 43 Vipera latastei 0.989 0.004 0.007 VL VL 44 Vipera latastei 0.985 0.008 0.007 VL VL 45 Vipera latastei 0.978 0.018 0.004 VL VL 46 Vipera latastei 0.987 0.01 0.003 VA VA 47 Vipera latastei 0.964 0.027 0.009 VL VL 48 0.991 0.004 0.006 VL VL 49 Vipera latastei 0.955 0.041 0.005 VL VL * 50 Vipera latastei 0.991 0.005 0.004 VL VL 51 Vipera latastei 0.991 0.004 0.005 VL VL * 52 Vipera latastei 0.989 0.005 0.006 VL VL * 53 Vipera latastei 0.968 0.007 0.025 VL VL 54 Vipera latastei 0.98 0.015 0.005 VL VL 55 Vipera latastei 0.975 0.008 0.018 VL VL 56 Vipera latastei 0.988 0.005 0.007 VL VL * 57 Vipera latastei 0.748 0.01 0.242 VL VL * 58 0.008 0.005 0.987 VS VS * 59 Vipera sp. 0.988 0.006 0.006 VL VL 60 Vipera sp. 0.088 0.903 0.009 VL VL 61 Vipera sp. 0.865 0.111 0.023 VL VL 62 Vipera sp. 0.429 0.566 0.005 VL VL 63 Vipera sp. 0.928 0.027 0.045 VS VS 64 Vipera sp. 0.017 0.978 0.005 VA VA 65 0.078 0.007 0.915 VL VL 66 Vipera sp. 0.033 0.962 0.005 VL VL 67 Vipera aspis 0.456 0.535 0.009 VL VL 68 Vipera aspis 0.829 0.165 0.006 VL VL 69 Vipera aspis 0.6 0.392 0.008 VL VL 70 Vipera aspis 0.005 0.992 0.003 VA VA 71 Vipera aspis 0.287 0.706 0.007 VA VA 72 Vipera aspis 0.007 0.976 0.017 VA VA * 73 Vipera aspis 0.008 0.97 0.021 VA VA 74 Vipera aspis 0.052 0.911 0.037 VS VS 75 Vipera aspis 0.013 0.983 0.004 VA VA 76 Vipera aspis 0.006 0.99 0.004 VA VA * 77 Vipera aspis 0.014 0.98 0.006 VA VA 78 Vipera aspis 0.006 0.988 0.006 VA VA 79 Vipera aspis 0.008 0.988 0.005 VA VA 80 Vipera aspis 0.01 0.985 0.005 VA VA 81 Vipera aspis 0.006 0.987 0.006 VA VA Continued on next page FCUP 235 Appendix E Table E.1 – Continued from previous page Number morphology pop1 pop2 pop3 NlaIII XspI HW 82 Vipera aspis 0.009 0.983 0.008 VA VA 83 Vipera aspis 0.041 0.95 0.01 VL VL 84 Vipera aspis 0.017 0.969 0.014 VA VA 85 Vipera aspis 0.012 0.981 0.007 VA VA 86 Vipera aspis 0.006 0.972 0.022 VA VA 87 Vipera aspis 0.036 0.953 0.011 VA VA 88 Vipera aspis 0.004 0.991 0.004 VA VA 89 Vipera aspis 0.005 0.99 0.005 VA VA 90 Vipera aspis 0.035 0.96 0.005 VA VA 91 Vipera aspis 0.032 0.959 0.008 VA VA 92 Vipera aspis 0.007 0.987 0.005 VA VA 93 Vipera aspis 0.014 0.981 0.005 VA VA 94 Vipera aspis 0.012 0.981 0.007 VA VA 95 Vipera aspis 0.006 0.99 0.004 VA VA 96 Vipera aspis 0.013 0.98 0.007 97 Vipera aspis VA VA 98 Vipera aspis 0.056 0.935 0.009 VA VA * 99 Vipera aspis 0.007 0.988 0.005 VA VA * 100 Vipera aspis 0.03 0.959 0.011 * 101 Vipera aspis 0.005 0.99 0.005 * 102 Vipera aspis 0.003 0.993 0.003 VA VA 103 Vipera aspis 0.121 0.873 0.006 VA VA 104 Vipera aspis 0.005 0.989 0.006 VA VA 105 Vipera aspis 0.021 0.975 0.004 VA VA 106 Vipera aspis 0.004 0.949 0.047 VA VA 107 Vipera aspis 0.006 0.99 0.005 VA VA 108 Vipera aspis 0.006 0.988 0.006 VA VA 109 Vipera aspis 0.004 0.992 0.004 VA VA 110 Vipera aspis 0.874 0.122 0.004 VA VA 111 Vipera aspis 0.004 0.992 0.004 VA VA 112 Vipera aspis 0.007 0.99 0.003 VA VA 113 Vipera aspis 0.004 0.989 0.007 VA VA * 114 Vipera aspis 0.088 0.908 0.004 VA VA 115 Vipera latastei VL 116 Vipera latastei 0.981 0.01 0.009 VL VL 117 Vipera latastei 0.992 0.004 0.004 VL VL 118 Vipera latastei 0.982 0.007 0.011 VL VL 119 Vipera latastei VL VL 120 Vipera latastei 0.988 0.008 0.004 VL VL 121 Vipera latastei 0.983 0.009 0.008 VL VL 122 Vipera latastei 0.982 0.014 0.004 123 Vipera latastei 0.138 0.853 0.01 VL VL 124 Vipera latastei 0.648 0.345 0.008 VL VL 125 Vipera latastei 0.99 0.006 0.005 VL VL 126 Vipera latastei 0.963 0.015 0.022 VL VL 127 Vipera latastei 0.987 0.008 0.006 VL VL 128 Vipera latastei 0.991 0.005 0.004 VL VL 129 Vipera latastei 0.946 0.037 0.017 VA VA 130 Vipera latastei 0.989 0.005 0.006 VL VL Continued on next page