Figure 5). Comparison of fit scores summed across the landscape (overall) and for eight individual sites (see Farrell et al., unpublished for details) for scenarios where each coring point is modelled as an area of mire surrounded by an alder carr belt, and mire size, carr width and proportion of clearings in the wider landscape are all varied. Colour shows relative fit score value for each dataset (single column, 2,816 scores), where green is low and red is high. The lowest 2% of values (identifying the scenarios which produce simulated pollen assemblages most similar to the actual pollen recorded) are highlighted with blue overlay.

In the case studies presented here, each time slice was treated individually, in order that a problem in one time slice would not be propagated into others, and could more easily be detected as an anomaly. An alternative way to use the MSA is to work on timeslices sequentially, beginning each timeslice with one or more grids from the best-fit family for the previous timeslice, and modifying the vegetation until suitable reconstructions are identified. This would ensure a plausible sequence, but would need to be carefully planned and carried out to avoid missing out on equifinal solutions along the way.

As an alternative to maps, results can be presented quantitatively, as landscape metrics. Examples here showed metrics taken from the output file, based on the model parameters of the best fit scenarios (e.g., the mean and standard deviation of the probability that a pixel in a defined landscape setting could be occupied by oak-pine woodland), which are not the same as land cover metrics (due to the random nature of land cover map creation). The trends they show can be interpreted with confidence, but the absolute values will not be accurate. Values of area of landscape occupied by a given community can be extracted directly from the land cover grids, and are a more attractive metric for end users, but at present are also far more labor-intensive to obtain. Improving the extraction of information from grids is a development priority for the MSA.

# The MSA and the Palynological Workflow

Whilst the LandPolFlow package is complete and documented, the use of the MSA is in its infancy and not all package features are considered here. A suggested workflow for future studies is shown in **Figure 8** based on our experience with these case studies, and **Figure 9** shows how the methods could be built into palynological studies more generally. For these case studies we used R to speed up the manipulation of output files, and see considerable potential to use automation to increase the capacity of the method. More detailed validation of the MSA using modern pollen and vegetation datasets, currently under way, will also help understand its capacity and limitations better, and improve its deployment.

As more sets of pollen dispersal and deposition parameters are published for different taxa and geographic areas, and as increasing processor speed reduces the time and computational resource needed to run LandPolFlow, we believe that modelling approaches should become a routine part of the pollen analysts' tool kit, as zonation using a constrained clustering approach such as CONISS rather than merely visual inspection has become routine since its introduction (Grimm, 1987).

# CONCLUSION

Modelling approaches such as the MSA offer a clear route to getting more information from existing published pollen records, which represent a huge investment of time and money. This paper presents two case studies of the MSA in use at a landscape level, filling the gap between the regional land cover reconstructions offered by REVEALS (Sugita, 2007a; 100 × 100 km squares in e.g., Trondman et al., 2015) and the single-site reconstruction of distanceweighted plant abundance made possible by LOVE (Sugita,

2007b; distance of c. 500 m−2.5 km depending on site-specific Relevant Source Area of Pollen), and generating reconstructions in forms which are both transparent and useful for a wide range of end users, whether paleoecologists, environmental archaeologists, conservation biologists, neoecologists, or science communicators. The land cover hypotheses produced by the MSA can serve as the starting point for more structured, hypothesis-testing uses of paleoecology and environmental archaeology methods, including not only pollen analysis but also other techniques recording different aspects of land cover.

Our reconstructions of Neolithisation in both a heavily wooded and a largely treeless landscape show that, although the extent of land cover modification was limited and did not alter the essential character of each landscape, it was persistent. In Orkney, the timing of land cover change occurred at the same time as the start of settlement use as determined from an independent set of radiocarbon dates from archaeological sites, strongly supporting our interpretation of land cover change as a proxy for human activity in the wider landscape in the Somerset case study where archaeological evidence for Neolithic settlement is minimal. Increased levels of clearance or disturbance in both locations persisted even when archaeological evidence implied a reduction in intensity of landscape use and occupation, suggesting a significant transition in landscape character and ecology occurred with Neolithisisation. Modelling, and its ability to support targeted selection of sites for pollen analysis in specific environmental contexts and to address specific archaeological questions, clearly has a useful role to play in improving our understanding of landscape impacts and development, and in fostering better collaboration between

paleoecologists and other disciplines concerned with past landscapes.

We argue that the MSA deserves a place within the pollen analyst's standard tool kit (**Figure 9**). There is still work to do in improving the use of the approach and the communication of outputs, and as with all models, outputs are dependent on assumptions and calibration and therefore "wrong" in an absolute sense. The case studies presented here show how LandPolFlow makes environmentally informed reconstruction of land cover accessible to pollen analysts, and therefore has the capacity to improve research design, prompt more overt, and thoughtful exploration of assumptions and site taphonomy, and improve communication between subject specialists and with wider stake holders.

# AUTHOR CONTRIBUTIONS

MB and MF: worked together on designing model strategy, iterative analysis of output of each run, interpretation of findings, and led the authoring of paper. MF: prepared input GIS grids, collated pollen targets, and carried out most model runs. PM: collated datasets and modelled radiocarbon chronologies for allocation of samples to timeslices with AB. PM, AB, and AW: selected sites in light of archaeological relevance and contributed to hypothesis development and output interpretation and representation.

# FUNDING

This study was supported as part of The Times of Their Lives project (www.totl.eu), funded by the European Research Council (Advanced Investigator Grant 295412) and led by AW and AB.

# ACKNOWLEDGMENTS

