Glacial refugia and the prediction of future habitat coverage of the South American lichen species Ochrolechia austroamericana

The biogeographic history of lichenized fungi remains unrevealed because those organisms rarely fossilize due to their delicate, often tiny and quickly rotting thalli. Also the ecology and factors limiting occurrence of numerous taxa, especially those restricted in their distribution to tropical areas are poorly recognized. The aim of this study was to determine localization of glacial refugia of South American Ochrolechia austroamericana and to estimate the future changes in the coverage of its habitats using ecological niche modeling tools. The general glacial potential range of the studied species was wider than it is nowadays and its niches coverage decreased by almost 25% since last glacial maximum. The refugial areas were covered by cool and dry grasslands and scrubs and suitable niches in South America were located near the glacier limit. According to our analyses the further climate changes will not significantly influence the distribution of the suitable niches of O. austroamericana.

and habitat requirements 26 . Moreover, species distribution modeling methods were used in more or less local conservation planning [27][28][29] . Our work is the first one, which presents data on the potential distribution from the past to the future as well the first such dealing with lichenized fungi in South America. Because so far MaxEnt application was rarely used in lichen species distribution modeling the additional aim of the present study was to evaluate differences in the models created using this application with various input data.

Results
Ecological niche modeling evaluation. The area under the ROC curve (AUC) for each model created based on different dataset was calculated to estimate the trustworthiness of the analysis. AUC values were calculated by MaxEnt application automatically and it describes the probability that the model scores a random presence site higher than a random background site 30,31 (Phillips et al.). All projected niche models received high AUC scores ( Table 1) that indicates a high reliable performance of the analysis. Models were then projected onto the spatial layers representing the current, past and future climatic conditions. The niche overlap test was also used to evaluate differences between models created with all available climatic and altitudinal data and those based on reduced "bioclims" dataset. The obtained results indicated some incongruities between the created models -for LGM simulation the niche overlap statistics were calculated as I = 0.899, D = 0.699 and for the present time I = 0.899, D = 0.621. The internal tests of GARP models showed high reliability of these analysis (Supplementary Table S1). Current potential range. According to the most reliable model (Fig. 1), the niches of the studied species are distributed from Venezuelan Cordillera de Merida along the Andean range to southern Argentina with the highest concentration of the suitable habitats in Central Andes. Somewhat less appropriate conditions are found in the area of Altiplano. In Northern America the proper niches are found along the Pacific coast, west of Coast Ranges. Less suitable habitats are found along Sierra Nevada. In Mexico some potentially available niches are located in the southern part of Sierra Madre. Moreover, the model indicate the existence of suitable niches outside the Neotropics -in French Massif Central, Ethiopian Highlands, African Great Karoo as well as in southern New Zealand and Australian Great Dividing Range. Two other models which received lower AUC scores are presented in Supplementary Figure S2.
Glacial refugia. The general distribution of the suitable niches of O. austroamericana in LGM in the most reliable "All" model is congruent with the present potential range of this species (Fig. 2). There are, however, several regions which current bioclimatic conditions are not appropriate for this lichen. In LGM potentially available habitats were located in Falkland Islands, Atlantic coastal regions of Argentina and Uruguay as well as in south-western Iceland and Ireland. The general glacial potential range of O. austroamericana was apparently wider than it is nowadays.
In South America the areas indicated in the models as potentially available for O. austroamericana were covered by cool and dry grasslands and scrubs 32 . Interestingly, the suitable niches in North America were located near the glacier limit.
Two other models which received lower AUC scores are presented in Supplementary Figure Table 2). The most optimal altitude for the species occurrence is about 4500 m. In the models with reduced climatic and altitudinal datasets the temperature seems to be crucial for the studied species distribution. In "SelArea" models the most important variables contributing in the models were mean temperature of the wettest quarter (bio8) and annual mean temperature (bio1). The latter factor together with isothermality (bio3) were critical for "SelLay" models. The niche overlap statistics in "SelLay" is very similar to that calculated for the most reliable models (D = 0.629, I = 0.875), but when considering exclusively South America ("SelArea") the niche overlap is higher receiving the values of D = 0.725 and I = 0.922.
Further changes. Apparently, the further climate changes will not significantly influence the distribution of the suitable niches of O. austroamericana. The most reliable models created for three different climate change scenarios are presented in Figs 3 and 4. Two other models which received lower AUC scores are presented in Supplementary Figure S4.
The most insignificant global changes in the geographical distribution of the suitable niches in comparison to the present time will be observed in A2a scenario (Table 3) based on the conducted niche overlap test. The coverage of the most proper habitats, with the suitability above 0.7, will decrease in their surface in A1b and B2a models while their coverage in A2a scenario will increase in both analysed areas -whole World and South America only. The coverage of the areas characterized by the different suitability for the studied species is compiled in Tables 4 and 5. The similar tendency was observed based on GARP models analysis (Supplementary Figure S5, Supplementary Table S6).

