Coastal Pine-Oak Glacial Refugia in the Mediterranean Basin: A Biogeographic Approach Based on Charcoal Analysis and Spatial Modelling

During the glacial episodes of the Quaternary, European forests were restricted to small favourable spots, namely refugia, acting as biodiversity reservoirs. the Iberian, Italian and Balkan peninsulas have been considered as the main glacial refugia of trees in Europe. In this study, we estimate the composition of the last glacial forest in a coastal cave of the Cilento area (SW Italy) in seven time frames, spanning from the last Pleniglacial to the Late Glacial. Charcoal analyses were performed in seven archaeological layers. Furthermore, a paleoclimate modelling (Maxent) approach was used to complement the taxonomic identification of charcoal fragments to estimate the past potential distribution of tree species in Europe. Our results showed that the mesothermophilous forest survived in this region in the core of the Mediterranean basin during the Last Glacial Period (LGP, since ~36 ka cal BP), indicating that this area played an important role as a reservoir of woodland biodiversity. Here, Quercus pubescens was the most abundant component, followed by a wide variety of deciduous trees and Pinus nigra. Charcoal data also pointed at the crucial role of this coastal area, acting as a reservoir for warm temperate trees of genera Tilia, Carpinus and Sambucus, in LGP, in the Mediterranean region. Our modelling results showed that P. nigra might be the main candidate as a “Pinus sylvestris type” in the study site in the Last Glacial Maximum (LGM). Furthermore, we found that P. nigra might coexist with Q. pubescens in several European territories both currently and in the LGM. All models showed high levels of predictive performances. Our results highlight the advantage of combining different approaches such as charcoal analysis and ecological niche models to explore biogeographic questions about past and current forest distribution, with important implications to inform today’s forest management and conservation.


Introduction
Present-day Mediterranean vegetation is the result of the interaction of several factors and processes including past glaciation, location of refuge areas, biogeographic barriers and from the middle Holocene onwards, also anthropogenic influence [1][2][3].
During the Pleistocene, the advance and retreat of the ice sheets, due to climatic oscillations, had severe impacts on the distributions of many animal and plant species, which were able to survive only in suitable unglaciated habitat available during at least part of the ice ages [4].

Study Area
The Cilento region ( Figure 1) presents several cave human settlements embracing a wide period from the Paleolithic up to the Bronze Age [29]. Serratura Cave (hereafter abbreviated as SC, Figure 1, Figure S1), located in Camerota Bay (40 • 00 N; 15 • 22 E), belongs to the geological unit of Mt Bulgheria (1,225 m a.s.l.). The SC is located ca. 20 m inland from the coastline, at 2 m a.s.l. and has an area of 70 m 2 . The cave preserves long-term Late Pleistocene-Holocene stratigraphic successions, which have a prominent position in the prehistoric framework of the southern Italian mainland [30,31].
Forests 2020, 11, x FOR PEER REVIEW 3 of 18 As species identification by charcoal analysis alone is problematic due to the absence of specific diagnostic key features, this hypothesis was explored using an ENM for all the species belonging to the P. sylvestris type found in the study site through charcoal analyses and developing models for the LGM, assessing past potential environmental suitability for the different species [15].

Study Area
The Cilento region ( Figure 1) presents several cave human settlements embracing a wide period from the Paleolithic up to the Bronze Age [29]. Serratura Cave (hereafter abbreviated as SC, Figure 1, Figure S1), located in Camerota Bay (40°00′ N; 15°22′ E), belongs to the geological unit of Mt Bulgheria (1,225 m a.s.l.). The SC is located ca. 20 m inland from the coastline, at 2 m a.s.l. and has an area of ~70 m 2 . The cave preserves long-term Late Pleistocene-Holocene stratigraphic successions, which have a prominent position in the prehistoric framework of the southern Italian mainland [30,31]. The current climate of Camerota Bay is subhumid thermo-Mediterranean. Temperatures never fall below 0 °C and the mean annual temperature is 14.8 °C (Capo Palinuro, 185 m a.s.l.). Precipitation has a mean annual value of 762 mm (referred to years 1958-1999) and is irregularly distributed throughout the year, with only 4.7% falling during the summer (July and August).
In the coastal sector, Mediterranean maquis dominates vegetation, with Pistacia lentiscus, Phillyrea latifolia, Juniperus phoenicea, Euphorbia dendroides, Calicotome villosa, Spartium junceum Myrtus communis and several Cistus species. In xeric sites degraded by recurrent wildfires, the tussock grass Ampelodesmos mauritanicus is the dominant species. Scattered large, old Quercus pubescens trees occur along the coast, whereas further inland, the vegetation is dominated by Q. ilex, accompanied by P. latifolia and M. communis and some deciduous species such as Q. pubescens, Fraxinus ornus and Ostrya carpinifolia. In the submontane and montane sector small Q. pubescens woods together with Castanea sativa and Q. cerris woods are present. From 1,000 m a.s.l. up to the treeline, large Fagus sylvatica forests are found [32]. The current climate of Camerota Bay is subhumid thermo-Mediterranean. Temperatures never fall below 0 • C and the mean annual temperature is 14.8 • C (Capo Palinuro, 185 m a.s.l.). Precipitation has a mean annual value of 762 mm (referred to years 1958-1999) and is irregularly distributed throughout the year, with only 4.7% falling during the summer (July and August).
In the coastal sector, Mediterranean maquis dominates vegetation, with Pistacia lentiscus, Phillyrea latifolia, Juniperus phoenicea, Euphorbia dendroides, Calicotome villosa, Spartium junceum Myrtus communis and several Cistus species. In xeric sites degraded by recurrent wildfires, the tussock grass Ampelodesmos mauritanicus is the dominant species. Scattered large, old Quercus pubescens trees occur along the coast, whereas further inland, the vegetation is dominated by Q. ilex, accompanied by P. latifolia and M. communis and some deciduous species such as Q. pubescens, Fraxinus ornus and Ostrya carpinifolia. In the submontane and montane sector small Q. pubescens woods together with Castanea sativa and Q. cerris woods are present. From 1000 m a.s.l. up to the treeline, large Fagus sylvatica forests are found [32].

