Background and Aims

Climate change is expected to alter the geographic range of many plant species dramatically. Predicting this response will be critical to managing the conservation of plant resources and the effects of invasive species. The aim of this study was to predict the response of temperate homosporous ferns to climate change.

Methods

Genetic diversity and changes in distribution range were inferred for the diploid rock fern Asplenium fontanum along a South–North transect, extending from its putative last glacial maximum (LGM) refugia in southern France towards southern Germany and eastern-central France. This study reconciles observations from distribution models and phylogeographic analyses derived from plastid and nuclear diversity.

Key Results

Genetic diversity distribution and niche modelling propose that genetic diversity accumulates in the LGM climate refugium in southern France with the formation of a diversity gradient reflecting a slow, post-LGM range expansion towards the current distribution range. Evidence supports the fern's preference for outcrossing, contradicting the expectation that homosporous ferns would populate new sites by single-spore colonization. Prediction of climate and distribution range change suggests that a dramatic loss of range and genetic diversity in this fern is possible. The observed migration is best described by the phalanx expansion model.

Conclusions

The results suggest that homosporous ferns reproducing preferentially by outcrossing accumulate genetic diversity primarily in LGM climate refugia and may be threatened if these areas disappear due to global climate change.

INTRODUCTION

Understanding and forecasting the response of plant species to climatic fluctuation is one of the top priorities for current biodiversity research because of the critical need to conserve and manage natural resources and biodiversity (Bellard et al., 2012; Jay et al., 2012; Maiorano et al., 2012; D'Amen et al., 2013). Climate fluctuations are not a new phenomenon. Plants have responded to global, regional and local climate change via migration and/or adaptation since their origin (Hewitt, 2004; Corlett and Westcott, 2013). In turn, slow and/or little response to climate change (e.g. slow migration rate) increases the probability of local or global extinction. By constructing the spatio-temporal dynamics of plant response to climate change from the past, it may be possible to improve our ability to predict future changes in the range and distribution of species and their genetic diversity. Low diversity coincides with high climate change velocity – defined as a measure of the local rate of displacement of climatic conditions over the Earth's surface (Sandel et al., 2011) – that contributed to extinctions in the Quaternary period (Hewitt, 2004). Recently researchers have tried to untangle the response of plants to changing climates at the microevolutionary scale, by integrating species distribution models (SDMs) and statistical phylogeography (SPG) (Richards et al., 2007; Hickerson et al., 2010; Knowles and Alvarado-Serrano, 2010). Combining these two techniques will not only overcome their individual limitations, but will also improve our understanding of the spatio-temporal population dynamics involved (Cordellier and Pfenninger, 2009; Marske et al., 2011; Schorr et al., 2012; Schurr et al., 2012).

Recent studies using this integrative approach have emphasized the importance of migration capacity – including the discrete components of establishment and dispersal of propagules – as a major factor in the assembly of plant diversity under climate change (Cobben et al., 2011; Gugger et al., 2011; Normand et al., 2011). By studying the genetic diversity of seed plants (mainly trees) in Central and Northern Europe, these studies have recovered evidence for the complex interactions between the availability of a suitable environment through time and migration capacity (Normand et al., 2011; Alsos et al., 2012). The prediction of sensitivity of plants to climate change cannot be inferred just by using current distribution range assessments and associated climate (Thuiller et al., 2005). It is necessary to include other factors that control the spatio-temporal dynamics of the distribution of genetic diversity. The majority of studies have focused on angiosperms (in particular trees), highlighting important traits that influence colonization strategy and the putative impact of warming climates on the maintenance of genetic diversity (Normand et al., 2011; Alsos et al., 2012). Migration capacities vary considerably among plants. To date, little attention has been given to plants that are assumed to have very high migration capacities, such as homosporous ferns and bryophytes. In these plants, newly available habitats are likely to be occupied rapidly via long-distance dispersal with numerous wind-borne propagules. For this reason, the range expansion of ferns is expected to fit best with the pioneer expansion model, to which both short- and long-distance dispersal events contribute. In contrast to that, short-distance dispersal is the major process in the phalanx and stepping stone models (Nichols and Hewitt, 1994; Ibrahim et al., 1996).

In this study, we explore the putative impact of warming climates on plant species with high migration capacity, by studying the European rock fern Asplenium fontanum subsp. fontanum. Ferns are generally recognized for their high migration capacity (Wolf et al., 2001) through production of large quantities of small, wind-dispersed spores. Utilizing this mechanism, fern migration may be less restricted by dispersal limitation. This allows ferns to maintain their genetic diversity despite range shifts. However, recent studies have shown there is considerable variation in the traits that influence the migration capacity in ferns (Lloyd, 1974; de Groot et al., 2012a, b). Traits associated with fern reproductive biology include breeding systems and polyploidy (Lloyd, 1974; Wolf et al., 2001; Wubs et al., 2010). In turn, intragametophytic mating and homozygosity has been proposed to take place preferably in taxa occurring in less competitive pioneer habitats, whereas intergametophytic mating and heterozygosity is favoured in more mature (i.e. stable) habitats (Lloyd, 1974).

The focal species for this study, A. fontanum, was selected as it exemplifies the most common combinations of those traits which influence the migration capacity of ferns (Ranker and Geiger, 2008). It is a diploid (Herrero et al., 1997) derived fern, which appears to reproduce preferably via outcrossing (Hunt et al., 2009). In addition, the taxon lacks the ability to reproduce asexually. Phylogenetic analyses have demonstrated that A. fontanum is distantly related to other Asplenium species occurring in Europe (Schneider et al., 2004). Asplenium fontanum subsp. fontanum is found on shaded limestone outcrops in the mountains of South-west Europe, expanding in a more or less continuous distribution to Central-west Europe (Tutin et al., 1993). Previous studies have demonstrated that the mountains of north-east Spain and the Maritime and Western limestone Alps were putative refugia during the last glacial maximum (LGM). On the other hand, the more northern sites of A. fontanum along the Rhone valley, the Jura mountains in France, Switzerland and Germany, are likely to be the result of post-glacial colonizations (Hunt et al., 2009). The focus of this study is the latitudinal gradient in a continuous distribution from the putative LGM refugium in the Maritime Alps to the lower altitudes of the Western limestone Alps in the north, including occurrences in eastern France, southern Germany, and eastern to northern Switzerland.

To explore the potential impact of global warming on this fern, we carried out an integrative study reconciling SDMs and SPG. Species distribution models were projected on two distinct LGM scenarios and four scenarios of future climates under the A2 Intergovernmental Panel on Climate Change (IPCC) emission scenario. To enable the integration, we recognized three categories of area. Populations in the LGM refugia were predicted to remain stable under both the current and the LGM climatic conditions, whereas populations in the ‘expansion zone’ were predicted to be extinct during the LGM. Populations in these two zones should form a rather continuous range. Geographically isolated populations along the northern range of the current distribution are treated as ‘outliers’. This distinction allows the inference of the assembly of genetic diversity in the context of range expansion models incorporating or excluding the contribution of pioneer populations (Nichols and Hewitt, 1994). The population genetics of the species was explored using allozyme and chloroplast DNA sequence data.

Specifically, we ask two general questions. (1) Is there evidence for the influence of past climate change on the current geographic distribution and genetic diversity of A. fontanum? (2) Is there a threat from future climate warming to the survival of the historically accumulated genetic diversity of A. fontanum? The first specific research question (RQ1) assumes that survival populations from localities with relatively stable climate (i.e. climate refugia) retain higher genetic diversity due to long-term persistence and/or a minimal demographic change. The second research question (RQ2) assumes that range expansion/migration into newly available habitat creates a gradient of genetic diversity. To explore the second question, we relate the predicted geographic location and size of the niche to the current location and size taking into consideration the observed ability for range expansion. These research questions are linked to another two specific research questions concerning the process that enables ferns to occupy new niches. The third research question (RQ3) infers the contribution of long- and short-distance dispersal to the process of range expansion in ferns. Frequent long-distance dispersal results in genetic spatial distribution as expected under the pioneer model, whereas low frequency of successful long-distance dispersal generates genetic patterns expected under the phalanx model (Nichols and Hewitt, 1994). It has been shown that in homosporous ferns most spores are deposited close to the mother plant (Ranker and Geiger, 2008) and a higher rate of short-distance events is predicted. The fourth question (RQ4) concerns the reproductive system of this fern; here we specifically infer evidence of changing inbreeding/outbreeding ratios.

MATERIALS AND METHODS

Statistical distribution modelling

Asplenium fontanum subsp. fontanum distribution data were retrieved from the GBIF database (http://www.gbif.org, last accessed 20 February 2012) and visually examined for potentially incorrect data points. The data set was expanded by the addition of 101 of our own field records. We confirmed most of the local occurrences by consulting local floras and checklists or by communication with local experts. The frequency of misidentifications is expected to be low given the unique morphological characters of A. fontanum. For this reason, many of the issues associated with the GBIF data do not apply to this study. Some of the data points used were from specimens used to generate molecular data from Hunt et al. (2009) or newly obtained for this study (Table 1). In total, 651 unique collection localities were used for modelling, 40 of which had associated genetic information. We ignored the distribution of A. fontanum subsp. pseudofontanum (= A. pseudofontanum), which occurs in the Western Himalaya (Reichstein and Schneller, 1982).

Table 1.

Summary of the Maxent logistic output values for populations in the studied South–North transect and their genetic diversity observed based on allozymes (genetic diversity) and chloroplast DNA sequence variation

POP LAT LON BLC LGM GRO DIS n pPL Ar Ae He Ho f SR K1 K2 
F-80 49·57 7·59 0·31 0·000 OUT 634 39 0·00 1·019 1·004 0·003 0·003 –0·013 0·61 0·03 0·97 – – – 
F-128 48·62 9·785 0·27 0·000 OUT 558 15 0·07 1·043 1·010 0·009 0·009 –0·037 0·57 0·02 0·98 – – – 
F-129 48·62 9·85 0·29 0·000 OUT 558 15 0·00 1·000 1·000 0·000 0·000 0·000 0·59 0·02 0·98 – – – – 
F-130 48·63 9·85 0·35 0·000 OUT 560 14 0·07 1·066 1·040 0·026 0·014 0·458 0·56 0·12 0·88 – – – – 
F-51 47·42 8·00 0·55 0·000 OUT 396 17 0·27 1·278 1·180 0·093 0·082 0·117 0·46 0·12 0·87 – – 
F-84 47·38 7·76 0·28 0·000 OUT 390 17 0·20 1·197 1·141 0·083 0·071 0·158 0·45 0·13 0·87 – – – – 
F-34* 47·29 7·71 0·15 0·000 OUT 381 50 0·27 1·212 1·114 0·073 0·055 0·249 0·48 0·19 0·81 – – 
F-134 47·43 4·93 0·24 0·000 OUT 421 14 0·07 1·065 1·034 0·023 0·019 0·188 0·57 0·04 0·95 – – – – 
F-135 47·24 4·73 0·36 0·000 OUT 405 20 0·13 1·153 1·135 0·069 0·067 0·033 0·43 0·23 0·77 – – – – 
F-114 46·31 6·97 0·15 0·000 EXZ 273 40 0·40 1·459 1·224 0·126 0·112 0·117 0·41 0·21 0·79 – – 
F-115 46·34 6·91 0·23 0·000 EXZ 277 30 0·47 1·568 1·266 0·152 0·133 0·123 0·43 0·29 0·70 – – – – 
F-116 46·31 6·90 0·38 0·000 EXZ 273 30 0·33 1·454 1·264 0·143 0·140 0·018 0·38 0·32 0·68 – – – – 
F-113* 46·26 5·90 0·63 0·000 EXZ 244 30 0·33 1·389 1·182 0·105 0·084 0·198 0·46 0·26 0·74 – – – 
F-127* 46·36 5·90 0·55 0·000 EXZ 255 30 0·27 1·227 1·101 0·067 0·044 0·345 0·50 0·09 0·90 – – – – 
F-126 46·09 5·90 0·37 0·000 EXZ 225 30 0·33 1·360 1·144 0·101 0·093 0·080 0·44 0·14 0·86 – – – – 
F-125 45·93 5·86 0·41 0·000 EXZ 207 33 0·33 1·276 1·120 0·083 0·067 0·203 0·48 0·23 0·77 – – 
F-124 45·76 6·06 0·25 0·000 EXZ 221 30 0·47 1·542 1·369 0·199 0·140 0·299 0·44 0·75 0·25 – 
F-122 45·35 5·76 0·24 0·000 EXZ 181 30 0·40 1·496 1·227 0·138 0·147 -0·063 0·39 0·37 0·63 – – – – 
F-123 45·43 5·80 0·21 0·000 EXZ 189 30 0·33 1·442 1·269 0·138 0·124 0·103 0·43 0·34 0·65 – – – – 
F-121 45·37 5·66 0·20 0·000 EXZ 185 31 0·27 1·522 1·291 0·144 0·146 –0·015 0·40 0·30 0·70 – – – – 
F-120 45·16 5·43 0·24 0·000 EXZ 167 30 0·40 1·485 1·320 0·142 0·138 0·033 0·40 0·38 0·62 – 
F-119* 44·83 5·39 0·17 0·000 EXZ 132 30 0·47 1·502 1·286 0·161 0·107 0·343 0·46 0·49 0·50 – – – – 
F-117* 44·69 5·50 0·36 0·001 EXZ 114 42 0·53 1·607 1·261 0·163 0·152 0·066 0·40 0·70 0·29 – 
F-118 44·77 5·40 0·40 0·001 EXZ 125 40 0·33 1·413 1·234 0·122 0·102 0·172 0·44 0·41 0·59 – – – – 
F 53* 44·68 5·09 0·28 0·002 EXZ 120 24 0·46 1·434 1·202 0·165 0·128 0·223 0·38 0·79 0·21 – – 
F-76 44·14 5·87 0·22 0·008 EXZ 56 22 0·60 1·801 1·417 0·233 0·227 0·027 0·38 0·85 0·15 – – – – 
F-77 44·12 5·86 0·47 0·005 EXZ 54 33 0·60 1·650 1·328 0·193 0·176 0·089 0·40 0·86 0·14 – – – 
F-7 44·05 5·33 0·39 0·010 EXZ 49 44 0·67 1·772 1·364 0·216 0·205 0·055 0·39 0·88 0·11 – – – – 
F-82 44·62 9·42 0·16 0·010 EXZ 177 30 0·67 1·670 1·351 0·198 0·200 –0·012 0·38 0·93 0·07 – 
F-5* 43·96 5·33 0·38 0·021 EXZ 40 41 0·53 1·627 1·250 0·167 0·150 0·105 0·40 0·77 0·23 – – – – 
F-131 44·056 5·92 0·40 0·053 EXZ 46 30 0·67 1·769 1·425 0·236 0·211 0·107 0·38 0·85 0·15 – – – – 
F 48 44·11 5·33 0·41 0·078 EXZ 54 32 0·46 1·419 1·187 0·140 0·128 0·083 – – – – – 
F-132 43·53 5·63 0·25 0·813 REF – 30 0·67 1·687 1·348 0·218 0·233 –0·069 – – – – – – – 
F-12 43·88 7·45 0·45 0·751 ERF – 20 0·53 1·665 1·293 0·171 0·180 –0·057 0·38 0·94 0·06 – – 
F-50 43·79 6·43 0·36 0·190 REF – 28 0·60 1·669 1·288 0·181 0·174 0·042 0·39 0·81 0·19 – – 
F-49 43·79 6·40 0·36 0·141 REF – 36 0·53 1·693 1·344 0·199 0·183 0·081 0·40 0·81 0·19 – – – – 
F-9 43·76 6·32 0·46 0·166 REF – 35 0·60 1·713 1·339 0·199 0·215 –0·083 0·40 0·76 0·06 – – – – 
F-46 43·66 5·62 0·14 0·152 REF – 29 0·60 1·581 1·283 0·170 0·170 –0·003 0·40 0·76 0·24 – – – – 
F-45 43·45 5·55 0·36 0·297 REF – 25 0·47 1·643 1·321 0·177 0·160 0·096 0·41 0·66 0·34 – 
F-75 43·36 5·79 0·49 0·407 REF – 36 0·60 1·585 1·271 0·162 0·137 0·154 0·43 0·85 0·15 – – – – 
POP LAT LON BLC LGM GRO DIS n pPL Ar Ae He Ho f SR K1 K2 
F-80 49·57 7·59 0·31 0·000 OUT 634 39 0·00 1·019 1·004 0·003 0·003 –0·013 0·61 0·03 0·97 – – – 
F-128 48·62 9·785 0·27 0·000 OUT 558 15 0·07 1·043 1·010 0·009 0·009 –0·037 0·57 0·02 0·98 – – – 
F-129 48·62 9·85 0·29 0·000 OUT 558 15 0·00 1·000 1·000 0·000 0·000 0·000 0·59 0·02 0·98 – – – – 
F-130 48·63 9·85 0·35 0·000 OUT 560 14 0·07 1·066 1·040 0·026 0·014 0·458 0·56 0·12 0·88 – – – – 
F-51 47·42 8·00 0·55 0·000 OUT 396 17 0·27 1·278 1·180 0·093 0·082 0·117 0·46 0·12 0·87 – – 
F-84 47·38 7·76 0·28 0·000 OUT 390 17 0·20 1·197 1·141 0·083 0·071 0·158 0·45 0·13 0·87 – – – – 
F-34* 47·29 7·71 0·15 0·000 OUT 381 50 0·27 1·212 1·114 0·073 0·055 0·249 0·48 0·19 0·81 – – 
F-134 47·43 4·93 0·24 0·000 OUT 421 14 0·07 1·065 1·034 0·023 0·019 0·188 0·57 0·04 0·95 – – – – 
F-135 47·24 4·73 0·36 0·000 OUT 405 20 0·13 1·153 1·135 0·069 0·067 0·033 0·43 0·23 0·77 – – – – 
F-114 46·31 6·97 0·15 0·000 EXZ 273 40 0·40 1·459 1·224 0·126 0·112 0·117 0·41 0·21 0·79 – – 
F-115 46·34 6·91 0·23 0·000 EXZ 277 30 0·47 1·568 1·266 0·152 0·133 0·123 0·43 0·29 0·70 – – – – 
F-116 46·31 6·90 0·38 0·000 EXZ 273 30 0·33 1·454 1·264 0·143 0·140 0·018 0·38 0·32 0·68 – – – – 
F-113* 46·26 5·90 0·63 0·000 EXZ 244 30 0·33 1·389 1·182 0·105 0·084 0·198 0·46 0·26 0·74 – – – 
F-127* 46·36 5·90 0·55 0·000 EXZ 255 30 0·27 1·227 1·101 0·067 0·044 0·345 0·50 0·09 0·90 – – – – 
F-126 46·09 5·90 0·37 0·000 EXZ 225 30 0·33 1·360 1·144 0·101 0·093 0·080 0·44 0·14 0·86 – – – – 
F-125 45·93 5·86 0·41 0·000 EXZ 207 33 0·33 1·276 1·120 0·083 0·067 0·203 0·48 0·23 0·77 – – 
F-124 45·76 6·06 0·25 0·000 EXZ 221 30 0·47 1·542 1·369 0·199 0·140 0·299 0·44 0·75 0·25 – 
F-122 45·35 5·76 0·24 0·000 EXZ 181 30 0·40 1·496 1·227 0·138 0·147 -0·063 0·39 0·37 0·63 – – – – 
F-123 45·43 5·80 0·21 0·000 EXZ 189 30 0·33 1·442 1·269 0·138 0·124 0·103 0·43 0·34 0·65 – – – – 
F-121 45·37 5·66 0·20 0·000 EXZ 185 31 0·27 1·522 1·291 0·144 0·146 –0·015 0·40 0·30 0·70 – – – – 
F-120 45·16 5·43 0·24 0·000 EXZ 167 30 0·40 1·485 1·320 0·142 0·138 0·033 0·40 0·38 0·62 – 
F-119* 44·83 5·39 0·17 0·000 EXZ 132 30 0·47 1·502 1·286 0·161 0·107 0·343 0·46 0·49 0·50 – – – – 
F-117* 44·69 5·50 0·36 0·001 EXZ 114 42 0·53 1·607 1·261 0·163 0·152 0·066 0·40 0·70 0·29 – 
F-118 44·77 5·40 0·40 0·001 EXZ 125 40 0·33 1·413 1·234 0·122 0·102 0·172 0·44 0·41 0·59 – – – – 
F 53* 44·68 5·09 0·28 0·002 EXZ 120 24 0·46 1·434 1·202 0·165 0·128 0·223 0·38 0·79 0·21 – – 
F-76 44·14 5·87 0·22 0·008 EXZ 56 22 0·60 1·801 1·417 0·233 0·227 0·027 0·38 0·85 0·15 – – – – 
F-77 44·12 5·86 0·47 0·005 EXZ 54 33 0·60 1·650 1·328 0·193 0·176 0·089 0·40 0·86 0·14 – – – 
F-7 44·05 5·33 0·39 0·010 EXZ 49 44 0·67 1·772 1·364 0·216 0·205 0·055 0·39 0·88 0·11 – – – – 
F-82 44·62 9·42 0·16 0·010 EXZ 177 30 0·67 1·670 1·351 0·198 0·200 –0·012 0·38 0·93 0·07 – 
F-5* 43·96 5·33 0·38 0·021 EXZ 40 41 0·53 1·627 1·250 0·167 0·150 0·105 0·40 0·77 0·23 – – – – 
F-131 44·056 5·92 0·40 0·053 EXZ 46 30 0·67 1·769 1·425 0·236 0·211 0·107 0·38 0·85 0·15 – – – – 
F 48 44·11 5·33 0·41 0·078 EXZ 54 32 0·46 1·419 1·187 0·140 0·128 0·083 – – – – – 
F-132 43·53 5·63 0·25 0·813 REF – 30 0·67 1·687 1·348 0·218 0·233 –0·069 – – – – – – – 
F-12 43·88 7·45 0·45 0·751 ERF – 20 0·53 1·665 1·293 0·171 0·180 –0·057 0·38 0·94 0·06 – – 
F-50 43·79 6·43 0·36 0·190 REF – 28 0·60 1·669 1·288 0·181 0·174 0·042 0·39 0·81 0·19 – – 
F-49 43·79 6·40 0·36 0·141 REF – 36 0·53 1·693 1·344 0·199 0·183 0·081 0·40 0·81 0·19 – – – – 
F-9 43·76 6·32 0·46 0·166 REF – 35 0·60 1·713 1·339 0·199 0·215 –0·083 0·40 0·76 0·06 – – – – 
F-46 43·66 5·62 0·14 0·152 REF – 29 0·60 1·581 1·283 0·170 0·170 –0·003 0·40 0·76 0·24 – – – – 
F-45 43·45 5·55 0·36 0·297 REF – 25 0·47 1·643 1·321 0·177 0·160 0·096 0·41 0·66 0·34 – 
F-75 43·36 5·79 0·49 0·407 REF – 36 0·60 1·585 1·271 0·162 0·137 0·154 0·43 0·85 0·15 – – – – 

Information is provided in the table from left to right. The first three columns provide the general information; numbers of populations (POP) and their locations expressed as latitudes (LAT) and longitudes (LON). Populations are sorted according to their latitudinal co-ordinates from North to South. The next two columns provide the Maxent logistic output values in the studied North–South gradient for the baseline climate logistic output (BLC) and the last glacial maximum logistic output (LGM). The next two columns provide information on the phylogeographic hypothesis implemented: GRO summarizes the assignment of the populations as outlier populations (OUT) at the northern range limit; as part of the proposed post-LGM extension zone (EXZ), as part of the proposed LGM refugium (REF), and the distance of the populations in the extension zone to populations in the refugium (DIS). The next seven columns provide information on the genetic diversity observed with allozymes: n = number of individuals studied, pPL = polymorphic loci using 0·95 criterion, Ar = allelic richness, Ae = effective number of alleles per locus, He = expected heterozygosity, Ho = observed heterozygosity, and f = inbreeding coefficient; mSR = mean selfing rate, Bayesian assignment using INSTRUCT to K1 and K2. The last four columns report the genetic diversity observed in the chloroplast genome. Haplotype numbers (1, 2, 4 and 7) are given according to Hunt et al. (2009); we report here the number of individuals with a haplotype 1, 2, 4 and 7.

*Populations departing from Hardy–Weinberg equilibrium.

The environmental data describing baseline climate (BioClim layers for the period 1950–2000 at a spatial resolution of 30 arc s) and the LGM climate (BioClim layers derived from the Paleoclimate Modelling Intercomparison Project Phase II database at a spatial resolution of 2·5 arc min), were retrieved from the WorldClim database (Hijmans et al., 2005). Climate change downscaled data from the IPCC were downloaded from the International Centre for Tropical Agriculture (CIAT) data portal (http://www.ccafs-climate.org/). To explore variability of the LGM climate scenarios, we used two generalized climate models: the CCSM-GCM (http://www.ccsm.ucar.edu) and MIROC-GCM (http://www-pcmdi.llnl.gov/ipcc/model_documentation/MIROC3.2_hires.htm). To explore variability of the future climate change scenarios, we used the outputs of four general circulation models (CCSR/NIES, CGCM2, CSIRO-Mk2 and HadCM3) projected to 2080–2090 under the A2 IPCC emission scenario. Compared with the rest of the available scenarios, A2 IPCC predicts the second (after A1FI scenario) most dramatic globally averaged surface warming of 3·4 °C and a sea level rise of 0·23–0·51 m in 2090–2099 relative to 1980–1999 (Solomon et al., 2007). In contrast, the B1 IPCC, a low emission scenario, predicts temperature change of 1·8 °C and a sea level rise of 0·18–0·38 m, and would be expected to have the smallest (compared with other scenarios) influence with regard to habitat loss and range retreat of the target species. Thus, A2 IPCC was chosen with the aim of showing the projections of species distribution models under the assumption of a dramatic change in climate.

Our choice of environmental predictors was determined by the climate data available for the past and future climates. The final set included those predictors that characterize seasonality of climate (bio4, temperature seasonality; bio7, temperature annual range; bio15, precipitation seasonality), limits to climatic tolerance of the species (bio6, minimum temperature of the coldest month; bio18, precipitation of the warmest quarter; bio8, mean temperature of the wettest quarter) and water availability (bio19, precipitation of the coldest quarter). As an additional layer, we used the Generalised Geology map of Europe that describes the generalized geological age of the European surface outcrops of bedrock downloaded from the USGS website (http://energy.usgs.gov/). The geology of the study area was assumed to be important in shaping species' response to the environment, because the known distributions of A. fontanum are closely associated with limestone outcrops. For all selected variables, Spearman's ρ statistic was <0·85, implying low covariation between these environmental factors.

Maxent software v. 3·3.3 (Phillips et al., 2006) has been successfully used to predict potential distribution outside the known environmental and geographic range (Marske et al., 2011). For each climate change scenario, we made a projection of a logistic output within the geographic extent of the study area generated by Maxent, using 10-fold cross-validation with random seed and ‘regularization multiplier’ set to 0·01. Although setting the regularization multiplier very small was expected to result in a highly overfitted output that is characterized by a very close fit to the given presence data, it allowed us to reduce the areas of low reliability of predictions. To identify such areas, we used the multivariate environmental similarity surfaces (MESS maps) available in Maxent. The MESS maps measure the similarity of any given point to a reference set of points, with respect to the chosen predictor variables. We also analysed the effects of ‘clumping’, the process by which features are constrained to remain within the range of values in the training data (Elith et al., 2011). The resulting LGM projections were compared with the current understanding of the distribution of ice sheets in Europe (Buoncristiani and Campy, 2011). To assess the model's performance, a threshold-independent measure, the area under the curve (AUC) score, was used.

To transform continuous logistic output values from the models into binary presence–absence values, we used threshold values based on maximum training sensitivity plus specificity (MTSS). The MTSS values were averaged over ten runs of each model. Although the choice of the threshold method was shown to alter estimates of range changes under global warming (Nenzen and Araujo, 2011), there is little agreement on what is the most suitable threshold. Therefore, we used MTSS, which found a high degree of support when evaluated against other methods and is regarded as ‘stringent’ (Richmond et al., 2010).

To explore the variability of predictions by different modelling techniques, we used Biomod ensemble software (Thuiller et al., 2009). Following Schorr et al. (2012), we chose four models: GLM (generalized linear models), GBM (generalized boosting models), GAM (generalized additive models) and RF (Breiman and Cutler's Random Forest for classification and regression). The models were trained and tested on the A. fontanum presence data described in the first paragraph of this section. Absence data were generated in the Biomod environment by using 1000 points randomly selected from the pool of localities randomly generated in the GIS environment of Europe. The procedure was repeated twice, resulting in two sets of pseudo-absences. Model performance was accessed using Roc, TSS and Kappa statistics. The best model was projected onto the area of interest stretching from the coastline of the Maritime Alps to the North, over 6·2° of latitude, using current climate and two LGM climate scenarios, CCSM and MIROC.

Genetic diversity

Leaf material of A. fontanum subsp. fontanum was collected to generate allozyme and chloroplast DNA sequence markers. We compiled a data set from a total of 40 locations representing a South–North climatic transect running from the Maritime and Western limestone Alps, along the Jura mountains, towards outlying occurrences in southern Germany, and the Dijon area in France. Twenty-two locations were new collections for this study, in addition to 18 sites used previously (Hunt et al., 2009). Locations were selected to achieve a uniform density across the transect. Up to 30 individual plants were sampled per location (Table 1). Total population sampling comprised per country: France (n = 29), Switzerland (n = 6), Germany (n = 4) and Italy (n = 1) (Fig. 1).

Fig. 1.

Distribution of A. fontanum modelled by Maxent under the assumption of the four climate scenarios: (A) LGM climate, based on the output of CCSM-GCM; (B) LGM climate, based on the output of MIROC-GCM. In (A) and (B), the extent of the LGM ice sheets is shown as solid white with dots. (C) Baseline climate (=current); (D) future climate based on the outputs of CGCM2, HadCM3, CSIRO-Mk2 and CCSR/NIES; the decrease of the land surface area due to the rising sea levels is not shown. The logistic outputs of four Maxent models were transformed to binary presence–absence values using the MTSS threshold. In the LGM and baseline (=current) climate outputs, light grey indicate the area interpreted as presence; in the future climate, darker shades of grey indicate a higher number of the overlapping model's outputs interpreted as presence (see Supplementary Data Fig. S4 for individual maps). Thus, the darkest grey corresponds to the highest number of the overlapping models (four in our case). Black dots correspond to all available collection localities, and white dots circled black, to the populations with associated genetic information.

Fig. 1.

Distribution of A. fontanum modelled by Maxent under the assumption of the four climate scenarios: (A) LGM climate, based on the output of CCSM-GCM; (B) LGM climate, based on the output of MIROC-GCM. In (A) and (B), the extent of the LGM ice sheets is shown as solid white with dots. (C) Baseline climate (=current); (D) future climate based on the outputs of CGCM2, HadCM3, CSIRO-Mk2 and CCSR/NIES; the decrease of the land surface area due to the rising sea levels is not shown. The logistic outputs of four Maxent models were transformed to binary presence–absence values using the MTSS threshold. In the LGM and baseline (=current) climate outputs, light grey indicate the area interpreted as presence; in the future climate, darker shades of grey indicate a higher number of the overlapping model's outputs interpreted as presence (see Supplementary Data Fig. S4 for individual maps). Thus, the darkest grey corresponds to the highest number of the overlapping models (four in our case). Black dots correspond to all available collection localities, and white dots circled black, to the populations with associated genetic information.

Allozyme experimental conditions followed Hunt et al. (2009), with the exception that allelic variations at the aspartate aminotransferase locus 2 (AAT-2; EC 2·6.1·1) and diaphorase locus 1 (DIA-1; EC1·8.1·4) were included. Genotype diversity was recorded for 15 loci. The statistical independence and selective neutrality of the 15 allozyme loci were examined by exact tests between all pairs of polymorphic loci using GDA v. 1·16 (Lewis and Zaykin, 2001); and by the Ewens–Watterson test implemented in POPGENE v. 1·32 (Yeh and Boyle, 1997). In both cases, a liberal sequential Bonferroni correction (i.e. P = 0·005) was applied for multiple tests (Rice, 1989). Population structure and inbreeding were assessed by the analysis of variance procedure, which provides estimators (f, F, θ) of F-statistics, FIS, FIT and FST. Using POPGENE, the following diversity measures were estimated for each population over all loci: the proportion of polymorphic loci (pPL) (0·95 criterion), the effective number of alleles per locus (Ae), the expected heterozygosity (He) and the observed heterozygosity (Ho). The inbreeding coefficient (f) and allelic richness (Ar) was estimated to take the relative differences in sampling sizes into consideration, using FSTAT v. 2·9.3 (Goudet, 2001).

Chloroplast DNA sequence variation was analysed in the rps4-trnS intergenic spacer region, using previously reported primers and protocols (Hunt et al., 2009). DNA sequence variation was examined for 4–10 individuals from 21 out of the 40 locations. Haplotype data were used to infer effective population size through time. Such analyses were carried out using Bayesian skyline plots (Drummond et al., 2005) generated with the BEAST 1·7.5 software package (http://beast.bio.ed.ac.uk) under the assumption of a constant clock with a rather conservative estimated rate of 8 e−10.

Comparing patterns of genetic diversity with SDM outputs

Geographic patterns in diversity were assessed by Mann–Whitney U-tests and by linear regressions of genetic diversity against geographic distance (km) using INSTAT v. 3·036 (University of Reading, UK). These patterns were used to test the two (interlinked) hypotheses (see the Introduction). (RQ1) Refugial populations from localities with relatively stable ecological niche space under changing climate should retain higher genetic diversity due to long-term persistence. (RQ2) Range expansion/migration into newly available niche space should create a gradient of genetic diversity. RQ1 was examined by a one-sided Mann–Whitney U-test, whereby relative intrapopulation genetic diversity was compared between areas predicted to be of low or high ecological suitability under the modelled glacial conditions. RQ2 was examined by linear regression, whereby intrapopulation genetic diversity was compared with the shortest linear kilometre distance to the nearest sampling location, within the area of highest predicted ecological suitability under the modelled glacial conditions. This approach was chosen to take the broader modelled ‘area of survival’ into consideration and avoid a bias regarding the potential distance of colonization from a single reference location. RQ2 was further tested by examining the relationship between genetic similarity (logFST) and geographic distance (logkm) using all sampling locations; and separately for the samples derived from the area of low ecological suitability under the modelled glacial conditions. Mantel tests were performed using the R-package v. 4·00 (R Development Core Team, 2012), and the geographic distances were derived from the longitude and latitude co-ordinates (Table 1). Pairwise FST values were calculated using GDA. A partial Mantel test on genetic distance and categorized geographic distance classes (Table 1) was performed, determining over which geographic scales the genetic isolation by distance relationship remained statistically significant. Based on the Maxent logistic output values under the LGM scenario and genetic diversity statistics, populations were considered to be part of one (out of a maximum three) region(s) (see Table 1): (1) the LGM refugium; (2) the expansion zone (i.e. post-LGM colonization); and (3) the outlying populations of the northern range limit (Dijon, France; Schwaebische Alp, Germany; and north-eastern Switzerland). We further considered that the obtained Bayesian assignment could be distorted by allele frequencies changing under drift (Excoffier et al., 2009). This issue was investigated by inferring genetic isolation by genetic drift (IBD) using INSTRUCT v. 1·0 (Gao et al., 2007), which does not make explicit assumptions about Hardy–Weinberg equilibrium (HWE). Analyses were initially performed under mode 2 allowing for admixture (K = 2 clusters), then under mode 3 allowing for admixture (K = 3 clusters) and to estimate individual selfing rates. The estimated selfing rates are considered as relative values for allowing among-population comparisons. Both runs employed 600 × 104 Markov chain Monte Carlo (MCMC) replicates, including 100 × 104 for the burn-in period, and five chains. These assignment tests exploit the effect of allele frequency changes under drift/migration to result in assignment proportions, independently of historical population admixture.

In order to reconcile SDM and SPG, we inferred the relationships between SDM predictions and observed genetic diversity. To compare LGM climate models with patterns in population genetic data, we extracted the values of logistic output generated by SDM, at locations selected for collection of material for genetic analysis (40 in total). Spearman's ρ statistic was used to estimate the degree of association of expected heterozygosity (He) with the following environmental variables and SDM outputs: longitude (LON); latitude (LAT); present day annual mean temperature (ANMT); CCSM–ANMT anomaly calculated as the difference between the current climate ANMT and that at the LGM as estimated by CCSM (CCSM-ANO); ANMT anomaly by MIROC calculated as the difference between the current climate ANMT and that at the LGM as estimated by MIROC (MIROC-ANO); difference between output value produced by Random Forest (RF) for current climate and LGM climate based on MIROC (MIROC-RF); the same based on CCSM (CCSM-RF); the difference between logistic output produced by Maxent for current climate and LGM climate based on MIROC (MIROC-M); and the same based on CCSM (CCSM-M). ANMT is highly correlated with the minimum temperature of the coldest month used as an input variable in the models.

To estimate the impact of the spatio-temporal population dynamics on the maintenance of genetic diversity in the projected future habitat range, we inferred the total allelic diversity (ALD), the private allele fraction (PAF) and the common allele fraction (FAC) in populations assigned to the outlying populations, expansion zone and LGM refugium (Tables 1, 2).

RESULTS

Palaeo-distribution modelling

A suitable habitat for A. fontanum was predicted well beyond the available collection localities, in spite of the Maxent model based on the baseline climate (the time period 1950–2000) being deliberately overfitted, due to a low setting of the ‘regularization multiplier’ feature (Fig. 1). Projections based on CCSM-GCM predicted very different habitat suitability for A. fontanum across the area of interest, compared with those based on MIROC-GCM, with the major differences observed in France and Central Europe. According to Maxent predictions based on CCSM-GCM outputs, most of the current habitats potentially available in Central Europe today were climatically unsuitable for A. fontanum at the time of the LGM, with the refugium located on the Spanish Peninsula, in the Maritime Alps and along an expanded Mediterranean coastline. According to MIROC-GCM, reduction of the suitable habitat was observed in the northern part of the current range, while suitability values at the locations with the current species presence records remained relatively similar to those corresponding to the baseline climate (Supplementary Data Fig. S1). Despite apparent differences, both MIROC and CCSM were consistent in predicting the LGM distribution within the current range of A. fontanum. They suggest LGM contraction of the available habitat, particularly in its northern part, and consequently post-glacial expansion. The CCSM-GCM data suggested that the LGM range was much smaller compared with that of today and restricted to the southernmost part of the current range. The LGM predictions were free from clumping, but MESS analysis identified a number of variables outside the range present in the training data (Supplementary Data Fig. S2). The problematic areas identified by MESS analysis did not overlap with those locations of the populations with associated genetic information (white dots in Fig. S2). The distribution of ice sheets was consistent with CCSM predictions, but not with those by MIROC.

When the outputs of four different SDM algorithms were compared in the BIOMOD environment, the highest rated model by all statistics (Roc, TSS and Kappa) was that produced by RF. The LGM distributions produced by RF were largely similar to those generated by Maxent (Supplementary Data Fig. S3). As in the case of Maxent, two different LGM scenarios used by RF resulted in markedly different habitat suitability predictions; with the LGM distributions using MIROC climate being similar to distributions in the current climate, and with those using CCSM climate showing much lower habitat suitability values compared with the present day.

Genetic diversity

Polymorphisms were found in 13 out of 15 loci across the 13 enzyme systems (only ACP-1 and TPI-1 were monomorphic). In total, 52 alleles were discriminated. Four out of nine known European chloroplast haplotypes were found to occur in the regions studied here. The allozyme loci were inferred as statistically independent and selectively neutral for subsequent population genetic analysis based on the results of the test for linkage disequilibria (Supplementary Data Table S1), selective neutrality (Supplementary Data Table S2) and the Evans–Watterson test on the combined population level. An FST of 0·169 [95 % confidence interval (CI) 0·086–0·229)] was estimated over the entire data set, indicating modest differentiation and that the vast proportion of genetic variation was distributed within populations. A small but significant departure from HWE over the entire data set was evident by an overall FIS of 0·088 (95 % CI 0·044–0·128), and individual values of the inbreeding coefficient (f) varied from –0·069 to 0·458, although typically from –0·12 to 0·066 (Table 1). Populations departing from HWE were mainly found in the western or northern-most ranges (Table 1). The mean estimates of genetic variation, as measured by the proportion of polymorphic loci (pPL), allelic richness (Ar) and expected heterozygosity (He), were 0·400, 1·454 and 0·135, respectively. Similar estimates were obtained for populations from the Maritime and Western limestone Alps (Table 1). One-sided Mann–Whitney tests determined that populations of the LGM refugium had significantly higher values of pPL (P < 0·0001), Ar (P < 0·0001) and He (P < 0·0001), and lower value of f (P = 0·0226). Chloroplast sequence variation indicated there has been a relatively stable effective population size over a 750 000 year period, with a steady, but small increase from 1·08 × 10−7 (1·5 × 10−5 to 1·02 × 10−9) to 1·15 × 10−7 (1·5 × 10−5 to 1·04 × 10−8) towards the present.

Post-glacial range expansion

Samples were organized into two groups according to the LGM projections of the CCSM-GCM (Fig. 1, Table 1). Populations with LGM Maxent probabilistic output values of >0·14 (0·14–0·81, mean 0·36) for ecological suitability under LGM climate conditions were grouped as LGM refugium (n = 8), while scores <0·09 (0·00–0·09, mean <0·01) were grouped as ‘expansion zone’ (n = 31). Linear regression analysis determined that there was a significant positive relationship between genetic variation of individual populations within the expansion group and the shortest distance (km) from the nearest population sampled from within the refugial group (Fig. 2). The r2 values for measures pPL, Ar, He and Ho were, respectively, 0·843 (P < 0·0001), 0·798 (P < 0·0001), 0·795 (P < 0·0001) and 0·7645 (P < 0·0001). There was no relationship between the expansion distance and within-population inbreeding levels (f in Table 1) where r2 was 0·001 (P = 0·867). However, the overall estimate of inbreeding was higher in the expansion zone plus outlier populations FIS = 0·111 (95 % CI 0·065–0·158) compared with the LGM refugium FIS = 0·020, which did not differ significantly from zero (95 % CI –0·031 to 0·080). Similarly, population differentiation was higher in the expansion zone, FST = 0·174 (95 % CI 0·092–0·224), compared with the LGM refugium, FST = 0·065 (95 % CI 0·027–0·104). Separation of populations in the expansion zone and the outlier populations resulted in further differentiation between these regions. Populations in the expansion zone only showed an FIS = 0·104 (95 % CI 0·062–0·062) and FST = 0·126 (95 % CI 0·074–0·171), whereas outlier populations showed an FIS = 0·178 (95 % CI 0·026–0·342) and FST = 0·412 (95 % CI 0·088–0·544). A Mantel test on matrices of genetic distance (logFST) and physical distance (logkm) among all pairwise combinations of samples showed an r2 value of 0·336 (P < 0·001), indicating a significant pattern of genetic isolation by distance. A lower r2 value of 0·277 (P < 0·001) was obtained for the reduced data set comprising samples from the expansion area alone.

Fig. 2.

Decay of genetic diversity in the expansion zone as a function of the distance from the stable LGM area: estimated as distance of populations within the expansion zone to the nearest population within the inferred stable LGM area; expressed as the proportion of polymorphic loci (pPL), allelic richness (Ar), expected heterozygosity (He) and inbreeding coefficient (f). Only the inbreeding coefficient did not show correlation with distance.

Fig. 2.

Decay of genetic diversity in the expansion zone as a function of the distance from the stable LGM area: estimated as distance of populations within the expansion zone to the nearest population within the inferred stable LGM area; expressed as the proportion of polymorphic loci (pPL), allelic richness (Ar), expected heterozygosity (He) and inbreeding coefficient (f). Only the inbreeding coefficient did not show correlation with distance.

Along the gradient of genetic diversity, the degree of association between expected heterozygosity (He in Table 1), a set of environmental variables and SDM outputs was the highest regarding latitude (0·85), followed by annual mean temperature anomaly at the LGM estimated by CCSM (CCSM-ANO; 0·76). All remaining scores were below 0·70 (Fig. 3). Within the model outputs, the difference between the logistic output of the Maxent model based on the current climate and that based on the LGM climate estimated by CCSM (CCSM-M) had the highest value of association (0·59) with He across the distribution range.

Fig. 3.

Smoothed regression (below diagonal) and rank-based measure of association estimated by Spearman's ρ statistic (above diagonal), which are more robust and have been recommended if the data do not necessarily come from a bivariate normal distribution, for longitude (LON), for latitude (LAT), for the difference between the output value produced by Random Forest for current climate and LGM climate based on MIROC (MIROC-RF), for the same but based on CCSM (CCSM-RF), for the difference between logistic output produced by Maxent for current climate and LGM climate based on MIROC (MIROC-M), for the same but based on CCSM (CCSM-M), for annual mean temperature (ANMT), for ANMT anomaly by CCSM (CCSM-ANO), for ANMT anomaly by MIROC (MIROC-ANO) and for expected heterozygosity (He). ANMT is highly correlated with the minimum temperature of the coldest month used as an input variable in the models (0·95). The relationships between He and LAT, CCSM-M and CCSM-ANO discussed in the text are highlighted. For example, the value of Spearman's ρ corresponding to the relationship between He and LAT (0·85) can be found at the crossing point between the respective column (He) and row (LAT); the corresponding smoothed regression is diagonally opposite this value.

Fig. 3.

Smoothed regression (below diagonal) and rank-based measure of association estimated by Spearman's ρ statistic (above diagonal), which are more robust and have been recommended if the data do not necessarily come from a bivariate normal distribution, for longitude (LON), for latitude (LAT), for the difference between the output value produced by Random Forest for current climate and LGM climate based on MIROC (MIROC-RF), for the same but based on CCSM (CCSM-RF), for the difference between logistic output produced by Maxent for current climate and LGM climate based on MIROC (MIROC-M), for the same but based on CCSM (CCSM-M), for annual mean temperature (ANMT), for ANMT anomaly by CCSM (CCSM-ANO), for ANMT anomaly by MIROC (MIROC-ANO) and for expected heterozygosity (He). ANMT is highly correlated with the minimum temperature of the coldest month used as an input variable in the models (0·95). The relationships between He and LAT, CCSM-M and CCSM-ANO discussed in the text are highlighted. For example, the value of Spearman's ρ corresponding to the relationship between He and LAT (0·85) can be found at the crossing point between the respective column (He) and row (LAT); the corresponding smoothed regression is diagonally opposite this value.

Spatial changes in the genetic diversity were associated with changes in the rate of selfing as supported by the Mann–Whitney test and the IBD. Populations in the outlying areas had significantly higher selfing rates (mean 0·525, n = 9, P < 0·001) compared with populations in the refugium (mean 0·397, n = 7, P = 0·001) or expansion areas (mean 0·417, n = 22, P = 0·001). These changes are associated with the results of the Bayesian assignment test with K = 2. Populations in the refugium were assigned to one cluster (K1, mean value 0·825, n = 7), whereas populations of the outlying area were assigned to the second cluster (K2, mean value 0·898, n = 9). Populations in the expansion zone had mixed assignment proportion with a progressive increase of the K2 proportion along the South–North gradient. This gradient was supported in the regression plot by r2 = 0·618 (P = 0·001).

Future climate change

The outputs of the Maxent models based on four scenarios for future climate change (CGCM2, HadCM3, CSIRO-Mk2 and CCSR/NIES) varied greatly over the area of study (Fig. 1; Supplementary Data Fig. S4). When logistic outputs of Maxent models were transformed to binary presence–absence values, three out of four models were in agreement, predicting a northern shift of the potentially suitable area by 2080–2090 (Fig. 1; Supplementary Data Fig. S4). The majority of the models also predicted the area of Maritime and Western limestone Alps as potentially suitable sites for A. fontanum in the warmer climate, while the current habitats in the Iberian Peninsula, France, Germany and Switzerland were indicated as potentially becoming climatically unsuitable. These predictions assume a drastic climate change (A2 IPCC SRES scenario), with a species distribution model closely fit to the available presence records and a ‘stringent’ presence–absence transformation threshold.

DISCUSSION

As predicted in RQ1, we found evidence for higher genetic diversity in populations located in the putative LGM refugium compared with populations occurring in the expansion zone or outlier zone. This is probably caused by the longer persistence of these populations. As expected (RQ2), evidence was found for range expansion into an area predicted to be unoccupied during the LGM. In the ‘expansion zone’, the genetic diversity was arranged along a gradient, with the highest diversity in the population closest to the refugium. Thus, the predictions of RQ2 were sustained for populations in this zone. However, populations located in the ‘outlier zone’ were not part of this gradient even though they are also located outside of the LGM refugium. These populations are the most distant from the LGM refugium and have higher inbreeding rates. The distinction between ‘expansion zone’ and ‘outlier zone’ gives an insight into RQ3 and RQ4. As predicted in RQ3, both long- and short-distance dispersal contribute to the migration of this fern. Overall, the recovered pattern suggests a lower contribution of long-distance dispersal events to the colonization of previously unoccupied niches, compared with the contribution of short-distance dispersal events. However, populations in the ‘outlier zone’ carry the genetic signature of pioneers founded by long-distance dispersal. Thus, the contribution of long-distance spore dispersal to the range expansion may not be controlled by this dispersal event alone. These findings relate to RQ4 in that populations in the ‘outlier zone’ carry the signature of an increased inbreeding rate. Populations in the other zones show a rather low inbreeding rate, suggesting a preferred reproduction via outbreeding. The high frequency of inbreeding may diminish the probability of successful establishment of pioneer populations as a result of the accumulation of genetic load. This argument is consistent with previous research on the migration process in ferns that showed differences in genetic load in pioneer and non-pioneer ferns (Lloyd, 1974). The answers to the four research questions provided the framework to address the two main questions (see the Introduction), which are discussed in the following sections.

Climate change and the spatial genetic structure of Asplenium fontanum

The recovered LGM distributions suggested two alternative scenarios. In the CCSM-GCM-based scenario, the range was highly reduced compared with the current range, whereas MIROC-GCM indicated a rather limited range reduction, restricted to the most north-eastern parts of the current range. The genetic diversity distribution was found to be consistent with the CCSM-GCM scenario. The one-sided Mann–Whitney tests found significant differences in the genetic diversity between populations in the predicted LGM refugia and in the predicted expansion. The area predicted as the LGM refugium under the CCSM-GCM scenario is widely considered as a glacial refugium for both animals and plants (Hewitt, 2004) and is consistent with the current understanding of the extent of the LGM ice sheets in this area of Europe (Buoncristiani and Campy, 2011).

In support of the MIROC-GCM scenario, we found evidence for a constant effective population size throughout the studied time interval. This result may be seen as a conflict with the assumption in RQ1. However, constant effective population sizes may also be explained by other factors such as limited contributions of individuals located in the expansion area, or extinction of populations in the southern range. Of concern is the overpredicted current range in the Maxent analyses. The discrepancy between prediction and observed data may be the result of limitations of the SDM approach. The overpredicted areas may have a suitable climate, but lack the continuous forests with suitable limestone outcrops as the result of habitat destruction. The sporadic appearance of isolated populations, such as at man-made walls, e.g. in the port of Amsterdam (Bremer, 2003), are consistent with the Maxent predictions indicating the possibility of colonizing isolated available niches via pioneer populations. Observations based on these populations support the existence of suitable habitats at the northern range but also show the limited contribution of pioneer populations to range expansion in this fern.

The allozyme data provided strong evidence for a latitudinal gradient of genetic diversity associated with climatic variables, in particular annual mean temperature. The distribution of intrapopulation genetic diversity was consistent with the assumption of an LGM refugium in the unglaciated parts of the Maritime and Western Alps and a post-LGM range expansion into the previously glaciated parts of the Western Alps towards the east, along the Rhone valley and Jura Mountains to the north. It is noteworthy that the genetic diversity within the expansion zone does not show a drop as expected under the pioneer model of range expansion (Nichols and Hewitt, 1994). Instead, the observations are consistent with the expectations of the phalanx model, which has been reported for several alpine species (Hewitt, 2004). In this model, a high proportion of individuals colonize the new area via short-distance colonization events instead of processes that result in expanding patches of homozygosity (Nichols and Hewitt, 1994). The distribution of the haplotypes and allozyme data support the phalanx range expansion model (RQ3) because the diversity in the expansion zone maintains a large proportion of the genetic diversity (Table 1). A further important aspect of the phalanx-type expansion is its slower colonization of newly available niche space. For such reasons, A. fontanum may show a delayed response to post-glacial climate warming similar to that observed in other calcicolous species (Dullinger et al., 2012).

The delayed response to climate change in the case of A. fontanum may appear counterintuitive given the generally assumed high migration capacity of ferns (Wolf et al., 2001). However, this argument only takes into account spore dispersal and ignores other critical aspects of the reproductive cycle of ferns when establishing new populations in previously unoccupied but suitable habitats (de Groot et al., 2012a, b). As mentioned, A. fontanum reproduce preferably via outcrossing despite a low frequency of inbreeding occurring throughout the study area. The observed low level of inbreeding is expected to occur as a reproductive reassurance in diploid, sexually reproducing ferns that reproduce preferably via outcrossing (Ranker and Geiger, 2008; Wubs et al., 2010). In ferns, the rate of successful outcrossing is constrained by the need for two gametophytes to be growing very close together and being in the reproductive state at the same time. The distance is constrained by the swimming range of fern spermatozoids, whereas the time is restricted by the short life expectation of gametophytes (Ranker and Geiger, 2008). Taking these factors into account, the probability of successful outcrossing depends on the low probability of spore arrival from two different sporophytes. Finally, we need to consider that founder populations of diploid ferns may carry a genetic load (Ranker and Geiger, 2008), which limits the successful establishment of populations founded via long-distance migration. These limitations may explain the recovered pattern of a continuous distribution for the genetic diversity over the studied latitudinal gradient.

Survival of Asplenium fontanum in a warming world

Asplenium fontanum together with other ferns is expected to migrate northwards into the UK and Ireland under the regime of a warmer climate (Robinson, 2009). This suggestion is consistent with the projection of the SDM under 2080–2090 climate predictions. Interestingly, the predicted potential range of A. fontanum not only includes new highly suitable areas such as the limestone regions in south-west England, but also suggests a high extinction risk for populations in the area predicted as the LGM refugia as well as the described post-LGM expansion zone. Only one out of the four employed climatic models did not result in the prediction of a major range loss from the current distribution range. From this, we conclude that A. fontanum may require migration to newly available niches in the north of its current range to avoid increased extinction risk.

Simulations demonstrated correlations of the dynamics of range shifts (i.e. displacement) and the maintenance of genetic diversity, e.g. that there is a high probability of extinction for populations under a fast range shift scenario (Arenas et al., 2012). As demonstrated by the delayed response to the post-LGM warming, the phalanx expansion mode may limit the response of this fern to future changes, especially if they are rapid. Based on observations of the current expansion zone, we predict that the vast majority of the genetic diversity will be maintained, as long as the range expansion is achieved mainly via short-distance dispersal. The current expansion zone contains 81 % of the total allelic diversity, whereas the LGM stable region contains only 77 % (Table 2). Given these observations, we predict a substantial reduction of the total allelic diversity if the southern populations were lost. However, a large proportion of this diversity is comprised of rare genotypes that only occur in one or two populations (Table 2). If these rare genotypes are excluded, our analyses predict the survival of 100 % of the genetic diversity (CAF in Table 2).

Table 2.

Estimates of the allele frequency in the areas that may function as sources of genetic diversity, to the projected range under the global warming models

 N-AL N-AL% P-AL P-AL% N-CAF N-CAF% MPS 
OUT 24 45·3 0·00 24 77·4 17·22 
EXP 43 81·1 12 54·5 100 25·70 
REF 41 77·4 10 45·5 31 100 27·75 
TOT 53 100 22 100 31 100 24·08 
 N-AL N-AL% P-AL P-AL% N-CAF N-CAF% MPS 
OUT 24 45·3 0·00 24 77·4 17·22 
EXP 43 81·1 12 54·5 100 25·70 
REF 41 77·4 10 45·5 31 100 27·75 
TOT 53 100 22 100 31 100 24·08 

The zones recognized are outlier populations at the northern edge of the current range (OUT), post-LGM expansion zone (EXP) and LGM refugium (REF), and all populations (TOT).

Distribution of genetic diversity considering all observed alleles: number of all alleles observed per zone (N-AL), percentage of alleles per zone (N-AL%), number of private alleles per zone (P-AL), percentage of private alleles (P-AL), number of alleles belonging to the common allele fraction (N-CAF), percentage of the common allele fraction (N-CAF%) and mean population score (MPS). Private alleles were defined as occurring in one or two populations, whereas common alleles occurred in at least three populations.

These estimates assume that the range expansion follows the model of a phalanx expansion (Nichols and Hewitt, 1994). A much more pessimistic prediction is likely under the assumption that the future range will be occupied by spores that are mainly contributed by the outlier populations at the northern edge of the distribution range in France (Dijon), Germany (Schwaebische Alp) and Switzerland (northern Swiss Jura and eastern Swiss Mitteland). These populations show reduced genetic diversity and may have been founded via a single spore colonization event (Table 2). If these populations are the major source of new northward colonizations, then only 45·3 % of the total genetic diversity (or 77 % of the diversity when excluding rare genotypes) may be present in the new expanded areas in northern France, Germany and south-west England. This more positive prediction assumes that the assemblage of the outlier populations occurs via multiple colonizations from areas holding considerable amounts of genetic diversity. Calculating the common allele proportion (alleles found in >3 populations) for single populations in the outlier regions reduces the average to 62 % of the CAF score of the stable zone.

The reproductive strategy of this fern sustains the conservation of genetic diversity in a world of climate change as long as the dynamics of range migration coincides with the colonization of new areas by more than a single spore. Should single spore colonization and subsequent intragametophytic selfing start to prevail, this would result in a substantial loss of the genetic diversity of this species, and subsequently increase the extinction risk of the European A. fontanum.

Dispersal mechanisms and migration capacity

Our study contributes to the increasing understanding of the complexities of range expansion in plant lineages that are generally assumed to have a high migration capacity, such as homosporous ferns and bryophytes. Based on spatial resolution analyses, some authors have argued that ferns do not maintain historic-spatial signals as a result of spore-based dispersal systems (Linares-Palomino and Kessler, 2009). However, this assumption has been repeatedly refuted (Hunt et al., 2009; Shepherd et al., 2009; Wang et al., 2011, 2012). This discrepancy is likely to be caused by reproductive traits, e.g. ploidy level and breeding system, influencing the capacity to shift ranges in response to changing climates. Available studies support the existence of different reproductive strategies, pioneer vs. non-pioneer, among ferns (Lloyd, 1974; Wubs et al., 2010) that have substantial consequences for the migration rate. Some ferns may be able to occupy newly available niche spaces rapidly via single spore colonization. However, many ferns show a different pattern, indicating a lower capability to colonize newly available niche space (Loyd, 1974). As argued above, these ferns may share several critical traits in that they are diploid, sexually reproducing with preference for outcrossing, and have limited vegetative reproduction or totally lack it. According to available estimates (Ranker and Geiger, 2008), this condition is found in the majority of ferns and cannot be regarded as an exception. In a world of rapid climate change, ferns with these traits may carry a higher extinction risk than ferns with traits suited to support long-distance dispersal via a single spore.

This conclusion requires confirmation from further studies, although substantial support is provided by studies applying very different research approaches and geographic/ecological settings. Differences in the reported results are probably due to differences in the reproductive systems and ploidy levels. Furthermore, our study documented that the migration limitations are of importance not only to plants dispersing with large propagules such as seeds but also for spore-bearing plants such as ferns. Comparable results were previously recorded for the moss Homalothecium sericeum (Desamore et al., 2012). Subsequently, the study contributed to the rapidly growing body of studies detangling the role of climate and migration in shaping species ranges (Normand et al., 2011; Corlett and Westcott, 2013). Research on the fate of the populations in the ‘outlier zone’ may provide further insight into the contribution of mating system shifts on the trailing edge (Levin, 2012).

More broadly, the results of this study underpin the need to focus on species-specific traits and how they are affected by rapid climate change, on the distribution of plant diversity, the maintenance of genetic diversity and/or conservation of species diversity (Cunze et al., 2013). The dynamics of climate change is very important given its putative correlation with the dynamics of range change via expansion, local extinction and migration. The dynamics depends strongly on the contribution of pioneer populations to the colonization of newly available niches. Populations at the leading edge were found to show increased frequency of selfing that subsequently leads to the reduction of genetic diversity (Levin, 2012). The reported results are consistent with other studies that raised awareness of the enhanced extinction risk for plants in a rapidly changing global climate (Corlett and Westcott, 2013; Cunze et al., 2013). We have demonstrated that similar mechanisms influence the persistence/extinction effects on plants dispersing via either wind-dispersed spores or seeds (Zhu et al., 2012).

SUPPLEMENTARY DATA

Supplementary data are available online at www.aob.oxfordjournals.org and consist of the following. Figure S1: disparity maps showing the differences between Maxent predictions assuming baseline climate and LGM climate. Figure S2: distribution map of A. fontanum modelled by Maxent assuming an LGM climate. Figure S3: distribution maps of A. fontanum modelled by RF. Figure S4: distribution of A. fontanum modelled assuming a warmer climate, for the period 2080–2090. Table S1: populations and loci with significant departures from linkage. Table S2: loci with allele frequencies significantly different from neutral expectations under Ewans–Watterson test.

ACKNOWLEDGEMENTS

This work was supported by The Natural History Museum (NHM), London. We thank the team of the DNA sequencing laboratories at the NHM and the curation team of the NHM, in particular Alison Paul, for their support. We acknowledge the support of various colleagues in providing us with locality information and/or sharing materials.

LITERATURE CITED

Alsos
IG
Ehrich
D
Thuiller
W
et al.  
Genetic consequences of climate change for northern plants
Proceedings of the Royal Society B: Biological Sciences
 
2012
279
2042
2051
Arenas
M
Ray
N
Currat
M
Excoffier
L
Consequences of range contractions and range shifts on molecular diversity
Molecular Biology and Evolution
 
2012
29
207
218
Bellard
C
Bertelsmeier
C
Leadley
P
Thuiller
W
Courchamp
F
Impacts of climate change on the future of biodiversity
Ecology Letters
 
2012
15
365
377
Bremer
P
Chandra
S
Srivastava
M
Some aspects of the fern flora (Filicopsida) of the Netherlands
Pteridology in the new millenium
 
2003
Dodrecht
Kluwer
327
340
Buoncristiani
JF
Campy
M
Ehler
J
Gibbard
PL
Hughes
PD
Quaternary glaciations in the French Alps and Jura
Quaternary glaciations: extent and chronology: a closer look
 
2011
Amsterdam
Elsevier Science Ltd
Cobben
MMP
Verboom
J
Opdam
PFM
et al.  
Projected climate change causes loss and redistribution of genetic diversity in a model metapopulation of a medium-good disperser
Ecography
 
2011
34
920
932
Cordellier
M
Pfenninger
M
Inferring the past to predict the future: climate modelling predictions and phylogeography for the freshwater gastropod Radix balthica (Pulmonata, Basommatophora)
Molecular Ecology
 
2009
18
534
544
Corlett
RT
Westcott
DA
Will plant movements keep up with climate change?
Trends in Ecology and Evolution
 
2013
28
482
488
Cunze
S
Heydel
F
Tackenberg
O
Are plant species able to keep pace with the rapidly changing climate?
PLoS One
 
2013
8
e67909
D'Amen
M
Zimmermann
NE
Pearman
PB
Conservation of phylogeographic lineages under climate change
Global Ecology and Biogeography
 
2013
22
93
104
Desamore
A
Laenen
B
Stech
M
et al.  
How do temperate bryophytes face the challenge of a changing environment? Lessons from the past and predictions for the future
Global Change Biology
 
2012
18
2915
2924
Drummond
AJ
Rambaut
A
Shapiro
B
Pybus
OG
Bayesian coalescent inference of past population dynamics from molecular sequences
Molecular Biology and Evolution
 
2005
22
1185
1192
Dullinger
S
Willner
W
Plutzar
C
et al.  
Post-glacial migration lag restricts range filling of plants in the European Alps
Global Ecology and Biogeography
 
2012
21
829
840
Elith
J
Phillips
SJ
Hastie
T
Dudik
M
Chee
YE
Yates
CJ
A statistical explanantion of MaxEnt for ecologists
Diversity and Distributions
 
2011
17
43
57
Excoffier
L
Foll
M
Pettit
RJ
Genetic consequences of range expansions
Annual Review of Ecology, Evolution and Systematics
 
2009
40
481
501
Gao
H
Williamson
S
Bustamante
CD
An MCMC approach for joint inferences of populaiton structure and inbreeeding rates from multi-locus genotype data
Genetics
 
2007
176
1635
1651
Goudet
J
FSTAT, a program to estimate and test gene diversities and fixation indices. V. 2·9.3
 
2001
de Groot
GA
Verduyn
B
Wubs
ERJ
Erkens
RHJ
During
HJ
Inter-and intraspecific variation in fern mating systems after long-distance colonization: the importance of selfing
BMC Plant Biology
 
2012a
12
3
de Groot
GA
During
HJ
Ansell
SW
et al.  
Diverse spore rains and limited local exchange shape fern genetic diversity in a recently created habitat colonized by long-distance dispersal
Annals of Botany
 
2012b
109
965
978
Gugger
PF
Gonzalez-Rodrıguez
A
Rodrıguez-Correa
H
Sugita
S
Cavender-Bares
J
Southward Pleistocene migration of douglas-fir into Mexico: phylogeography, ecological niche modeling, and conservation of ‘rear edge’ populations
New Phytologist
 
2011
189
1185
1199
Herrero
A
Prada
C
Pajaron
S
Pangua
E
Numeros cromosoaticos de plantas occidentales, 723–726
Anales Jardin Botanico de Madrid
 
1997
55
135
Hewitt
GM
Genetic consequences of climatic oscillations in the Quaternary
Philosophical Transactions of the Royal Society B: Biological Sciences
 
2004
359
183
195
Hickerson
MJ
Carstens
BC
Cavender-Bares
J
et al.  
Phylogeography's past, present, and future: 10 years after Avise, 2000
Molecular Phylogenetics and Evolution
 
2010
54
291
301
Hijmans
RJ
Cameron
SE
Parra
JL
Jones
PG
Jarvis
A
Very high resolution interpolated climate surfaces for global land areas
International Journal of Climatology
 
2005
25
1965
1978
Hunt
HV
Ansell
SW
Russell
S
Schneider
H
Vogel
JC
Genetic diversity and phylogeography in two diploid ferns, Asplenium fontanum subsp. fontanum and A. petrarchae subsp. bivalens, in the western Mediterranean
Molecular Ecology
 
2009
18
4940
4954
Ibrahim
KM
Nichols
RA
Hewitt
GM
Spatial patterns of genetic variation generated by different forms of dispersal during range expansion
Heredity
 
1996
77
282
291
Jay
F
Manel
S
Alvarez
N
et al.  
Forecasting changes in population genetic structure of alpine plants in response to global warming
Molecular Ecology
 
2012
21
2354
2368
Knowles
LL
Alvarado-Serrano
DF
Exploring the population genetic consequences of the colonization process with spatio-temporally explicit models: insights from coupled ecological, demographic and genetic models in montane grasshoppers
Molecular Ecology
 
2010
19
3727
3745
Levin
DA
Mating system shifts on the trailing edge
Annals of Botany
 
2012
109
613
620
Lewis
PO
Zaykin
D
Genetic data analysis: computer program for the analysis of allelic data. V. 1·0
 
2001
Linares-Palomino
R
Kessler
M
The role of dispersal ability, climate and spatial separation in shaping biogeographical patterns of phylogenetically distant plant groups in seasonally dry Andean forests of Bolivia
Journal of Biogeography
 
2009
36
280
290
Lloyd
RM
Mating systems and genetic load in pioneer and non-pioneer Hawaiian Pteridophyta
Botanical Journal of the Linnean Society
 
1974
69
23
35
Maiorano
L
Cheddadi
R
Zimmermann
NE
et al.  
Building the niche through time: using 13,000 years of data to predict the effects of climate change on three tree species in Europe
Global Ecology and Biogeography
 
2012
22
302
317
Marske
KA
Leschen
RAB
Buckley
TR
Reconciling phylogeography and ecological niche models for New Zealand beetles: looking beyond glacial refugia
Molecular Phylogenetics and Evolution
 
2011
59
89
102
Nenzen
HK
Araujo
MB
Choice of threshold alters projections of species range shifts under climate change
Ecological Modelling
 
2011
222
3346
3354
Nichols
R
Hewitt
G
The genetic consequences of long distance dispersal during colonization
Heredity
 
1994
72
312
317
Normand
S
Ricklefs
RE
Skov
F
Bladt
J
Tackenberg
O
Svenning
JC
Postglacial migration supplements climate in determining plant species ranges in Europe
Proceedings of the Royal Society B: Biological Sciences
 
2011
278
3644
3653
Phillips
SJ
Anderson
RP
Schapire
RE
Maximum entropy modeling of species geographic distributions
Ecological Modelling
 
2006
190
231
259
R Development Core Team
R Foundation for Statistical Computing
R: A language and environment for statistical computing
 
2012
Vienna
Austria
Ranker
TA
Geiger
JMO
Ranker
TA
Haufler
CH
Population genetics
Biology and evolution of ferns and lycophytes
 
2008
Cambridge
Cambridge University Press
107
133
Reichstein
T
Schneller
JJ
Asplenium pseudofontanum Kossinosky (Aspleniaceae, Pteridophyta). Studies in Asplenium for Flora Iranica
Candollea
 
1982
37
117
128
Rice
WR
Analyzing tables of statistical tests
Evolution
 
1989
43
223
225
Richards
CL
Carstens
BC
Knowles
LL
Distribution modelling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses
Journal of Biogeography
 
2007
34
1833
1845
Richmond
OMW
McEntee
JP
Hijmans
RJ
Brashares
JS
Is the climate right for Pleistocene rewilding? Using species distribution models to extrapolate climatic suitability for mammals across continents
PLoS One
 
2010
5
e12899
Robinson
RC
Invasive and problem ferns: a European perspective
International Urban Ecology
 
2009
4
83
90
Sandel
B
Arge
L
Dalsgaard
B
et al.  
The influence of Late Quaternary climate-change velocity on species endemism
Science
 
2011
334
660
663
Schneider
H
Russell
SJ
Cox
CJ
et al.  
Chloroplast phylogeny of asplenoid ferns based on rbcL and trnL-F spacer sequences (Polypodiidae, Aspleniaceae) and its implications for biogeography
Systematic Botany
 
2004
29
260
274
Schorr
G
Holstein
N
Pearman
PB
Guisan
A
Kadereit
JW
Integrating species distribution models (SDMs) and phylogeography for two species of Alpine primula
Ecology and Evolution
 
2012
2
1260
1277
Schurr
FM
Pagel
J
Cabral
JS
et al.  
How to understand species' niches and range dynamics: a demographic research agenda for biogeography
Journal of Biogeography
 
2012
39
2146
2162
Shepherd
LD
De Lange
PJ
Perrie
LR
Multiple colonizations of a remote oceanic archipelago by one species: how common is long-distance dispersal?
Journal of Biogeography
 
2009
36
1972
1977
Solomon
S
Qin
D
Manning
M
et al.  
Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change
 
2007
Cambridge
Cambridge University Press
Thuiller
W
Lavorel
S
Araujo
MB
Niche properties and geographical extent as predictors of species sensitivity to climate change
Global Ecology and Biogeography
 
2005
14
347
357
Thuiller
W
Lafourcade
B
Engler
R
Araujo
MB
BIOMOD – a platform for ensemble forecasting of species distributions
Ecography
 
2009
32
369
373
Tutin
TG
Burges
NA
Schater
AO
et al.  
Flora Europaeae. Volume I: Psilotaceae to Platanaceae
 
1993
Cambridge
Cambridge University Press
Wang
L
Schneider
H
Wu
Z
He
L
Zhang
X
Xiang
Q
Indehiscent sporangia enable the accumulation of local fern diversity at the Qinghai-Tibetan Plateau
BMC Evolutionary Biology
 
2012
12
158
Wang
L
Wu
ZQ
Bystriakova
N
et al.  
Phylogeography of the Sino-Himalayan fern Lepisorus clathratus on ‘the roof of the world
PLoS One
 
2011
6
e25896
Wolf
PG
Schneider
H
Ranker
TA
Geographic distributions of homosporous ferns: does dispersal obscure evidence of vicariance?
Journal of Biogeography
 
2001
28
263
270
Wubs
ERJ
de Groot
GA
During
HJ
et al.  
Mixed mating system in the fern Asplenium scolopendrium: implications for colonization potential
Annals of Botany
 
2010
106
583
590
Yeh
FC
Boyle
TJB
Population genetic analysis of co-dominant and dominant markers and quantitative traits
Belgian Journal of Botany
 
1997
129
157
Zhu
K
Woodall
CW
Clark
JS
Failure to migrate: lack of tree range expansion in response to climate change
Global Change Biology
 
2012
18
1042
1052

Author notes

Present address: Centre for Biotechnology, University Bielefeld, D-33615 Bielefeld, Germany.
Present address: Museum für Naturkunde, D-10115 Berlin, Germany.

Comments

0 Comments