CREST (Climate REconstruction SofTware): a probability density function (PDF)-based quantitative climate reconstruction method

. Several methods currently exist to quantitatively reconstruct palaeoclimatic variables from fossil botanical data. Of these, probability density function (PDF)-based methods have proven valuable as they can be applied to a wide range of plant assemblages. Most commonly applied to fossil pollen data, their performance, however, can be limited by the taxonomic resolution of the pollen data, as many species may belong to a given pollen type. Consequently, the climate information associated with different species cannot always be precisely identiﬁed, resulting in less-accurate reconstructions. This can become particularly problematic in regions of high biodiversity. In this paper, we propose a novel PDF-based method that takes into account the different climatic requirements of each species constituting the broader pollen type. PDFs are ﬁtted in two successive steps, with parametric PDFs ﬁtted ﬁrst for each species and then a combination of those individual species PDFs into a broader sin-gle PDF to represent the pollen type as a unit. A climate value for the pollen assemblage is estimated from the likelihood function obtained after the multiplication of the pollen-type PDFs, with each being weighted according to its pollen per-centage.Totest its performance, we have applied the method to southern Africa as a regional case study and reconstructed a suite of climatic variables (e.g. winter and summer temperature and precipitation, mean annual aridity, rainfall sea-sonality). The reconstructions are shown to be accurate for both temperature and precipitation. Predictable exceptions were areas that experience conditions at the extremes of the regional climatic spectra. Importantly, the accuracy of the reconstructed values is independent of the vegetation type where the method is applied or the number of species used.


Introduction
Reconstructing past climates, while being an important element in the global effort to understand climate system dynamics and their potential future structure and characteristics, is often limited to qualitative assessments of past conditions. This limits the potential for comparisons with the general circulation model (GCM) simulations, and the integration of palaeoenvironmental information in modelling initiatives (Braconnot et al., 2012). As a result, while inconsistencies exist both between GCM simulations and between GCM simulations and fossil records, it is difficult to use the bulk of the palaeodata available to evaluate GCM simulations in an efficient and effective way.
Many techniques have been developed to quantitatively reconstruct past climates from palaeobotanical data (Guiot et al., 1993;Huntley et al., 1995;Overpeck, 1985;Kühl et al., 2002). They rely on the fundamental hypothesis that a causal relationship exists between the modern distributions of plants and the associated climates (Jackson and Williams, 2004, and references therein). These techniques can be divided into two types: (1) those based on plant assemblages -

