Latest Holocene paleoenvironmental and paleoclimate reconstruction from an alpine bog in the Western Mediterranean region: The Borreguil de los Lavaderos de la Reina record (Sierra Nevada)

Several organic and inorganic geochemical analyses have been carried out in the sedimentary record of the Borreguil de los Lavaderos de la Reina (BdlR-03), an alpine peat bog located on the north face of the Sierra Nevada (southern Iberian Peninsula). This study permitted a high-resolution reconstruction of paleoenvironmental evolution for the last ~2700 cal yr BP in the highest mountain range of southern Iberian Peninsula. An overall trend towards a climatic aridification and a reduction of aquatic environments is observed in this record. Insolation and long-term positive North Atlantic Oscillation trends were the most important factors controlling this aridification, forcing regional and local environmental changes. Four phases are differentiated within the paleoenvironmental evolution of BdlR-03: (1) a pre-bog environment from 2700 to 2600 cal yr BP with important siliciclastic sedimentation and low organic content; (2) a bog environment with important presence of terrestrial vascular plant, water availability and maximum humid conditions between ~2600 and ~1870 cal yr BP, coinciding with the Iberian-Roman Humid Period; (3) a subsequent drier bog environment, between ~1870 and ~300 cal yr BP during the Dark Ages, the Medieval Climate Anomaly and the first stages of the Little Ice Age, characterized by lower productivity that affected the development of terrestrial vascular plants and aquatic environments; and (4) a wetland environment with an increase of aquatic algae, development of ephemeral pools and unstable climate conditions between ~300 cal yr BP and the present, coinciding with the final stages of the Little Ice Age and the Modern Global Warming. This recent environmental pattern, which is opposite to the general aridification trend in the western Mediterranean, is likely explained by the gradual melting of the perennial ice and/or snow packs at higher elevations in the last centuries. Aeolian inputs would have continuously contributed nutrients to these nutrient-impoverished alpine environments. High dust inputs are especially noticed during the last ~200 years, which can be explained by human-induced enhanced aridification and the development of the commercial agriculture in some North African regions. In addition, the environmental signal in the last century seems to be significantly affected by human activities.


Introduction
Climate change is one of the greatest challenges for modern societies and natural environments since it represents a serious global threat for food provisioning, security, economy, biodiversity and ecosystems (IPCC, 2014). The Intergovernmental Panel on Climate Change shows that the last decades have been the warmest of the Common Era in the Northern Hemisphere based on the increase of ~ 0.85 • C in the globally averaged temperature between 1880 and 2012 (IPCC, 2014). Future scenarios for the Northern Hemisphere are not optimistic, since significant increases in mean temperatures are expected (Pauling et al., 2006;IPCC, 2014). Furthermore, projections for future precipitation patterns and their impact on specific regions are more uncertain (IPCC, 2014). These abrupt climatic changes will especially affect environmentallyvulnerable areas such as the Mediterranean region, which is one of the most prominent "Hot-Spots" in future climate change projections (Giorgi and Lionello, 2008) due to, among other factors, the intensification of seasonal droughts (Smith et al., 2016). Additionally, a recent study based on estimated changes in sea surface temperature (SST) showed a mean warming trend of 0.041 ± 0.006 • C/year over the whole Mediterranean Sea from 1982 to 2018 (Pisano et al., 2020) coupled with changes in the amplitude of the Mediterranean seasonal signal for this period (increase and decrease of summer and winter mean temperatures, respectively) (Pisano et al., 2020). These effects along with a precipitation decrease and an increase in the frequency and severity of extreme events will contribute to the enhancement of summer droughts, boosting desertification patterns, and affecting the soil water availability and thus, the sustainability of agriculture (Carsan et al., 2014;Smith et al., 2016;Muñoz-Rojas et al., 2017;Cox et al., 2018;Abd-Elmabod et al., 2020). Nevertheless, the study of current and past climate change fails to establish a general pattern for the whole Mediterranean region  because the complex topography, coastline, and vegetation cover of the Mediterranean region modulate the regional climate signal at smaller spatial scales (e.g., Lionello et al., 2006), making it difficult to generalize in future predictions.
From a temporal point of view, the western Mediterranean area has also registered large past climate and environmental shifts throughout the Holocene, being very sensitive to intense short-term climatic fluctuations . Therefore, paleoclimatic evidence from high resolution paleoenvironmental records can help understand the current warming and drying scenario in the western Mediterranean basin from a long-term perspective of natural climate variability (Ljungqvist et al., 2016). In this regard, environmentally-vulnerable alpine areas in the Mediterranean region such as the alpine wetlands from the Sierra Nevada range (southern Iberian Peninsula) are of particular interest since their sedimentary records register fast ecological and biogeochemical responses to different environmental forcings (Catalan et al., 2013;García-Alix et al., 2017).
In this work we perform a high-resolution (sub-centennial scale) multi-proxy analysis in a peat record covering the last ~2700 cal yr BP from the Borreguil de los Lavaderos de la Reina peat bog with the goal of assessing the environmental impacts caused by past climate change or human impact on the Sierra Nevada alpine areas. The study of the sedimentary record of BdlR-03, located in a glacial cirque on the north face of the Sierra Nevada, improves the time resolution of previous studies in this area, allowing for a better understanding of the environmental and climatic fluctuations during the latest Holocene in alpine Mediterranean mountains.

Sierra Nevada: Climate and vegetation
Sierra Nevada is the highest mountain range in southwestern Europe. This mountain range has an east -west orientation and includes the highest peak of the Iberian Peninsula (Mulhacén -3479 masl). Metamorphic rocks, primarily mica schists, are the main lithologies characterizing the high-elevation outcrops in Sierra Nevada. Soils are scarcely developed in those alpine rocky areas where the geomorphology is strongly influenced by past glacial activity (Castillo Martín, 2009). During the Late Pleistocene valley glaciers reached down to 2300-2400 masl on the north faces and 2400-2500 masl on the south faces of the Sierra Nevada (Gómez-Ortiz et al., 2015). The subsequent melting of cirque glaciers during the latest Pleistocene-earliest Holocene resulted in the formation of small lakes and peat bogs (locally known as borreguiles) in the formerly glaciated basins. These wetlands are mostly located between 2450 and 3230 masl (Castillo Martín, 2009;Palma et al., 2017) within the total protection perimeter of the Sierra Nevada National Park. Although large glaciers in Sierra Nevada disappeared in the earliest Holocene, small glaciers with limited extension reappeared on the north faces of the highest peaks in certain intense cold periods such as during the Little Ice Age (LIA), disappearing in the first part of the 20th Century (Gómez-Ortiz, 1987;Gómez-Ortiz et al., 2009;Palma et al., 2017;Oliva et al., 2018).
The climate in the Sierra Nevada area is Mediterranean, characterized by dry and warm summers and humid and mild winters. The mean annual temperature at 2500 masl is ~4.5 • C, while the mean annual rainfall is ~700-750 mm/year, mainly concentrated as snow between October and April (Spanish National Weather Agency -AEMet Open Data, 2020). The North Atlantic Oscillation (NAO)  largely influences precipitation in the study area, since the interaction between the anticyclone of the Azores and the low pressures of Iceland determine the latitudinal position of the North Atlantic storms entering in the Mediterranean area . Rainfall and thermal gradients determine the occurrence of four natural vegetation belts in Sierra Nevada situated in an altitudinal gradient (Valle, 2003). The cryoromediterranean vegetation belt is located above 2800 masl and includes grasslands and plants with basal rosettes. The oromediterranean belt, located between 2800 and 1900 masl, is principally characterized by Pinus sylvestris and P. nigra, Juniperus spp., Genista and Cytisus shrubs. Between 1900 and 1400 masl, the supramediterranean belt is mainly dominated by deciduous and evergreen Quercus, with Acer, Fraxinus, Sorbus and various shrubs such as Adenocarpus, Helleborus and Daphne. Finally, between 1400 and 600 masl, the mesomediterranean belt includes sclerophyllous shrublands and evergreen Quercus forests (El Aallali et al., 1998;Valle, 2003). However, in last centuries several anthropogenic activities like Pinus reforestation (Valbuena-Carabaña et al., 2010) or Olea cultivation (Anderson et al., 2011;Jiménez-Moreno and Anderson, 2012;Ramos-Román et al., 2019) have altered the natural spatial distributions of some plant species.