Discussion
Species distribution models have become increasingly noticeable in ecological and biogeographical research 33 . Mostly because ecologists need ways of rapidly assessing the impacts of climate change on large numbers of species for which the occurrence data are often the only source of information 34 . While several critical opinion on ENM analyses were presented in the last years 35 , MaxEnt seems to be the most reliable application for modeling species distribution. Its usefulness was also tested in case of rare organisms 36,37 . Noteworthy, some of the studies indicating inappreciable usefulness of ENM tools 38 (Beale et al.) were called into question by the subsequent researchers 39 .
In our research we tested three approaches to evaluate the past, present and future distribution of the suitable niches of poorly known lichen species. Some of the previous studies indicated that the correlations between the environmental data used in the modeling should be reduced and that the correlated variables should be excluded from the analysis 40 . Also it was suggested that using restricted area in ENM analysis is more reliable than calculating habitat suitability in the global scale 41 . According to the received AUC scores the most reliable model was created based on all available climatic and altitudinal data and it was constructed for the whole globe. We therefore      would postulate that all potentially useful climatic variables and altitudinal data should be used in ENM studies, especially when ecological information about the studied taxon are poor and they do not allow to discriminate any climatic data as irrelevant. Noteworthy, despite slight differences in the trustworthiness of the three conducted analyses all created models indicated similar areas that could be occupied by Ochrolechia austroamericana. In our opinion lichen species distribution modeling with MaxEnt may be extended into new fields and it would be especially useful in reconstructing their past distribution and potential migration routes. While future habitat loss of the lichens became emerging question in lichenology 42,43 , so far MaxEnt was not implemented in any of those studies. While the limitations of such predictions are well-known and their verification is not possible in the present time, we believe that the vulnerability of specific areas indicated in the modeling should be taken into consideration prior to planning conservation actions.
Despite Ochrolechia austroamericana has been reported from rather limited number of localities 23,24 , the ENM method has shown its potential distribution range can be much wider.
The analyses have shown the suitability of habitats on almost every continent; however its occurrence is highly improbable in most regions as the genus Ochrolechia A.Massal. was a subject of several taxonomic treatments over the past c. 30 years 24,25,44-48 and O. austroamericana was never found outside South America. In general, most Ochrolechia species have distributions rather restricted to one continent or region, and that can be possibly related to the limited dispersal of diaspores (ascospores), which are relatively large 49 . Only very few species, e.g. the tropical O. africana Vain., are known to have wider range 25,45 but if they truly represent one evolutionary lineage or several cryptic species, has not been settled yet. However, in the light of recent studies, which showed several lichenized fungi to represent numerous phylogenetically distinct clades 50,51 , we suspect this scenario to be more adequate in this case.
Concerning the most probable potential distribution range we consider that O. austroamericana occurs only South America from Venezuela to southern Argentina with, as our modeling has shown, the highest concentration of the suitable habitats in Central Andes. As it appears to be restricted to cooler climate conditions, it could not spread more northward due to the lack of appropriate spreading passages in Central America.
Apparently, the current distribution of the suitable habitats of O. austroamericana results from location of its glacial refugia and no long-distance dispersal of this lichen is observed. Its niches coverage decreased by almost 25% since LGM. Most of the loss is observed within the Pampas and in the high regions of Central and Southern Andes. We interpret the first issue as related with the coast line regression after Late Glacial that significantly affected local climatic conditions 52 . The loss of the habitats along the Andes may be caused by the warming of the high-Andean regions that was documented for numerous South American regions and which has been intensified in the last three decades, and, as the consequence, the vertical shift of vegetation belts in the altitudinal gradient, i.e. uppering the forest line [53][54][55][56][57][58][59][60] .
Surprisingly, the future distribution of O. austroamericana will not change much from the present range of the species. The difference between the three created models for 2080 differ in the suitability of the Atacama Desert and lower parts of the eastern slopes of Central Andes for the studied species occurrence. The studies of Boulanger et al. 61 indicated that in twenty-first century the amplitude of the seasonal cycle will tend to increase in southern South America, while in northern South America the amplitude of the seasonal cycle would be reduced; that explains the relatively little changes in the general distribution of the species. However, while all emission paths tend to show the same pattern of warming, the highest warming is predicted in A2 storyline. This scenario is the only one which predicts increase of the suitable niches cover for O. austroamericana. We interpret it as a result of a tree-line location change 62 as the lowering of the forest limit due to the climate change will reduce vegetation cover and may expose additional areas appropriate for O. austroamericana.

