Studies exploring past changes in species realized niches have done so by relating historical species distributions to climate variables obtained from General Circulation Models (GCMs) (e.g., Pearman et al., 2007; Nogués-Bravo et al., 2008; Maiorano et al., 2013). Changes in these relationships have been studied either in time snapshots [usually 6000 and 21,000 years before present (BP)] or in transient simulations since the last glacial period. Notably, the goal of the simulations that produced the employed paleoclimate data was not to determine past climate but rather to examine distinct GCMs under different forcing conditions through comparison with reconstructed paleoclimate from observed data (i.e., Masson et al., 1999). There remains substantial disagreement among paleoclimate simulations and significant biases compared to observations exist, especially for rainfall (Braconnot et al., 2012), that may limit the suitability of these simulations for determining past climatic niches of species.

Here, we use January temperature (Tjan), quantitatively reconstructed over the past 10,000 years in Europe, in recognition of its ecological importance to fulfillment of chilling requirements for budburst and growth (Nienstaedt, 1967; Heide, 1993; Kramer, 1994, 1995; Harrington and Gould, 2015) and the influence of winter temperature on soil temperature, an important factor in determining treeline globally (Körner and Paulsen, 2004). Winter frost may damage needles and wood and, therefore, restrict the range of A. alba and F. sylvatica in areas where temperature descends below −10°C. P. abies is more frost tolerant than F. sylvatica and may withstand temperatures down to −17°C. Cold temperatures are a limiting factor for expansion of many tree species in the temperate latitudes, explaining about 80% of the variation in their range sizes (Pither, 2003). Cold temperatures also limit species distributions and diversity along elevation gradients worldwide (Körner and Paulsen, 2004). Furthermore, cold thermal limits of European Holarctic plants may be more conserved than warm limits (Pellissier et al., 2013), thus providing a conservative test for niche changes. Another important abiotic variable for plant species is water availability as influence by precipitation. Palaeoecological studies show that annual precipitation and its seasonal distribution may have complex interactions with temperature in shaping ecosystems over the long term (Cheddadi and Khater, 2016). In the European temperate zone where P. abies, A. alba, and F. sylvatica occur, temperature is a limiting factor for growth (see i.e., Kramer, 1994 for F. sylvatica; Heide (1974) for P. abies, and Maxime and Hendrik, 2011 for both A. alba and F. sylvatica), but precipitation is not. In the present study we have focused on January temperature (Tjan) rather than on precipitation as it has a direct impact on the three focal species and its future increase will probably affect these species more than the expected up to 20% increase of annual precipitation (IPCC, 2014).

Tjan was calculated using fossil pollen records from a network of sediment cores from across Europe, combined with empirical probability density functions of plant taxa identified for the fossil samples, following the method described by Chevalier et al. (2014). Compared to other pollen-based reconstruction methods, this technique avoids the problem related to the lack of analogy between past and modern ecosystems (Jackson and Williams, 2004).

We express taxa-climate relationships (Figure 1) as probability density functions (pdf, Kühl et al., 2002) that are obtained from modern species distributions (Figure 1A) and values of the focal climate variable (Tjan, Figures 1B,C) in these distributions. Fossil pollen grains are often identified to a genus or a family level which does not allow an accurate identification of the originating species. In order to use the modern species distributions and their related climate, the pdf method requires an assignment of each taxon to a plant species. In the present study we assigned fossil taxa to modern plant species that have the widest climate range and occur in the area. For instance, Pinus pollen grains have been assigned to the species P. sylvestris whose geographical range encompasses a wide temperature range. Pinus pollen grains cannot be identified to the species level which prevents the distinction between temperate and Mediterranean pines. Assigning a pollen grain to a species that has a restricted range (i.e., P. pinaster instead of P. sylvestris) would bias the climate reconstruction due to an unjustified reduction of the temperature range of the unknown species. Taxa pdfs (Figure 1D) for Tjan were built from a database of 270 georeferenced European plant species from published maps (Jalas and Suominen, 1973; Hultèn and Fries, 1986) and from WorldClim 10' gridded, interpolated weather station data for modern climate (1950–2000 period; Hijmans et al., 2005). The instrumental modern climate data-set used in the present study is set within a warming trend that is related to ongoing anthropogenic climate change (IPCC, 2014). Between years 1861 and 2000 global temperature increased by about 0.6°C (Folland et al., 2001). Over a comparable time period (1901–2000), the winter temperature (DJF) in Europe recorded a warming trend of ~0.08°C per decade (Luterbacher et al., 2004), which is slightly higher than the global value. The human-induced temperature increase recorded over the last century (< 1°C) has already affected the niche of many plant and animal species (Walther et al., 2002). However, this recent temperature increase is lower than the reconstructed amplitude of temperature change over the Holocene (Davis et al., 2003; Cheddadi and Bar-Hen, 2009).