First and most significantly, we wish to acknowledge the work of Richard Middleton. Now retired and busy with other projects (e.g., as Vice County Recorder of Botany for VC 61), he did not choose to be a co-author of this manuscript, but without his talent, enthusiasm, and willingness to explore unfamiliar intellectual territory this paper, and the wider MSA project, would never have come to fruition. We next acknowledge the many researchers who contributed pollen count data to the Somerset case study. We also wish to acknowledge the many members of the PolLandCal and LandClim networks, and of the PAGES Landcover 6k Working Group, for fruitful discussion of land cover reconstruction. We consider this work to be in part a contribution to the work of the PAGES Landcover 6k Working Group (http://www.pages-igbp.org/ini/ wg/landcover6k/intro), included as part of the Past Global Changes (PAGES) project, which in turn received support from the US National Science Foundation and the Swiss Academy of Sciences. We thank two reviewers for comments on an earlier version of the manuscript.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00036/full#supplementary-material

# REFERENCES

models using pollen, charcoal and non-pollen palynomorph data. Veg. Hist. Archaeobot. 19, 545–558. doi: 10.1007/s00334-010-0263-1

**Conflict of Interest Statement:** 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.

Copyright © 2018 Bunting, Farrell, Bayliss, Marshall and Whittle. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Quantitative Palynology Informing Conservation Ecology in the Bohemian/Bavarian Forests of Central Europe

Vachel A. Carter 1 \*, Richard C. Chiverrell 2 , Jennifer L. Clear 3, 4, Niina Kuosmanen4 , Alice Moravcová1 , Miroslav Svoboda4 , Helena Svobodová-Svitavská5 , Jacqueline F. N. van Leeuwen6 , Willem O. van der Knaap6 and Petr Kuneš 1

1 Department of Botany, Faculty of Science, Charles University, Prague, Czechia, 2 Department of Geography and Planning, University of Liverpool, Liverpool, United Kingdom, 3 Department of Geography and Environmental Science, Liverpool Hope University, Liverpool, United Kingdom, 4 Department of Forest Ecology, Faculty of Forestry and Wood Sciences, Czech University of Life Sciences, Prague, Czechia, 5 Institute of Botany, v.v.i., Czech Academy of Sciences, Pr ◦ uhonice, Czechia, 6 Institute of Plant Sciences and Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland

## Edited by:

Jesse L. Morris, University of Utah, United States

### Reviewed by:

Heikki Tapani Seppä, University of Helsinki, Finland Ranae Dietzel, Iowa State University, United States

> \*Correspondence: Vachel A. Carter<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Agroecology and Land Use Systems, a section of the journal Frontiers in Plant Science

> Received: 28 July 2017 Accepted: 27 December 2017 Published: 17 January 2018

### Citation:

Carter VA, Chiverrell RC, Clear JL, Kuosmanen N, Moravcová A, Svoboda M, Svobodová-Svitavská H, van Leeuwen JFN, van der Knaap WO and Kuneš P (2018) Quantitative Palynology Informing Conservation Ecology in the Bohemian/Bavarian Forests of Central Europe. Front. Plant Sci. 8:2268. doi: 10.3389/fpls.2017.02268 In 1927, the first pollen diagram was published from the Bohemian/Bavarian Forest region of Central Europe, providing one of the first qualitative views of the long-term vegetation development in the region. Since then significant methodological advances in quantitative approaches such as pollen influx and pollen-based vegetation models (e.g., Landscape Reconstruction Algorithm, LRA) have contributed to enhance our understanding of temporal and spatial ecology. These types of quantitative reconstructions are fundamental for conservation and restoration ecology because they provide longterm perspectives on ecosystem functioning. In the Bohemian/Bavarian Forests, forest managers have a goal to restore the original forest composition at mid-elevation forests, yet they rely on natural potential vegetation maps that do not take into account long-term vegetation dynamics. Here we reconstruct the Holocene history of forest composition and discuss the implications the LRA has for regional forest management and conservation. Two newly analyzed pollen records from Prášilské jezero and Rachelsee were compared to 10 regional peat bogs/mires and two other regional lakes to reconstruct total land-cover abundance at both the regional- and local-scales. The results demonstrate that spruce has been the dominant canopy cover across the region for the past 9,000 years at both high- (>900 m) and mid-elevations (>700–900 m). At the regional-scale inferred from lake records, spruce has comprised an average of ∼50% of the total forest canopy; whereas at the more local-scale at mid-elevations, spruce formed ∼59%. Beech established ∼6,000 cal. years BP while fir established later around 5,500 cal. years BP. Beech and fir growing at mid-elevations reached a maximum land-cover abundance of 24% and 13% roughly 1,000 years ago. Over the past 500 years spruce has comprised ∼47% land-cover, while beech and fir comprised ∼8% and <5% at mid-elevations. This approach argues for the "natural" development of spruce and fir locally in zones where the paleoecology indicates the

**84**

persistence of these species for millennia. Contrasting local and regional reconstructions of forest canopy cover points to a patchwork mosaic with local variability in the dominant taxa. Incorporation of paleoecological data in dialogues about biodiversity and ecosystem management is an approach that has wider utility.

Keywords: biodiversity, Holocene, land-cover, palynology, pollen, REVEALS

# INTRODUCTION

Quantitative reconstructions are fundamental for providing long-term perspectives of ecosystem processes because they can be used to develop baselines for conservation and restoration ecology (National Research Council, 2005; Froyd and Willis, 2008). Paleoecological records have utility in conservation strategies related to biodiversity maintenance, ecosystem naturalness, conservation evaluation, habitat alteration, changing disturbance regimes, and species invasions (e.g., Birks, 1996; Jackson, 1997; Landres et al., 1999; Swetnam et al., 1999; Foster et al., 2003; Gillson and Willis, 2004). Unfortunately, paleoecological research is still largely ignored by conservation biologists and conservationists (Willis and Birks, 2006; Birks, 2012). In 1916, Lennart von Post published the first pollen diagram, providing one of the first qualitative reconstructions of vegetation change that extended over millennial timescales, as well as benchmarking the foundation of palynology (Manten, 1967). Müller (1927) published the first pollen diagram from the Bohemian Forest more than a decade later. Over the past 100 years, palynology has developed from a qualitative tool to a more quantitative analysis of vegetation dynamics that is increasingly well constrained in time (Davis, 2000).

Pollen diagrams produced from mires, peat bogs, and lake sediment profiles, with the data expressed as percentages against time or depth, have been the main way to present data (e.g., Stalling, 1987; Knipping, 1989; Svobodová et al., 2001, 2002; Jankovská, 2006). While percentage pollen diagrams identify changes in vegetation through time, changes in one taxon can affect the percentage proportions of all other taxa. This effect, termed data-closure (Birks and Birks, 1980), can produce disconnections between the trends of pollen and the actual vegetation land-cover. Several quantitative methodologies have been developed to address these issues. Data closure artifacts in pollen percentages can be identified by comparison with pollen influxes (i.e., pollen accumulation rates; grains cm−2 yr−1 ), but this requires robust age-to-depth relationships for the sediments. Even then, both percentage and influx data fail to take account of differences in the ecological and environmental factors that influence pollen production and dispersal, nor do they define a spatial scale in reconstructions of vegetation cover (Loidi et al., 2010; Loidi and Fernández-González, 2012). Additionally, because pollen influxes are dependent upon age-to-depth relationships, changes in sedimentation rates can influence pollen concentrations and the subsequent accuracy of influx calculations. Therefore, Sugita (2007a,b) developed the Landscape Reconstruction Algorithm (LRA) which incorporates pollen productivity and pollen dispersal capacity, as well as factors that influence pollen dispersal, such as the size and type of the sedimentation basin. Thus, the LRA provides a more quantitative approach of estimating past vegetation abundance in a defined space. The LRA has two steps: (1) the REgional Vegetation Estimates from Large Sites (REVEALS) model (Sugita, 2007a) estimates pollen-derived regional vegetation cover from large sites (>100 ha) or alternatively many small sites across an area of 106 km2 (e.g., Abraham et al., 2014, 2016); and (2) the LOcal Vegetation Estimates (LOVE) model (Sugita, 2007b) uses the regional estimates from REVEALS to estimate pollen-derived local vegetation composition in smaller areas (5–104 ha).

Long-term perspectives on forest dynamics are beneficial for conservation and restoration, and here we demonstrate the utility of quantitative palynological data in the management of Bohemian/Bavarian Forests, including the use of pollen influx data and pollen-landscape reconstruction models (e.g., REVEALS: Sugita, 2007a). Abraham et al. (2016) used the REVEALS model to estimate vegetation cover by integrating pollen sequences from the Šumava region and estimate that spruce has comprised over 50% of the regional forest canopy for the past 7,000 years. Here, we synthesize two newly analyzed, well-dated pollen records from Prášilské jezero (Šumava, Bohemian Forest; Czech Republic) and Rachelsee (Bavarian Forest; Germany) and several previously analyzed pollen records from peat bogs/mires and lakes using multiple quantitative methods to determine the ecological development and vegetation change during the Holocene. The objectives are: (1) to critically evaluate the multiple quantitative approaches used in palynology to reconstruct the history of forest composition; and (2) to explore the utility of quantitative land-cover reconstructions in regional conservation and restoration.

# BACKGROUND: THE BOHEMIAN/BAVARIAN FOREST REGION OF CENTRAL EUROPE

The Bohemian/Bavarian Forest region of central Europe includes two national parks located along the border of the Czech Republic (Šumava National Park, Bohemian Forest), and Germany (Bavarian Forest National Park1 , Bavarian Forest; **Figure 1**). Together, these two parks with their surrounding area comprise one of the largest forested landscapes in central Europe, which is home to many endangered flora (see Kˇrenová and Hruška, 2012) and fauna, and provides a close-to-nature ecosystem of ecological value (Meyer et al., 2009). The Bavarian Forest National Park

1Bavarian Forest National Park Authority: http://www.nationalpark-bayerischerwald.de/english/nationalpark/management/index.htm.

was established in 1970 with part of Šumava declared a UNESCO Biosphere Reserve in 1990, and Šumava National Park established in 1991.

Subject to long-term human influence, management practices in the Bohemian/Bavarian Forests have ranged from natural development to intensive forest management which has modified both the structure and composition of the forest. Currently, the two parks are divided into wilderness (i.e., not managed), and non-wilderness (i.e., managed). Within the non-wilderness areas, management focuses on sanitary logging of spruce trees killed by previous disturbance events, specifically bark beetle outbreaks. Within Šumava National Park, forest management restoration practices involve reducing spruce populations from the current 84% to a target "natural" representation of 30– 40%, and increase beech (ca. 6%) and fir (ca. 1%) up to 35% (Šumava NP Authority). Forest management policies in the Czech Republic direct that the lands between 400 and 900 m altitude should be dominated by beech forests with spruce as a secondary component (Vacek and Mayová, 2000; Pr ◦ uša, 2001). This management plan is likely based on the Geobotanical Map (Mikyška et al., 1968–1972) and the Map of Potential Natural Vegetation (Neuhäuslová et al., 1998) with beech considered the natural dominant species at mid-elevations in central Europe (Ellenberg and Leuschner, 1996). Within the Bavarian Forest National Park, land managers are concerned with the rapid decline of fir populations in valley bottoms and have focused their attention to factors that could facilitate the natural development of this species (Heurich and Englmaier, 2010), rather than focusing on the removal of spruce which has doubled its range in the Bavarian Forest since the nineteenth century. This strategy is in response to previous paleoecological work documenting that the Bavarian Forests consisted of ∼32% fir ∼3,000 years ago (Stalling, 1987).

Palynological research has a long history in the Bohemian/Bavarian region of central Europe with Müller (1927) publishing one of the first qualitative glimpses of longterm vegetation development in pollen diagrams from several moors in the region (**Figure 2**). Müller (1927) documented the long-term presence of spruce in the Bohemian/Bavarian Forests, but was unable to discuss these changes relative to time. With widespread application of radiocarbon dating, pollen data from peat sequences (i.e., small peaty wetlands within forests) in Šumava National Park confirm that spruce has dominated the forest canopy cover at elevations >700 m a.s.l. for the past 7,000–8,000 years (Svobodová et al., 2001, 2002). Notwithstanding the better chronological control, the reconstructions presented by Svobodová et al. (2001, 2002) are still limited in that; (1), pollen percentage diagrams offer

a qualitative reconstruction of vegetation change through time; and (2) peat bogs/mires and wetlands capture changes in vegetation at various spatial scales [i.e., a forest-stand scale (forest mire) to the extra-local scale (open, wetland mire)]. Reconstructing extra-local to regional-scale vegetation changes is possible by analyzing sediments from lakes as well as open mires, which typically have a much larger pollen catchment area than forest mires. However, pollen-based land-cover models (i.e., REVEALS) are the best approach for reconstructing vegetation abundances at the regional-scale. Currently, the only quantitative land-cover reconstruction of vegetation change in the region using the REVEALS model (Sugita, 2007a) and 14 sediment sequences suggests that spruce forests have comprised ∼50% of the forest canopy throughout the past ∼7,000 years (Abraham et al., 2016). However, to develop more detailed descriptions of regional-scale vegetation dynamics, further quantitative data are needed. This study provides a more detailed description of land-cover change by analyzing lakes and peat bogs/mires that vary spatially across mid- (700–900 m a.s.l.) and high-elevations (>900 m a.s.l.) from the Bohemian/Bavarian Forests (**Table 1**), and addresses whether beech has been the dominant forest canopy type at mid-elevations over millennial time-scales.

# METHODS

# Core Retrieval, Sediment Limnology, and Radiocarbon Dating

A 2.18 m sediment profile was collected in August 2015 from the deepest (14.8 m) part of Prášilské jezero (49◦ 04′ N, 13◦ 24′ E, 1,079 m a.s.l.) and is comprised of two parallel and overlapping cores (PRA 15-2-1 and PRA 15-2-2). The sediment profile was sampled from a floating platform using a handpercussion Russian-style corer (1.5 × 0.075 m). The sediment– water interface was collected using a 0.1 m diameter gravity corer (PRA15-2GC) (Boyle, 1995). At Rachelsee (48◦ 58′ 29" N, 13◦ 24′ 7" E, 1,071 m a.s.l.), 11.8 m of sediment was collected in August 2012 at 12.5 m water depth using a Livingstone piston corer (Wright, 1967) (RAA1-4). Sediment was not recovered between depths 0–57 cm, therefore, the mud-water interface depth begins at 57 cm. Sediment age-depth relationships were established using 14C radiocarbon dating at Prášilské jezero (n = 10) and Rachelsee (n = 7), with an additional 210Pb series (Appleby, 1978) at Prášilské jezero (**Table 1**). For all sites, the geochronological data including the sediment surface were compiled within the Bayesian routine "BACON" (Blaauw and Christen, 2011). This analysis partitioned both cores into 0.05 m thick sections and estimated the accumulation rates for each segment using a Markov Chain Monte Carlo (MCMC) approach (Christen and Pérez, 2009). The analyses were constrained by a prior model of sediment accumulation rate (a gamma distribution with mean 50-year cm−1 and shape 1.5) and its variability (memory, a beta distribution with mean 0.5 and shape 16) for both sites. All 14C ages were calibrated and modeled in "BACON" using the IntCal13 curve (Reimer et al., 2013), with a Student-t distribution to account scatter in the 14C measurements and to allow for statistical outliers (Blaauw and Christen, 2011). The weighted mean modeled ages against depth were smoothed using a 21 point moving average (**Figure 3**).

TABLE 1 | Summary of age-depth relationships for Prášilské jezero (Core IDS, PRA-15GC2, PRA 15-2-1, and PRA 15-2-2), Czech Republic and Rachelsee (Core ID, RAA1-4), Germany.

# Pollen Analysis

At Prášilské jezero, pollen analysis was conducted at 1–2 cm resolution throughout the core by V. A. Carter. For each sample, 0.5 cm−3 was processed using standard pollen procedures (Faegri et al., 1989). At Rachelsee, pollen analysis was conducted at 1 cm resolution for 72–124 cm, 2 cm for 126–184 cm, and 4 cm for 188–540 cm by J. F. N. van Leeuwen. The sediment between depths 57–71 cm was too watery and could not be subsampled for pollen. Here, we present data from depths 72–368 cm. At Rachelsee, a minimum of 300 tree pollen grains, and at Prášilské jezero, a minimum of 500 pollen grains were counted in each sample. A Lycopodium tablet with a known amount of spores were added to each sample as an exotic tracer in order to calculate pollen concentration and influx rates (Stockmarr, 1972). Total pollen counts were converted into pollen percentages and plotted using Tilia software (Grimm, 1987) based on the abundance of each pollen type relative to the sum of all terrestrial identified pollen. The pollen profiles were divided into pollen assemblage zones based on optimal splitting using the broken-stick model (Bennett, 1996). Pollen counts were converted to concentrations (grains cm−3 ) in Tilia using the counts of Lycopodium tracers and sample sediment volumes (Grimm, 1987). The smoothed Bayesian age-depth model for both sites was used to calculate and plot pollen influx rates from the concentration data using Tilia software.

# REVEALS

The REVEALS model (Sugita, 2007a) was applied to estimate pollen-derived regional vegetation cover during the Holocene. Previous studies from the region (e.g., Abraham et al., 2014) mainly used records from peat bogs to calculate a REVEALS output, with the model showing the robustness of peat sequences in portraying present-day vegetation in a realistic manner. Here, the REVEALS model was used to estimate land-cover abundance from both lakes and non-lakes (i.e., peat bogs and mires) using new (this paper) and previously published and unpublished pollen data obtained from the Czech Quaternary Pollen Database (Kuneš et al., 2009; **Table 2**). Pollen counts were aggregated into 500-year intervals for the entire Holocene, and then all available sites in each 500-year time window were used to estimate mean vegetation abundances within a 60 km radius. The following parameters were selected: 4 m s−1 wind speed, different dispersal-deposition models for bogs and mires (Prentice, 1985) and lakes (Sugita, 1994), and taxon-specific relative pollen productivities for 28 selected pollen-equivalent taxa (see Supplementary Table 1; Abraham and Kozáková, 2012; Mazier et al., 2012; Abraham et al., 2014). For model calculations, we adjusted the basin sizes of non-lakes to 4.5 ha. Larger basin sizes produce unrealistic estimates for some taxa as some trees may still grow within the basin, violating the model's assumptions (Abraham et al., 2014). Therefore, assuming smaller non-lake basins was more appropriate for our calculations. Openness was calculated by the summation of all herbs included in the model. All calculations were conducted in R (R Core Team, 2016) using functions to calculate mean vegetation abundances and their error estimates based on bootstrap methods (https://github.com/ petrkunes/LRA).

# RESULTS

# Regional Vegetation Development

Müller (1927) originally presented data for 12 tree taxa (**Figure 2**), however, the discussion of this research will focus on the three dominant canopy species growing in the Bohemian/Bavarian Forests at modern times; spruce (Picea abies), beech (Fagus sylvatica), and fir (Abies alba). Based on

TABLE 2 | Regional sites uded to create the REVEALS pollen-based vegetation model.

Sites with a "\*" symbol indicate new, but unpublished data that was used to refine the REVEALS model (see Supplementary Material).

the broken-stick model (Bennett, 1996), the pollen percentages profiles from Prášilské jezero and Rachelsee were divided into seven and eight pollen assemblage zones, but for simplicity, the pollen reconstructions were grouped into three main periods; (1) Early- to Mid-Holocene, (2) Mid- to Late-Holocene, and (3) Last-Millennium.

# Early- to Mid-Holocene (12,000–∼6,800 cal. Years BP)

Spruce first arrived at the high-elevation lakes Prášilské jezero and Rachelsee between 10,500 and 10,000 years ago, but was present at mid-elevation peat bogs/mires by ∼11,500 cal. years BP (**Figure 4**). Spruce rises to dominance between 10,000–8,500 cal. years BP across all elevations (**Figure 4**) replacing pine (Pinus) and hazel (Corylus) as the dominant canopy species (see Supplementary Figures 3–5). Thereafter, spruce percentages fluctuated between 35% and >20% pollen at both Prášilské jezero and Rachelsee, and around 30% at mid-elevation peat bogs/mires, respectively. Pollen influx data for spruce conversely differ between the two sites with a peak of 4,000 grains cm−2 yr−1 at Prášilské jezero, whereas the increases at Rachelsee are more subdued at <500 grains cm−2 yr−1 (**Figure 5**). Quantitative landcover reconstructions of spruce percentages calculated using REVEALS (Sugita, 2007a) rise around 10,000 cal. years BP from ∼5% to >60% calculated on a local (using peat bogs/mire records) basis (**Figure 6**). On a regional basis (using lake records), spruce comprised >50% total land-cover. The REVEALS landcover reconstructions suggest that spruce forms close to double the land-cover percentages compared to pollen percentages. Beech begins to expand at both sites from 7,000 cal. years BP onwards, but is at low percentages in the Early- to Mid-Holocene. The dynamics of fir are restricted to after 5,000 cal. years BP. Total herb percentages decrease from ∼20% and oscillate around 10% at both Prášilské jezero and Rachelsee. Pollen influx data for total herbs decrease ∼10,500 cal. years BP from values >2,000 grains cm−2 yr−1 . The REVEALS reconstruction also suggests that the percent of landscape openness is close to double that of pollen percentage data at both the regional and local-scale (**Figure 6**).

# Mid- to Late-Holocene (∼6,800–1,000 cal. Years BP)

Between 6,800 and 6,000 cal. years BP beech increases gradually in pollen percentages at both Prášilské jezero and Rachelsee, and then oscillates around 25% (**Figure 4**). Beech percentages from mid-elevation peat bogs/mires begin to increase ∼500 years earlier than in the high-elevation lakes. At the localscale (peat bogs/mires), beech percentages gradually increase throughout the Mid- to Late-Holocene peaking at ∼25% roughly 1,000 cal. years BP. Fir expansion commences 5,500–5,000 cal. years BP to around 20% at both Prášilské jezero and Rachelsee. These increases in pollen percentages are associated with no substantial reduction in spruce pollen percentages at Prášilské jezero and Rachelsee, but spruce percentages decrease from ∼40% to ∼20% <900 m a.s.l. (**Figure 4**). Viewed from a perspective of pollen influx, the down core pattern of increase and stability in species percentage data are reflected in the influx data (**Figure 5**). Though as with spruce, influx rates for beech and fir are much lower at Rachelsee. The REVEALS land-cover reconstructions for spruce and fir are lower at the regional scale (i.e., lake sites) relative to the more local-scale (i.e., peat bog/mire sequences), whereas beech is an equally abundant component in both reconstructions comprising ∼25% land-cover (**Figure 6**). The land-cover contribution of spruce remains double (>40%) the values indicated by pollen percentage data (>15%), but these differences are not present for beech and fir. Land-cover percentages of spruce remained relatively stable at the lake sites however, spruce decreases from >60% to ∼40% at mid-elevation peat bog/mire sequences during this period. The REVEALS reconstruction also suggests that the percent of openness is <5% at both the regional and local-scales (**Figure 6**).

well as the average percentages and maximum/minimum percentages of all peat bog and mire sites grouped into 500 year bins.

Bohemian/Bavarian Forests of central Europe.

# Last-Millennium (1,000 cal. Years BP-Present)

At both lakes, spruce percentages remain unchanged initially and then increase ∼5% in the last 300 years (**Figure 4**). Spruce percentages remain unchanged at mid-elevation peat bogs/mires. Both beech and fir pollen percentages decline at both lakes and mires sequences through the last 1,000 years (**Figure 4**). Openness percentages increased dramatically over the past 1,000 years to ∼12% (**Figure 4**). Pollen influx rates for all three taxa increase during the last 1,000 years at both lakes, and sharply in the last 300 years which likely reflects focusing of pollen toward the lake center (i.e., in-lake lateral redistribution of materials toward the deeper waters) rather than any change in species abundance (**Figure 5**). The REVEALS reconstructions broadly confirm the indications seen in the percentage pollen data, with declines in beech and fir occupied by open-ground indicators rather than any expansion in spruce. This is also shown in the pollen influx data, showing a dramatic rise in the influx of open-ground taxa (**Figure 5**).

# DISCUSSION

# Developments in Understanding Dominant Ecology of the Bohemian/Bavarian Forests

Over the past 100 years, methodological improvements in the field of paleoecology such as radiocarbon dating and pollen-based quantitative land-cover modeling (i.e., REVEALS; Sugita, 2007a) have led to more realistic reconstructions of past vegetation change across temporal and spatial scales. Applying these approaches to the palaeoenvironmental records from Prášilské jezero and Rachelsee illustrates the long-term presence and dominance of spruce in the Bohemian/Bavarian Forests. First appearing in the pollen record ∼10,000 cal. years BP, the local presence of spruce in the region is supported by Picea stomata at Prášilské jezero and Rachelsee, which agrees with previous paleoecological research in the region (Müller, 1927; Svobodová et al., 2001, 2002; Jankovská, 2006; Abraham et al., 2016). At the regional scale (i.e., lake sites), the REVEALS model suggests that spruce has comprised an average of ∼50% of the total forest canopy, whereas at the more local scale (i.e., peat bogs/mires) spruce has comprised ∼59%. Thus, spruce has been a dominant component of the forest canopy across a range of elevations >700 m a.s.l. in the Bohemian/Bavarian Forests.

Modern quantitative approaches also illustrate a unique history of forest composition throughout the region with beech and fir contributing as secondary canopy species relative to spruce at mid-elevations found between 700 and 900 m a.s.l. Beech established around 7,000 cal yr BP and became an abundant canopy species around 6,000 cal. years BP, with the REVEALS model indicating that it contributed ∼20% of the total land-cover at the regional scale. However, at the more local-scale between 700 and 900 m a.s.l. beech comprised just 16% of the total land-cover during the Mid-Holocene. After the establishment of beech, Müller (1927) found an increase in fir percentages. Svobodová et al. (2001, 2002) suggested that fir had spread to elevations >700 m a.s.l. across the entire region between 6,300 and 3,400 cal. years BP. However, these records were constrained by few radiocarbon dates. At Prášilské jezero and Rachelsee, fir first appeared between 7,000–6,300 cal. years BP but did not become an important component of the canopy until 5,000–4,500 cal. years BP. The REVEALS model output indicates that for the last 4,000 years, fir has comprised an average of ∼7% of the total land-cover at the regional-scale, but slightly higher (10%) at the more local-scale between 700 and 900 m a.s.l. As beech and fir expanded around 4,000 cal. years BP, spruce percentages declined to some extent but still comprised ∼40% of the regional total land-cover, and >60% of the local-scale (i.e., mid-elevation) forest composition in the Bohemian/Bavarian Forests (**Figure 6**). This further demonstrates that beech has been a secondary canopy species relative to spruce at mid-elevations, specifically between 700 and 900 m a.s.l. Total land-cover of beech and fir populations found at mid-elevations began to decrease from their maximum values of 24% and 13% total landcover ∼1,000 years ago at both the regional and more localscale (**Figure 6**). The decrease in beech and fir forests are likely the result of increasing anthropogenic landscape modifications. However, climate could also explain the decline in beech and fir as temperatures were likely unfavorable for these species during the Little Ice Age (Grove, 2001; Brázdil et al., 2017).

Declines in beech and fir characterize the last millennium with an expansion in open ground taxa in both the percentage pollen data and the quantitative REVEALS land-cover reconstructions, though spruce is slightly more abundant in terms of total land-cover (**Figure 6**). Svobodová et al. (2002) also indicate that beech-spruce-fir forests have decreased over the past three centuries, and are being replaced by spruce and Pinus rotundata in forest bogs at lower elevations. The regional reconstructions presented here point to a change in forest composition toward more open ground within the forest, and slightly increased dominance by spruce. Over the past 500 years spruce has comprised ∼47% land-cover, while beech and fir has comprised ∼8% and <5% at mid-elevations. The pollen influx values show a radically different story with the records from Prášilské jezero and Rachelsee illustrating an increase in all pollen taxa through the last 500 years. This increase in pollen influx likely reflects increased in-lake focusing of pollen to the lake center assisted by wind currents from a more open landscape (see **Figure 5** for percent openness), rather than increased pollen rain from spruce plantations found at low-elevations. Openings within highelevation forests are generally attributed to human clearance, fire, windthrow, or bark beetle outbreaks (Ips typographus) which are key disturbance drivers influencing Norway spruce forest dynamics (Holeksa and Cybulski, 2001; Fischer et al., 2002; Holeksa et al., 2007; Svoboda et al., 2012, 2014; Cada et al., ˇ 2016; Janda et al., 2017; Kulakowski, 2017). During the 1990s, an I. typographus bark beetle outbreak occurred in the region creating gaps (i.e., openings) in the canopy. Severe windstorms occurred in the region in 1986, 1999, 2007, and 2008, which also caused widespread forest damage (Svoboda et al., 2010). Yet, Zeppenfeld et al. (2015) found that spruce seedlings and young spruce trees already occur in areas affected by these largescale disturbances, encouraging rapid recovery of spruce and not limiting flux of pollen to the lakes in the region. Therefore, increase in pollen fluxes especially at Prášilské jezero is likely a result of natural processes and not a result of spruce plantations found at lower elevations.

While the REVEALS model demonstrates that spruce has been the dominant canopy cover at both high- and mid-elevations, the model only provides an average percentage of land-cover for the entire Bohemian/Bavarian Forest. Local site conditions such as slope and aspect affect forest composition at the more local-scale (i.e., stand-scale). For example, spruce generally dominates on north-facing mesic sites with poorly-drained soils, whereas beech typically dominates on south-facing xeric sites where soils drain well. These local conditions can lead to stands of monospecific spruce or beech at mid-elevations, specifically in the Bavarian Forest National Park. However, the REVEALS model used in this study was used to inform managers of the natural average forest composition for elevations >700 m a.sl. for the area through time.

# Applications to Conservation and Restoration in the Bohemian/Bavarian Forests

The overall management goal of the parks is to restore the original forest composition (IUCN, 1994). Within the Šumava National Park forest management has likely based their target "natural" forest composition from Geobotanical Maps (Mikyška et al., 1968–1972), and the Map of Potential Natural Vegetation (Neuhäuslová et al., 1998), which suggest that beech dominate over spruce at mid-elevations in central Europe (Ellenberg and Leuschner, 1996). More recently, Vacek and Mayová, 2000, suggest that beech is dominant over spruce at mid-elevations between 650 and 900 m. Together, these previous maps and research have influenced the management strategy of removing spruce from mid-elevation forests within Šumava National Park so that spruce will comprise 30–40% instead of its current representation of "84%" (Šumava National Park Authority2 ). However, these maps and previous research do not consider paleoecological studies which provide long-term dynamics of vegetation abundance through time. Paleoecological studies have utility in conservation and restoration in that they can provide baseline targets (i.e., more accurate representations of forest composition), as well as provide the natural range of variability within vegetation composition through time. While the REVEALS model shows that beech and fir forests reached a combined maximum of 37% at mid-elevations (between 700 and 900m a.s.l.) ∼2,000 years ago, spruce forests still comprised 45% of the total forest composition at these elevations (**Figure 6**). A regional paleoecological perspective of the past 500 years shows that the forest composition at mid-elevations has been a mixed-forest cover of spruce (∼50%), beech (∼22%), and fir (∼10%), which counters the findings of Vacek and Mayová (2000) Our findings support previous investigations that highlight that mixed-forests consisting of spruce-beech-fir typically occur between 650 and 1,150 m in the Šumava National Park (Röder et al., 2010), as well as showing the long-term dominance of coniferous forests in the Bohemian-Moravian

2 Šumava National Park Authority. http://www.npsumava.cz/en/3291/sekce/thespecies-composition-of-the-forests/.

highlands over the past 9,000 years (Szabó et al., 2017). The long-term presence of spruce throughout the Holocene illustrates that forest management practices involving the formation of spruce plantations at mid-elevations have not suppressed beech (Szabó et al., 2017). During the late Holocene the results from the regional-scale sites show some suppression of beech (∼15% land-cover) that is not related to spruce dominance, but rather to the extent of open ground. This supports previous research from Šumava National Park with beech rarely reaching >20% from mid-elevation sites (Hrubý et al., 2014).

Within the Bavarian Forest National Park, soil conditions and south-facing slopes create more favorable conditions for beech and fir forests. While spruce has doubled its range since the nineteenth century in the region, foresters have been more concerned with the recent decline of fir. At the time of establishment of the Bavarian Forest National Park, fir comprised only 3.2% of the total forest composition (Heurich and Englmaier, 2010). This is supported by the regionalscale REVEALS reconstruction which shows the decline of fir over the past 500 years to 3.8%. However, at the more local-scale (i.e., peat bogs/mires) at mid-elevations fir has declined to just 2.9% (**Figure 6**). Stalling (1987) found that fir pollen frequencies reached a peak of ∼32% fir roughly 3,000 years ago on the basis of 15 sites with a mean elevation of 837 m a.s.l. in the Bavarian Forest. However, quantitative landcover reconstructions presented here indicate that fir reached a maximum regional land-cover of just 13% between 3,500 and 1,000 cal. years BP, and a maximum of ∼20% at the more local-scale at mid-elevations (**Figure 6**). Methodological improvements such as quantitative land-cover modeling offer a more robust reconstruction of vegetation abundance from the Bohemian/Bavarian Forests of central Europe. Thus, the results presented here provide a more accurate and updated description of fir abundances than the previous work by Stalling (1987). However, the REVEALS model is limited in that all of the mid-elevation sites used in the model are located within the park boundaries of Šumava National Park. Mid-elevation paleoecological reconstructions do exist from the Bavarian Forest National Park, yet, the sites were not incorporated into the model because they were either not available publicly or they lacked radiocarbon dating. Therefore, the results interpreted here may be more reflective of vegetation abundances from midelevation in Šumava National Park. Regardless, the REVEALS model reconstructed the regional decline of fir, and it is therefore a concern that should be addressed by foresters in the Bavarian Forests. Heurich and Englmaier (2010) suggest that the trend of declining fir populations relative to increasing spruce populations in the Bavarian Forests is likely not solely related to human activity but also to natural causes such as the Little Ice Age and bark beetle outbreaks. Because spruce is a cold adapted species, conditions during the Little Ice Age would have favored spruce over fir. Additionally, downed logs from windthrow and bark beetle outbreaks are the preferred natural substrate that gives spruce seedlings a competitive advantage over other tree species (Svoboda et al., 2010). While several studies suggest that fir is sensitive to human activities and anthropogenically caused air pollution (Tinner et al., 1999; van der Knaap et al., 2004; Tinner and Lotter, 2005; Feurdean and Willis, 2008), human impacts in the form of pasturing, selective logging and litter raking have been shown to be beneficial for fir growth (Nožicka, ˇ 1957; Málek, 1981; Kozáková et al., 2011). Tinner et al. (2013) suggest that the range of fir has the potential to increase with climate change as long as precipitation does not decrease below 700–800 mm/yr, or anthropogenic disturbances such as fire and grazing do not become excessive.

One of the benefits of the Šumava and Bavarian Forest National Parks is that it serves to protect montane spruce forests (Svoboda et al., 2010). Spruce is likely to be the species most impacted by climate change, and so is under threat in a warmer world (Bolte et al., 2009). Norway spruce is considered a cold-adapted species growing on mesic sites. Therefore, expected changes in temperature and precipitation in the future specifically the increased likelihood of more frequent and intense summer droughts across central and southern Europe (Gao and Giorgi, 2008; Feyen and Dankers, 2009), will likely result in the loss of available habitat for spruce (Hanewinkel et al., 2012). In the Bohemian/Bavarian Forests, Norway spruce has already occupied the highest possible altitudes and ostensibly the species has nowhere to go with climate change (Spathef et al., 2014). Conversely, beech and fir forests, the main competitors of spruce in the region, could continue to impinge on and expand in the lower altitude habitat of spruce. However, areas where beech predominates likely will also be threatened with increasing temperatures (Cheddadi et al., 2016). Migration of the optimum climate windows for these dominant forest taxa is a concern for forest managers, yet our paleoecological data show these plant communities have comprised a mixed-forest for >5,000 years in the region surviving previous climatic fluctuations during the Holocene. Mixed-forests with multiple tree species provide higher levels of ecosystem services (Gamfeldt et al., 2013) and should therefore be incorporated into forest management in the Bohemian Forests. We suggest that forest managers in the Šumava National Park consider utilizing information gained by paleoecological reconstructions that use pollen-based land-cover models that more accurately represent the natural vegetation abundances found at mid-elevations.

# CONCLUSIONS

Methodological advances in the field of palynology and paleoecology, including quantitative land-cover models and robust chronologies have allowed for more accurate reconstructions of vegetation dynamics through time. Using these methodologies the results of this study demonstrate that spruce has been the dominant canopy cover in the Bohemian/Bavarian Forests for the past 9,000 years across a range of elevations including mid-elevations (between 700 and 900 m a.s.l.) where beech and fir forests were previously thought to dominant. When beech and fir forests peaked around 2,000 years ago, together they comprised a total land-cover of 37%. However, spruce contributed 40% of the total land-cover, further documenting the dominance of this species in mixed beech/fir forests. Over the past 500 years, spruce has comprised ∼47% land-cover, while beech and fir comprised ∼8% and <5% at mid-elevations.

Paleoecological results have considerable utility in forest management and ecosystem conservation (e.g., Birks, 1996; Jackson, 1997; Landres et al., 1999; Swetnam et al., 1999; Foster et al., 2003; Gillson and Willis, 2004). For example, forest management in the Šumava National Park base current management strategies on outdated information that states that beech and fir forests are the natural, dominant vegetation at mid-elevations. However, the results presented here demonstrate that mid-elevations forests have been a mixed-forest with spruce dominating the forest canopy. The current management plan is to remove spruce to a natural representation of 30–40%, while increasing beech and fir up to a total of 35% (Šumava National Park Authority2 ). The longer perspective afforded by paleoecology suggests that the Šumava National Park should aim for 45% spruce as that is within the natural range of total land-cover for the species, but ecosystem managers should take into account the potential for climate change to significantly impact the range of spruce within the park and take steps toward preserving the species. As for the Bavarian Forest National Park, the REVEALS model was able to capture a decline of fir cover to levels previously recorded (see Heurich and Englmaier, 2010), yet the results may be more indicative of vegetation abundances within Šumava National Park. Additional paleoecological sites with robust chronologies are needed from the region, specifically from between 400 and 700 m a.s.l. in order to compare local-scale vegetation dynamics through time.

# REFERENCES

# AUTHOR CONTRIBUTIONS

PK and JC: obtained funding to support the project; PK, JC, RC, and AM: collected the data; PK, JvL, RC, and VC: analyzed the data; VC and PK: conceived ideas and led the manuscript writing; VC, RC, JC, NK, AM, MS, HS-S, JvL, and WvdK: edited the manuscript during all phases.

# FUNDING

This research was supported by the Czech Science Foundation EUROPIA project no. 16-06915S, and PEDECO project no. 16- 23183Y.

# ACKNOWLEDGMENTS

We would like to thank Vojtech Abraham for his suggestions ˇ on how to improve the REVEALS model, and Daniel Vondrák and Günther Kletetschka for their help with radiocarbon dating. We also thank Vlasta Jankovská and Maurice Reille for providing their published data through the Czech Quaternary Pollen Database.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2017. 02268/full#supplementary-material

responses of silver fir forests to global warming. Ecol. Monogr. 83, 419–439. doi: 10.1890/12-2231.1

**Conflict of Interest Statement:** 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.

Copyright © 2018 Carter, Chiverrell, Clear, Kuosmanen, Moravcová, Svoboda, Svobodová-Svitavská, van Leeuwen, van der Knaap and Kuneš. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Temperature Range Shifts for Three European Tree Species over the Last 10,000 Years

Rachid Cheddadi 1 \*, Miguel B. Araújo2, 3, 4, Luigi Maiorano5 , Mary Edwards 6, 7 , Antoine Guisan8, 9, Matthieu Carré1 , Manuel Chevalier 1 and Peter B. Pearman10, 11

1 Centre Nationnal de la Recherche Scientifique, Institut des Sciences de l'Evolution, University Montpellier II, Montpellier, France, 2 Departamento de Biogeografía y Cambio Global, Museo Nacional de Ciencias Naturales, CSIC, Madrid, Spain, 3 Center for Macroecology, Evolution and Climate, Natural History Museum of Denmark, University of Copenhagen, Copenhagen, Denmark, 4 InBIO-CIBIO, University of Évora, Évora, Portugal, 5 Dipartimento di Biologia e Biotecnologie "Charles Darwin," Università di Roma "La Sapienza," Roma, Italy, 6 Geography and Environment, University of Southampton, Southampton, UK, 7 College of Natural Science and Mathematics, University of Alaska, Fairbanks, AK, USA, 8 Department of Ecology and Evolution, University of Lausanne, Lausanne, Switzerland, 9 Institute of Earth Science Dynamics, University of Lausanne, Lausanne, Switzerland, 10 Department of Plant Biology and Ecology, Faculty of Sciences and Technology, University of the Basque Country, UPV/EHU, Leioa, Spain, 11 IKERBASQUE, Basque Foundation for Science, Bilbao, Spain

### Edited by:

Valentí Rull, Institute of Earth Sciences Jaume Almera (ICTJA-CSIC), Spain

### Reviewed by:

Donatella Magri, Sapienza University of Rome, Italy Rüdiger Grote, Karlsruhe Institute of Technology, Germany

\*Correspondence: Rachid Cheddadi<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Agroecology and Land Use Systems, a section of the journal Frontiers in Plant Science

> Received: 09 June 2016 Accepted: 06 October 2016 Published: 25 October 2016

### Citation:

Cheddadi R, Araújo MB, Maiorano L, Edwards M, Guisan A, Carré M, Chevalier M and Pearman PB (2016) Temperature Range Shifts for Three European Tree Species over the Last 10,000 Years. Front. Plant Sci. 7:1581. doi: 10.3389/fpls.2016.01581 We quantified the degree to which the relationship between the geographic distribution of three major European tree species, Abies alba, Fagus sylvatica and Picea abies and January temperature (Tjan) has remained stable over the past 10,000 years. We used an extended data-set of fossil pollen records over Europe to reconstruct spatial variation in Tjan values for each 1000-year time slice between 10,000 and 3000 years BP (before present). We evaluated the relationships between the occurrences of the three species at each time slice and the spatially interpolated Tjan values, and compared these to their modern temperature ranges. Our results reveal that F. sylvatica and P. abies experienced Tjan ranges during the Holocene that differ from those of the present, while A. alba occurred over a Tjan range that is comparable to its modern one. Our data suggest the need for re-evaluation of the assumption of stable climate tolerances at a scale of several thousand years. The temperature range instability in our observed data independently validates similar results based exclusively on modeled Holocene temperatures. Our study complements previous studies that used modeled data by identifying variation in frequencies of occurrence of populations within the limits of suitable climate. However, substantial changes that were observed in the realized thermal niches over the Holocene tend to suggest that predicting future species distributions should not solely be based on modern realized niches, and needs to account for the past variation in the climate variables that drive species ranges.

Keywords: Holocene, past climate reconstruction, niche conservatism, Abies, Fagus, Picea

# INTRODUCTION

Changes in climate may cause shifts in the geographic distribution of species (Parmesan and Yohe, 2003; Root et al., 2003), increase extinction rates (e.g., Thuiller et al., 2005), and alter provision of ecosystem services (Schröter et al., 2005). Adaptation of human societies to these global changes requires accurate predictions of the future potential distributions of key species, such as endemic forest trees. Prediction of potential distributions of trees generally involves developing models of the relationship between current climate, distribution, and range dynamics (Guisan and Zimmermann, 2000; García-Valdés et al., 2013), then applying these models to data derived from models of potential future climate. These predictions rely on key aspects of the species climate requirements. In this context, one assumption for predicting species distributions, using nichebased models, is that species occupy all or most suitable geographic areas (Svenning and Skov, 2004). A second necessary assumption for these statistical models is that the species realized environmental niche (Hutchinson, 1957), as determined from empirical observations, remains stable over time, and space, a phenomenon called ecological niche conservatism (Wiens and Graham, 2005; Pearman et al., 2008). The validity of these assumptions remains a matter of contention.

One reason for uncertainty regarding niche conservatism is that species ecological flexibility and/or genetic adaptation can impact the observed degree of conservatism and the accuracy of predicted species distributions. Veloz et al. (2012) find at the generic level that between the late glacial period and today, the realized niches of Fraxinus, Ostrya/Carpinus, and Ulmus shifted substantially while the niches of other taxa, such as Quercus, Picea, and Pinus strobus, remained relatively stable. On the other hand, studies using current climate and species distributions to compare species realized niches between native and invasive ranges indicate that niches tend to remain stable during species invasions, but many exceptions (niche shifts) were also reported (Petitpierre et al., 2012; Guisan et al., 2014; Early and Sax, 2014). Additional studies suggest that over spans of decades or a 100 years, species environmental niches rarely change, but as time scales increase to millions of years, niche lability increases (Peterson, 2011).

The temporal niche lability of tree species may be partly due to adaptive evolutionary processes. Genetic modifications may be observed spatially (Heywood, 1991), however, they seldom develop over relatively short periods such as decades while they easily happen throughout several millennia, such as the Quaternary, with distinct and contrasting climate periods (Davis and Shaw, 2001; Davis et al., 2005). Molecular markers from the nuclear and organelle genomes have allowed reconstruction of the evolutionary history of plant species since the last glacial period (Taberlet et al., 1998; Petit et al., 2003). Using neutral markers, phylogeographic studies show that the modern spatial genetic structure over the range of many temperate tree species originated during the past few millennia (Hewitt, 2000, 2004; Petit et al., 2003). These phylogenetic studies show that there are strong interactions between past climate change, changing species ranges, and genetic constitution. While climate variability during the Holocene period was much reduced compared to the climate transitions between glacials and interglacials, species had to migrate and adapt nonetheless during the post-glacial recolonization process, which may have had an impact on their climate requirements and tolerances. The Holocene seems to be an ideal time period for investigating potential temperature shifts of tree species because (1) its reduced climate variability and relative favorability compared to the glacial period allowed taxa to migrate to either track climate changes and potentially to adapt and (2) there are more fossil data-sets for species occurrences and past climates than for older Quaternary time periods.

This study uses observational data at the Holocene time scale to investigate whether or not some European tree species remained within their modern temperature range over the last 10,000 years. Obviously, temperature is only one component of the species niche. However, any shift of the modern range/climate relationship would challenge the assumption of niche conservatism. We examine whether January temperatures experienced by Abies alba, Fagus sylvatica, and Picea abies are similar to those experienced over approximately 10,000 years during the Holocene. Substantial changes in realized thermal niches over the last several millennia would demonstrate the difficulty of predicting future range shifts solely based on contemporary realized climatic niches and further reflect on the generality of the assumption of niche conservatism.

# MATERIALS AND METHODS

# Quantifying Past January Temperature

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:

$$T\text{jan}\ (s) = \operatorname\*{argmax}\left[\left(\prod\_{t \ge \iota\_1}^{\text{tax}\_{\iota}} \mathbf{p} \mathbf{d} f^{W\_i(s)}\right)^{\left(\sum \mathcal{W}\_i(s)\right)^{-1}}\right]$$

FIGURE 1 | Schema summarizing the approach used for quantifying a past climate variable from a set of plant species that has been matched to a fossil pollen assemblage. The modern species ranges are obtained from a plant database (A). The climate variable (January temperature) is obtained from the WorldClim database (Hijmans et al., 2005; B). A probability density function (pdf) is used to infer the species-climate relationship (C) and the combination of all pdfs provides the intersecting value (D).

where:

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

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

Wi(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.

 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

using their overall Holocene Tjan range (blue histograms) based on reconstructed Tjan between 10 and 3 ka (see Figure 3). Numbers over the histograms correspond to the percentages of pollen samples that deviate from 0. Between 83 and 89% of the reconstructed Tjan deviate by less than 1◦C to the observed Tjan values from WorldClim (Hijmans et al., 2005). The map in the right panel shows the location of the modern pollen samples used for the dissimilarity analysis.

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.

# Recovering Past Species Ranges

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.

# Recovering January Temperature Range for Each Species

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).

# RESULTS

# Reconstructed January Temperature

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 TABLE 1 | Comparison of the medians of the reconstructed Tjan at different time slices in the past (10 to 3 ka and the full range as well) for each species with their modern range.

We present Wilcoxon tests (p-value) and estimation of the standard error (SE) related to the number of sites where occurrences (OCC) are identified for the three studied species. Since the Tjan medians of all time slices (n = 9) for the three species have been compared to the same one (present), the p-value (α = 0.05) has been adjusted to α/n = 5.6e-3 in accordance with the Bonferroni criterion. The bold p-values are significant at this corrected threshold. P-value and p-value (w) correspond to the Wilcoxon test, with unweighted and weighted Tjan values respectively (see Figure 4).

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.

# January Temperature Range of Focal Species

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**).