Borreguil de los Lavaderos de la Reina (BdlR)
BdlR is a small peat bog located in the depression of a glacial cirque on the north face of Sierra Nevada (37 • 07 ′ 48"N, 3 • 16 ′ 21"W) at 2421 masl ( Fig. 1). There is a set of waterfalls upstream of the site at 2730 masl called the Lavaderos de la Reina that give name and provide water to the BdlR peat bog. The catchment basin bedrock consists of lowdegree metamorphic rocks, mainly mica schists with feldspar and amphibole, of the Nevado Filábride complex  with scarce soil development and typical vegetation of the oromediterranean (mostly Juniperus shrubs) and cryoromediterranean tundralike belt (Valle, 2003).

Core sampling, lithology and magnetic susceptibility
Three sediment cores (BdlR-01, BdlR-02 and BdlR-03) were recovered in September 2006 from the BdlR wetland with a Livingstone piston corer. They were wrapped with tin foil and cling film and initially stored at 4 • C at the Northern Arizona University and subsequently at the University of Granada. The longest core, BdlR-03 (97 cm), was selected for this multi-proxy study.
Core BdlR-03 was split longitudinally and described in the laboratory with respect to lithology and colour. Magnetic susceptibility (MS) was also measured with a Bartington MS2E meter in SI units (x 10 − 5 ) every 0.5 cm, with a measuring time period of 10 s, for the totality of the record (Fig. 2). and a 239+240 Pu activity profile in the topmost 34.5 cm. Radiocarbon dates were converted into calendar year before present (cal yr BP) using the IntCal20.14c curve (Reimer et al., 2020) (Table 1). 239+240 Pu activity profile was obtained from ten sediment samples taken every ~3.5 cm. Samples were prepared following chemical procedures from Ketterer et al. (2004) and analysed using a Thermo X2 quadrupole ICPMS system located at Northern Arizona University. The detection limit for 239+240 Pu, based on 0.5 g samples, was 0.1 Bq/kg (Table S1). A bayesian age-depth model for BdlR-03 was constructed using 14 C and 239+240 Pu data and the R modeling package "BACON (v2.4.3-August 2020)" (Blaauw and Christen, 2011). The better result, in terms of model smoothness without over-parameterization, was obtained using 24 sections of 4 cm and a prior accumulation rate of 25 yr/cm with a gamma distribution of 0.75. The prior memory and shape were 0.15 and 25, respectively, to accommodate a shift in sediment accumulation rates in the upper portion of the core (Fig. 2).

Inorganic geochemistry
The inorganic geochemical composition of the sediments was obtained using continuous X-Ray Fluorescence (XRF) analysis by means of an Avaatech XRF core scanner, equipped with an Oxford Rhodium X-ray source (operating voltage range between 4 and 50 kV), at the CORELAB laboratory at the University of Barcelona. Measurements were taken using a prism flushed with helium in contact with the surface of the sediment at a resolution step of 0.5 cm and under two different conditions depending on the atomic weights of the elements: (a) 10 s count time, 700 μA of X-ray current, and 10 kV of X-ray voltage for elements with atomic weights between Al and Fe (Al, Si, P, S, Cl, Ar, K, Ca, Ti, V, Cr, Mn and Fe); and (b) 35 s count time, 1300 μA of X-ray current, and 30 kV of X-ray voltage for elements with atomic weights between Ni and Pb (Ni, Cu, Zn, Ga, Br, Rb, Sr, Y, Zr, Nb, Mo, Rh, Ag, Au, Hg and Pb). Elemental variations from XRF core scanners are expressed as counts per second (CPS = Peak area/Count time) since they depend on the excitation conditions used in the measurements (Frigola et al., 2015).