Materials and Methods
The list of localities used in the modeling was compiled based on the examined samples exclusively. As the source of data we used specimens, including type material, from herbaria B, BM, GZU, H, KRAM, LPB, PRA, and UGDA, which were revised by the first author. All samples were investigated by thin-layer chromatography to detect diagnostic secondary lichen metabolites, crucial for the identification of Ochrolechia species. Methods were followed Orange et al. 63 . Most of those records have been already published by Messuti and Lumbsch 24 and Kukwa et al. 23 (data from Bolivia and Ecuador). Only the locations that could be precisely placed on the map were applied   in the ENM analysis. The database of a total of 44 specimens was created from 21 different localities (Fig. 5). One specimen was excluded from the study as most probably it was mislabelled (the specimen was supposed to be collected in the middle of town in area without rock outcrops, the only substrate which O. austroamericana inhabits). The ecological niche modeling was conducted using maximum entropy method implemented in MaxEnt version 3.3.2 18,64,65 based on the species presence-only observations. In total, 20 different locations of O. austroamericana were used (Supplementary Table S7), which is more than the minimum number required to obtain reliable predictions in MaxEnt application 66,67 .
Three different approaches were used to conduct analysis. In the first one ("All") as input data, all 19 climatic variables ("bioclims", Table 6) in 2.5 arc minutes (± 21.62 km 2 at the equator) developed by Hijmans et al. 68 were used together with the altitudinal data (Alt) and the models were created for the whole globe. From the second dataset ("SelLay"), we removed altitudinal data and seven "bioclims" due to their significant correlation (above 0.9) as evaluated by the Pearsons' correlation coefficient calculation computed using ENMTools v1.3. The following variables were excluded from the dataset: alt, bio6, bio7, bio9, bio10, bio11, bio16 and bio17. This analysis was also made for the whole globe. The last set of models ("SelArea") was made using the same, reduced number of variables and the area was restricted to longitude of 95°-30°W and latitude of 13°N-60°S.
In all analysis the maximum iterations was set to 10000 and convergence threshold to 0.00001. For each run, 15% of the data were used to be set aside as test points 69 . The "random seed" option which provided random test partition and background subset for each run was applied. The run was performed as a bootstrap with 1000 replicates, and the output was set to logistic. All operations on GIS data were carried out on ArcGis 9.3 (ESRI). The bioclimatic data for the LGM were developed and mapped by Paleoclimate Modeling Intercomparison Project Phase II 70 based on an atmosphere-ocean coupled general circulation model (AOGCM).
To assess potential range of O. austroamericana in 2080 three different climate change scenarios were applied into modeling 71,72 . A1b (CCCMA-CGCM3 simulation), A2a (CCCMA-CGCM2 simulation) and B2a (CCCMA-CGCM2 simulation). A1b scenario is characterized by the balance across all energy sources (where balanced is defined as not relying too heavily on one particular energy source, on the assumption that similar improvement rates apply to all energy supply and end-use technologies). The A2 storyline describes a highly heterogeneous future world with regionally oriented economies. The main driving forces are a high rate of population growth, increased energy use, land-use changes and slow technological change. The B2 is scenario with a general evolution towards environmental protection and social equity. The datasets used in the analysis are available on CIAS website (http://ccafs-climate.org). The coverage of the suitable niches were calculated in order to measure further changes in the habitat availability for both whole world and South America only.
The corresponding analyses based on dataset that received the highest AUC scores using GARP algorithm applying the best subsets 73 were conducted to verify MaxEnt models for the present and future time. These models were created using openModeller 74 . The niche overlap tests implemented in ENMTools application were used to evaluate the overlap of the potential ranges modeled for LGM and present time as well as to estimate differences between tested models "All" and "SelLay". The overlap was measured using Schoener's D (D) 75 and I statistic (I) 76,77 . Schoener's D was developed initially to compare diet and microhabitats and here it is used with assumption that direct measures of local species density are compared with each other. The I statistic is basing on Hellinger distance and measures the ability of the model to estimate the true suitability of habitat. Both metrics range from 0 (niches are completely different) to 1 (overlap).