# DISCUSSION

The use of observed data from fossil records instead of data from model simulations to evaluate the relationship between past species ranges and climate is a challenging issue. First, data on past climate and on species distributions should originate from two independent data-sets and both should have a spatial coverage that minimizes errors of spatial interpretation. Second, the methods for identifying species, recovering their past occurrences and reconstructing their contemporaneous climate values from different proxies should be accurate, reproducible, and have little error. Over the past decades, pollen data have proven to be excellent proxies for the reconstruction of past species distributions (Huntley and Birks, 1983) as well as past climates (Prentice et al., 1991). However, using fossil pollen data to evaluate the relationship between a species occurrence and any climate variable through time requires two important assumptions: First, fossil pollen data represent an image of an ecosystem that is in equilibrium with contemporaneous climate. Second, the observed changes in the past ecosystems are a result of individualistic responses of taxa to climate (Webb, 1986).

Paleoecologists also assume that species evolution is a slow process over the late-Quaternary because fossil data do not provide information on species genetic adaptation to changing climate but rather on species persistence in situ, their range expansion/contraction through migration and their extinction (Davis et al., 2005). Fossil pollen assemblages represent a static image of their originating ecosystem, one that prevents evaluating the level of competition between taxa. In the present study, we considered these fossil data as an instantaneous record of that species that occurred in various past ecosystems, each of which occupied portions of a common climate space.