Charcoal Analysis
Sediment samples were collected during the archaeological excavations [30,31] and then sieved in situ by water through a sieving column. All charcoal fragments ranging between 2 and 4 mm mesh sizes were sorted under a dissection microscope and then analysed using a reflected light microscope (100X-1000X). Taxonomic identification relied on the reference collection of plant and wood anatomy, and wood anatomy atlases [33][34][35][36][37]. Relevant specific literature was used to reach the species level identification in the taxonomic group of deciduous Quercus [36].
We analysed charcoal from the Late Pleistocene layers, whose stratigraphy, the corresponding 14 C dating and cultural facies, are shown in Table 1. Sampling layers were selected to collect only scattered charcoal (sensu Chabal) [38], because these fragments, resulting from long-term burning activities, can be considered representative of local vegetation and thus suitable for paleoecological studies [10,[39][40][41][42]. A minimum of 200 charcoal fragments were examined for each layer except Layer 11, where the available number of fragments was limited to 100. Charcoal frequencies were calculated for all layers. While quantitative data were not shown for Layer 12, a fireplace where charcoals probably represent the remnants of a single burning event of collected wood [38].

Training and Projection Area
As training and projection areas we considered the European and North Africa territories comprised between latitudes 30 • N-70 • N and longitudes-15 • W-35 • E.

Data Collection
To collect occurrence records of the taxa of interest identified through charcoal analysis, we used several sources as: (1) public access databases, including the European forest genetic resources programme (http://www.euforgen.org) and the European Information System on Forest Genetic Resources (http://portal.eufgis.org/) and (2) shapefiles obtained by Caudullo et al. [43]. We screened all the records in ArcGis (version 10.2.2; http://www.esri.com/software/arcgis) for spatial autocorrelation using average nearest neighbour analyses and Moran's I measure of spatial autocorrelation to remove spatially correlated data points and guarantee independence [44,45]. After spatial autocorrelation analysis of our dataset, to generate ENMs we used only fully independent presence records falling within the native range as described by Caudullo et al. [43] (Supplementary Materials, Figures S2-S5).

Environmental Variables
To build the ENMs for our taxa, we started from a set of 19 bioclimatic variables obtained from the WorldClim database vers. 2.0. (www.worldclim.org/current) [21]. We downloaded the bioclimatic variables in a consistent format (ESRI grid file) and resolution (2.5 min resolution, approximately Forests 2020, 11, 673 5 of 18 5 km grid cell sizes at the equator). Then, we converted all the bioclimatic variables in ASCII files and generated a Pearson's correlation matrix with SDMtoolbox (version 2.2) [46] in ArcGis (version 10.2.2; http://www.esri.com/software/arcgis). From the matrix, we selected only the variables for which r < 0.70 [47][48][49] in order to remove any variables that could be highly correlated with one another before developing the models. Pearson's correlation matrix led to a final set of 5 climatic variables: temperature seasonality (%), mean temperature of the wettest quarter ( • C), mean temperature of the driest quarter ( • C), precipitation of driest month (mm) and precipitation of coldest quarter (mm). We used these variables to carry out ENMs in current and in Last Glacial Maximum (LGM) scenarios.

Maxent Models
We built ENMs using the maximum entropy modelling approach, Maxent ver. 3.4.1 (http://biodiversityinformatics.amnh.org/open_source/maxent/) [50]. This algorithm usually provides an excellent predictive approach compared with other modelling methods and is especially suited to deal with scarce presence-only data [51][52][53][54]. Since this technique relies on a generative rather than a discriminative approach, it performs well when the amount of training data is limited. Because our study area was small, we trained models using data from the entire European territory to account for a more comprehensive niche representation.
We carried out an ENM for each of the taxa identified in the cave by charcoal analysis using three steps: (1) we ran current models using presence records and the bioclimatic variables selected as described above; (2) we carried out paleoclimate models projecting current distribution in the LGM scenarios; and (3) we again ran paleoclimate models projecting current distribution in the LGM scenarios, but only for some taxa selected in Step 2. Unlike the previous model, here, we added the potential distribution map of the best-represented taxon in the study area using it as an "bioclimatic variables".
From the Maxent's setting panel, we selected the following options: random seed; remove duplicate presence records; write plot data; regularisation multiplier (fixed at 1) [50]; 1000 maximum iterations, 10,000 background points, cloglog format (this output appears to be most appropriate for estimating the probability of presence) [54,55]; and, finally, we used a 20-replicate effect with cross-validation run type. This run type makes it possible to replicate n-sample sets removing one locality at a time [22,56]. We fixed to 1 the default regularisation value as this is based on the different performances recorded across a range of taxonomic groups [57]. The remaining model values were set to default values [22,58].
For each species, the average final map had a cloglog output format with suitability values from 0 (unsuitable habitat) to 1 (suitable habitat). The 10th percentile (the value above which the model classifies correctly 90% of the training locations) was selected as the threshold value for defining the species' presence. This is a conservative value commonly applied to ecological niche modelling studies, particularly those relying on datasets collected over a long time by different observers and methods [22,56,58,59]. We used this threshold to reclassify our model into binary presence/absence maps.
We generated the paleoclimate models using the same climatic variables above described. These models were trained with all occurrences collected, and projected to Europe in the LGM (23,000-18,000 years BP). We developed the paleoclimate models using the most used LGM scenarios: CCSM4 [19,60,61]. Projecting ENMs to regions other than those where models were calibrated, or to past or future times is a common approach to make inferences such as forecasting the spreading of alien organisms, providing paleo-reconstructions or predicting distributional patterns in future epochs [62,63]. In order to project to new area models calibrated elsewhere, whether in the current epoch or in the LGM, variables in the projection area must meet a condition of environmental similarity to the environmental data used for training the model. Therefore, we first ascertained that this condition occurred by inspecting the Multivariate Environmental Similarity Surfaces (MESS) generated by Maxent [64,65].

Model Validation
We tested the predictive performance of models using the receiver operating characteristics, for which we analysed the Area Under Curve (AUC) [66], and the minimum difference between training and testing AUC data (AUC diff ) [67]. These model evaluation statistics range between 0 and 1. Excellent model performances are expressed respectively by AUC values close to 1 and AUC diff close to 0.

Resolving Taxonomic Ambiguity by Means of ENM
In this study, several charcoal fragments were identified as "Pinus sylvestris type" (see Results). Materials classified as Pinus sylvestris type might correspond to Pinus sylvestris, Pinus nigra and Pinus mugo/uncinata, but species identification by charcoal analysis is somewhat problematic due to the absence of specific diagnostic key features.
Therefore, we first generated ENMs for the three species mentioned above, along with Q. pubescens, because this was the most common taxon identified at species level at the site (see Results). In a second step, we used the Q. pubescens potential distribution as an environmental layer for the above-mentioned Pinus models, assuming that the species found in the same charcoal assemblage should have similar ecological requirements since, these prehistoric communities collected wood fuel in the immediate surroundings of human settlements [41].

Charcoal Analysis
We identified 1677 charcoal fragments and established the occurrence of 12 taxa ( Figure 2; Supplementary Materials, Figure S6).
Forests 2020, 11, x FOR PEER REVIEW 6 of 18 environmental data used for training the model. Therefore, we first ascertained that this condition occurred by inspecting the Multivariate Environmental Similarity Surfaces (MESS) generated by Maxent [64,65].

Model Validation
We tested the predictive performance of models using the receiver operating characteristics, for which we analysed the Area Under Curve (AUC) [66], and the minimum difference between training and testing AUC data (AUCdiff) [67]. These model evaluation statistics range between 0 and 1. Excellent model performances are expressed respectively by AUC values close to 1 and AUCdiff close to 0.

Resolving Taxonomic Ambiguity by Means of ENM
In this study, several charcoal fragments were identified as "Pinus sylvestris type" (see Results). Materials classified as Pinus sylvestris type might correspond to Pinus sylvestris, Pinus nigra and Pinus mugo/uncinata, but species identification by charcoal analysis is somewhat problematic due to the absence of specific diagnostic key features.
Therefore, we first generated ENMs for the three species mentioned above, along with Q. pubescens, because this was the most common taxon identified at species level at the site (see Results). In a second step, we used the Q. pubescens potential distribution as an environmental layer for the above-mentioned Pinus models, assuming that the species found in the same charcoal assemblage should have similar ecological requirements since, these prehistoric communities collected wood fuel in the immediate surroundings of human settlements [41].

Charcoal Analysis
We identified 1677 charcoal fragments and established the occurrence of 12 taxa (Figure 2; Supplementary Materials, Figure S6).  Quercus deciduous type, namely Q. pubescens, was the best-represented taxon in almost all investigated layers, ranging from 61% (Layer 8) to 21% (Layer 11, Figure 2). Although deciduous Quercus can be hard to distinguish at the species level, the porous ring with only one to two rows of large vessels in the earlywood and their diameter never exceeding 250 µm, allowed us to identify Q. pubescens. Fraxinus ornus-oxycarpa was represented in all the SC layers, apart from Layer 11. Based on the comparison of autecological features, both species are conceivable since F. ornus could belong to the mesoxerophilous forest community while F. oxycarpa could be restricted to lowland and riparian areas forming the mesohygrophilous azonal forest together with Populus and/or Salix attested in Layers 9, 8F and 8E. Acer (excluded A. pseudoplatanus and A. platanoides) was present in all samples; the maximum value (28%) is attested in Layer 11 and the minimum (4%) in Layer 8E. Tilia, Carpinus and Sambucus were attested with low values (0.5-1.3%) especially in Layer 9. Coniferous wood was represented by Pinus sylvestris type and Juniperus. Pinus sylvestris type was present in all layers, with its maximum value in the Layer 11 (30%), while Juniperus was attested by a few wood charcoals in Layer 9. Rosaceae Maloideae were always attested in all the charcoal assemblages, their value ranging from 3% (Layer 10) to 13% (Layer 11). Prunus was present in SC Layers 11, 10C, 8G and 8E with a maximum value of 7.9%.

Ecological
Niche Models of Pinus spp. and Q. pubescens

Current Models
The analysis of single variable contribution showed that temperature seasonality, mean temperature of wettest quarter and precipitation of driest month were the main factors influencing the model performance for all the tree species. We found that the accumulated contribution of these three variables was 88%, 83%, 93% and 81% for P. mugo/uncinata, P. nigra, P. sylvestris and Q. pubescens, respectively. The current model predicted a high probability of potential distribution of P. mugo/uncinata especially in the Alps and other European mountain areas, whereas P. nigra was more likely to occur in central and southern Europe and in western Asia, P. sylvestris in central and northern Europe and in western Asia and Q. pubescens in central and southern Europe (Figures 3 and 4). Quercus deciduous type, namely Q. pubescens, was the best-represented taxon in almost all investigated layers, ranging from 61% (Layer 8) to 21% (Layer 11, Figure 2). Although deciduous Quercus can be hard to distinguish at the species level, the porous ring with only one to two rows of large vessels in the earlywood and their diameter never exceeding 250 μm, allowed us to identify Q. pubescens. Fraxinus ornus-oxycarpa was represented in all the SC layers, apart from Layer 11. Based on the comparison of autecological features, both species are conceivable since F. ornus could belong to the mesoxerophilous forest community while F. oxycarpa could be restricted to lowland and riparian areas forming the mesohygrophilous azonal forest together with Populus and/or Salix attested in Layers 9, 8F and 8E. Acer (excluded A. pseudoplatanus and A. platanoides) was present in all samples; the maximum value (28%) is attested in Layer 11 and the minimum (4%) in Layer 8E. Tilia, Carpinus and Sambucus were attested with low values (0.5-1.3%) especially in Layer 9. Coniferous wood was represented by Pinus sylvestris type and Juniperus. Pinus sylvestris type was present in all layers, with its maximum value in the Layer 11 (30%), while Juniperus was attested by a few wood charcoals in Layer 9. Rosaceae Maloideae were always attested in all the charcoal assemblages, their value ranging from 3% (Layer 10) to 13% (Layer 11). Prunus was present in SC Layers 11, 10C, 8G and 8E with a maximum value of 7.9%.

Current Models
The analysis of single variable contribution showed that temperature seasonality, mean temperature of wettest quarter and precipitation of driest month were the main factors influencing the model performance for all the tree species. We found that the accumulated contribution of these three variables was 88%, 83%, 93% and 81% for P. mugo/uncinata, P. nigra, P. sylvestris and Q. pubescens, respectively. The current model predicted a high probability of potential distribution of P. mugo/uncinata especially in the Alps and other European mountain areas, whereas P. nigra was more likely to occur in central and southern Europe and in western Asia, P. sylvestris in central and northern Europe and in western Asia and Q. pubescens in central and southern Europe (Figures 3 and 4).

Last Glacial Maximum Projection Models
The LGM model predicted a high probability of potential distribution of P. mugo/uncinata especially in central and western Europe, whereas P. nigra was more likely to occur in southern Europe and in western Asia, P. sylvestris in western, central and northern Europe and in western Asia and Q. pubescens in southern Europe (Figures 5 and 6).

Last Glacial Maximum Projection Models
The LGM model predicted a high probability of potential distribution of P. mugo/uncinata especially in central and western Europe, whereas P. nigra was more likely to occur in southern Europe and in western Asia, P. sylvestris in western, central and northern Europe and in western Asia and Q. pubescens in southern Europe (Figures 5 and 6).

Last Glacial Maximum Projection Models
The LGM model predicted a high probability of potential distribution of P. mugo/uncinata especially in central and western Europe, whereas P. nigra was more likely to occur in southern Europe and in western Asia, P. sylvestris in western, central and northern Europe and in western Asia and Q. pubescens in southern Europe (Figures 5 and 6).    Figures S7-S10). Maxent models for P. mugo/uncinata, P. nigra, P. sylvestris and Q. pubescens showed AUC of 0.968 ± 0.032, 0.908 ± 0.082, 0.816 ± 0.076 and 0.859 ± 0.073, respectively. AUCdiff mean and standard deviation's values for all the Maxent models were <0.1.

LGM Projection Models Adding the Q. Pubescens Distribution
The analysis of single variable contribution showed that temperature seasonality, mean temperature of wettest quarter and precipitation of driest month were the main factors influencing the model performance for P. mugo/uncinata and P. sylvestris. We found that the accumulated contribution of these three variables was 80% and 92% for P. mugo/uncinata and P. sylvestris, respectively. Instead, the potential distribution of Q. pubescens was the main layer influencing the model performance of P. nigra with a contribution of ca. 70% versus 30% and barely 8% of influence in P. sylvestris and P. mugo/uncinata model, respectively. The MESS analysis and the potential distributions of P. mugo/uncinata, P. nigra and P. sylvestris were not affected by the presence of the variable Q. pubescens' potential distribution. In fact, both the trend of MESS results and of the potential distributions of the three Pinus were similar to those obtained from the previous first models.

Pinus Type Sylvestris: Ecological Niche Models and Ecological Considerations
Our ENMs showed considerable performances in estimating the current and past distributions of P. mugo/uncinata, P. nigra, P. sylvestris and Q. pubescens in Europe, as also shown by model validation. AUC values such as the ones that we obtained (>0.8) are among the highest reported for  Figures S7-S10). Maxent models for P. mugo/uncinata, P. nigra, P. sylvestris and Q. pubescens showed AUC of 0.968 ± 0.032, 0.908 ± 0.082, 0.816 ± 0.076 and 0.859 ± 0.073, respectively. AUC diff mean and standard deviation's values for all the Maxent models were <0.1.

LGM Projection Models Adding the Q. pubescens Distribution
The analysis of single variable contribution showed that temperature seasonality, mean temperature of wettest quarter and precipitation of driest month were the main factors influencing the model performance for P. mugo/uncinata and P. sylvestris. We found that the accumulated contribution of these three variables was 80% and 92% for P. mugo/uncinata and P. sylvestris, respectively. Instead, the potential distribution of Q. pubescens was the main layer influencing the model performance of P. nigra with a contribution of ca. 70% versus 30% and barely 8% of influence in P. sylvestris and P. mugo/uncinata model, respectively. The MESS analysis and the potential distributions of P. mugo/uncinata, P. nigra and P. sylvestris were not affected by the presence of the variable Q. pubescens' potential distribution. In fact, both the trend of MESS results and of the potential distributions of the three Pinus were similar to those obtained from the previous first models.

Pinus Type Sylvestris: Ecological Niche Models and Ecological Considerations
Our ENMs showed considerable performances in estimating the current and past distributions of P. mugo/uncinata, P. nigra, P. sylvestris and Q. pubescens in Europe, as also shown by model validation.
AUC values such as the ones that we obtained (>0.8) are among the highest reported for published models (e.g., [44,53]) and documented a high predictive power of habitat suitability [68]. Our study was further supported by AUC diff values (e.g., [47]).
In agreement with Caudullo et al. [43], our current models for Europe matched well the observed distribution of P. mugo/uncinata, P. nigra P. sylvestris, while, with regard to Q. pubescens, the models identified a larger area than that shown in Caudullo et al. [43], in our case comprising the whole of the Iberian Peninsula. P. mugo/uncinata occurs in the mountains of Central and Eastern Europe, but is especially abundant in the subalpine belt of the Eastern Alps and the Carpathians [42]. Disjunct ranges occur in the lower mountains of the Jura and the Vosges, and at high altitudes in the Mediterranean and Balkan Mountains, such as the Apennines, the Albanian Alps and the Rila-Pirin-Rhodopes in Bulgaria [43]. Indeed, in Italy, P. mugo/uncinata belongs strictly to the subalpine belt, occurring above the forest treeline in the Alps and locally in the northern-central Apennines [69]. P. sylvestris ranges from Scotland, Ireland and Portugal in the west, east to eastern Siberia, south to the Caucasus Mountains and north to the Arctic Circle in Scandinavia [7] while, in Italy, it spreads to the Alps and occurs as a relic in the northern Apennines [69]. The P. nigra group is a widely distributed Mediterranean mountain conifer with a discontinuous range extending from North Africa (35 • N), through the northern Mediterranean, eastwards to the Black sea, finally in the western Mediterranean on the islands of Corsica and Sicily (both as P. nigra subsp. laricio). In southern Italy, P. nigra forests are today confined to a few relic carbonatic rocky mountains where they form open vegetation on the steep slopes between the grasslands and the broadleaved forest [27]. Q. pubescens ranges from the Atlantic coast of France to the shores of the Mediterranean Sea, and across peninsular Italy, the Balkan Peninsula and the Aegean regions, to the coasts of the Black Sea and most of Anatolia. Although our models showed that Q. pubescens might potentially occur in the entire Iberian Peninsula, the western extreme of its geographic range, to the best of our knowledge, pubescent oak lives only in northern Spain [43]. This is not surprising since ENMs do not take into account biotic interactions such as competition, nor they may include historical factors which might play a role in influencing actual distribution [70]. In Italy, according to our models, this species occurs on almost the entire territory [69].
The substantial matching between our models and the maps shown in Caudullo et al. [43] represent an encouraging piece of evidence supporting the reliability of the maps of past predicted distribution of P. mugo/uncinata, P. nigra, P. sylvestris and Q. pubescens in Italy. The past dynamics of Pinus species, and consequently their paleobiogeography, are poorly documented because pollen and charcoal studies do not allow confident species identification. Regarding the Iberian Peninsula, paleoecologists mention the occurrence of P. sylvestris and P. nigra probably because both species are currently present in the mountain areas of this region [71]. An integrated paleobotanic, genetic and modelling approach pointed at the existence of western Europe of potential glacial refugia of P. sylvestris up to 40 • N [3]. The few preceding modelling studies of P. sylvestris past potential distribution (e.g., [3,16,26]) showed contrasting results. Habitat suitability modelling of the boreal P. sylvestris [3] indicated a potentially wide LGM range in southern Central Europe and eastwards, but those models covered a smaller area than the one we predicted and indicated that potential glacial refugia of P. sylvestris were located between ca. 40 • N and 50 • N with a patchy geographical distribution. Svenning et al.'s Maxent models [26], instead, showed a potential distribution that extended over southern Italy more than ours. Such discrepancies may be due to the different modelling approaches or climatic variables used in such studies. Based on our models, we predicted that in our study area and more generally over southern Italy, P. nigra represents the most likely candidate. In fact, neither past potential distribution of P. mugo/uncinata nor that of P. sylvestris reached our study region in the LGM.
Besides, based on the assumption that the ecological requirements of P. nigra match closely those of Q. pubescens, we found that the potential distribution of the latter species contributed only 30% to the potential occurrence of P. sylvestris vs. 70% to that of P. nigra, further corroborating our findings.
Two subspecies of P. nigra are described for Italy [72], the subsp. nigra and the subsp. laricio (P. nigra and P. laricio sensu Pignatti [69]), respectively. The former has several scattered populations from 200 to 1500 m a.s.l. in the northeastern Alps, while small patches grow on calcareous slopes in the central-southern Apennines between 100 and 1350 m a.s.l. [73]. P. n. laricio mostly grows on siliceous soils, in the Sila region (Calabria) between 900 and 1600 m a.s.l. and in Mt. Etna (Sicily) between 1200 and 2000 m a.s.l. [69,73]. In our area, the dominance of carbonatic substrate is a strong argument against the occurrence of P. nigra subsp. laricio which is strictly linked to siliceous substrates.
Today, P. nigra forests are extremely rare in southern Italy. It is an early successional species, and pure self-replacing forests are constrained to few mountainous Mediterranean areas where they can be considered an edaphic climax limited to thin soils. More often this species is part of precursor or transitional associations towards deciduous broadleaved forests. However, P. nigra communities are still present, with a very small population, just ca. 80 km north to our study site on a subcoastal Mesozoic limestone ridge (Mts. Picentini, Vallone della Caccia), where these stands are part of the xerophilic open vegetation that occurs on the steep slopes as transitional vegetation between the grasslands and the broadleaved forest vegetation [27]. The Pinus nigra-Q. pubescens association which characterises the Pleniglacial could indicate at the local scale a warm-cool bioclimate sensu Finlayson et al. [74].

Vegetation Cover of the Pleniglacial
The vegetation, inferred at three temporally distinct episodes, was located close to the coastal caves, in a wide coastal plain including currently submerged sectors. Indeed, between 30 and 19 ka cal BP the coastline would have been about 8 km away from the present one, due to sea level lowering [75].
In both the oldest layers, namely at~36 ka cal BP and at~29 ka cal BP vegetation cover at the site can be envisaged as woodland with coexisting P. nigra and deciduous mesoxerophilous trees with Q. pubescens and Acer, which summed together, account for the 50% of the charcoal assemblage. Interestingly, the youngest layer falls in the cold arid period detected in the Salerno Gulf,~80 km to the north, spanning between 34 and 27 ka cal BP [76]. In this phase, despite the occurrence of the lowest values of annual temperatures [TANN (7 • C)], July temperatures [TJUL (20 • C)] and annual precipitation [PANN (400 mm year −1 )], it seems that in Camerota Bay the amount of precipitation may still support the development of forest vegetation.
At~19 ka cal BP, vegetation is still dominated by Q. pubescens and Fraxinus (up to~70%), but the number of pine charcoals appears significantly reduced. The age of this layer falls in the LGM, as defined by EPILOG [77] during which cold arid conditions are inferred for the whole Mediterranean area [78]. However, at the regional scale, this period coincides with the end of a relatively warm-humid phase (25.5 to 18.5 ka cal BP [76]) also recorded at Monticchio (~100 km NE) between 25 and 20 ka BP [79] which might explain the high amount of Q. pubescens and Fraxinus. This may well have disadvantaged the light-demanding P. nigra [80] which probably was restricted to the rocky slopes with thin soils. Overall, during the Pleniglacial, our data suggest a stable forest cover with temporal phases of dominance variation between Pinus and deciduous trees.
The presence of a continuous forest cover in the Camerota Bay is also confirmed by the presence of woodland mammals such as Gliridae, Muridae and Clethrionomys [81,82], while species related to open and steep slope environments are very rare [83]. The broad diversity of micromammal assemblages throughout the Pleniglacial indicates the contemporaneous presence of probably both forested and open environments that periodically expanded and retreated, offering habitats for many taxa [81].
In a wider spatiotemporal perspective, our data are consistent with the pollen record from the Gulf of Salerno [25], where arboreal pollen values are almost always over 30% and deciduous oaks were stable at around 10% (Pinus excluded) during the entire Last Glacial Period (LGP). Our data, combined with the evidence of the almost treeless grass-dominated landscape inland and the much higher elevation sites of the Italian peninsula [79,84,85], suggest that most of the arboreal pollen content of these pollen spectra should be ascribed to the coastal sector.
The only Italian pollen record showing a vegetation cover comparable to that shown by our data concerns the north Tyrrhenian coast where until~28 ka cal BP deciduous Quercus and Pinus with Abies alba combined to form a dense forest [86].
Roughly at the same latitude as Camerota Bay, the southern Adriatic coast of Italy, at a direct distance of~200 km, at~28.5 ka cal BP (24,410 ± 320 BP), was characterised by evergreen vegetation with P. halepensis, Juniperus and Pistacia [87], testifying warmer and drier conditions in this eastern coastal sector.
Such evidence highlights the remarkable vegetation cover heterogeneity in the Italian peninsula both in terms of latitude and longitude during the late Pleistocene, probably due also to a west-east precipitation gradient due to the longitudinal split by the Apennines which intercept and block westerly humid air masses. Our data also highlight the role of local topography on climate: indeed, here, the proximity of the Bulgheria massif probably played a preeminent role in trapping the clouds of the western weather systems spreading from the Tyrrhenian sea.
Interestingly, a good match with our data has also been found along the west Atlantic coast of Portugal, where wood remains of Pinus nigra/sylvestris, deciduous Quercus and Fraxinus dated between 34 to 20 ka BP suggest a very similar forest cover [88]; also, in this case, it seems that oceanic humid air masses played a major role.
In our reconstruction of the forest structure, the spatial position of the pioneer and shade-intolerant black pine with low-density canopy overtops a matrix of mesothermophilous winter deciduous broadleaved species. This structure is today detectable in P. nigra and P. leucodermis relict forests of southern Italy. Indeed, Rauh's canopy architecture model of these pines exhibits in the mature-old ontogenetic stage a tabular canopy with the green crown restricted to the upper third of the stem and large branches similar to the stem [89][90][91]. This feature favours direct light transmission and thus permits the establishment of relative shade-tolerating trees, grasses and shrubs. We should also speculate on the engineering capability of the pine canopies, which positively modifies the microclimate by affecting near-ground temperatures, soil moisture and wind speed. In such circumstances, these trees, by acting as nurse plants, should facilitate both establishment and survival of the broadleaved tree species.
In Europe, today, mixed P. nigra deciduous forests can be seen in the Eastern Alps where it occurs between 200 to 1200 m a.s.l. together with O. carpinifolia and F. ornus; in southern Bulgaria P. nigra grows at low altitudes mixed with Q. frainetto and Q. pubescens; further, in southern France mixed forests with Q. pubescens and P. nigra can also be found [92].
Forest vegetation probably that is very similar to the one suggested by our charcoal assemblage is the Beynam forest, located in the Kuyrukçu Mountains, not far from Ankara, in central Turkey [93]. According to Emberger's climate classification, this region is characterised by semi-arid and very cold Mediterranean climate [94]. Thus, the summer temperatures and the coldness of winter are the main factors characterising the climate. This forest is entirely surrounded by steppe and it is dominated by P. nigra subsp. pallasiana, Q. pubescens and Juniperus oxycedrus. Interestingly, in this region holly oak is not a tree but a shrub, and characterises the lowest layer of the forest cover. This present landscape seems to be the most similar vegetation to the Pleniglacial cover in Camerota Bay.

Vegetation Cover of the Late Glacial
During the last glacial-interglacial transition (Late Glacial sensu Orombelli et al. [95]), Q. pubescens is still the dominant species, accompanied by other deciduous taxa; P. nigra is also still present. It should be pointed out that around~16 ka cal BP new temperate warm taxa such as Tilia, Carpinus and Sambucus appear. This evidence can be interpreted as a clear consequence of the temperature rise and increased soil moisture supply following the beginning of deglaciation. The expansion of Tilia recorded in several Late Glacial pollen sequences of southern [96] and central Italy [86,[97][98][99][100][101] might reflect an early expansion of that tree from nearby refugia. This hypothesis agrees with the previous ones of shelter areas for Tilia in the lower thalwegs of the Mediterranean coastal rivers [102]. The presence of Populus agrees with the riparian forest evidence inferred by paleo-shell analysis carried out in the cave [103].
In Layers 8G, 8F, 8E (13.9-13.4 ka cal BP) P. nigra declines when Fraxinus increases. At this time, in the Gulf of Salerno, a rapid climatic change, culminating at 13.8 ka cal BP, marks the Bølling-Allerød chronozone characterised by the increase in atmospheric temperatures, especially the summer values [TJUL (24 • C)], and precipitation [PANN 900 mm)] [76]. Additionally, stable isotopes of land snail shells from the SC layers dated to 14-13.4 ka cal BP suggest moisture conditions quite similar to those of the present day must have occurred. These climatic conditions could have reduced the competitive advantage of P. nigra over the broadleaf species, especially with respect to pioneer species such as F. ornus. In this respect, it is interesting to note that, in medium and high belt wooded landscapes of Iberian Peninsula, cryophilous pines, that were the main species during the Late Glacial, declined in favour of broadleaved species following Holocene climate amelioration [104].
To sum up, our data suggest that the Cilento coast acted as a refugium for temperate deciduous tree species and P. nigra, confirming the coastal environment as a potential reservoir of biodiversity [105,106] and contrasting with the mid-altitude theory based on the assumption that precipitation in this region would have been higher than on the plains [1,2,96,[107][108][109].
Our results remark the usefulness of combining different approaches to explore biogeographic questions about past and current forest distribution-a fundamental step to informing forest management and conservation.

Conclusions
Our results give a very clear picture of bioclimatic conditions in the surroundings of the Camerota caves during the LGP and as late as the Lateglacial. They indicate that the climatic conditions were always able to sustain forest cover. The data show the presence of mesothermophilous forest during the LGP (from~36 ka cal BP), proving that this area played an important role as a reservoir of woodland biodiversity in which Q. pubescens was the most abundant component, followed by a wide variety of deciduous trees and mountain pines, most likely P. nigra. ENM projections provided a useful complement to our paleoecological studies, refining charcoal evidence and offering a less subjective picture of past geographic distributions of Pinus species in the LGM. Ours is the first study that, by using paleoclimate model and charcoal analysis, suggests the potential presence of a glacial refugium of P. nigra on the Tyrrhenian coastal sector of southern Italy, confirming our initial hypothesis. Finally, this work provides punctual evidence of the crucial role of coastal areas as reservoirs for temperate tree taxa in the Mediterranean basin.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/6/673/s1, Figure S1: Serratura Cave entrance; Figure S2: Presence records of P. mugo/uncinata after spatial autocorrelation analysis; Figure S3: Presence records of P. nigra after spatial autocorrelation analysis; Figure S4: Presence records of P. sylvestris after spatial autocorrelation analysis; Figure S5: Presence records of Q. pubescens after spatial autocorrelation analysis; Figure S6: Microscopic anatomical features of main identified taxa; Figure S7: Multivariate Environmental Similarity Surfaces (MESS) maps for P. mugo/uncinata obtained with CCSM4 models for LGM scenarios; Figure S8: Multivariate Environmental Similarity Surfaces (MESS) maps for P. nigra obtained with CCSM4 models for LGM scenarios; Figure S9: Multivariate Environmental Similarity Surfaces (MESS) maps for P. sylvestris obtained with CCSM4 models for LGM scenarios; Figure S10: Multivariate Environmental Similarity Surfaces (MESS) maps for Q. pubescens obtained with CCSM4 models for LGM scenarios.