Organic geochemistry
A total of 119 samples with a sampling interval of 0.5 cm were taken from the BdlR-03 core for CHNS analyses. One cubic centimeter of each sample was decalcified overnight with HCl 1 N, and rinsed several times with MilliQ water until a neutral pH. Two aliquots were obtained from each sample: one for organic elemental analyses and another for carbon and nitrogen isotope analyses. About 5 mg of the first aliquot was encapsulated in tin capsules and analysed by means of a CHNS Elemental Analyzer Thermo Scientific Flash 2000 at the Centro de Instrumentación Científica of the University of Granada (CIC-UGR). Helium was used as gas carrier and the flash combustion was produced at 1000 • C. The equipment was calibrated every day using a certified sulfanilamide standard and the precision of the measurements was better than ±0.1%. Atomic C/N ratio was obtained from the division of total organic carbon (%TOC) and total nitrogen (%TN) normalized by their respective atomic weights. The methodology is explained thoroughly in García-Alix et al. (2018).
The other decalcified aliquot was also encapsulated tin capsules for carbon and nitrogen stable isotope analyses (δ 13 C and δ 15 N). The sample weight ranged from 1.8 to 5.4 mg, depending on the organic content of each sample (previously estimated with the elemental analyses). Samples were measured in duplicate by means of isotope-ratio mass spectrometry (IRMS-Thermo-Fisher DeltaV advantage) with a coupled elemental analyzer (Thermo-Fisher EA IsoLink) and helium as gas carrier at the CIC-UGR. Carbon and nitrogen isotope analyses were performed simultaneously at the same run. The isotopic measurements were calibrated using internal and international standards. An international standard Cafeine-IAEA 600 (δ 15 N +1.00‰ and δ 13 C − 27.771‰) along with two certified standards Sorghum -OAS-(δ 15 N +1.58‰ and δ 13 C − 13.68‰) and EBD23 (δ 15 N +9.94‰ and δ 13 C − 18.6‰) were used for data correction. One certified standard Wheat -OAS-(δ 15 N +2.85 and δ 13 C − 27.21‰) and another internal standard PEAT (δ 15 N +4.41‰ and δ 13 C − 27.93‰) were measured as internal control. Isotopic results are expressed in δ notation, using the standard VPDB (δ 13 C) and AIR (δ 15 N). After correction for mass spectrometer daily drift, the calculated precision was better than ±0.15‰ for δ 13 C and δ 15 N.
Lipid biomarkers were extracted from a total of 47 samples taken every ~2 cm from the BdlR-03 sedimentary record. The initial dry sample weight ranged from 40 to 260 mg. The total lipid extract (TLE) was obtained from freeze-dried samples by means of sonication (20 min) and temperature (38 • C for 1 h) using a DCM:MeOH (3:1) solution in closed glass vials. The neutral fraction of this TLE was separated by means of aminopropyl-silica gel chromatography (~4-5 cm aminopropyl-silica gel in a pre-combusted glass pasteur pipette) and a solution of 1:1 DCM:isopropanol (~4 ml). Later, the aliphatic hydrocarbon fraction with the n-alkanes was eluted from the neutral fraction using ~4 ml of hexane trough a 230-400 mesh / 35-70 μm silica-gel chromatographic column (~4-5 cm silica-gel in a pre-combusted glass pasteur pipette). This fraction was analysed at the BECS laboratory (University of Glasgow, UK) by means of gas chromatography coupled to a mass spectrometer (GC--MS Shimadzu QP2010-Plus Mass Spectrometer interfaced with a Shimadzu 2010 GC) in order to identify the different compound in some pilot samples and to a flame ionization detector (GC-FID Shimadzu 2010) in order to quantify the n-alkanes following Moossen et al. (2013) settings. A BP1 (SGE Analytical Science) column (60 m, diameter: 0.25 mm, film thickness: 0.25 μm; coating: 100% dimethyl-polysiloxane) and splitless injection mode were used. The initial GC oven temperature was set at 60 • C for two minutes. Afterwards, the temperature increased up to 120 • C (at 30 • C min − 1 ) and subsequently to 350 • C (at 3 • C min − 1 ) for 20 min. The carrier gas used was hydrogen for GC-FID analyses and helium for GC--MS analyses. An external standard-mix, consisting of C 16 , C 18 , C 19 , C 20 , C 23 , C 25 , C 26 , C 28 , C 30 , C 32 and C 37 n-alkanes, was measured every five samples in order to correctly identify the retention times of the different n-alkane peaks in  Table 1 Radiometric ages from BdlR-03. All ages were calibrated using IntCal20.14c curve (Reimer et al., 2020)  the GC-FID. The n-alkanes were quantified by comparing their peak areas with the peak areas of the C 26 n-alkane (concentration 13.3 μg/ml) of the standard injection. The analytical error (standard deviation) of the n-alkanes integration, obtained from duplicate analyses of 15 n-alkane containing fractions, was better than 5%. An accurate n-alkane quantification would have needed an internal standard in each sample; thus, n-alkane ratios will be used to summarize the n-alkane distribution, preventing from any quantification bias between the different samples. The ACL is the weighted average chain length of the odd chain n-alkanes (Poynter and Eglinton, 1990). We will use the ACL 25-33 that considers the odd chain lengths ranging from 25 to 33 carbon atoms (Eq. 1). The CPI is the carbon preference index that compares the distribution of odd-over-even chain lengths (Bray and Evans, 1961). We will use the CPI 25-33 chain lengths (odd and even) ranging from 25 to 35 carbon atoms (Eq. 2). The P aq is the proportion of (sub)aquatic plants according to Ficken et al. (2000), which assesses the contribution of C 23 -C 25 and C 29 -C 31 n-alkanes in one sample (Eq. 3).

Lithology and magnetic susceptibility
The sedimentary record of the BdlR-03 core is mainly composed of peat (Fig. 2). The bottom of the sedimentary record, specifically between 97 and 94 cm (2690-2565 cal yr BP), consists of dark grey coarse peat with silts, sands and small stones (2.5Y 4/1). The highest MS values are recorded in this bottommost part of the core. The rest of the sedimentary record is composed of peat and registered negative values of MS. Brown (10YR 4/3) and black (5YR 2.5/1) coarse peat dominated from 94 to 82 cm (2565-1965 cal yr BP) and from 82 to 64 cm (1965-1340 cal yr BP), respectively. The section between 64 and 29 cm (1340-190 cal yr BP) consisted of dark reddish-brown (5YR 2.5/2) coarse peat. The lowest MS values are registered between 29 and 22 cm (190-130 cal yr BP) in a black coarse peat. Finally, the section from 22 cm to the core top (from 130 to − 56 cal yr BP) consisted of dark reddish-brown coarse peat with very coarse roots in the uppermost 3 cm along with some sand and small rock fragments (from − 25 to − 56 cal yr BP).

Chronology and sedimentation rates
The bayesian age-depth model for the BdlR-03 record indicates that the sedimentary record covered approximately the last 2700 cal yr BP (Table 1; Fig. 2). 239+240 Pu values, above the detection limit, were registered in the uppermost ~16 cm of the record. Plutonium maximum activity was obtained between 5.5 and 3.6 cm (6.07 Bq/Kg) (Table S1). Sediment accumulation rates (SARs) were calculated from the linear interpolation of the AMS radiocarbon dates. SARs below 46 cm varied between 0.16 and 0.58 mm/yr, with the lowest SARs recorded between 46 and 60.5 cm (0.16 mm/yr). SARs above 46 cm fluctuated between 0.84 and 1.2 mm/yr, with the highest SARs occurring in the uppermost 27.3 cm (1.2 mm/yr).

Inorganic geochemistry
The XRF analyses allowed the identification of relative variations of certain chemical elements in the sedimentary record. Nonetheless, the high organic matter content in the BdlR-03 sediment record, which consist mainly in C, N, and H, can dilute the counts of the other elements; for this reason, we applied a normalization to the raw data (Löwemark et al., 2011) using a conservative element such as Al or Ti (Calvert and Pedersen, 1993;Martinez-Ruiz et al., 2015). In order to elucidate fluctuations of K and Ca in the lithogenic fraction we normalized them by the Ti, since the Al detection in the XRF corescanner might be biased due to the low intensity of this element (Tjallingii et al., 2007). Values of Ca/Ti ratio range from 6.2 to 0.6 with a mean value of 2.6 ± 0.9. The lowest and highest values are concentrated from 2680 to 2600 cal yr BP and from 60 cal yr BP to the present, respectively. On the other hand, K/Ti ratio range from 1.82 to 0.59 with a mean value of 1.01 ± 0.27. In this case, the lowest and highest values are concentrated from 820 to 420 cal yr BP and from 2680 to 2600 cal yr BP, respectively. Both ratios recorded a general increasing trend, being more pronounced in the Ca/Ti ratio (Fig. 3).