Tjan was calculated from a pollen assemblage of n taxa in a fossil sample “s” as the temperature of maximum probability in the intersection of the n pdfs (Figure 1D), as follows: Tjan(s)=argmax[(∏tax1taxnpdf taxi,TjanWi(s))(∑Wi(s))-1] where:

Tjan(s) is the reconstructed climate for sample s;

pdf is the probability density function of plant taxon i (tax_i) to taxon n (tax_n) for the variable Tjan; n is the number of taxa identified in a sample;

W_i(s) is the weight associated with taxon i for sample s determined by the pollen percentage of taxon i. Pollen percentages were calculated after removing the pollen counts of the three taxa from the total pollen sum. argmax is the function that returns the value for which a function is maximum.

This phenomenological, statistical method does not account for, or depend on, variation in ecological or environmental processes (competition, phenology, dispersal, or evolutionary adaptation, etc.) or their relative importance in structuring paleocommunities. Similarly, estimated variable values do not depend on or account for values that are estimated for preceding periods.

Average Tjan values were calculated for these records to represent fixed time slices every 1000 years (±250 years). Fossil pollen records were selected from the European Pollen Database (EPD, http://europeanpollendatabase.net) for their time resolution and the quality of the associated calibrated age model. Pollen records with less than three radiometric datings or reversed dates were excluded as well as those with a sampling resolution that does not provide a series of at least one sample every 500 years over the covered period of time.

A conceptual difficulty of our climate reconstruction technique lies in its assumption that plant climate envelope is stable, which is what we aim to test. We partially circumvent this issue by excluding observations of the three focal tree species from the data-set used for paleoclimate reconstruction, so that any change in their thermal requirements would not affect the paleoclimate data. In order to evaluate the effect of the exclusion of the three focal taxa, we computed the difference between Tjan reconstructed from a modern pollen data-set and the same pollen data-set excluding the focal three species (Figure 2). In order to estimate the impact of a shift on the reconstructed Tjan, we used the Holocene Tjan range of the three species in the reconstruction procedure and calculated the difference with Tjan based on the full set of species using a modern pollen data-set (Figure 2).

In this preliminary analysis we find that the reconstructed Tjan values remain robust to the exclusion of the three focal tree species from the suite of species used to make the climate reconstruction. Other geological or biological indicators of temperature, such as isotope ratios or chironomid flies, might be considered independent climate proxies, their sparse geographical distributions unfortunately do not allow for spatial climate gridding at the extent of the European continent. Moreover, complete independence of these potential proxies with terrestrial vegetation is unlikely. Similarly, data from GCM simulations of past climates have cryptic dependencies with plant distributions because of the need to specify boundary conditions that depend on vegetation albedo, biomass distribution, and the carbon cycle.

The Tjan values obtained in fossil records were spatially interpolated for each time slice (10 to 3 ka) over Europe onto a 0.5° longitude × 0.5° latitude grid using universal kriging (R gstat package, Pebesma, 2004). The interpolated Tjan was then mapped as a smoothed surface using a spherical surface spline in tension (Smith and Wessel, 1990) within the generic mapping tools (GMT, Wessel et al., 2013; Figures 3A,B).

The number of time series from which Tjan is obtained increases from 10 to 3 ka, which may potentially affect the comparison of the estimated Tjan range through time. A similar issue has been raised and statistically tackled by weighting the observed temperature values by the availability of climate over a gridded geographical space (Broennimann et al., 2012). The method in the latter work is not directly applicable to our data since these are non-gridded temperature values that are determined by the location of pollen cores and the age of samples. Therefore, we cannot project species range between 10 and 3 ka to a standardized climate space. However, we were able to generate an analogous normalization of the distribution of temperature by weighting each observation by a proportion that represented occupation of similar climate space, which we present for Picea at 10 ka, 3 ka, and today as an example (Figure 4). We weighted the reconstructed Tjan values by the relative frequencies of observations in relation to all cores with similar temperature values. Similar temperatures were determined by binning temperature values corresponding to all cores that represented a particular time slice. This normalization still integrates an unavoidable potential bias that is related to the availability of pollen cores in the Tjan gridded space in each bin.

We reconstructed the past geographic distributions of A. alba, F. sylvatica, and P. abies at each time period using macrofossil data and pollen data. These species are ideal case studies because (1) they represent dominant arboreal species in Europe, and (2) they can be identified with a high degree of confidence to species in the fossil record. Unlike many other genera, Abies, Fagus, and Picea each has one dominant species in Europe, substantially reducing potential errors in the identification of fossil pollen taxa.

We restricted this study to the period prior to 3 ka to avoid human disturbances that strongly impacted European forests in the recent millennia (Kaplan et al., 2009). Therefore, changes observed in distributions of taxa in the past are primarily the result of climate change and ecological processes. The macrofossil data used for identifying species occurrences were obtained from published data-sets (F. sylvatica, Magri et al., 2006; A. alba, Terhürne-Berson et al., 2004 and P. abies, Latałowa and Van Der Knaap, 2006). Occurrences from pollen data were obtained from the EPD using a threshold of pollen percentage of 1% for detection of the three taxa. Since neither macroremains nor pollen of these three taxa were used for paleoclimate reconstructions, we assumed that the reconstructed geographic distributions are independent of the climate variables.

The extent to which the distribution of Tjan values obtained from fossil data reliably represents the thermal range of a species depends on sample size and geographical distribution. The coring sites are reasonably well distributed over the whole of Europe, but the number of sites where occurrences are observed decreases between 3 and 10 ka. The sample size ranges from 375 at 3 ka for P. abies, to 12 at 10 ka for A. alba (Table 1) because the geographic distributions of European trees, including the focal three species, were dramatically reduced during glacial times and expanded progressively throughout Europe during the Holocene. As a result, species thermal ranges may be less accurately represented at 10 ka, which corresponds to early post-glacial recolonization. We estimated the uncertainty related to random spatial sampling and to sample size using a Monte Carlo simulation. The modern spatial distribution of the focal species was randomly sampled to extract a Tjan distribution with the same sample sizes as available in the fossil data. Standard errors for Tjan medians were obtained after 1000 iterations and represented as error (Figure 5). The sampling uncertainty associated with Tjan medians was lower than 1°C for the three species at all time periods.

Several issues exist regarding the reconstruction of species past occurrences from fossil data. In general, fossil records allow the identification of in situ presence of a species, but tell less about whether the species was effectively absent. Undetected (thus probably small or sparse) populations may occur at a given site but be “silent” in the fossil record (e.g., trees not producing pollen; Hicks, 2006). Also, since tree plant species have different dispersal capacities, long generation time and their propagules may travel considerable distances from parent trees, it is generally accepted that there is “an inevitable time lag between the establishment of seedlings and the maturation of the trees to the stage that they are able to produce pollen” (Hicks, 2006). Reconstructions of broad-scale Holocene vegetation changes, including Fagus (Huntley et al., 1989), suggest that the bioclimatic response exhibits only a minor time lag. Although Svenning et al. (2008) argue that the expansion in Europe was delayed by several centuries, Tinner and Lotter (2006) conclude that, due to efficient animal dispersers, the lag was probably of minor importance in the overall postglacial migration process. To address both these issues, we considered that a species was present when it occurred within a broad time window of ± 250 years. In such a window, the potential time lag between the first arrival of many European tree species at a site and the establishment of a substantial population should be barely detectable (Tinner and Lotter, 2006).

The Tjan values reconstructed in the occupied grid cells yield frequency distributions of Tjan for each species at each time period, allowing us to estimate changes in the thermal range of the three species (Figure 5). Tjan distributions (Figures 5A–C) are generally not normal (Shapiro test, p < 0.05) and, thus, they are represented by their medians and the 0.25 and 0.75 quartiles. This interval corresponds to the core range of the species environmental distribution. The significance of differences between past and modern median Tjan values are derived from non-parametric Wilcoxon tests (Table 1).

Discarding the three taxa from the modern plant data-set to quantify Tjan does not affect significantly the reconstructed temperature values (Figure 2). About 83% of the values deviate with less than 1°C, while about 95% deviate less than 2°C in the negative values (Figure 2). More than 99% of the reconstructed positive values deviate by less than 1°C. In a test of temperature reconstruction in the reciprocal direction, use of the Tjan range from the Holocene to quantify modern Tjan values results in 89% of the values deviating less than 1°C and 99% of values less than 2°C when negative (Figure 2). These two dissimilarity tests show (1) that the reconstruction method is robust and, therefore, that the reconstructed Tjan values are reliable and (2) that using a Tjan range of P. abies, F. sylvatica, and A. alba, as obtained for each species between 10 and 3 ka (Figure 5), to quantify Tjan has a very minor effect on the reconstructed values. This effect impacts mostly, if not exclusively, the reconstructed negative Tjan values.

The interpolated reconstructed Tjan values show Tjan ubiquitously lower than 10°C until 7 ka. Early warming (Tjan >0°C) in Europe first took place in the SW part of Europe (Iberian Peninsula), with a gradient of more than 10°C developing toward NE Europe (Figures 3A,B). Subsequently, warmer winters developed over southern Europe, a change corresponding to the onset of “Mediterranean” climate, with dry summers and frost-free winters. Mountainous and other areas lacking pollen data may exhibit stronger interpolation deviation, potentially resulting in low Tjan values, such as those in Sicily and the Peloponnese area at 9 ka. The longitudinal gradient of the Early Holocene switched to a latitudinal gradient after 7 ka and remained so until 3 ka. The Tjan difference between the Mediterranean region and northern Scandinavia after 6 ka was as high as 15–20°C. This latitudinal distribution of Tjan across Europe seems to have persisted until today.

The overall Tjan range within all the areas recolonized during the Holocene by P. abies lies between approximately −16°C and 7°C while the modern range is between −16°C and slightly less than 6°C (Figure 5A). The Tjan at which P. abies populations are frequent today (the inter-quartile range) are about 3.5°C colder than at 5 ka, when P. abies occupied areas that tended to be relatively warm (Figure 5A). The recolonization of most of Scandinavia by P. abies took place after 4 ka (Tollefsrud et al., 2008), which explains the strong shift of the Tjan range toward lower values.

Unlike Picea, the main European Fagus populations today occupy, on average, warmer areas than at any time during the Holocene (Figure 5B). Nonetheless, half of the frequently occupied temperature range of the modern period overlaps with temperatures also occupied areas in the Holocene. We observe that during the early recolonization process (at 10 and 9 ka), Fagus populations occurred in areas with Tjan values that were generally lower than at other times during the Holocene. In contrast to P. abies and F. sylvatica, A. alba occurs in areas that are within its Holocene Tjan range (Figure 5C). Today the geographical range of the species spans a less extended temperature range than during the Holocene.

The Wilcoxon tests (Table 1) show that the medians of P. abies and F. sylvatica Tjan range differ substantially between modern values and those of the period 10 and 3 ka. This indicates shifts of the thermal range of theses two species between the Holocene and today. The thermal range for A. alba was significantly different from that of today during some periods but not in the overall Holocene range (Table 1, Figure 5C).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.