Published by Copernicus Publications on behalf of the European Geosciences Union.
modern analogue technique (MAT) (Overpeck, 1985;Guiot, 1990), weighted averaging (WA) (ter Braak and Looman, 1986;ter Braak and van Dame, 1989;Birks et al., 1990), weighted averaging-partial least-squares regressions (WA-PLS) (ter Braak and Juggins, 1993), artificial neural networks (ANNs) (Malmgren et al, 2001) or regression trees (Salonen et al., 2012) -and (2) those based on plant distributions -mutual climatic range (MCR) (Atkinson et al., 1987;Sinka and Atkinson, 1999;Elias, 1997), the coexistence approach (CA) (Mosbrugger and Utescher, 1997;Utescher et al., 2014) or probability density functions (PDFs) (Kühl et al., 2002;Truc et al., 2013). These methods are fully detailed in Birks et al. (2010). Amongst these methods, MAT, WA and WA-PLS are the most commonly used. To date, no consensus has been reached as to which performs the best, and since the paper of Telford and Birks (2005) the debate has been focussed on the sensitivity of different methods to spatial autocorrelation in the training data set, which can cause an overestimation of the performance of a method (Telford and Birks, 2005Guiot and de Vernal, 2011a, b).
Beyond these conceptual issues, the calibration data set is another key concern. Methods based on plant assemblages need to be calibrated on robust modern data sets covering various environmental and biotic conditions, data sets that are not available from many parts of the world (e.g. drylands where pollen is often poorly preserved in surface samples). In such situations, the flexibility of methods based on plant distributions (MCR, CA, PDFs) becomes more evident, expanding the range and scope of "reconstructible" environments. They also offer the possibility to reconstruct climate from non-analogue assemblages, provided that most species from that palaeoassemblage still exist. Conceptually, PDF-based methods evolved from MCR techniques as a way to model the strength of the relationship between plants and climate. Indeed, MCR (which considers a rectangular envelope defined by minimum and maximum values for a given climate variable) can be seen as the simplest PDF-based method. These methods are based on the correlation between plants' modern geographical distribution and climate gradients, with the climate value that is the most common in the plant distribution being its "optimum". Among the approaches that have already been proposed within the last decade (Kühl et al., 2002;Gebhardt et al., 2007;Truc et al., 2013), a recurrent issue concerns the assumptions made about the morphological characteristics of the envelope (width, skewness, central tendencies). Kühl et al. (2002) fitted a multidimensional Gaussian surface that excluded both multimodality and asymmetry, which are however common features when dealing with botanical assemblages. Later, Gebhardt et al. (2007) proposed to fit mixture models (combination of several Gaussian surfaces) to relax the constraints of a unimodal Gaussian shape, and more recently Truc et al. (2013) proposed the application of non-parametric PDFs to improve the fit between PDF and data.
In addition to the issue of the shape, the accuracy of such models is also a function of the taxonomic resolution at which pollen can be identified (usually family to generic level) and the number of species making up a given pollen type. Pollen types often become climatically non-informative due to a saturation effect wherein too many species result in the climatic information conveyed by each species being averaged and lost. Contrary to the problem of the shape of the climate envelope, the problem of low taxonomic resolution has rarely been discussed as its effects are usually not significant when plant diversity is relatively low. However, in areas where pollen types can comprise a high number of plant species ( > 30), it becomes increasingly significant and can result in saturated PDFs. Truc et al. (2013) proposed a species selection method (SSM) that recursively alters the taxonomic composition of a pollen type by taking into account the coexistence with other pollen types. In order to minimize PDF saturation, the SSM removes species that have climate requirements that are different from that of the assemblage.
However, the SSM only removes species with optima at the extremes of climatic gradients, leaving a certain number of climatically undifferentiated species around the median climate. We believe that the problems of PDF shape and plant diversity are in fact intimately related to the strategy used for fitting PDFs. A pollen type is not a homogeneous ecological unit in the sense that many species with different climate requirements can be classified in the same pollen type. From this point of view, fitting a density function directly to a pollen type is questionable. On the basis that species are the ecological units that respond to climate gradients, we propose a two-step procedure to define the PDF of the pollen types: (1) univariate and unimodal parametric PDFs are fitted for the species (PDF sp 's), and (2) those parametric PDF sp are combined to produce the PDF of the pollen-type (PDF pol ). The PDF pol reflects the diversity that exists among its species by considering independently each species. To reconstruct a climate value, we propose to combine the PDF pol with a weighted geometric mean. The multiplication of PDF pol ensures the conservation of the mutual climatic range.
To quantify the method's capability to reconstruct different variables in different environments, we have reconstructed a set of modern climatic conditions (20 variables) over a large area (3389 quarter-degree grid cells representing southern Africa). Southern Africa -composed of South Africa, Botswana, Lesotho, Swaziland and Namibia -is well suited for such a test as it is characterized by a strong topographic, geological and climatic heterogeneity (Tyson, 1986;Partridge and Maud, 2000;Chase and Meadows, 2007), leading ultimately to a great diversity of plant species (Goldblatt and Manning, 2002). Statistical tests were performed on climate anomalies (1) to analyse where and why the model was reliable, and (2) to measure the effects of parameters, such as the type of variable, the number of taxa used and/or the vegetation type.
Clim. Past, 10, 2081-2098, 2014 www.clim-past.net/10/2081/2014/ The method presented here has been implemented in a software package entitled CREST (Climate REconstruction SofTware). CREST -presented in Appendix A -is intended to make quantitative climate reconstructions more accessible to the wider community. Our hope is that a proliferation of quantitative reconstructions of past climate conditions will facilitate the consideration of palaeoenvironmental data in the assessment of GCM performance and ultimately allow for an improved understanding of both past and potential future climate change.

Methodology
The climate reconstruction method we propose is based on univariate PDFs. Here, a PDF represents the probability of a species existing along a climate gradient and is a surrogate for the species' realized niche (see for example Kearney, 2006). The process follows three general steps: (1) the plant-climate relationship is quantified (i.e. PDFs are fitted), (2) information conveyed by each taxon is combined and finally (3) a climate value from the resulting climate likelihood function is extracted. This method relies on the assumption that the plant-climate relationship has remained relatively constant since the deposition of the fossil assemblage.

Fitting of the PDFs
This step is critical for all PDF-based methods. Many different strategies have been proposed (Kühl et al., 2002;Gebhardt et al., 2007;Truc et al., 2013), all of them fitting a PDF to the pollen types identified in the fossil record. This strategy leads to a loss of certain information because (1) individual signals are mixed and (2) rare species are masked by the most extended ones.
Here we propose a two-step procedure to fit PDFs that better integrates the diversity that can exist within some pollen types. First, we fit a PDF to each species (noted PDF sp ), and secondly we combine the PDF sp into PDF pol . The latter considers more clearly the pollen type's diversity.

Creating PDF sp
Based on observations, we propose that distributions of climatic values where a species is found -its niche -can be classified into two shapes: a log-normal shape (Fig. 1a) or a normal shape (Fig. 1b) (Austin, 1987;Austin and Gaywood, 1994;Hirzel and Le Lay, 2008). The normal shape is symmetric, while the log-normal shape is markedly right-skewed (left-skewed distributions have been observed but are uncommon). In addition, the log-normal function is null for negative values, which is of interest when modelling variables such as rainfall amounts. Both curves are defined by two parameters: the mean x sp (Eq. 1) and the variance s 2 x,sp (Eq. 2) of the species niche, with x being the studied climatic gradient.
The estimation of x sp and s 2 x,sp can be biased if the abundance of each climate value is not considered (Telford and Birks, 2011a). Following the model of Kühl et al. (2002) on that particular point, we propose dividing the climatic space into bins of equal width. All the climate values (a total of N) are sorted into J bins. The number and/or length of the bins is variable and depends on the dispersion of the climate values. A weight k j is defined for each bin as the ratio of N with the number of pixels n j in the bin j (Eq. 3). Each climate value from bin j will have the same weight k j .
The shape and the position of the PDF sp along a gradient can be calculated with Eqs. (4) and (5) representing the normal law and the log-normal law ( Fig. 1b and a), respectively.

Creating PDF pol
To create the PDF pol , all the PDF sp 's are added with a weight determined by their geographical extent (represented by the number of grid cells occupied n i , Eq. 6). Due to the absence of detailed, functional information regarding the pollen production of the different plant species, we are forced to consider that it is a constant among all the species of a pollen type and thus that each species is likely to have equally contributed to the observed pollen biomass. While PDF sp 's have an imposed shape, PDF pol 's are not constrained since no assumptions are being made. A PDF pol can thus be multimodal when the diversity within the pollen type results in two or more climatically separated groups of species. Figure 1c and d highlight the advantage of that method: for instance, the climatic signals conveyed by the three species Tribulus cristatus, T. pterophorus and T. zeyheri are not masked by the signals of the most extended one, T. terrestris. . An example of PDF fitting for two variables (Prec dry Q and Tmean ann) for the pollen type Tribulus, which is composed of four species in our database. Four PDF sp 's are then fitted for each variable (a and b) and combined to create the PDF pol (c and d). The dashed lines on (c) and (d) are the PDFs obtained by Truc et al. (2013). The difference between the two methods is more marked for prec dry Q, where (i) the PDF pol is null for negative precipitation values (more realistic) and (ii) the optimum is more marked and reflects the optima of the different species.

Combination of the PDF pol to create the PDF var
We propose the combination of different PDF pol 's with a weighted geometrical mean (Eq. 7). The multiplication of PDF pol 's ensures that the reconstructed climate value will be in the mutual climate range of the taxa considered. In addition, since plants may pollinate more when they live close to their climate optimum (Birks and Seppä, 2004;Jackson and Williams, 2004), the PDF pol 's are weighted according to a monotonically increasing function of their pollen percentage ω pol (s), with s representing a sample (Eq. 8).
Using pollen percentages p pol (s) to weight taxa is problematic, as it is with traditional interpretive techniques, because pollen production can vary substantially from one plant to another (Jackson and Williams, 2004). The pollen production is unknown for the vast majority of plant species, and one cannot effectively employ this information to fit the PDF pol 's. To address this issue, we follow the method developed by Truc et al. (2013), wherein percentages were rescaledbetween 0 and 1, with 1 corresponding to the highest percentage observed for the pollen type. This normalization is, however, very sensitive to outliers. In CREST, percentages are rescaled by the mean percentage of the pollen type when it is present in a sample. In other words, to calculate the mean, only strictly positive values are considered (Eq. 8). For a given pollen type, our weights have a correlation of 1 with those of Truc et al. (2013); the difference lies in the relative weights between taxa.

Climate reconstruction
The reconstructed climate corresponds to the abscissa x(s) of the optimum of PDF var (Eq. 9). PDF var describes the likelihood of any climatic value to be the target value when considering the presence of many pollen types.

Error estimations
PDF var 's provide access to the complete distribution of errors. They can be estimated at different thresholds (noted α).
The α% confidence interval (CI) is more appropriate than a standard deviation because PDF var 's are rarely symmetrical (Fig. 2).

Validation
As a case study, we have used a modern botanical database to reconstruct a set of contemporary climate values to highlight and explore the strengths and weaknesses of the approach, and to quantify its accuracy and robustness. Using southern Africa as a study area, we consider five countries:  (2007): the winter rainfall zone (WRZ; > 66 % winter rain), the summer rainfall zone (SRZ; < 33 % of winter rain) and the year-round rainfall zone (YRZ) in between.

Climate system
Most of southern Africa is dominated by summer rainfall related to the seasonal dynamics of the Intertropical Convergence Zone (ITCZ) and the advection of moist tropical air masses off the Indian Ocean. Annual rainfall is highest along the eastern escarpment (∼ 1200 mm yr −1 ; Hijmans et al., 2005;Mitchell and Jones, 2005) and decreases westward. Conversely, in the Cape region (southern tip of Africa), most of the rain falls during the winter months as a result of frontal systems embedded in the southern westerlies (Tyson, 1986) and can reach annual totals of more than 900 mm yr −1 . A complex mosaic of rainfall regimes are found at the boundary between those two systems: from year-round rainfall along the south coast of South Africa (> 900 mm yr −1 distributed in more than 100 rain events per year) to the super-arid Namib Desert (< 20 mm, < 10 rain events). The orographic effects of the Drakensberg escarpment and the Cape Fold Belt are very marked, creating a strong rain shadow effect in their lee.
The west coast is cooled by upwelling associated with the northward-flowing Benguela Current, whereas the south and east coasts are warmed by the southward-flowing Agulhas and Mozambique currents, respectively. At a given latitude, the difference in temperature between the two coasts can exceed 6 • C. The greatest diurnal temperature ranges are found in the interior, especially in the Kalahari and the Karoo region, where the altitude is greater than 1000 m a.s.l. in many areas (Hijmans et al., 2005).
The study region currently supports four primary biomes: deserts and xeric shrublands (54.7 %); montane grasslands and shrublands (16.8 %); tropical and subtropical grasslands, savannas and shrublands (25.3 %); and Mediterranean forests, woodlands and scrub (3.2 %) (Olson et al., 2001). The latter is better known as the Cape Floristic Region, which is dominated by the Fynbos Biome. Each biome is divided into ecoregions (Fig. 3), which will be used to describe the model's properties.

Data
We have extracted botanical data for all grid cells where at least one plant with more than 25 pixels in its distribution had been recorded, leading to a total of 3389 "samples" (Fig. 4). We have then selected 20 climatic variables of interest: 9 temperature-like and 11 moisture-like variables (Table 1). A total of 4969 species distributions have been used.

Botanical data
Botanical data were extracted from a series of databases held by the South African National Biodiversity Institute (SANBI, 2003;Rutherford et al., 2003Rutherford et al., , 2012. The data from these sources, which are derived mainly from herbarium collections and documented observations, are most commonly available as "presence" within a particular 0.25 • ×0.25 • grid square. We have used this resolution for our analyses, upscaling more precisely located data to this common resolution. Data were obtained from field surveys performed during the late twentieth century, between 1970 and 2000 with a peak in the 1980s. In this study, we only consider species with at least 25 occurrences, resulting in a number of species (n sp ) available per pixel between 1 and 1371 (median = 47). This strong heterogeneity is mainly due to both the range of environments found in our study area (Fig. 3) and the strong difference that exists between the different countries ( Fig. 4), with South Africa providing by far the most extensive data set.

Climatic data
To define PDFs, the species distributions have to be associated with climate data. For this study we have reconstructed 20 climate variables from our data set. While some of those variables represent climatic features playing a strong role in a plant life cycle (e.g. number of frost days or different precipitation variables), the impact of some others is much more indirect (e.g. mean diurnal range or temp ann range). The purpose of this strategy is to assess the extent to which different climate variables can be reliably reconstructed with CREST.
Climate variables for this study were obtained from WORLDCLIM1.4 (Hijmans et al., 2005), which, along with monthly precipitation and temperature data, provides a data set of 19 bioclimatic variables that are consideredimportant Clim. Past, 10, 2081-2098, 2014 www.clim-past.net/10/2081/2014/  Mitchell and Jones (2005) elements in studying the eco-physiological tolerance of plants species. These data were then upscaled to match the resolution of the botanical data (0.25 • × 0.25 • ). Additional variables of interest have also been derived from WORLD-CLIM's data, including the soil water content (SWC; Trabucco and Zomer, 2010) for both summer and winter, the mean annual aridity (Trabucco and Zomer, 2009) and winter rainfall percentage (WRP). We have also used two variables from the CRU 2.10 (Climate Research Unit) data (Mitchell and Jones, 2005): the number of frost and wet days during the year. Those data (0.5 • ×0.5 • grid cells) were downscaled to meet our resolution. All the climatic values used in this study are representative of the period 1960-1990; the period is extended to 1950-2000 in some situations.
The description of all variables as well as their original reference is summarized in Table 1.

Accuracy of the model
We have measured the climate anomalies δ v (s) for each sample s and each variable v between the reconstructed climate Recon v (s) and the instrumental value Instru v (s) according to Eq. (10). A positive/negative anomaly is equivalent to an under/overestimation of the targeted climate.
The dispersion of the anomalies for each variable was calculated with the R software (as all the statistics presented here; R core team, 2014) and has been compiled in Table 2.
The distributions of anomalies are all centred around 0. Even if all the median values are statistically different from 0 (sign test, p < 10 −5 ), the relatively low value of these medians indicates that the model is not subject to undue bias. A major dichotomy can be observed between the two types of variables (χ 2 df =1 test with Yates' continuity correction, p = 1.8 × 10 −3 ): for the temperature-like variables, the median is positive for eight out of nine variables (general underestimation), while for moisture-like variables the opposite is observed, with negative medians for 10 out of 11 variables (general overestimation, Table 2). The different percentiles we have calculated provide insight regarding the dispersion of the reconstructed values, as do the histograms in Fig. 10. The skewness is most often negative (for 14 variables), meaning that when errors are negative (overestimation) their absolute value is higher than when they are positive (75 and 95 % percentiles respectively higher than the 25 and 5 %).
The root mean square deviation (RMSD v ; Eq. 11) is an index that reflects the mean error of a model, but it is sensitive to outliers. It does, however, allow for a good evaluation of the performance of the model. All the values are compiled in Table 2.
The amplitude of δ v (s) and RMSD v are functions of the variable range. Direct comparisons between variables cannot be performed -except for those with a similar range www.clim-past.net/10/2081/2014/ Clim. Past, 10, 2081-2098, 2014 . Four variables present a high NRMSD: mean diurnal range (0.75), temp ann range (0.72), prec seasonality (0.70) or temp seasonality (0.68). The climatic signal of these four variables does not seem to be well captured by the botanical data, and plant distribution is apparently not directly driven by those variables. They represent annual climatic variability, and a range of climatic scenarios could result in the same values. For example, the variable prec seasonality takes identical values for seasonal rainfalls whether they occur mainly in winter or summer. This major difference is incorporated into WRP, which has been reconstructed with a much better accuracy (NRMSD = 0.44).

Geographical analysis of the errors
Generally, southern African climates are accurately reconstructed with CREST. The anomalies that do exist are not randomly dispersed throughout the study area. On the contrary, regions of enhanced or diminished error are observed for each variable (Figs. 5 and 6). On these figures, anomalies have been normalized by the RMSD (Eq. 13) to make all the maps comparable.
This observation is validated with the measure of the spatial autocorrelation of the anomalies with Moran's I (Moran, 1950) (Eq. 14 and Fig. 7). This index measures the (dis)similarity of nearby locations in space. To compute this index, a neighbourhood matrix of weights w is defined. We have measured the spatial autocorrelation at different distances: from a local perspective, where only adjacent grid cells are neighbours, to the continental scale, where all the grid cells are considered neighbours. Under the null hypothesis (no spatial autocorrelation), Moran's I is normally distributed. However, δ v (s) is not, and mean and standard deviations were estimated empirically with 999 permutations for this vector. We ran 50 tests for each variables (one for each distance) and applied a Bonferroni correction α = 0.05 50 . Most of the tests were highly significant (blue and red dots in Fig. 7).
The main feature is that the anomalies are highly correlated at local scales and that correlation decreases with distance. No large-scale structure is observed in the data (Figs. 5, 6 and 7). These results show that anomalies are spatially clustered: in some areas the model performed very well, while it was less reliable in others. Temperature anomalies are clustered more  on the local scale than precipitation anomalies (higher values at distances lower than seven-eight grid cells), but the correlation decreases faster with distance (no distinction beyond 15 grid cells). Seasonality of temperature and precipitation have a distinct pattern, being correlated at longer distances, highlighting -in association with high NRMSDs (Table 2) that their errors are not limited to distinct regions. Four areas present a group of outliers for several variables: (1) the Namibian coast (for temperature and precipitation), (2) the high mountains of Lesotho (for temperature and humidity variables), (3) the eastern part of the Great Escarpment (precipitation) and (4) the southern coast of South Africa (precipitation) (Figs. 5 and 6).

Factors impacting the reconstructions
As the errors are spatially clustered, we have looked for factors that could explain this distribution. There is no clear linear relation between the anomalies absolute values and n sp . The slopes of the linear models we fitted were statistically significant at the 5 % threshold but the R 2 were always low (3.3 % of variance explained on average, Table 3). The models can, however, be biased by the uneven distribution of n sp ; half of the grid cells were reconstructed with 47 or fewer species, while some others were reconstructed with more than 1000 (Fig. 4). Some of our clusters of errors are found in mountainous regions, and we have hypothesized that the errors may arise from a mix of low-and high-altitude plants, with the anomalies observed being proportional to the degree of mixing. Thus, we have calculated the intra-pixel variation of altitude (the standard deviation of all the 30 arc-second altitude values in each quarter-degree grid cell, later called Alt). We fitted linear models to explain the anomalies as a function of n sp and Alt. However, the gain of explained variance was relatively small (+0.9 % on average). These results indicate that the anomalies are not a result of the number of species used for the reconstruction.
We also considered the impact of vegetation type on the anomalies. We used the Olson et al. (2001) classification to assign a biome and an ecoregion to each grid cell (Fig. 3). We used an ordination technique called between-groups principal component analyses (PCAs) (Thioulouse et al., 1997), available in the R package ade-4 (Dray and Dufour, 2007), to reveal the differences that may exist between vegetation types. With all the variables considered in the same analysis, we measured whether the type of vegetation impacted the reconstructions. At the biome level (seven levels; Fig. 8), the between-groups variance only explained 9 % of the total variance, meaning that more than 90 % of the variance was not explained by the differences between the biomes. The length of the boxes in Fig. 8 highlights that there is more variance within each group than between them. The between-groups PCA run on ecoregions explains 25 % of the total variance, but this is low relative to the number of groups (25). Again, more variance remained within the groups than between them. Figure 9 summarizes the mean dispersion of errors within each ecoregion. Some ecoregions appear to concentrate outliers, but these are always composed of 25 or fewer samples (low geographical extension and/or low amount of botanical data). Thus, despite the high botanical diversity that exists in southern Africa, we were not able to demonstrate the type of vegetation (forests, grasslands, savannas, etc.) having any effect on the quality of the reconstructions.
The only factor that explains a significant part of the dispersion is the distance of the expected value from the most represented value of the variable over the study area (Table 3). We fitted linear models to explain the anomalies as a function of the expected value (Fig. 10). All were signifi-cant (p value < 0.001) with positive slopes. A noticeable difference between temperature-like (R 2 = 42 % on average) and moisture-like variables (R 2 = 22 % on average) is observed, indicating that values that lie far from the most represented climate exhibit the highest anomalies (on the left and/or right-hand side(s) on the x axes in Fig. 10).

Discussion
Our results indicate that the PDF-based method of CREST performs well (Table 2, Figs. 5 and 6), even if some differences in terms of reconstruction quality exist between variables. The variables that were best reconstructed were those that have a direct impact on the physiology of plants, and thus strongly constrain their distribution (e.g. Tmean wet Q, Frost days, Prec dry Q or Prec wet Q) (referred to as direct gradients by Guisan and Zimmermann, 2000). The impact of other variables such as mean diurnal range or temp seasonality on the plant life cycle is indirect. Thus, they are less likely to be accurately reconstructed from pollen data. Variables that are surrogates for direct gradient may show strong ability to describe modern data but poor predictive power to describe past conditions (see, e.g., the palaeolimnological example of Juggins, 2013). Selection of variable(s) of interest should always be conditioned to an appropriate analysis of the data. Statistics could help in that process (Mac Nally, 2000;Telford and Birks, 2011c), but the final decision about the variables to reconstruct should always derive from an enlightened choice based on both statistical and ecological/environmental con-siderations. In the semi-arid to arid environments of southern Africa, precipitation and/or water availability strongly constrain plants distributions, which probably explains why we get lower NRMSDs for moisture-related variables in our case study.
The method performed well regardless of vegetation type. We were not able to show any differences in accuracy between the different biomes and/or ecoregions, provided that the distribution of the biome and/or ecoregion was sufficiently spatially extensive. Based on these results, we have found that the method works best for vegetation types represented by at least ∼ 25 to 50 quarter-degree grid cells   (estimation based on Fig. 9) in order to adequately determine the plant-climate relationship.
While our expectation was that a high number of species would result in more precise reconstructions, we were not able to observe any relationship between anomalies and the number of species. Anomalies do decrease when the number of species begins to increase (from 1 to ∼ 20-30), but then the tendency is reversed, and the large anomalies were observed in samples with the largest number of species. This may be related to a saturation problem, wherein more is not necessarily better. As we used a presence/absence weighting strategy, species far from their climate optimum have the same importance as those living in their optimal climate. The increase in the number of species could increase these marginal elements, biasing the reconstructions. The role of the number of taxa on the accuracy is not yet fully understood and is the subject of ongoing studies.
Other studies (Kühl et al., 2002;Scott et al., 2003;Truc et al., 2013) have shown that selecting a subset of the There is more dispersion within each biome (length of the boxes) than between, confirming the results of the between-groups PCA (90 % of variance not explained by the groups).
recorded taxa was sometimes more appropriate when attempting to capture a given climate signal. In order to improve the quality of the reconstructed variables, consideration should be given to reconstructing each variable with a different subset of the total of the available species list. Reducing this list to a shorter list of responsive species reduces the noise and consequently leads to better reconstructions. These choices, however, are not wholly objective. A first approach consists in observing directly the PDF pol 's since flat and multimodal PDFs may indicate insensitivity to a given climate variable. Scott et al. (2003) made choices based on considerations of the ecology of the given species, while Kühl et al. (2002) opted for a more statistical approach. To avoid using redundant information, only species that were statistically different (based on the Mahalanobis distance between the PDFs) were conserved. Truc et al. (2013) combined those two approaches and based their choices both on niche-based modelling and ecological considerations. The process of selection is thus not straightforward, and, while it may improve reconstructions, care needs to be taken to avoid undue bias in the results. The selection of the different sets of pollen types probably lowers the impact of spatial autocorrelation Birks, 2005, 2009). However, it should be kept in mind that the "supposed" reconstructed variable may in fact be a composite of primary and secondary variables. Secondary variables have the power to strongly bias reconstructions when spatially correlated with the variable 7 -S o u t h e r n A f r ic a m a n g r o v e s 2 a -D r a k e n s b e r g a lt i-m o n t a n e g r a s s la n d s a n d w o o d la n d s ** ** ** * ** * ** ** ** ** Figure 9. Box plots representing the dispersion of the normalized anomalies (Eq. 13) for each ecoregion. There is globally more dispersion within each ecoregion (length of the boxes) than between, confirming the results of the between-groups PCA (75 % of variance not explained by the groups). * means the ecoregion is composed of fewer than 50 grid cells; * * means fewer than 25. The numbers match those of Fig. 3. of interest (Juggins, 2013), so that only part of the reconstructed palaeovariability can be related directly to it. CREST (Appendix 1) provides a range of outputs that indicate the sensitivity of different taxa to given climatic parameters and allow the user to assess the data being considered and make informed choices in the selection of such subsets. When plotted on a map (Figs. 5 and 6), the reconstruction anomalies appear to be spatially clustered. Those patches of large anomalies can be explained by the position of the local climate along the climate gradients (Fig. 10) and are a direct consequence of the hypotheses underlying the model. The method is correlative, and consequently it is biased towards the best-represented climate values. This uneven sampling of the environmental gradients biases the estimation of plants' optima by shifting the "real" optimum towards the portion of the gradient with the most observations (Telford and Birks, 2011a). In most cases, lowest/highest values along the studied climate gradient are rare, but there are exceptions. For example, low rainfall amounts are common in southern Africa, and as a result they are well represented and the signal is easily captured by the model.
To offset the impact of the climate distribution's heterogeneity, we upweighted rare climate values as proposed by Kühl et al. (2002) and Truc et al. (2013). This method shifts PDF optima towards the rarest climate values. The climate abundance weighting did decrease the errors for the extreme climates but also increased them for the most common ones (data not shown). The overall impact is nevertheless positive since it decreased the RMSDs for all the variables. It also reduced the clustering of errors. Despite its advantages, the strategy has the drawback that artificial geographical limits must be selected (e.g. mountain ranges or country borders; Kühl et al., 2002) to compute the weights. A finite number of grid cells must be selected and sorted into bins. Any change in the boundaries would affect -potentially significantlythe weights, and thus the reconstructions. It is also possible that the climate abundance weighting may be the cause of the small but significant bias observed between temperature and moisture variables, which are, respectively, under-and overestimated for rare climates in the region.
Even with the climate abundance weighting, it is apparent that reconstructing the rarest climates is extremely complicated with models such as those described here. This is why, for example, the Cape region is poorly reconstructed for the precipitation-like variables but not for the temperature-like variables. The temperature of the area is common in southern Africa -so its signal is well captured -but it is an outlier in terms of quantity and seasonality of rainfall. Other areas of notable climatic rarity in southern Africa include (1) the eastern portion of the Great Escarpment (high precipitation); (2) www.clim-past.net/10/2081/2014/ Clim. Past, 10, 2081-2098, 2014 Figure 10. Anomalies plotted against their expected values. The density of points is heterogeneous: being very dense around the median climate and sparser at the extremes. This is illustrated by the histograms that represent the marginal distributions. The anomalies are smaller for the best-represented climate values and increase with distance from the median climate. The blue line represents the linear model fitted, with its associated R 2 .
the high mountains of Lesotho, where temperatures are very low and precipitation is high; (3) the thin coastal band along the southern coast of South Africa, where moist forests can develop as a result of significant aseasonal rainfall; and (4) along the Namibian coast (stable temperature and extremely low precipitation). All these areas lie at an extreme (relative to the study area) of one or several climatic gradients, giving rise to clusters of high anomalies.  Figure 11. Scatterplot representing a 2-D projection of the climatic space of southern Africa for the two variables Tmean ann and prec ann. In green and red are the modern positions of two fictious palaeoarchives. Those two points represent two very different situations relative to the climatic space: well-represented (green) vs. rare (red) climate. Reconstructing climate changes for the green palaeoarchive should be more accurate because it can "move" in several directions around its modern climate. However, the only major direction in which the red sample can move is towards warmer and drier conditions. Colder temperatures should be "reconstructible" but with an amplitude that may not reflect actual variability.
In terms of reconstructing quantitatively long-term climate variations it should be kept in mind that PDF sp is defined by the modern climatic space. Climatic space varies over time, and certain elements of some past climate regimes may be more or less abundant and/or more or less accessible to some species than in the modern climatic space (Veloz et al., 2012). Depending on the location of the site vis-à-vis the climatic space, the potential to estimate the amplitude of climate change varies. As shown schematically in Fig. 11, samples located in the mean climate space have greater potential to "move" in several directions and with greater amplitude than samples that are already at the margin of the climatic space. In the latter case, the exact amplitude of change may be underestimated, but the overall trends and direction of change may still be accurate. It is expected that, even under a different climate, the relative position of the different taxa along a climatic gradient would stay the same, so that the replacement in the past of a taxon by another that currently lives in colder environments will effectively indicate colder conditions with -possibly large -uncertainties regarding the amplitude of change (Veloz et al., 2012).

Conclusions
The PDF-based method we have presented in this paper provides robust results across a range of climates and vegetation types. We have demonstrated that the accuracy does not vary significantly as a function of vegetation type or the number of species considered, and it is thus a useful tool for reconstructing climates in many regions and biomes. The accuracy of the reconstructions is, however, strongly impacted by the climate variable being reconstructed (direct or indirect gradients) and primarily by the position of the targeted climate on the climate gradient of the study area. To ensure a robust reconstruction, one should 1. select climate variables that directly impact the distribution of the species and, inversely, use only species whose distributions are significantly defined by the climatic variable; 2. where possible, work with samples collected in widespread vegetation types to fit the most reliable PDFs; 3. define a climatically coherent study area to take advantage of the climate abundance weighting.
The results presented in this paper highlight our current understanding of the potential and limitations of the CREST method for reconstructing climates from botanical data. Recent work has shown the potential of the models upon which CREST has been based, particularly in regards to long-term climate reconstructions (Chase et al., 2015;Truc et al., 2013). Our goal with CREST is to make these techniques more accessible to the wider scientific community, and it is our hope that this tool will be applied to study other areas where long-term climate variations still need to be quantitatively described. CREST is freely available from the authors as well as at www.hyrax.univ-montp2.fr.