Organic geochemistry
TN ranges from 0.3% to 2.8% with a mean value of 2.2 ± 0.4%. The lowest and highest TN values are recorded from 2680 to 2600 cal yr BP and from 820 to 420 cal yr BP, respectively. TN data show a general increasing trend. TOC ranges from 6% to 50% with a mean value of 40 ± 8%. The lowest and highest TOC values are recorded from 2680 to 2600 cal yr BP and from 2500 to 1870 cal yr BP, respectively. TOC data show a general decreasing trend, although not as pronounced as in the case of C/N ratio. C/N ratio ranges from 30.7 to 14.4 with a mean value of 22.1 ± 4.1. The highest values are registered between 2600 and 1870 cal yr BP, being extremely high from 2500 to 2200 cal yr BP. The lowest values are registered between 420 cal yr BP and the present, being extremely low from 300 to 60 cal yr BP (Fig. 3).
Carbon and nitrogen stable isotopes (δ 13 C and δ 15 N) values range from − 28.7‰ to − 27.3‰ with an average of − 28.1 ± 0.3‰, and from − 0.6‰ to +5.4‰ with an average of +1.6 ± 1.1‰, respectively. Both isotopic records registered the highest values between 2680 and 2600 cal yr BP. Their lowest values were registered from 2500 to 1870 cal yr BP and from 60 cal yr BP to the present, respectively (Fig. 3).
The n-alkane record of BdlR-03 shows high concentrations of chain lengths ranging from 21 to 35 carbon atoms (C 21 to C 35 ) with a strong predominance of odd-over-even chain lengths (Fig. S1). Although C 19 , C 20 , C 36 , and C 37 chain lengths are also present in some samples, their low concentration prevents from quantification. The total n-alkane concentration in BdlR-03 record ranges from 9.5 to 125.9 μg gsed − 1 (microgram per gram of dry sediment). The dominant n-alkane through most part of the record is the C 31 (from ~2580 to 2300, 1780 to 440, from 310 to 220, and from 180 to 80 cal yr BP) with concentrations ranging from 1.7 from 34.6 μg gsed − 1 . The C 29 also predominates at some specific periods (from ~2220 to 1870, from 400 to 330, at 200, and from 60 cal yr BP to the present). The concentration of this n-alkane exhibits a range of variation from 2.0 to 38.5 μg gsed − 1 . The C 27 n-alkane only predominates at the beginning of the record (~2700 cal yr BP) with a concentration of 4.5 μg gsed − 1 . Its concentration throughout the record varies from 0.8 to 16.5 μg gsed − 1 . The ranges of variation of the other odd n-alkanes are: from 0.4 to 3.6 μg gsed − 1 for the C 21 n-alkane, from 0.6 to 6.3 μg gsed − 1 for the C 23 n-alkane, from 0.4 to 7.1 μg gsed − 1 for the C 25 n-alkane, from 0.6 to 17.3 μg gsed − 1 for the C 33 n-alkane, and from 0.2 to 3.4 μg gsed − 1 for the C 35 n-alkane (Fig. S1). This n-alkane distribution can be summarized throughout the record by means of the following n-alkane indices. Values of ACL 25-33 and CPI 25-33 ranges from 30.6 to 28.7 with an average of 29.7 ± 0.4, and from 7.3 to 3.5 with an average of 4.9 ± 0.9, respectively. The lowest values of ACL 25-33 and CPI 25-33 are recorded from 2500 to 1870 cal yr BP and from 420 to 300 cal yr BP, respectively. Both proxies registered the highest values between 300 and 60 cal yr BP. A general increasing trend is depicted by the ACL 25-33 and CPI 25-33 records. On the contrary, a general decreasing trend is observed in the P aq throughout the record. Values of P aq range from 0.38 to 0.09 with an average of 0.2 ± 0.05. This proxy recorded the highest values between 1300 and 900 cal yr BP, although registered a maximum peak at ~100 cal yr BP. The lowest values were recorded between 420 and 300 cal yr BP (Fig. 3).
Principal component analysis (PCA) was performed on the concentrations of the odd n-alkanes to assess the distribution and associations of the studied n-alkanes in this record. Previously, data were normalized using a Z-score function ((x-mean)/standard deviation). The first two PCs explain more than the 85% of the variance (PC1 = 70.1% and PC2 = 15.6%). The biplot diagram helps to differentiate two n-alkane groups: C 21 , C 23 , C 25 and C 27 , and C 29 , C 31 , C 33 , and C 35 (Fig. S1). The same grouping is obtained with the loadings of PC2. Since the PC1 clustered all the n-alkanes together, it does not provide specific information about the relationships among them (Fig. S1). Thus, we chose the PC2 scores to illustrate the variations between these two n-alkane groups during the latest Holocene, that are ultimately related to environmental fluctuations. Values of PC2 range from − 2.5 to 2.3 with an average of − 0.2 ± Fig. 3. Geochemical proxies studied in BdlR-03: (a) K/Ti ratio, (b) Ca/Ti ratio, (c) PC2, (d) CPI 25-33 according to Eq. 2, (e) ACL 25-33 according to Eq. 1, (f) P aq according to Eq. 3, (g) C/N ratio, (h) TOC (%), (i) TN (%), (j) δ 13 C (‰), and (k) δ 15 N (‰). Scales of δ 13 C and δ 15 N values are reversed. All proxies are represented with respect to age and zonation done using a Euclidean cluster analyses: Zone 0: Pre-bog environment (~2700-2600 cal yr BP), Zone 1: Bog environment (2600-1870 cal yr BP), Zone 2: Bog environment (Subzones A, B, C and D; 1870-300 cal yr BP), and Zone 3: Wetland environment (Subzones A and B; 300 cal yr BP-present).
1.1. The lowest and highest values are recorded from 2500 to 1870 cal yr BP and from 300 to 60 cal yr BP, respectively. A general increasing trend is depicted by the PC2 record (Fig. 3).