FIGURE 5 | Thermal amplitude of (A) Picea abies, (B) Fagus sylvatica, and (C) Abies alba for each 1000-year time slice between 10,000 and 3000 years cal. BP (gray boxplots), the Holocene overall (white boxplots), and the modern (orange boxplots) ranges, respectively. The black line is the median, the boxes represent the first and third quartiles (25 and 75th quartiles, respectively), and the whiskers represent the minimum and the maximum Tjan range. We consider the boxes and the whiskers as the Tjan ranges where populations are frequent and sparse, respectively. Red arrows indicate differences between current limits of Tjan distribution and those over the Holocene. Maps show today's areas of Picea abies (A′ ), Fagus sylvatica (B′ ), and Abies alba (C′ ) where Tjan values correspond to the range between the first and third quartiles (orange, corresponding to areas where the species are most abundant) and more extreme values (yellow). The red bars inside the boxes correspond to the estimated uncertainty related to random spatial sampling and to sample size in Monte Carlo simulations.

Our results demonstrate that within Holocene ecosystems, the thermal ranges of P. abies and F. sylvatica differ significantly from the corresponding modern ones. These two species currently occur in areas where Tjan is either warmer (F. sylvatica) or cooler (P. abies) than the occupied regions of the Holocene. Conversely, A. alba is roughly within the same thermal range as during most of the Holocene. The finding that the thermal ranges of two species have varied during the Holocene is partially consistent with previous studies that used modeled Holocene temperatures (Maiorano et al., 2013).