Proxy integration
A stratigraphically constrain cluster Euclidean analyses of inorganic and organic proxies has allowed the identification of four zones and six subzones in the sedimentary record of BdlR-03 (Figs. 3 and 4). These geochemically different periods will help understand the environmental change in the BdlR area during the last ~2700 cal yr BP (see below).  This zone is characterized by the highest mean values of K/Ti (1.8), CPI 25-33 (5.9 ± 0.3), δ 15 N (+3.9 ± 0.6‰) and δ 13 C (− 27.7 ± 0.3‰). Also, this section records the lowest TOC, TN, ACL 25-33 and Ca/Ti with averages of 9 ± 4%, 0.4 ± 0.1%, 29 ± 0.3 and 0.6 ± 0.1. The C 27 nalkane predominates in this part of the record.
Subzone 2B [From 1370 to 820 cal yr BP (66-53.4 cm depth)]. This subzone is characterized by C/N below 20 (average of 19.2 ± 0.7), moderate TOC (average of 39 ± 3%) and ACL 25-33 (average of 29.7 ± 0.1) values as well as the highest P aq values in the zone (average of 0.24 ± 0.03). Mean values of δ 13 C and δ 15 N slightly increase with respect to the previous subzone. Minor changes are observed in the inorganic ratios, except for the K/Ti ratio that registers a mean value of 1.22 ± 0.27 with a major peak at ~1270 cal yr BP. The dominant n-alkane is the C 31 .
Subzone 2C [From 820 to 420 cal yr BP (53.4-47.1 cm depth)]. This subzone is characterized by high values of ACL 25-33 , TOC and TN with averages of 30.2 ± 0.2, 43 ± 2% and 2.5 ± 0.1%, respectively. However, P aq values drop (mean, 0.14 ± 0.01) and the C/N ratio oscillate around 20 (20.5 ± 0.7). Mean values of δ 13 C and δ 15 N slightly decrease with respect to the previous subzone. On the other hand, inorganic proxies show a stable pattern during this interval. The dominant n-alkane is the C 31 .
Subzone 2D [From 420 to 300 cal yr BP (47.1-39.4 cm depth)]. This subzone registers lower C/N values than subzone 2C (average of 17.9 ± 1.1). During this period TOC, P aq , and Ca/Ti values show abrupt decreases reaching the lowest mean values for the Late Holocene. An opposite trend displayed ACL 25-33 , δ 13 C and δ 15 N values with relatively high averages. A maximum peak of ACL 25-33 and δ 15 N are reached at ~300 cal yr BP, in the boundary with the next stage of the bog. The dominant n-alkane is the C 29 .

Zone 3 [between 300 cal yr BP and the present (39.4-0 cm depth)]
This zone registers the lowest C/N (generally below 20) and P aq values of the record with averages of 17.9 ± 2 and 0.17 ± 0.08 respectively. On the other hand, this section is characterized by the highest ACL 25-33 , PC2 and Ca/Ti values with averages of 30 ± 0.4, 0.7 ± 1 and 4.1 ± 1, respectively. Furthermore, δ 15 N and δ 13 C mean values slightly increases with respect to the previous zone.
Subzone 3-A [From 300 to 60 cal yr BP (39.4-12.9 cm depth)]. This subzone registers the lowest averaged C/N values of the record (mean, 17.7 ± 0.9) and relatively low TOC content. Furthermore, TN content, δ 13 C and δ 15 N show abrupt increases and δ 13 C reaches the highest values for the entire sedimentary record (− 27.7 ± 0.2‰). In addition, the ACL 25-33 (mean, 30.2 ± 0.3) and PC2 (mean, 1.2 ± 0.8) register the highest values of the record at ~300 and ~ 200 cal yr BP, respectively, followed by gradual decreasing trends. A maximum peak of P aq was reached at ~100 cal yr BP. Fluctuating Ca/Ti and K/Ti ratios with overall increasing trends are observed. The dominant n-alkane is the C 31 .
Subzone 3-B [From 60 cal yr BP to the present (12.9-0 cm depth)]. This subzone is characterized by a progressive increase in the C/N ratio (average of 18.5 ± 3.3) after reaching the lowest value of the record at ~40 cal yr BP, and exhibiting values above 20 from ~0 cal yr BP until present. On the other hand, ACL 25-33 values were very low (around 29.5), like δ 13 C (− 28.1 ± 0.4‰) and δ 15 N (+0.4 ± 0.9‰). These are the lowest δ 15 N of the whole record. K/Ti and Ca/Ti ratios increase remarkably and Ca/Ti ratio reaches the highest values for the entire Late Holocene. The C 29 n-alkane predominates in this part of the record.

Proxy interpretation
Three major elements (K, Ca and Ti) will be used in our paleoenvironmental interpretations. The catchment basin of BdlR-03 is mainly composed of mica schists enriched in K-Feldspar, but depleted in Ca . Consequently, the K content, expressed as the K/Ti ratio, can be used to identify erosion/runoff in these alpine wetlands (Mesa-Fernández et al., 2018), whereas high concentrations of Ca might represent allochthonous aeolian inputs. Aeolian deposition in the alpine areas of the Sierra Nevada is mostly related to North African dust, which usually travels at elevations between 3000 and 5000 masl (Mladenov et al., 2010). This is the primary origin of aeolian dust input in the southeastern Iberian Peninsula, with 85% of the total aeolian content and the main source of allochthonous inputs at the alpine areas of the Sierra Nevada (Moreno et al., 2006;Morales-Baquero and Pérez-Martínez, 2016;Jiménez et al., 2018;Jiménez et al., 2019). Particle geochemistry of dust from North African can vary depending on the geological source. Concretely, western Sahara soil is a carbonate-rich geological depocentre which regularly supplies mineral particulates enriched in detrital carbonate (both calcite and dolomite) to the southeastern Iberian Peninsula (Moreno et al., 2006). In addition, an increasing algal content in nutrient-impoverished alpine Sierra Nevada wetlands could be strongly related to the input of aeolian nutrients (Pulido-Villena et al., 2006;Morales-Baquero and Pérez-Martínez, 2016). Thus, the Ca/Ti ratio would indicate variation in the North African dust input in the peat bog.
Bulk organic data can provide general information about the characteristics of the organic matter and productivity in the wetland where the organics are produced/deposited. In this regard, variations in the C/ N ratio have previously been used to decipher the organic matter source in lacustrine sediments (Meyers, 1994). Freshwater algae are rich in Fig. 4. Comparison of the last 2700 cal yr BP between selected proxies from BdlR-03, regional proxies of sedimentary records from Sierra Nevada and other climatic proxies: (a) Total Solar Irradiance (TSI; Steinhilber et al., 2009), (b) NAO index reconstruction (Trouet et al., 2009;Olsen et al., 2012) protein and poor in cellulose (enriched in nitrogen and impoverished in carbon) respect to vascular plants; therefore, low atomic C/N ratios (<10) suggest algal organic matter, whereas high C/N ratios (>20) indicate the predominance of a vascular plant source (enriched in carbon and impoverished in nitrogen) in the sediments. C/N values between 10 and 20 involve a mixed contribution (Meyers, 1994;Meyers and Teranes, 2001). High TOC values could indicate either aquatic or terrestrial enhanced primary productivity but are usually related to increases in land-derived organic matter (Meyers and Horie, 1993).
Isotopic values of δ 13 C and δ 15 N can help specify the organic matter provenance and to understand biogeochemical processes like productivity levels (Meyers and Horie, 1993;Meyers, 1994;Herczeg et al., 2001;Talbot, 2001). Aquatic plants and algae preferentially incorporate the lighter isotopes of dissolved inorganic carbon (DIC) and dissolved inorganic nitrogen (DIN) from the water masses, which mainly depends on allochthonous aeolian inputs in Sierra Nevada (Pulido-Villena et al., 2005;Morales-Baquero and Pérez-Martínez, 2016). Therefore, a productivity increase in Sierra Nevada alpine aquatic systems would generate an isotopic enrichment of these elements in the water pools as well as in algae or aquatic plants (O'Leary, 1988;Hodell and Schelske, 1998;Wolfe et al., 2002), unless a continuous input of these allochthonous elements occur. Therefore, high values of δ 13 C and δ 15 N in the studied sediments might suggest high primary productivity in the aquatic systems (Brenner et al., 1999;Talbot and Laerdal, 2000;Teranes and Bernasconi, 2000). δ 13 C is also used as an indicator of photosynthetic pathways in land plants, and oscillating isotopic carbon values in the C3-dominated-vegetation BdlR-03 record, might be related to changes in their water use efficiency (Farquhar et al., 1982). Even though isotopic data from current vegetal organic matter are unavailable from the study site, carbon isotopic values of plants from other alpine wetlands in the Sierra Nevada exhibit values from − 30.3‰ to − 25.8‰ . Unpublished plant isotopic data from other alpine wetlands in the area show higher values, up to − 23.8‰, and semiaquatic plants such as bryophytes and Ranunculaceae usually display values higher that − 25.5‰. Algal organic matter registers more isotopically-enriched values ranging from − 22.8‰ to − 20.2‰ . These later values agree with the δ 13 C data obtained in the BdlR-03 record and suggest algal occurrence, which is also supported by C/N values <20, such as from ~400 to 60 cal yr BP, usually related to higher carbon isotope values. Although we do not know the precise source of the signal contributed by algae, we assume it represents either periodic flooding of the bog surface by the stream, perhaps fed by increased productivity emanating from the pools upstream, or possible from algae living within the bog deposit itself.
Leaf wax indices (n-alkanes) help to characterize more specifically the source of organic matter through time. The length of the carbon chains in n-alkanes can be related to different kinds of vegetation in the catchment basins and/or different ranges of aridity/temperatures (Bai et al., 2009;Luo et al., 2012;Bush and McInerney, 2015). Three indices have been calculated in order to summarize all the information obtained from the n-alkanes. The ACL is a useful index to characterize the length of the carbon chains in a specific sample (Poynter and Eglinton, 1990). The CPI informs about the predominance of even or odd carbon chains lengths in one sample. Values lower than 2 (even n-alkanes preference) suggesting diagenetic alteration or algae/bacteria influence and values higher than 2 (odd n-alkanes preference) suggesting terrestrial plant source and/or thermal immaturity of the source rock (Bush and McInerney, 2013). The P aq represents the ratio between typical aquatic and terrestrial n-alkanes (Ficken et al., 2000) in one sample. The analyses of these indices in sedimentary records help understand the local ecological responses to environmental changes throughout time . The environmental responses registered in the n-alkanes vary depending on the study areas (Vonk and Gustafsson, 2009;Bai et al., 2009;Luo et al., 2012;Bush and McInerney, 2015) and a correct interpretation of these biomarkers must consider the study of present-day vegetation around the study sites. Previous studies based on a plant survey in Sierra Nevada wetlands proposed that the distance between plants and the main water sources controls the length of the nalkane carbon chains, and thus, plants near water pools displayed strong predominance of the shorter carbon chains, registering low ACL 25-33 and high P aq values (García-Alix et al., 2017. Therefore, fluctuation of these indices in the local sedimentary record can help interpret changes in either the local water availability or the contribution of aquatic plants to the sediments (García-Alix et al., 2017). This interpretation can also be applied to the studied alpine record, where PCA analyses (Fig. S1) highlighted the differences between mid-long and long chain lengths (C 21-27 vs C 29-35 ), which would ultimately reflect temporal changes in the local humid conditions (as recorded by PC2). Although both proxies (PC2 and P aq ) exhibit the same trends during the latest Holocene (Pearson correlation = − 0.74, p < 0.0001), we will focus the discussion on the scores of the PC2 since it integrates the entire odd n-alkane distribution. The P aq will also be referred occasionally.

Origin and evolution of the BdlR within the regional paleoclimatic context during the last 2700 years
The record begins with a significant deposition of detrital sediments at the bottom of the sedimentary record between ~2700 and ~2600 cal yr BP, probably related to significant erosion rates due to runoff on the steep slopes of the glacial cirque where the BdlR-03 is located (maximum values of MS and K/Ti; Figs. 2 and 3a). The record evolved from the eroded bedrock material with scarce organic matter (minimum values of TOC and TN; Fig. 3h and i; Zone 0) into a bog environment.
Organic matter from the BdlR-03 record lacks any potential evidence of degradation. This is indicated by the CPI 25-33 values (>3.5) that would point towards a major vascular plant source (Niedermeyer et al., 2010;Bush and McInerney, 2013), in agreement with the predominant peaty lithology described in the record. The good preservation of organic matter allows for a correct interpretation of both organic biomarkers and bulk organic proxies analysed in this study.

Humid period and bog environment between ~2600 and ~1870 cal yr BP
Maximum humid conditions and water availability in the BdlR-03 record are observed between ~2600 and 1870 cal yr BP ( Fig. 3; Zone 1). Bog peaty deposits accumulated between ~2600 and 2200 cal yr BP, during the most humid period interpreted for this record, agreeing with an increase in water availability and a greater development of terrestrial vascular plants deduced by a decrease in PC2 and ACL 25-33 , high C/N values and C 31 like dominant n-alkane (Fig. 3c, e, g and Fig. S1a). Subsequently, a gradual reduction of the runoff and a moderate decrease in the occurrence of terrestrial vascular plants occurred in the bog environment between ~2200 and ~1870 cal yr BP (decrease of K/Ti and C/N values and C 29 like dominant n-alkane; Fig. 3a, g and Fig. S1a). These evidences with the moderate increase of PC2 and ACL 25-33 values could also point towards a transition to more arid environments. Nevertheless, low values of PC2 are registered suggesting important water availability in the organic deposits (Fig. 3c). After ~1870 cal yr BP there is an aridification trend indicated principally by a general reduction of water availability (increase in PC2 and ACL 25-33 ; Fig. 3 c and e) and vascular plants (C/N decrease; Fig. 3g).
Similar climatic variability to that recorded in the BdlR-03 between 2600 and 1870 cal yr BP is shown in other previous paleoclimatic studies in Iberian lake records (Martín-Puertas et al., 2009;Currás et al., 2012). For example, the humid bog phase identified in BdlR-03 coincides with the timing of the Iberian Roman Humid Period (IRHP) when the highest humidity of the Late Holocene was reached in the western Mediterranean. This humidity maximum is also observed in different paleoclimatic marine and lacustrine records (Fletcher and Sánchez Goñi, 2008;Combourieu Nebout et al., 2009;Martín-Puertas et al., 2009;Jiménez-Moreno et al., 2013;Ramos-Román et al., 2016;Mesa-Fernández et al., 2018;among others). Overall humid conditions in the Mediterranean between ~2800 and ~2100 cal yr BP have previously been explained by high winter rainfall triggered by persistent negative (or low) NAO conditions (Olsen et al., 2012) (Fig. 4b). This period also registered significant high values in the total solar irradiance (TSI) related to a solar maximum (Solanki et al., 2004) (Fig. 4a) and coeval with warm temperatures in Sierra Nevada, according to a lipid-derived temperature reconstruction obtained from a neighbor Sierra Nevada alpine wetland, Laguna de Río Seco (LdRS; Toney et al., 2020) (Fig. 4c). The warm and humid climate contributed to the important development of terrestrial vascular plant cover observed in the BdlR-03 until ~1870 cal yr BP (maximum values of C/N and TOC), in agreement with García-Alix et al. (2017) and Toney et al. (2020) that suggested that biomass development in these extreme environments of Sierra Nevada was favored by warmer temperatures or shorter cold seasons that provide more ice-free surfaces for longer periods and greater soil and vegetation development. The abundance of Mediterranean forest species (such as Pinus and deciduous Quercus, among others) observed in pollen records from surrounding alpine wetlands suggested high relative humidity (Fig. 4e) and significant forest development at lower elevations during this period (Anderson et al., 2011;Jiménez-Moreno et al., 2013;Ramos-Román et al., 2016;Mesa-Fernández et al., 2018;Manzano et al., 2019).

Arid period in the bog between ~1870 and ~300 cal yr BP
The BdlR-03 record shows a general trend towards climatic aridification linked with a reduction in the aquatic environments and a lower presence of terrestrial vascular plants between ~1870 and ~300 cal yr BP ( Fig. 3; Zone 2). Similar aridification have been depicted in different western Mediterranean records (Carrión, 2002;Fletcher and Sánchez Goñi, 2008;Jalut et al., 2009;Carrión et al., 2010). Concretely, in the case of the Sierra Nevada wetlands, this aridification is coeval with a decrease in natural forest species such as Pinus and a progressive increase in xerophytes such as Artemisia (Anderson et al., 2011;Ramos-Román et al., 2016;Mesa-Fernández et al., 2018;Manzano et al., 2019). In addition, Jiménez-Moreno and Anderson (2012) interpreted the decline of Pinus in this region as a movement of tree species towards lower elevations due to a climatic cooling produced by a decrease in summer insolation at these latitudes during the Late Holocene (Laskar et al., 2004).
The high-resolution record from BdlR-03 also permits to identify three phases with contrasted environmental conditions within this period: a first relatively arid phase between ~1870 and ~1370 cal yr BP (Subzone 2A), a second phase with a moderate increase in water availability between ~1370 and ~820 cal yr BP (Subzone 2B) and finally a third phase with a general arid and progressive cold conditions between ~820 and 300 cal yr BP (Subzones 2C and 2D).
The BdlR-03 record shows a significant aridification phase between ~1870 and ~1370 cal yr BP deduced by the increasing trend of PC2 and ACL 25-33 and the decreasing trend of C/N ratio (Fig. 3c, e and g), coinciding with the transition between the end of the IRHP and the "Dark Ages" (DA) (Fig. 4). The general decrease in terrestrial vascular plants (TOC decrease, C/N value slightly above 20 and C 31 like dominant n-alkane; Fig. 3h, g and Fig. S1a) was also probably linked to an overall climate cooling according to the LdRS estimated temperature record . A lower biomass along with lower substrate fixation would have probably triggered an increase in erosion and runoff (increasing trend in K/Ti; Fig. 3a) and further detrital sediment accumulation in the bog. Despite the regional aridity trend, a moderate increase in water availability in the BdlR-03 (low PC2 values; Fig. 3c) was observed between ~1370 and ~820 cal yr BP. This could have been due to an occasional local wetland development most likely caused by enhanced runoff (high but decreasing values of K/Ti after ~1290 cal yr BP, when the highest value of K/Ti in the record was registered; Fig. 3a) as a consequence of a significant and gradual melting of perennial ice and/or snow packs at higher elevations boosted by gradual increase in Sierra Nevada temperatures Toney et al., 2020) during the DA-MCA (Medieval Climate Anomaly) transition (Fig. 4c).
Finally, generally arid and progressive cold conditions are recorded in the BdlR-03 record from ~820 to 300 cal yr BP coinciding with the end of the MCA and the first stages of the LIA. During this period, water availability decreased, confirmed by the PC2 (Fig. 3c), and a significant climatic cooling took place in the area (García-Alix et al., 2020) (Fig. 4c). A greater occurrence of terrestrial vascular plants coincides with warmer temperatures at the end of the MCA (high TOC, C/N > 20 and C 31 like dominant n-alkane; Figs. 3h, 4g and Fig. S1a). Although overall dry conditions are observed for the period, a slight humidity recovery might be interpreted by the decrease in PC2 and ACL 25-33 values ( Fig. 3c and e), a small increase in runoff (K/Ti; Fig. 3a), and a moderate decrease in aeolian inputs (Ca/Ti; Fig. 3b) between 600 and 300 cal yr BP in agreement with persistent episodes of NAO negative conditions (Trouet et al., 2009) (Fig. 4b).
More arid conditions and lower temperatures between ~820 and 300 cal yr BP probably led to less biomass development (low TOC and drop in % of Pinus) and an increase in herbaceous and xerophytic species (increases ACL 25-33 and % of Artemisia) in Sierra Nevada, as observed in the BdlR-03 and other previously published records such as Laguna de la Mula (LM; Jiménez-Moreno et al., 2013), Borreguil de la Caldera (BdlC; Ramos-Román et al., 2016) and Laguna Hondera (LH; Mesa-Fernández et al., 2018) (Fig. 4c, d, e and f). These arid conditions were significant in other areas of southern Iberia, causing for example low lake levels in Zoñar Lake (Martín-Puertas et al., 2010). Climate in the western Mediterranean during this period alternated between arid and cold conditions with a warmer period during the MCA and wetter stages at the beginning of the LIA (Martín-Puertas et al., 2010;Nieto-Moreno et al., 2013. Previous studies in the area related droughty and cold conditions (Ramos-Román et al., 2016;García-Alix et al., 2017) with persistent positive NAO atmospheric mode (Trouet et al., 2009;Olsen et al., 2012) and overall decrease in temperatures Toney et al., 2020) probably linked with, among other factors, a decreasing solar output (Steinhilber et al., 2009) (Fig. 4a).

Unstable climate period and wetland environment between ~300 cal yr BP and the present
An increase in aquatic algae, contributed either from streamflow or in situ development, is interpreted in the BdlR-03 between ~300 and ~0 cal yr BP ( Fig. 3; Zone 3). The increase in aquatic algae is deduced by high values of δ 15 N and TN ( Fig. 3k and i), and the lowest and highest values of C/N ratio and δ 13 C, respectively ( Fig. 3g and j). Furthermore, the development of ephemeral pools is deduced by decreasing trends of PCA and ACL 25-33 , respectively. After the coldest period of the LIA in Sierra Nevada (~1690 CE) there was an important increase in temperatures (LdRS lipid-derived temperatures; Fig. 4c) and a significant increase in solar irradiance that could have probably influenced the gradual melting of the ice from higher altitudes along with the aerosol deposition that could have reduced the snow albedo ). An important melting event was recorded in the BdlR-03 record at ~1850 CE (absolute maximum of P aq ; Fig. 3f), slightly before to other melting events described in Sierra Nevada during the 19th-20th century transition (Gómez-Ortiz et al., 2009;García-Alix et al., 2017;Toney et al., 2020). This early melting is probably related to the lower elevation of the studied area and its catchment respect to other local sedimentary alpine wetlands. As a consequence, the increasing runoff (high K/Ti ratio; Fig. 3a) and high-water availability (low PCA and ACL 25-33 ; Fig. 3c and e) may have created a more humid wetland environment. This opposite local humidity pattern respect to the general aridification pattern in the western Mediterranean region have also been also observed in other alpine peat bog in Sierra Nevada, BdlC (García-Alix et al., 2017). The gradual isotopic depletion in carbon and nitrogen isotopes during the last ~200 years, when moderate C/N values are recorded, could be explained by the high aerosol inputs registered during this period that would have continuously contributed with nutrients to these nutrient-impoverished alpine environments, renewing the C and N pools and preventing from isotopic enrichment (Fig. 3b, j and k). In addition, this carbon isotope depletion could also be related to high water availability that would reduce the water stress in terrestrial plants (Farquhar et al., 1989;O'Leary, 1988). This increase in aerosol deposition is observed in other records from Sierra Nevada such as in LdRS (Jiménez-Espejo et al., 2014), and is probably related to the beginning of commercial agriculture in the Sahel region in the last 200 years  and the enhanced arid conditions (Jiménez-Espejo et al., 2014). This climatically-unstable period during the last centuries, with a persistent negative NAO phase and many short-term abrupt oscillations (Trouet et al., 2009), shows similar environmental responses in other alpine areas of the Sierra Nevada. Concretely, decreasing trends, although with fluctuations, of Artemisia and Pinus as well as a moderate increase in the water availability deduced were observed in the surrounding BdlC, LdRS or LH wetlands (Anderson et al., 2011;Ramos-Román et al., 2016;García-Alix et al., 2017Toney et al., 2020).
Finally, the last ~100 years in the BdlR-03 record are characterized by an increase in biomass deposition, mainly as terrestrial vascular plants (indicated by high C/N values), and significant runoff. However, significant eutrophication has not been detected in the BdlR ephemeral pools, according to the depletion in the nitrogen isotopes. The environmental signal in other Sierra Nevada alpine wetlands seems to be significantly affected by anthropogenic forcing (García-Alix et al., 2017; among others) with clear evidence of reforestation (Valbuena-Carabaña et al., 2010;Anderson et al., 2011). As a consequence, some natural species have altered their natural distributions as it is the case of Pinus, whose current forced-treeline is at ~2550 mas, while natural populations occur at ~2100 masl in the area (Anderson et al., 2011). This anthropically-forced shift in the natural distribution of plants could have produced the C/N increase in the studied record during the last ~60 years. Other human activities that are recorded at that time in Sierra Nevada are industrial activities such as mining  and agricultural activities and cattle (Ramos-Román et al., 2019).

Conclusions
The results obtained by multiproxy analyses in the BdlR-03 sedimentary record show: 1. An overall trend towards climatic aridification and reduction of the water availability in the BdlR area during the last 2700 years, interpreted from the progressive increase in PC2 and ACL 25-33 values and decrease C/N values. Oscillations in the C/N ratio and isotopic values of δ 13 C and δ 15 N throughout the record suggest a contribution of organic matter from terrestrial vascular plants in the early stages of the bog and a mixed contribution from terrestrial vascular plants and algae in the final stages. The progressive decrease in terrestrial vascular plants is consistent with the proposed climatic aridification, since in this extreme alpine environment water availability is one of the main controls on the vegetation development.
2. Solar forcing and NAO fluctuations could be the most important factors controlling this aridification trend, and thus, the regional and local environmental changes. The gradual insolation decrease together with the water availability reduction in the study area would likely cause the decrease on primary production in the catchment (vascular plant decline) deduced by the C/N trends.
3. Four phases are differentiated within the paleoenvironmental evolution of BdlR-03. After a first pre-bog stage (2700-2600 cal yr BP) characterized by a thin layer of eroded bedrock material and scarce organic matter, the first bog environment was established with important vascular plant content, water availability and maximum humid conditions between ~2600 and ~1870 cal yr BP. A second drier bog environment linked with a productivity decrease (lower TOC than in the previous stage) that gave rise to lower presence of terrestrial vascular plants and reduction of aquatic environments between ~1870 and ~300 cal yr BP. And finally, a wetland environment with an increase of aquatic algae and the development of ephemeral pools between ~300 cal yr BP and the present. This scenario can likely differ during the last ~100 years since the environmental signal in the BdlR-03, deduced by an increase in C/N values, seems to be significantly affected by anthropogenic forcing with clear evidence of reforestation as in other Sierra Nevada alpine wetlands.
4. The potential proliferation of wetland environments in the final stages of the record (especially after 1850 CE), interpreted by a decrease in PCA as well as low values of C/N (<20) and TOC, could be contrary to the regional aridification trend in the last millennium. The increase in algal/aquatic organic content can be explained by a gradual increase in runoff (detritic input as observed in the K/Ti ratio), probably due to the melting of perennial ice and/or snow packs that transformed the peat bog into a wetland environment with ephemeral pools. Another significant trigger for this aquatic primary production development would be the progressive increase in North African aeolian inputs (increase in the Ca/Ti ratio, reaching the highest values in the record) that would have supplied nutrients for the development of local biogeochemical cycles in aquatic environments during the last ~200 years. Similar enhanced water availability is observed in other alpine records of Sierra Nevada.
5. High North African dust inputs would have contributed to renew the isotopic signature of the local nutrient pool of these nutrientimpoverished alpine environments. This fact is specially marked during the last ~200 years agreeing with an aerosol deposition increase in the studied record, likely linked with the beginning of commercial agriculture in the Sahel region and the regional enhanced arid conditions.

Declaration of Competing Interest
None.