Discrepancies between modeled distributions and ones observed in data on fossil material (e.g., Pearman et al., 2008) may be partially explained by the relationship between the numerical abundance and geographic distributions of species, which is a matter of continuing debate. The traditional view is that species are most abundant near their geographic range center and less so toward their range edge (Brown, 1984). However, when considering the geographic distribution of species based on their environmental niche (i.e., the niche-biotope duality, see Guisan et al., 2014), findings are less conclusive for the hypothesis that species should be most abundant at the center of their environmental distribution (Sexton et al., 2009). Here, in close analogy to abundance, we consider frequency of occurrence along a climate gradient, rather than in geographical space. Thus, we treat the Tjan median value of each tree species as the point of greatest frequency of population abundance.

Today P. abies is at its greatest abundance in areas where Tjan is between ca. −12 and −6 ◦C (1st and 3rd quartiles, **Figures 5A,A'**). Assuming adequate plasticity and genetic variation, populations in these areas may potentially withstand a Tjan increase of about 4.5◦C in situ (**Figure 5A**), while a similar temperature increase may have greater effects on populations that grow at the upper Tjan limit (mainly from Romania to the Dalmatian coast). The 1◦C difference between the Holocene ranges and the modern one is not significant if we take into account the standard error related to the potential bias linked to the distribution of past species occurrences (**Table 1** and **Figure 5**). However, an increase of Tjan may have greater effects on Picea populations at the warm (southern) edge of the species range than on climatically central populations.

Fagus populations are frequent nowadays in areas where Tjan is between ca. −2.5 and 1.5◦C (inter-quartile range, **Figures 5B,B'**), but currently occur infrequently in areas with Tjan values that are 5◦C lower than those observed during the Holocene (**Figure 5B**). Unlike P. abies, F. sylvatica populations frequently occur today in areas that have Tjan values that are as high as those of the Holocene. Such thermal range occupancy suggests that, where it is frequent today, Fagus will likely be threatened by increasing temperature, while the sparse populations that occur at the cold edge of the distribution may withstand an increase of Tjan of up to 5◦C. F. sylvatica is more sensitive to cold winter temperature and spring frost than P. abies, so climate warming may promote its expansion into areas where the range of the two taxa overlap, such as southern Sweden and/or potentially in the mid to high elevations in Poland and the Carpathians. Simulations for the next century show that F. sylvatica may have to migrate further north while experiencing a loss of up to 29% of its current habitat (Kramer et al., 2010). This projection may well be valid since populations are most frequently located in areas that are at the highest Tjan values observed in the Holocene thermal range (**Figure 5B**). Furthermore, the species migration rate (100–250 m.yr−1 ) estimated by Tinner and Lotter (2006) is lower than the velocity of the expected temperature change (∼800 m.yr-1) for the temperate deciduous biome (Loarie et al., 2009). Nonetheless, the magnitudes of these velocities are highly dependent on the resolution of the underlying climate data-sets (Dobrowski et al., 2013).

The patterns shown by A. alba differ from those of P. abies and F. sylvatica. Unlike the other two species, the full modern Tjan range of A. alba is at its greatest extent since the early Holocene (**Figure 5C**). However, the most frequent populations occupy a much narrower Tjan range (ca. −3.5 to −1.5◦C, **Figures 5C,C'**) than during the Holocene (ca. −4 to +0.5◦C), which suggests these environmentally central populations may be able to withstand a somewhat wider temperature range than they currently experience (**Figure 5C**). Conversely, there are sparse A. alba populations today that experience Tjan as high as +10◦C, which is about 2◦C higher than for any known Holocene population (**Figure 5C**). Thus, the modern populations at the warmer end of the species distribution (Mediterranean mountains) may be threatened by further climate warming. While ecological niche models suggest that A. alba currently occupies less than 50% of its potential range (Svenning and Skov, 2004), our results suggest, in contrast, that this species currently occupies its full thermal range, and that the distribution of the populations within that thermal range may change as Tjan increases. This may occur through latitudinal migration, simple expansion of population density within other areas of the modern distribution (but without any further geographical expansion; see Figure 1E in Guisan et al., 2014), or expansion to fill sparsely and unevenly populated areas. The discrepancy with Svenning and Skov (2004) may be related to potential deviations during the interpolation of their climate data set, or the particular modeling algorithm they employ.

Ecological modelers and palaeoecologists both aim to address a crucial question: will major tree species withstand ongoing global warming in situ (because they can tolerate a wider temperature range than what they currently encounter), migrate to areas that acquire suitable conditions, or simply decline. Species modeling approaches suggest that the impact of future climate change will be substantial in Europe (Thuiller et al., 2005). On the other hand, some fossil-based distribution data, such as for A. alba, suggest little risk will arise with climate warming (Tinner et al., 2013). In contrast, a recent palaeoecological study (Seppä et al., 2015) shows that the range of hazel (Corylus avellana) previously expanded about 4◦ northward from its current geographical limit in Fennoscandia as a response to a temperature increase of 2.5◦C during the mid-Holocene, constituting an increase similar to some projections for response to on-going climate change. Species responses to climate warming may vary substantially and, thus, require focused study. While we concentrated on three tree species because of the direct correspondence between fossil pollen

taxa (identified to genus) and species-level identification (e.g., compared to the 22 European species Quercus that produce indistinguishable pollen), one might expect that the thermal ranges of additional tree species may have been different in the past. Additional gridded proxy data, other than pollen, as well as identification of fossil taxa using ancient DNA (Parducci and Petit, 2004) will be needed to address this issue.

Evaluating the climatic limits of species distributions through time is a key issue to understanding species responses to future climate change. Extended data-sets, such as those arising from fossil pollen, are sufficient for reconstructing past geographical distributions of many taxa and for estimating important climate variables. However, they do not contain any information about the adaptive evolution of species over the period of interest (i.e., the Holocene). When assessing the relationship between past species distributions and climate change, paleoecologists tend to assume that ecological responses, such as range shifts, migration rates, and recolonization, are the dominant processes that impact observed species distributions (Huntley et al., 1989; Huntley and Webb, 1989). However, other biotic responses such as phenotypic plasticity and genetic adaptation were also important during the Quaternary (Davis and Shaw, 2001). While evolutionary rates are likely slower than rates of ecological response to climate change, at least during the late Quaternary, promising developments in ancient DNA technology (Giguet-Covex et al., 2014) may help to determine whether or not species genetically adapted to past climate change.

Fossil records may indicate on-going and successful accommodation to climate when a species range remains in dynamic equilibrium with temperature over time (at least during the late-Quaternary period; Davis and Shaw, 2001). The assumption that species distributions in the past (or today) are at equilibrium with their contemporary climate has yet to be demonstrated. If species are not at equilibrium with climate, then there may be little justification to expect stasis in realized niches through space and time (Araújo and Peterson, 2012). Based on observed data, Araújo and Pearson (2005) suggest that assemblages of plants and breeding birds are closer to equilibrium with climate than other organisms and, therefore, their future ranges may be well predicted by niche-based models. Recent work based on fossil data, however, suggests significant potential disequilibrium between vegetation and climate at both leading and trailing edges of species ranges (Svenning and Sandel, 2013). This apparent disequilibrium arises in part from the estimated time lag of the response of individual species ranges to post-glacial warming (Webb, 1986). Thus, strong differences among species in the degree of distribution equilibrium may occur due to varying dispersal abilities, competition with other species, soil conditions, and other factors that may affect the ability of species distributions to track changing climate.

# CONCLUSIONS

Modeling approaches are necessary to project how species distributions could change in the future. Here we show that fossil data provide valuable information on how temperature changes might affect different parts of the range of a species in different ways. We demonstrate how fossil data can complement data on contemporary species distributions to provide valuable information on species tolerances to climate variation. Our results reveal changes in the thermal ranges of the three species with respect to an influential climate variable over the past 10,000 years. The results further suggest that differing impacts of future climate warming, depending on the density (sparsity) of populations. In parallel with approaches based on climate models, data on past species occurrences and reconstructed climates derived from them can provide independent estimates of whether species can cope with climate change and, potentially, what parts of their distributions will be the most threatened.

Nevertheless, the past reconstructions need to be constrained by an evaluation of the infra-specific adaptive ability of species, their dispersal capacity, their interspecific competitive ability, and the ecophysiological relationship of each species to different climate variables. Moreover, one climatic variable such as Tjan, may be considered as restrictive for evaluating the ranges of other species than our focal species. Therefore, other additional climate variables are necessary to better define the climatic niche of most plant species and their potential adaptive capacity through time. Water availability, through the annual amount of precipitation and its seasonal distribution, is a key variable for species spread, their persistence and populations expansion. Without an evaluation of these variables, the question remains as to whether a species is capable of withstanding a wider climatic range than that occupied currently and, thus, has portions of the full thermal range that remain unoccupied.

# AUTHOR CONTRIBUTIONS

RC performed all data analysis, programming, computational work, figures, and tables and has written the original manuscript. PBP has improved several versions of the manuscript and contributed substantially, in a decisive way to its final version. MCa contributed to improving many statistical aspects of the study and particularly to estimating the errors through the Monte-Carlo simulations.  All co-authors have contributed to designing the study and refining its objectives in the frame of the EU ECOCHANGE project. MA, LM, ME, AG, MCa, and PBP have contributed to improving the manuscript at different stages.

# ACKNOWLEDGMENTS

This work was funded by the EU ECOCHANGE project (Contract No. 066866 GOCE granted to AG, MA, ME, and RC). The original pollen data used in the present study were obtained from the European Pollen Database (http:// www.europeanpollendatabase.net). We thank N. Zimmermann for helpful discussion throughout the project and especially during the Bonn meeting. This is ISEM contribution number 2016–206.

# REFERENCES

**Conflict of Interest Statement:** 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.

Copyright © 2016 Cheddadi, Araújo, Maiorano, Edwards, Guisan, Carré, Chevalier and Pearman. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Filling a Geographical Gap: New Paleoecological Reconstructions From the Desert Southwest, USA

Andrea Brunelle1 \*, Thomas A. Minckley 2 , Jacqueline J. Shinker 2 and Josh Heyer 1

*1 Department of Geography, University of Utah, Salt Lake City, UT, United States, 2 Department of Geography, University of Wyoming, Laramie, WY, United States*

In 1916 the time stamp for quantitative palynology was set with Lennart von Post's initial paper on pollen analysis and environmental change in the Scandinavian peat bogs. In the 1930s, Von Post provided a map of the known palynological reconstructions. This map showed many conspicuous gaps of geographic coverage that have endured to this day. In particular, environmental reconstruction in arid lands remain much less known, largely due to the paucity of depositional environments in these areas as well as challenges with preservation in regions either uniformly dry or subjected to strong wet/dry phases. Over the last decade we have examined linkages between desert wetland development and episodes of wet and arid conditions. Desert wetland, or ciénegas, are recharged by groundwater and appear to be sensitive to climate-driven groundwater fluctuations. These systems appear to "grow" during wet periods potentially associated with enhanced El Niño activity, suggesting an important linkage with groundwater dynamics and the quantity and frequency of winter precipitation delivery. Hydrologic conditions in ciénegas are also important controls on the preservation of pollen, where episodes of aridity coincide with periods of poor pollen preservation. We assess modern El Niño events, as analogs of past wet conditions, to provide context on the atmospheric controls for delivery of moisture into the desert southwest during winter. Our analysis shows that anomalously high and persistent moisture delivery into the region during El Niño events enables the growth of ciénegas, improves preservation of pollen and promotes the growth of fuels necessary to support wildfire. This paper examines ciénega sites located in the southwestern region of North America at the US/Mexico Border and discusses results that addresses a geographical gap identified by Von Post's original work in paleoenvironmental research.

Keywords: pollen, charcoal, ciénega, desert wetland, ENSO, El Niño

# INTRODUCTION

The legacy of Lennart Von Post on quantitative paleoecology and palynology cannot be overstated. While Von Post did not invent palynology, he helped to cement this field of study through his foundational research conducted in Sweden and Norway during the late 1800s and early 1900s. Von Post's approach to paleoenvironmental reconstruction made its way to the American Deserts via the geochronological work of Ernst Antevs and American researcher Kirk Bryan (Manten, 1967; Smiley, 1974). While Antevs did not focus on pollen analysis himself, he was interested enough in the technique to opine on the naming of the methodology (Antevs, 1944) and his work was

### Edited by:

*Jesse L. Morris, Weber State University, United States*

### Reviewed by:

*Michael-Shawn Fletcher, University of Melbourne, Australia Alberto Saez, University of Barcelona, Spain*

\*Correspondence: *Andrea Brunelle<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Quaternary Science, Geomorphology and Paleoenvironment, a section of the journal Frontiers in Earth Science*

> Received: *12 December 2017* Accepted: *05 July 2018* Published: *31 July 2018*

### Citation:

*Brunelle A, Minckley TA, Shinker JJ and Heyer J (2018) Filling a Geographical Gap: New Paleoecological Reconstructions From the Desert Southwest, USA. Front. Earth Sci. 6:106. doi: 10.3389/feart.2018.00106*

critical to subsequent paleoenvironmental work in the desert landscapes of southwestern North America (hereinafter "southwestern deserts").

Early pollen work in the southwestern deserts focused on woodland expansion (Clisby and Sears, 1956; Martin and Gray, 1962; Martin, 1963; Mehringer, 1965; Meyer, 1973). This initial research reported observations of periods where no pollen or low pollen abundances were recovered. These results altered the relative abundances of pollen types, which was interpreted to be caused by differential preservation and the individual researchers ability to identify degraded pollen grains. For example, through careful examination, the crumbled or degraded pollen grains of Pinus, Artemisia, Poaceae, and Amaranthaceae are easily distinguishable by the presence of distinct wall or annulated pore characteristic of the particular type. However, degradation and preservation issues reveal the difficult nature of arid land palynology. In short, the lack of perennial lacustrine environments results in discontinuous records of environmental change. Further, simple abundance in arid lands seems lower than one might experience in temperate environments. Subsequent research concentrated on macrobotanical records within woodrat middens (see Betancourt et al., 1990), instead of the ephemeral sedimentary nature of desert wetlands (however see Davis et al., 2002; Mensing et al., 2008).

In arid regions, environments suitable for pollen preservation are limited to riparian deposits, ciénegas or wetlands (Davis et al., 2002; Minckley and Brunelle, 2007), or on the deflated surfaces of former lake beds (Clisby and Sears, 1956; Martin, 1963; Meyer, 1973). Of these environments, ciénegas are proven archives of past environments and climate conditions as they collect charcoal, pollen, and macrofossils over thousands of years (Martin, 1963; Hall, 1985; Hendrickson and Minckley, 1985; Davis et al., 2002; Minckley and Brunelle, 2007; Brunelle et al., 2010; Pigati et al., 2014). The global distribution of desert wetland deposits has not been mapped, but they are considered to be widespread, though in other literature they are referred to as ground water deposits associated with hot or cold springs, seeps, or other wetland complexes (Pigati et al., 2014; Sáez et al., 2016). However, systematic studies of these complex systems are relatively recent, variously focusing on lithological, isotopic, or biological indicators (Minckley and Brunelle, 2007; Quade et al., 2008; Minckley et al., 2013a; Brunelle et al., 2014; Pigati et al., 2014 and references therein, Sáez et al., 2016).

Arid and semi-arid environments occupy 40% of the terrestrial land area on earth and in many cases represent areas where recent human population growth rates are high relative to more mesic regions (18.5% between 1990 and 2000; Hassan et al., 2005). The consumption of natural resources associated with rapid population growth places intense pressure on water resources, such as ciénegas. In addition to being a ready source of water, ciénegas provide crucial ecosystem services including erosion mitigation and prevention of the associated lowering of the water table. Ciénegas have also been identified as regions of high conservation priority for many reasons including endemism (Minckley, 1969, 1992; Hendrickson and Minckley, 1985; Abell, 2000). In Arizona, USA where wetland environments occupy ∼2% of the land area, ciénegas are critical habitat for at least 19% of the threatened, candidate, or endangered species within the state (Baker et al., 2004). However, beyond species of concern, desert ciénegas and riparian corridors may increase regional biodiversity by up to 50% in some cases (Sabo et al., 2005).

Over the last decade, we examined sediment collected from ciénega systems located in the southwestern deserts to understand the hydrology, ecology, disturbance regimes, climatology, and environmental change from these critical systems. Our work builds on the legacies of scholars working in this region, including Antevs, Martin, and Mehringer (e.g., Antevs, 1962; Martin, 1963; Mehringer, 1967). Consideration of ciénegas as repositories of environmental history was largely untouched until Hendrickson and Minckley (1985) published "Ciénegas—Vanishing Climate Communities in the American Southwest." Paleoecological research of these sedimentary environments was largely unexplored until Davis et al. (2002) published a synthesis of several sites that examined changes in burning regime with European settlement. However, a largescale evaluation of ciénega form and function has not yet been undertaken (see Pigati et al., 2014).

In this paper, we present data from four ciénega sites in the Arizona/New Mexico-Mexico borderlands using multiproxy analyses. Herein, we present data from previously published research as well as new records to provide a regional perspective on vegetation and disturbance history from the borderlands region. We detail our approach to handling differential pollen preservation and explain how sedimentary indicators may also provide evidence for hydrology as a control on pollen preservation. Despite the analytical challenges associated with ciénega sediments, the environmental record we present here shows regional coherence of wet-dry episodes. We show that pollen records from ciénegas provide insight to past climatological states related to variations in winter precipitation. For the desert Southwest, this is currently associated with quasiperiodic El Niño Southern Oscillation (ENSO) variations. In the tradition of Von Post, our first objective is to identify conditions in the sedimentological proxy data (e.g., pollen and charcoal) that record changes in ciénega systems attributed to climate. For our second objective, we identify modern episodes that are analogous to conditions found in the paleoecological record (strong El Niño conditions) associated with ciénega growth (e.g., wet episodes) to assess the dominant atmospheric drivers related to enhanced winter precipitation in the arid southwest.

To address our first objective, we utilize indices of pollen preservation and charcoal abundance. We assess pollen preservation in ciénega systems of southwestern New Mexico and southeastern Arizona as well as the relationship to past ENSO variability based on our modern understanding of ENSO and southwestern North American hydroclimate teleconnections (Heyer et al., 2017). Further, we put forward past fire activity near the ciénegas was covariant with ENSO variability, similar to connections observed in historic southwestern forest fire research (Swetnam and Betancourt, 1990). Considering fuel availability produced during antecedent wet seasons is closely linked to wildfire in arid and semi-arid environments globally (Krawchuk and Moritz, 2011), we consider our linkage of past fire activity and ENSO variability viable. To address our second objective, on the impact that El Niño events have on Southwestern climate, we use high-resolution modern climate data to assess teleconnections between moisture delivery and southwestern precipitation anomalies during winter and spring of strong El Niño events as a mechanism for sustained wet conditions that are conducive for pollen preservation in ciénegas.

# REGIONAL DESCRIPTION

The southwestern deserts, particularly those in the Sonoran and Chihuahuan systems, have a high concentration of ciénegas (Hendrickson and Minckley, 1985; Minckley et al., 2013b; Cole and Cole, 2015). These wetlands are found in a variety of geomorphic contexts along riverine or stream channels, offchannel seeps, landslide backfill, and valley bottoms, to name a few. Often ciénegas are proximal to fault lines, which are common in the desert basin and range structures of western North America (Minckley, unpublished).

The modern climatic context of the arid southwest of North America is characterized as having bimodal precipitation (Mock, 1996; Shinker, 2010). Winter precipitation has broad spatialscale patterns with persistent light precipitation delivered by midlatitude frontal systems (Mock, 1996; Comrie and Glenn, 1998; Shinker, 2010). In contrast, summer precipitation is highly variable spatially and typically short in duration associated with convective thunderstorms delivered via the North American Monsoon System (NAMS; Mock, 1996; Adams and Comrie, 1997; Sheppard et al., 2002; Diem and Brown, 2006; Shinker, 2010). The modern NAMS developed after climate started to warm and the Laurentide ice sheet began to collapse at the end of the last ice age. This warming led to the formation of summer thermal lows that provide a mechanism for drawing moisture from the Gulf of Mexico and southeastern Pacific (Spaulding, 1991; Adams and Comrie, 1997; Minckley et al., 2004). The bimodal nature of precipitation in the southwest, especially the extra precipitation in the summer associated with the NAMS, provides effective moisture at the surface related to water- and energy-balance components important for ecosystem and water resources related to vegetation and the ciénega systems of the southwest (Shinker and Bartlein, 2010).

In addition to the bimodal climate of the southwest, largescale modes of variability such as ENSO impacts southwestern hydroclimate and ecosystems (Ropelewski and Halpert, 1986; Cole and Cook, 1998; Gershunov and Barnett, 1998a; Cayan et al., 1999; Dettinger et al., 2000). El Niño is the anomalous warming of the near-surface water off the western coast of South America that initiated about 5000 years ago (Enfield, 1992; Liu et al., 2014). Its atmospheric counterpart is the Southern Oscillation, which reflects changes in barometric pressure (Wallace and Gutzler, 1981; Enfield, 1992). Strong El Niño years in the modern record typically result in higher-thannormal winter precipitation in the southwestern United States; while strong La Niña typically causes lower-than-normal winter precipitation and is typically associated with regional drought (Shinker and Bartlein, 2009). Studies have identified a northto-south precipitation dipole correlated with ENSO in the western U.S., showing a strong relationship in the southwest between cool season (Oct-Mar) precipitation and El Niño events (Dettinger et al., 1998; Wise, 2010). Further, many studies have shown statistically significant correlations between ENSO and western North America precipitation, specifically in the southwest where spatially coherent correlations are found during the cool season (Dettinger et al., 1998; Heyer et al., 2017). However, various modes of ENSO are associated with different impacts on western NA hydroclimate (Ashok et al., 2007), and teleconnections can vary on decadal-to-centennial time scales (Gershunov and Barnett, 1998b), suggesting that southwest hydroclimate does not always respond to ENSO variability in a linear manner based on our modern understanding of southwest hydroclimate and ENSO correlations (Heyer et al., 2017). However, despite teleconnection variability (Gershunov and Barnett, 1998b), dipole variability (Wise, 2010) and different modes of ENSO (Ashok et al., 2007), paleo-sedimentary records suggest teleconnections between ENSO and southwest climate are similar to those observed in the modern record, and have persisted at century–millennial time scales (Barron and Anderson, 2010; Antinao and McDonald, 2013; Kirby et al., 2014; Hart et al., 2015). Further, paleoecological records have revealed the ENSO phenomenon is correlated with fire regimes of the southwest (Swetnam and Betancourt, 1990). Therefore, our research on San Bernardino Ciénega suggests that over centuries-to-millenial scales fire frequency trends are consistent with ENSO variability, with fire frequency increasing (decreasing) when El Niño (La Niña) was active (Brunelle et al., 2010).

Our work has been focused on ciénegas near the US/Mexico borderlands region in southern Arizona and New Mexico (**Figure 1**) in the monsoon-dominated ecotone of the Sonoran and Chihuahuan deserts and in a region well correlated with ENSO activity (Wise, 2010; Heyer et al., 2017). Total annual precipitation is around 254 mm (Western Regional Climate Center, 2017). Summer temperatures are generally around 37◦C while winter temperatures rarely fall below freezing.

# SITE DESCRIPTIONS

San Bernardino Ciénega (31.3◦N; 109.3◦W, 1,160 m asl; **Figure 1**) is located in the drainage of Black Draw Wash/Rio de San Bernardino (RSB) of southeastern Arizona, USA and northeastern Sonora, MX near the ecotone of Chihuahuan and Sonoran grassland and desert scrub (Minckley and Brunelle, 2007; Minckley et al., 2009, 2011; Brunelle et al., 2010). Currently the ciénega surface is dry except for a few artificial impoundments and small perennial springs.

Cloverdale Ciénega is located in the Peloncillo Mountains, Coronado National Forest, NM (31◦ 26.141′ N, 108◦ 58.517 W, 1,640 m asl; **Figure 1**). Cloverdale Ciénega is situated at the ecotone between Madrean woodland and grassland (Brunelle et al., 2014), making this location sensitive to changes in community composition over time. Currently, the Madrean woodlands form a mosaic of clustered tree stands interspersed with patches of shrub and grass (Brown, 1982).

Canelo Hills Ciénega (31◦ 33.833′N, 110◦ 31.711′W, 1,506 m asl; **Figure 1**) is a Nature Conservancy preserve located southeast

of Sonoita, Arizona on the east side of the Canelo Hills. The area is a Madrean grassland.

Babocamari Ciénega (31◦ 37.901′N, 110◦ 27.178′W, 1,391 m asl; **Figure 1**) is southeast of Sonoita, Arizona and located between the Whetstone and Huachuca Mountains. The site is considered Madrean grassland.

# METHODS

# Field

# San Bernardino

Sediments were collected from an active spring (Snail Spring), the incised channel wall of the Rio de San Bernardino arroyo (RSBA), and the ciénega surface of the San Bernardino National Wildlife Refuge (SBNWR). The RSBA samples were described and collected contiguously from a freshly exposed sediment face with a trowel. A 15-m long core was obtained from the SBNWR using a truck-mounted 5′′ (12.7 mm) hollow barrel auger (Minckley et al., 2007). The Snail Springs samples were collected using a modified Livingstone piston corer.

# Cloverdale, Babocamari, and Canelo Hills Ciénegas

Sediment cores were collected using a vibracore system where manual probes indicated maximum sediment accumulation. The vibracore apparatus uses a cement leveling motor attached to a 3-inch aluminum tube inserted into the sediments until refusal (bedrock). Sediment was retained in the aluminum tube by placing an airtight cap on the top of the tube to create a vacuum, allowing the sediment core to be removed from the ground without loss. The cores were transported in the aluminum tube to the University of Utah and split in the lab for description and sampling.

# Laboratory and Analytical

Once the cores were removed from the aluminum tube, the color and texture of the sediments were described. Every contiguous centimeter of the core was examined for charcoal analysis. For pollen analysis, samples were collected to acquire a minimum sample resolution of approximately every 200 years. The sediments were weakly stratified, massive, and organic.

Age-depth relationships were determined for the sediments based on AMS (accelerator mass spectrometry) dates on pollen concentrates, lead-210 chronology (when available), and the surface (collection year; **Table 1**). Ages are presented in calibrated years before present (cal yr BP) where 1950 CE is equal to year 0 cal yr BP. The age-depth models were generated in CLAM (Blaauw, 2010) using smoothing splines (**Figure 2**). The age model for San Bernardino was created using multiple linear regressions and was not rerun due to the multiple publications already based on the existing chronology (Brunelle et al., 2010 and citations therein).

Pollen extraction methods followed Faegri et al. (1989). For pollen analysis, we attempted to identify at least 300 terrestrial pollen grains to complete the rarefaction curve associated with pollen diversity (Maher, 1963). Pollen grains were identified at a magnification of 500X to the lowest possible taxonomic

classification using published atlases (e.g., Erdtman, 1969; Kapp et al., 2000). Counts were converted to percentages based on the sum of total terrestrial grains and influx using the sedimentation rate calculated using the age models. In desert sediments, a 300-grain terrestrial count is not always obtainable because of poor preservation (Minckley and Brunelle, 2007; Minckley et al., 2011). Our standard has been to attempt 300-grain terrestrial counts or 300–1000 Lycopodium spore (exotic spore added as a tracer) tallies for each depth. If pollen counts do not exceed at least 100 terrestrial pollen grains we plot pollen presence instead of abundance to avoid artificial inflation of identifiable pollen types as presenting ecological significance, an issue common in earlier studies as detailed in the Introduction.

A second way we avoid misinterpretation of desert pollen assemblages is by determining an index of pollen preservation ("preservation index") where we consider the number of Lycopodium tracers to the total number of identified grains. Studying the way pollen has been preserved can help to determine changes in the deposition or preservation environment. This type of examination is a relatively underused technique (Tipping, 2000) that could be more widely used as an interpretive tool (Tweddle and Edwards, 2010). In our study, the preservation index is used as an indicator of a decrease in water table height that would expose deposited grains to conditions of wetting and drying which are known to be consistent with deterioration (Campbell and Campbell, 1994; Tweddle and Edwards, 2010). The depth to water table in ciénega settings controls the water level at the surface. In turn, the water table is largely controlled by the amount of winter precipitation. High preservation indices are used to indicate wet (active) ciénega surfaces, while low indices indicate dry (inactive) surfaces.

Volumetric samples were taken from the sediments for charcoal analysis to reconstruct fire history. Sample volumes ranged from 1 to 10 cc based on the amount of charcoal found at various depths throughout the sections. Sample volume was adjusted to maintain a significant number of charcoal particles, but to reduce excessive counts (>1,000 particles). Methods and rationale follow Clark (1988) and Long et al. (1998). While taphonomic studies of large particle deposition on fluvial wetland surfaces have not been conducted, isotopic analysis of ciénega sediments indicates pooling and stagnation of surface waters (Minckley et al., 2009), suggesting sedimentation of charcoal fragments would be similar to those of lake and bog surfaces. Charcoal concentrations (particles/cm3 ) were converted to influx (particles/cm2 /yr) using the age models.

# Statistical Analyses

To compare ciénega CHAR record to regional charcoal variability, regional charcoal for Arizona, New Mexico, Southern Utah, and Southern Colorado (i.e., latitude >30, latitude <38, longitude ≤105, longitude ≥113) were collected from the paleofire database (Blarquez et al., 2014), and used in correlative analyses. A total of 13 sites were included, of which all sites are from mid-upper elevations (i.e., >2,300 meters) with the exception of the Montezuma Well site (near 1,100 meters). Prior to comparing regional charcoal to ciénega charcoal, regional charcoal data were averaged (i.e., long-term average 10,000 cal yr BP—present) and standardized to z-scores. Charcoal data were then composited using a locally weighted scatterplot smoother and bootstrapping techniques (Blarquez et al., 2014).

Prior to correlation analyses, the CHAR and pollen preservation ratio data were decomposed into pseudoannual values and then reaggregated into 100-year bins to be comparable to the Moy et al. (2002) El Niño event frequency reconstruction. The time series of all the ciénega sites were then averaged to provide a comprehensive trend of the series. To quantitatively describe the relationships between the ENSO index (Moy et al., 2002) and the preservation ratio and CHAR time series, a Pearson product-moment correlation was run on the time series providing correlation coefficients and p-values.

Correlations were performed between data from 7,700 cal yr BP to present for mean-pollen preservation, and from 8,100 cal yr BP to present for mean-charcoal accumulation. Correlations were between ciénega sites mean-pollen preservation, ciénega sites mean-charcoal accumulation, regional standardized charcoal accumulation, and El Niño events/100 years (Moy et al., 2002). From 8,100 to 5,400 cal yr BP, only the San Bernardino ciénega site was active, with pollen preservation and charcoal accumulation at or near zero, as expected considering El Niño event inactivity. Therefore, correlation analyses from 8,100 to 5,400 cal yr BP is not based on an actual average. However, we assume a lack of pollen preservation and charcoal accumulation at the inactive ciénega sites was due to a lack groundwater recharge, which may represent a reduction of persistent winter and spring precipitation potentially associated with reduced El Niño events (Scanlon et al., 2006; Heyer et al., 2017). We suggest if the other three ciénega sites had been active, then pollen preservation and charcoal accumulation would have been close to zero, similar to the San Bernardino ciénega. Thus, assuming the other ciénega sites would have also had pollen preservation and charcoal accumulation near zero, we only use pollen preservation and charcoal accumulation data from the San Bernardino ciénega to calculate the long term mean-pollen preservation (i.e., 7,700 cal yr BP—present) and mean-charcoal accumulation (i.e., 8,100 cal yr BP—present). From 5,300 to 0 cal yr BP at least two ciénega sites recorded data, allowing for a mean to be calculated.

# Modern Climate Analog Approach