High-Resolution Temperature Variability Reconstructed from Black Pine Tree Ring Densities in Southern Spain

: The presence of an ancient, high-elevation pine forest in the Natural Park of Sierras de Cazorla in southern Spain, including some trees reaching > 700 years, stimulated efforts to develop high-resolution temperature reconstructions in an otherwise drought-dominated region. Here, we present a reconstruction of spring and fall temperature variability derived from black pine tree ring maximum densities reaching back to 1350 Coefficient of Efficiency (CE). The reconstruction is accompanied by large uncertainties resulting from low interseries correlations among the single trees and a limited number of reliable instrumental stations in the study region. The reconstructed temperature history reveals warm conditions during the early 16th and 19th centuries that were of similar magnitude to the warm temperatures recorded since the late 20th century. A sharp transition from cold conditions in the late 18th century (t 1781–1810 = − 1.15 ◦ C ± 0.64 ◦ C) to warm conditions in the early 19th century (t 1818–1847 = − 0.06 ◦ C ± 0.49 ◦ C) is centered around the 1815 Tambora eruption (t 1816 = − 2.1 ◦ C ± 0.55 ◦ C). The new reconstruction from southern Spain correlates significantly with high-resolution temperature histories from the Pyrenees located ~600 km north of the Cazorla Natural Park, an association that is temporally stable over the past 650 years (r 1350–2005 > 0.3, p < 0.0001) and particularly strong in the high-frequency domain (r HF > 0.4). Yet, only a few of the reconstructed cold extremes (1453, 1601, 1816) coincide with large volcanic eruptions, suggesting that the severe cooling events in southern Spain are controlled by internal dynamics rather than external (volcanic) forcing.


Introduction
The climate of the Iberian Peninsula is characterized by a sustained warming trend of >1.0 • C since the late 20th century [1]. In southern Spain, this trend is accompanied by a precipitation decline of~50 mm since the 1960s [2], though these hydroclimatic changes are spatially and seasonally more

Tree-Ring Data and Detrending
Whereas most of the CNP is dominated by relatively young trees, there are a few high-elevation sites towards the southern end of the park covered by an ancient Pinus nigra forest including several individuals reaching ages >700 and even >800 years (Table 1). These sites have been subject to the development of large tree-ring datasets, including several hundred TRW series [25], though only a few MXD series have been produced so far [26]. Here, we combine the MXD series sampled in 2007, Here, we present a combined MXD data set from 51 high-elevation CNP black pines reaching back to 1196 CE. We detail the age trend inherent to these data and develop several detrended chronologies integrating differently old tree rings, so called age-band chronologies. The chronologies are calibrated against instrumental temperature data with particular emphasis on the effects of representative climate Atmosphere 2020, 11, 748 3 of 17 stations (or lack thereof) on MXD signal estimation and reconstruction uncertainty. We finally present a formal temperature reconstruction derived from Pinus nigra MXD data and compare this record with existing records from CNP TRW and Pyrenees MXD data.

Tree-Ring Data and Detrending
Whereas most of the CNP is dominated by relatively young trees, there are a few high-elevation sites towards the southern end of the park covered by an ancient Pinus nigra forest including several individuals reaching ages >700 and even >800 years (Table 1). These sites have been subject to the development of large tree-ring datasets, including several hundred TRW series [25], though only a few MXD series have been produced so far [26]. Here, we combine the MXD series sampled in 2007, produced in the tree-ring laboratory in Tharandt (Germany), with a larger compilation of the new MXD series sampled in 2015, produced in the laboratory in Mainz (Germany). The density profiles of these two campaigns have been developed using the same Walesch radiodensitometric setup [31][32][33], and the site MXD chronologies share high fractions of common variance over the past 450 years (r 1557-2006 = 0.80), justifying the amalgamation of MXD data into one pool of 88 series from 51 trees (see Supplementary Table S1 and Figures S1-S3). The combined MXD dataset covers the period from 1196-2014 CE and includes 18 series exceeding ages of 600 years ( Figure 2a). The data contain an age trend composed of a notable~0.2 g/cm 3 increase over the first 90 years of cambial age, followed by a persistent, though much shallower, negative trend over the subsequent 500 years (−0.12 g/cm 3 from 100-500 years). This Hugershoff-shape [34] age trend is characteristic for conifer MXD data [35,36] and represents noise from the perspective of a climate reconstruction. This noise is here removed by applying Regional Curve Standardization (RCS) to preserve high-to-low frequency variance in detrended index chronologies [37,38]. However, successful RCS-detrending typically requires the application to combined datasets integrating series from living and dead trees, so that younger and older tree-rings are spread throughout time, rather than being concentrated at the beginning (young rings) and end (old rings) of a chronology [39]. To produce such chronology variants, characterized by flat instead of monotonically increasing mean age curves, we removed all rings older than 300 (200) years from the MXD data and calculated 1-300  year age-band chronologies, here labeled ABC300 and ABC200 [40]. The differentiation by cambial age is illustrated in Figure 2, showing the replication and regional curves (RCs) of all MXD data in gray, and the 1-200 and 201-300 sections in blue and black, respectively. Figure 3 illustrates the effects of data truncation on the temporal distribution of MXD data throughout the past 800 years and Atmosphere 2020, 11,748 4 of 17 demonstrates that any chronology derived from the 1-200 year age band data is weakly replicated in the late 16th century, for instance (see the circle in Figure 3c).
The ALL, 1-300 and 1-200 year datasets were RCS-detrended by calculating ratios between the single MXD series and the smoothed mean of the age-aligned data, the so-called regional curves (RCs, shown in Figure 2b). The resulting index values were averaged using the arithmetic mean, and uncertainties estimated by calculating 95% bootstrap confidence limits derived from re-sampling and averaging the detrended MXD data with replacement 1000 times [37]. Mean tree age and segment length curves, as well as running inter-series correlations (Rbar) [41] and the Expressed Population Signal (EPS) [42] were calculated for each of the chronology, ALL, ABC300, and ABC200, to illustrate temporal changes in the structure and coherency of the underlying data over the past 700 years.

Climate Data, Calibration and Transfer
The ALL, ABC300 and ABC200 chronologies were calibrated against gridded monthly temperatures (CRU TS4.03) using the KNMI Explorer [43] to assess the seasonality and spatial extent of MXD climate signals. Station temperature records from five locations in the vicinity of the proxy site were used to evaluate the temporal robustness of climate signals as well as the fidelity of the observational data to represent the conditions at the high elevation CNP tree site. The station records are located in distances between 10 km (Cazorla) and 165 km (Murcia) from the MXD site, cover varying periods between 1905-2014 CE, and contain changing numbers of missing values ranging from 3% in Murcia to 24% in Jaen ( Table 2). Four of the station records (Jaen, Ciudad Real, Albacete, Murcia) are used in the gridded temperature products (CRU TS4.03), whereas the station closest to the tree site (Cazorla) is not included in the international climate databases [1,44,45].  The ALL, 1-300 and 1-200 year datasets were RCS-detrended by calculating ratios between the single MXD series and the smoothed mean of the age-aligned data, the so-called regional curves (RCs, shown in Figure 2b). The resulting index values were averaged using the arithmetic mean, and uncertainties estimated by calculating 95% bootstrap confidence limits derived from re-sampling and averaging the detrended MXD data with replacement 1000 times [37]. Mean tree age and segment length curves, as well as running inter-series correlations (Rbar) [41] and the Expressed Population Signal (EPS) [42] were calculated for each of the chronology, ALL, ABC300, and ABC200, to illustrate temporal changes in the structure and coherency of the underlying data over the past 700 years.

Climate Data, Calibration and Transfer
The ALL, ABC300 and ABC200 chronologies were calibrated against gridded monthly temperatures (CRU TS4.03) using the KNMI Explorer [43] to assess the seasonality and spatial extent of MXD climate signals. Station temperature records from five locations in the vicinity of the proxy site were used to evaluate the temporal robustness of climate signals as well as the fidelity of the observational data to represent the conditions at the high elevation CNP tree site. The station records are located in distances between 10 km (Cazorla) and 165 km (Murcia) from the MXD site, cover varying periods between 1905-2014 CE, and contain changing numbers of missing values ranging from 3% in Murcia to 24% in Jaen ( Table 2). Four of the station records (Jaen, Ciudad Real, Albacete, Murcia) are used in the gridded temperature products (CRU TS4.03), whereas the station closest to the tree site (Cazorla) is not included in the international climate databases [1,44,45]. The MXD chronologies were compared with the instrumental temperatures using the Pearson correlation coefficient calculated over a longer 1905-2014 period as well as over a shorter, but observationally better replicated, 1961-2014 period. Thirty-year running correlations were applied to emphasize temporal changes in signal strength throughout the 20th and early 21st centuries. In addition, we assessed the covariance of monthly and seasonal temperatures among the station records and evaluated the effects of deviating cold periods and extreme years on proxy data calibration.
The ABC300 chronology was finally transferred into estimates February-May and September-October (FMAM&SO) temperatures by scaling the MXD mean and variance against the mean of the five station records over the 1905-2014 period. The skill of the reconstruction was estimated considering the Reduction of Error (RE) and Coefficient of Efficiency (CE) statistics calculated over an early 1905-59 (55 years) and late 1960-2014 (55 years) calibration/verification period, after regressing the ABC300 chronology against the instrumental data [46]. The Durbin-Watson statistic (DW) was calculated to evaluate autocorrelation in the residuals between the instrumental temperatures and regressed MXD chronology [47]. Temporally changing uncertainties of the reconstruction were estimated using the 95% bootstrap confidence limits of the ABC300 chronology, after smoothing these estimates using a 30-year low-pass filter, to emphasize long-term confidence changes over the 1350-2014 reconstruction period. The coldest and warmest reconstructed 30-year periods, before and after 1700 CE, as well as the ten coldest and ten warmest years since the mid-14th century were highlighted for further discussion.
The FMAM&SO temperature reconstruction was compared with annually resolved temperature histories from a MXD network in the Pyrenees, located~600 km north of the CNP [19], as well as a previous-year September-October temperature reconstruction derived from TRW data of the CNP black pines [27]. These comparisons were conducted using the original as well as high-pass filtered versions of the reconstructions to evaluate covariance in the high-to-low frequency domains. Particular attention is paid to the reconstructed cold extremes, their coherency between southern (CNP) and northern Spain (Pyrenees), and the potential underlying forcing that exist as a result of large volcanic eruptions since 1350 CE [48][49][50][51][52].

Cazorla MXD Chronologies and Effects of Data Truncation
The RCS-detrended chronologies (ALL, ABC300, ABC200) share high fractions of high-to-low frequency variance, but also include substantial trend differences, e.g., since the mid-20th century, during the mid-19th century, and the mid-16th century (Figure 4a, Supplementary Figure S4). The original and 30-year low-pass filtered chronologies correlate at r Orig = 0.89 and r LF = 0.81 since 1300 CE, respectively, indicating that the truncation of tree-rings older than 300 and 200 years had a stronger effect on the low frequency variance. Over the most recent~70 years, the smoothed ABC200 chronology deviates substantially from the two other chronologies, including highest values in 2014. The ALL chronology, instead, shows a decline since the 1980s, whereas the ABC300 chronology remains in-between, i.e., shows high values in the 1980s but also a (minor) increase during the most recent years. Similar inter-chronology deviations are seen in the early and mid-19th century including high ABC200 and low ALL values, and the mid and late 16th century including low ABC200 and high ALL values.  Thirty-year smoothed Regional Curve Standardization (RCS) chronologies of the full MXD data (gray), the 1-300 age band data (black), and the 1-200 age band data (blue), shown together with their replication (b) and mean age curves (c). See Supplementary Figure S4 for the original (non-smoothed) chronologies and mean segment length curves.
The effects of data truncation are obvious in the chronology replication and mean age curves showing increasing inter-chronology differences towards present (Figure 4b,c). The sharp negative deviation of the ABC200 chronology in the mid-16th century occurs during a period of minimum replication (n1551-1560 ≤ 8 MXD series), indicating that this low frequency departure is potentially less reliable, compared to the better-replicated chronologies (nABC300 ≤ 17 and nALL ≤ 22 series). On the other hand, the truncation of rings older than 200 years in the ABC200 chronology, produced an almost horizontal mean age curve, revealing that the age-structure of this dataset is most suitable for the application of RCS detrending [38,39,53].
The ABC300 chronology is also derived from a relatively even-aged dataset characterized by an age range of only 103 years between the 1360s and 1540s (age1351-60 = 61 years, age1541-50 = 164 years). In comparison, the monotonically increasing mean age curve of the ALL chronology, from 60 to 331 years (see Supplementary Figure S4b for the segment length curves) biases the RCS detrending and produces an artificial shift towards the long-term mean, due to the dominance of old tree-rings from recent calendar years in the biologically old sections of the RC (the gray curve in Figure 2b; details in [38,54,55]). Mitigating this bias motivated the development of age-band detrending, a technique that was successfully applied to a large data set of Northern Hemisphere MXD sites [40].
The removal of data older than 300 and 200 years produced chronologies that are likely more reliable over the most recent centuries, during which the ALL chronology is characterized by a monotonic age increase. During earlier chronology periods, particularly before 1650 CE in ABC200, Thirty-year smoothed Regional Curve Standardization (RCS) chronologies of the full MXD data (gray), the 1-300 age band data (black), and the 1-200 age band data (blue), shown together with their replication (b) and mean age curves (c). See Supplementary Figure S4 for the original (non-smoothed) chronologies and mean segment length curves.
The effects of data truncation are obvious in the chronology replication and mean age curves showing increasing inter-chronology differences towards present (Figure 4b,c). The sharp negative deviation of the ABC200 chronology in the mid-16th century occurs during a period of minimum replication (n 1551-1560 ≤ 8 MXD series), indicating that this low frequency departure is potentially less reliable, compared to the better-replicated chronologies (n ABC300 ≤ 17 and n ALL ≤ 22 series). On the other hand, the truncation of rings older than 200 years in the ABC200 chronology, produced an almost horizontal mean age curve, revealing that the age-structure of this dataset is most suitable for the application of RCS detrending [38,39,53].
The ABC300 chronology is also derived from a relatively even-aged dataset characterized by an age range of only 103 years between the 1360s and 1540s (age 1351-60 = 61 years, age 1541-50 = 164 years). In comparison, the monotonically increasing mean age curve of the ALL chronology, from 60 to 331 years (see Supplementary Figure S4b for the segment length curves) biases the RCS detrending and produces an artificial shift towards the long-term mean, due to the dominance of old tree-rings from Atmosphere 2020, 11, 748 7 of 17 recent calendar years in the biologically old sections of the RC (the gray curve in Figure 2b; details in [38,54,55]). Mitigating this bias motivated the development of age-band detrending, a technique that was successfully applied to a large data set of Northern Hemisphere MXD sites [40].
The removal of data older than 300 and 200 years produced chronologies that are likely more reliable over the most recent centuries, during which the ALL chronology is characterized by a monotonic age increase. During earlier chronology periods, particularly before 1650 CE in ABC200, the data truncation likely weakened the age-band chronologies, as the already reduced replication of the ALL chronology is further lowered by removing old rings. Before 1500 CE, however, these effects become negligible, as the ALL chronology is composed of only a few young rings, i.e., no additional rings were removed in the ABC300 and ABC200 chronologies.

MXD Climate Signals and Uncertainties
The seasonality of MXD temperature signals is bimodal, including significant fields in February-March, May, and September-October surrounding the CNP ( Figure 5). February and September temperatures are most influential, particularly when focusing on the shorter, and observationally better replicated 1961-2014 period, during which correlations near the CNP tree site exceed r = 0.5 (Supplementary Figure S5). The bimodal nature of the signal, comprising a lack of forcing during the warm June-August summer months, is similar to the climate sensitivity reported from high-elevation Pinus uncinate MXD data in the Spanish Pyrenees [19]. The underlying physiological mechanisms are likely related to the insensitivity of cell wall formation to summer warmth, when temperatures do not fall below thresholds relevant to carbohydrate production and mobilization in high elevation black pines (see [19] for a detailed discussion).
Atmosphere 2020, 11, x FOR PEER REVIEW 8 of 18 warm June-August summer months, is similar to the climate sensitivity reported from high-elevation Pinus uncinate MXD data in the Spanish Pyrenees [19]. The underlying physiological mechanisms are likely related to the insensitivity of cell wall formation to summer warmth, when temperatures do not fall below thresholds relevant to carbohydrate production and mobilization in high elevation black pines (see [19] for a detailed discussion). The bimodal response characteristic for CNP (and Pyrenees) MXD data is unique if compared with other long MXD chronologies from cold environments in the Mediterranean, central and northern Europe, which are all unimodal, i.e., show a maxima during summer months [39,[56][57][58]. For instance, the millennium-length MXD chronologies from Greece [59], the Swiss Alps [60], and Fennoscandia [35,37,61] correlate best with July-September, June-September and June-August temperatures, respectively. The CNP MXD data, at the southern end of such a European transect, thereby reinforces a fading importance of the conditions during the warmest summer months, a tendency already demonstrated in the Pinus uncinata MXD data from the Spanish Pyrenees [19].
Besides the shift in seasonality towards bimodal, the correlations are also lower compared to the MXD counterparts from the Pyrenees (r = 0.72) and Greece (r = 0.73), but also the Alps (r = 0.69) and northern Europe (r = 0.77). The CNP black pine correlations against varying seasonal targets, including FMAM&SO, hardly exceed r = 0.4 over the 1905-2014 period ( Figure 6). The values increase over the shorter 1961-2014 period, though only the trials against the shorter February plus September (F&S), February-March plus September (FM&S), and February-March plus September-October (FM&SO) seasons benefit from constraining the calibration period to the most recent, and instrumentally best-replicated, decades (the blue, orange and gray bars in Figure 6a). On the other hand, the longer FMAM&S and FMAM&SO seasons correlate better, if both the instrumental and proxy data were filtered to emphasize high frequency variability in the timeseries. The bimodal response characteristic for CNP (and Pyrenees) MXD data is unique if compared with other long MXD chronologies from cold environments in the Mediterranean, central and northern Europe, which are all unimodal, i.e., show a maxima during summer months [39,[56][57][58]. For instance, the millennium-length MXD chronologies from Greece [59], the Swiss Alps [60], and Fennoscandia [35,37,61] correlate best with July-September, June-September and June-August temperatures, respectively. The CNP MXD data, at the southern end of such a European transect, thereby reinforces a fading importance of the conditions during the warmest summer months, a tendency already demonstrated in the Pinus uncinata MXD data from the Spanish Pyrenees [19].
Besides the shift in seasonality towards bimodal, the correlations are also lower compared to the MXD counterparts from the Pyrenees (r = 0.72) and Greece (r = 0.73), but also the Alps (r = 0.69) and northern Europe (r = 0.77). The CNP black pine correlations against varying seasonal targets, including FMAM&SO, hardly exceed r = 0.4 over the 1905-2014 period ( Figure 6). The values increase over the shorter 1961-2014 period, though only the trials against the shorter February plus September (F&S), February-March plus September (FM&S), and February-March plus September-October (FM&SO) seasons benefit from constraining the calibration period to the most recent, and instrumentally best-replicated, decades (the blue, orange and gray bars in Figure 6a). On the other hand, the longer FMAM&S and FMAM&SO seasons correlate better, if both the instrumental and proxy data were filtered to emphasize high frequency variability in the timeseries. Beyond these seasonal and frequency-dependent changes, the different chronologies-ALL, ABC300 and ABC200-all produce similar results, i.e., the calibration differences are statistically insignificant (not shown). Likely more important are the temporal changes in signal strength, as revealed by the running window correlations, demonstrating substantially higher proxy-target covariances since the 1980s for the shorter F&S, FM&S, and FM&SO seasons (the blue, orange and gray curves in Figure 6c). On the contrary, the longer FMAM&S and FMAM&SO seasons correlate better during the early period of overlap with instrumental data, whereas the shorter season correlations decline before the 1960s. Such early-calibration-period decreases are fairly common [62] and are typically related to increased observational temperature uncertainties due to changes in instrumentation, data gaps, station relocations, etc. [44,[63][64][65][66][67][68][69][70][71]. The recent, post-1980 correlation decline, seen in the longer seasonal means (FMAM&S and FMAM&SO), is more relevant, however, and deserves further attention before producing a formal temperature reconstruction based on the CNP MXD data.

Outlier Effects on Proxy Calibration
Considering the transferred ABC300 chronology, the fit with post-1960 FMAM&SO temperatures is characterized by (i) an offset during the 1970s and early 1980s, and (ii) two negative extremes in 1999 (−2.60 °C) and 2005 (−2.55 °C) that are not reflected in the observational data ( Figure  7). During the 1970s, the average reconstructed temperature is ~1 °C warmer than the average instrumental temperature (see the horizontal bars in Figure 7a), marking a substantial decadal scale proxy-target difference and questioning the reliability reconstructed temperatures at this frequency. However, this offset is substantiated by markedly cooler temperatures recorded at the Cazorla temperature station (dark green in Figure 7b  Beyond these seasonal and frequency-dependent changes, the different chronologies-ALL, ABC300 and ABC200-all produce similar results, i.e., the calibration differences are statistically insignificant (not shown). Likely more important are the temporal changes in signal strength, as revealed by the running window correlations, demonstrating substantially higher proxy-target covariances since the 1980s for the shorter F&S, FM&S, and FM&SO seasons (the blue, orange and gray curves in Figure 6c). On the contrary, the longer FMAM&S and FMAM&SO seasons correlate better during the early period of overlap with instrumental data, whereas the shorter season correlations decline before the 1960s. Such early-calibration-period decreases are fairly common [62] and are typically related to increased observational temperature uncertainties due to changes in instrumentation, data gaps, station relocations, etc. [44,[63][64][65][66][67][68][69][70][71]. The recent, post-1980 correlation decline, seen in the longer seasonal means (FMAM&S and FMAM&SO), is more relevant, however, and deserves further attention before producing a formal temperature reconstruction based on the CNP MXD data.

Outlier Effects on Proxy Calibration
Considering the transferred ABC300 chronology, the fit with post-1960 FMAM&SO temperatures is characterized by (i) an offset during the 1970s and early 1980s, and (ii) two negative extremes in 1999 (−2.60 • C) and 2005 (−2.55 • C) that are not reflected in the observational data (Figure 7). During the 1970s, the average reconstructed temperature is~1 • C warmer than the average instrumental temperature (see the horizontal bars in Figure 7a), marking a substantial decadal scale proxy-target difference and questioning the reliability reconstructed temperatures at this frequency. However, this offset is substantiated by markedly cooler temperatures recorded at the Cazorla temperature station (dark green in Figure 7b only station not included in the Global Historical Climatological Network (GHCN) [72], which could serve as an argument to exclude these data from calibration trials. On the other hand, Cazorla is the closest, and is therefore likely the most representative station for the calibration of high-elevation CNP MXD data (Table 2), which again supports its use. An assessment of the underlying reasons of the deviating 1970s temperatures, be it changes in the station's environment and instrumentation or (real) spatial variability, would require studying the station history and metadata, and monitoring current temperatures at historical sites [66,67], which is beyond the scope of this paper. From a tree-ring perspective, the 1970s proxy-target difference reported here adds uncertainty to any reconstructed lower frequency deviation derived from the CNP MXD data.
Atmosphere 2020, 11, x FOR PEER REVIEW 10 of 18 [72], which could serve as an argument to exclude these data from calibration trials. On the other hand, Cazorla is the closest, and is therefore likely the most representative station for the calibration of high-elevation CNP MXD data (Table 2), which again supports its use. An assessment of the underlying reasons of the deviating 1970s temperatures, be it changes in the station's environment and instrumentation or (real) spatial variability, would require studying the station history and metadata, and monitoring current temperatures at historical sites [66,67], which is beyond the scope of this paper. From a tree-ring perspective, the 1970s proxy-target difference reported here adds uncertainty to any reconstructed lower frequency deviation derived from the CNP MXD data.  (Figure 8b,8c). We do not propose that such data removal should be considered for the calibration and transfer of the CNP proxy data, but intend to emphasize the causes of the post 1980 correlation decline when considering the longer season FMAM&S and FMAM&SO temperature data for calibration. The temporal variability of proxy-target correlations, composed of (i) a pre-1960 decline in the shorter season means and (ii) a post-1980 decline in the longer season means, translates into large uncertainties when transferring the CNP MXD data into temperature estimates. In addition, it is not feasible to statistically validate the seasonality of such a reconstruction based on the calibration against heterogeneous instrumental temperatures. However, since the overall highest correlations were recorded using the ABC300 chronology and FMAM&SO season over the long 1905-2014 calibration period (Figure 6), and since the other chronologies either do not meet RCS-detrending requirements (ALL) or suffer from low sample replications (ABC200; Figure 4), we opt for this combination (ABC300 and FMAM&SO) when producing a formal reconstruction.  (Figure 8b,c). We do not propose that such data removal should be considered for the calibration and transfer of the CNP proxy data, but intend to emphasize the causes of the post 1980 correlation decline when considering the longer season FMAM&S and FMAM&SO temperature data for calibration. The temporal variability of proxy-target correlations, composed of (i) a pre-1960 decline in the shorter season means and (ii) a post-1980 decline in the longer season means, translates into large uncertainties when transferring the CNP MXD data into temperature estimates. In addition, it is not feasible to statistically validate the seasonality of such a reconstruction based on the calibration against heterogeneous instrumental temperatures. However, since the overall highest correlations were recorded using the ABC300 chronology and FMAM&SO season over the long 1905-2014 calibration period (Figure 6), and since the other chronologies either do not meet RCS-detrending requirements (ALL) or suffer from low sample replications (ABC200; Figure 4), we opt for this combination (ABC300 and FMAM&SO) when producing a formal reconstruction.

Temperature Reconstruction and Benchmarking
Given the overall weak interseries correlation among the MXD data, and the caveats of signal estimation against heterogenous regional climate data, one could argue that a proxy transfer into temperature estimates and publication of a formal reconstruction is not fully warranted. On the other hand, the CNP is perhaps the only location in southern Spain providing annually resolved estimates of temperature variability over the past several centuries. It therefore appeared reasonable to scale the dimensionless MXD index values against instrumental data [73] and produce a reconstruction of FMAM&SO temperature variability reaching back to 1350 CE, though we remind users of this reconstruction to consider the limitations und uncertainties detailed in the previous section. The ABC300 chronology correlates at r = 0.49 over the long 1905-2014 calibration period against mean FMAM&SO temperatures of five instrumental stations. Positive RE and CE values ranging from 0.18-0.39 and 0.15-0.29, respectively, indicate that the reconstruction has some statistical skill. This conclusion is supported by a high Durbin-Watson value of 1.4, indicating that the proxy-target residuals contain no substantial drift throughout the calibration period.
The FMAM&SO temperature reconstruction reveals warm conditions during the early 19th, early 17th, and 16th centuries that were of similar magnitude than warmth recorded since the second half of the 20th century ( Figure 9). The 95% bootstrap confidence limits demonstrate that none of these periods is significantly warmer than any other period, however, and it also remains unclear whether the reconstruction presented here captures the full spectrum of low frequency variance. The latter is concluded because the original CNP MXD data is composed of samples from only living trees, and we here "just" truncated the biologically older rings >300 years to produce a dataset characterized by relatively flat mean age and segment length curves (Figure 4, Supplementary Figure  S4) to meet requirements for RCS detrending [38]. The data truncation reduced the replication of the ABC300 chronology particularly after 1600 CE, but also the early chronology periods before 1584 CE are replicated by fewer than 20 MXD series. These values are overall not impressive if compared with the worldwide best-replicated MXD-based reconstruction from the Spanish Pyrenees integrating >100 MXD series throughout the 16th century [19].

Temperature Reconstruction and Benchmarking
Given the overall weak interseries correlation among the MXD data, and the caveats of signal estimation against heterogenous regional climate data, one could argue that a proxy transfer into temperature estimates and publication of a formal reconstruction is not fully warranted. On the other hand, the CNP is perhaps the only location in southern Spain providing annually resolved estimates of temperature variability over the past several centuries. It therefore appeared reasonable to scale the dimensionless MXD index values against instrumental data [73] and produce a reconstruction of FMAM&SO temperature variability reaching back to 1350 CE, though we remind users of this reconstruction to consider the limitations und uncertainties detailed in the previous section. The ABC300 chronology correlates at r = 0.49 over the long 1905-2014 calibration period against mean FMAM&SO temperatures of five instrumental stations. Positive RE and CE values ranging from 0.18-0.39 and 0.15-0.29, respectively, indicate that the reconstruction has some statistical skill. This conclusion is supported by a high Durbin-Watson value of 1.4, indicating that the proxy-target residuals contain no substantial drift throughout the calibration period.
The FMAM&SO temperature reconstruction reveals warm conditions during the early 19th, early 17th, and 16th centuries that were of similar magnitude than warmth recorded since the second half of the 20th century ( Figure 9). The 95% bootstrap confidence limits demonstrate that none of these periods is significantly warmer than any other period, however, and it also remains unclear whether the reconstruction presented here captures the full spectrum of low frequency variance. The latter is concluded because the original CNP MXD data is composed of samples from only living trees, and we here "just" truncated the biologically older rings >300 years to produce a dataset characterized by relatively flat mean age and segment length curves (Figure 4, Supplementary Figure S4) to meet requirements for RCS detrending [38]. The data truncation reduced the replication of the ABC300 chronology particularly after 1600 CE, but also the early chronology periods before 1584 CE are replicated by fewer than 20 MXD series. These values are overall not impressive if compared with the worldwide best-replicated MXD-based reconstruction from the Spanish Pyrenees integrating >100 MXD series throughout the 16th century [19]. The most striking, and statistically significant feature is the ~1 °C change from exceptionally cold late-18th-century conditions (t1781-1810 = −1.15 °C ± 0.64 °C) to record levels of warmth in the early 19th century (t1781-1810 = −1.15 °C ± 0.64 °C). This rapid temperature increase is centered around the Tambora eruption in 1815 that is marked by cold FMAM&SO temperatures in the 1816 "year without a summer" [74] (t1816 = −2.1 °C ± 0.55 °C, the 14th coldest year since 1350 CE). Whether this change from colder to warmer conditions is related to external climate forcings could perhaps be evaluated by comparison with climate model simulations [75]. The rapid transition from low-to-high temperatures throughout the Dalton Solar Minimum from 1790-1830 [76] suggests, however, that this lower frequency temperature change is not related to solar forcing. Before this period, and particularly before 1600 CE, the bootstrap confidence limits (red curves in Figure 9) markedly increase, demonstrating that the FMAM&SO temperature estimates are less reliable back in time (Supplementary Figure S7).
As with the early 19th century temperature shift, the annual cold extremes reconstructed over the past 650 years (blue dots in Figure 9) appear to be only partly related to external climate forcings, as only two large volcanic eruptions in 1600 (Huaynaputina in Peru) and 1452 (unknown) [77], recorded in bipolar ice core records [51], coincide with substantial cooling in subsequent years: t1601 = −2.4 °C ± 0.74 °C is seventh coldest year and t1453 = −1.8 °C ± 1.32 °C is the tenth coldest year since 1350 CE (ranks after high-pass filtering the data). Interestingly, the cold extremes are fairly evenly distributed throughout the past 650 years, whereas several of the warm extremes are clustered between the late 15th and early 16th centuries. This concentration of warm events is linked to an increase in interseries correlation (Rbar, Supplementary Figure S7a) enhancing the high-frequency variability of the mean chronologies during this period. We did not remove this variance surplus using adjustment techniques [78], as it appeared unrelated to structural changes in the underlying data. This variance could reflect a period during which the climate was potentially more variable. The underlying forcing for this high-variance period in the CNP reconstruction remains unknown, however.
Comparison of the CNP MXD record with warm season temperature reconstructions from a MXD network in the Spanish Pyrenees [7,16,19,20] and a previous-year September-October temperature reconstruction from CNP TRW data [27] reveals substantial covariances over the past 650 years (Figure 10). The correlations against the Pyrenees data, from sites located ~600 km north of The most striking, and statistically significant feature is the~1 • C change from exceptionally cold late-18th-century conditions (t 1781-1810 = −1.15 • C ± 0.64 • C) to record levels of warmth in the early 19th century (t 1781-1810 = −1.15 • C ± 0.64 • C). This rapid temperature increase is centered around the Tambora eruption in 1815 that is marked by cold FMAM&SO temperatures in the 1816 "year without a summer" [74] (t 1816 = −2.1 • C ± 0.55 • C, the 14th coldest year since 1350 CE). Whether this change from colder to warmer conditions is related to external climate forcings could perhaps be evaluated by comparison with climate model simulations [75]. The rapid transition from low-to-high temperatures throughout the Dalton Solar Minimum from 1790-1830 [76] suggests, however, that this lower frequency temperature change is not related to solar forcing. Before this period, and particularly before 1600 CE, the bootstrap confidence limits (red curves in Figure 9) markedly increase, demonstrating that the FMAM&SO temperature estimates are less reliable back in time (Supplementary Figure S7).
As with the early 19th century temperature shift, the annual cold extremes reconstructed over the past 650 years (blue dots in Figure 9) appear to be only partly related to external climate forcings, as only two large volcanic eruptions in 1600 (Huaynaputina in Peru) and 1452 (unknown) [77], recorded in bipolar ice core records [51], coincide with substantial cooling in subsequent years: t 1601 = −2.4 • C ± 0.74 • C is seventh coldest year and t 1453 = −1.8 • C ± 1.32 • C is the tenth coldest year since 1350 CE (ranks after high-pass filtering the data). Interestingly, the cold extremes are fairly evenly distributed throughout the past 650 years, whereas several of the warm extremes are clustered between the late 15th and early 16th centuries. This concentration of warm events is linked to an increase in interseries correlation (Rbar, Supplementary Figure S7a) enhancing the high-frequency variability of the mean chronologies during this period. We did not remove this variance surplus using adjustment techniques [78], as it appeared unrelated to structural changes in the underlying data. This variance could reflect a period during which the climate was potentially more variable. The underlying forcing for this high-variance period in the CNP reconstruction remains unknown, however.
Comparison of the CNP MXD record with warm season temperature reconstructions from a MXD network in the Spanish Pyrenees [7,16,19,20] and a previous-year September-October temperature reconstruction from CNP TRW data [27] reveals substantial covariances over the past 650 years ( Figure 10). The correlations against the Pyrenees data, from sites located~600 km north of the CNP, are all significant (p < 0.01) and increase when removing lower frequency variance from the timeseries using 10-year spline filters (the right panel in Figure 10a). The values are statistically indistinguishable among the four Pyrenees records and slightly increase from a minimum of r = 0.41 calculated over the full 1350-2005 period of overlap to a maximum of r = 0.46 calculated of the most recent period since 1901. Similarly, the correlations with the (reversed) CNP TRW-based reconstruction also increase towards present, though the change is much larger reaching from r 1350-2005 = 0.39 to r 1901-2005 = 0.66. Whereas the covariance between the CNP and Pyrenees MXD records declines when low-pass filtering the data (not shown), the CNP MXD and TRW records also correlate significantly at r 1350-2005 = 0.45 after smoothing the timeseries (Figure 10b). Note, however, that the TRW-based reconstruction was reversed and shifted by one year in these comparisons, as the correlations were otherwise close to zero (r < 0.1).
Atmosphere 2020, 11, x FOR PEER REVIEW  13 of 18 the CNP, are all significant (p<0.01) and increase when removing lower frequency variance from the timeseries using 10-year spline filters (the right panel in Figure 10a)  (Figure 10b). Note, however, that the TRW-based reconstruction was reversed and shifted by one year in these comparisons, as the correlations were otherwise close to zero (r < 0.1). Whereas the high covariance of the CNP FMAM&SO temperature reconstruction with the highly replicated MXD network from the Pyrenees is reconfirming that the reconstruction presented here contains useful climatic information (albeit the large uncertainty limits), the lacking coherence with the original TRW-based reconstruction from CNP Pinus nigra requires further explanation. After reversing the TRW-based reconstruction, we effectively demonstrate that MXD and TRW data from the same CNP pine forest correlate substantially. The correlation increases towards the present, likely affected by the replication changes of the underlying TRW and MXD data. This temporal change in co-variance is much smaller, however, when comparing the MXD-based reconstructions from the CNP and Pyrenees, suggesting a stronger temporal weakening of the signal in TRW. The lacking correlation between the MXD-based and original (non-reversed and shifted) TRW-based CNP reconstructions is in line with expectations based on the correlation of instrumental previous-year late summer temperatures (targeted by the TRW reconstruction) and current-year FMAM&SO temperatures (targeted by the MXD reconstruction), which are also insignificant. Disentangling this complicated association between previous-year and current-year climate signals recorded in CNP TRW and MXD data appears challenging, and we are currently in the process of developing a large earlywood/latewood width dataset from hundreds of the CNP black pines to address this conundrum. Whereas the high covariance of the CNP FMAM&SO temperature reconstruction with the highly replicated MXD network from the Pyrenees is reconfirming that the reconstruction presented here contains useful climatic information (albeit the large uncertainty limits), the lacking coherence with the original TRW-based reconstruction from CNP Pinus nigra requires further explanation. After reversing the TRW-based reconstruction, we effectively demonstrate that MXD and TRW data from the same CNP pine forest correlate substantially. The correlation increases towards the present, likely affected by the replication changes of the underlying TRW and MXD data. This temporal change in co-variance is much smaller, however, when comparing the MXD-based reconstructions from the CNP and Pyrenees, suggesting a stronger temporal weakening of the signal in TRW. The lacking correlation between the MXD-based and original (non-reversed and shifted) TRW-based CNP reconstructions is in line with expectations based on the correlation of instrumental previous-year late summer temperatures (targeted by the TRW reconstruction) and current-year FMAM&SO temperatures (targeted by the MXD reconstruction), which are also insignificant. Disentangling this complicated association between previous-year and current-year climate signals recorded in CNP TRW and MXD data appears challenging, and we are currently in the process of developing a large earlywood/latewood width dataset from hundreds of the CNP black pines to address this conundrum.

Conclusions
Ancient black pine trees from the Cazorla Natural Park (CNP) in southern Spain were used to develop one of the worldwide southernmost MXD-based temperature reconstructions produced so far. The underlying measurement series of this reconstruction were truncated at a biological age >300 years to support the application of a detrending method (RCS) capable of preserving low-frequency variability in the resulting index chronology. The reconstruction derived from this chronology is accompanied by large uncertainties arising from (i) relatively feeble calibration statistics (e.g., 25% explained FMAM&SO temperature variance), and (ii) a decline in sampled replication from >35 series in the 20th century to <18 series before 1500 CE. Additional assessments of the instrumental station records in the study region showed, however, that the calibration statistics were not only controlled by uncertainties inherent to the tree-ring proxy, but that the long distances of reliable instrumental data and inhomogeneities among these timeseries additionally affect the reconstruction skill estimates.
The final FMAM&SO temperature reconstruction shows relatively limited centennial scale variability, such as a prolonged Little Ice Age and transition into warmer 20th-century conditions, as reported from sites in central and northern Europe [35,59,60]. Instead, the most striking feature is a rapid warming trend from the late 18th to the early 19th century, subsequent to a gradual cooling trend throughout much of the 17th and 18th centuries. The influence of external (solar and volcanic) climate forcings appears to be of minor importance, as the high-to-low frequency variance of the new reconstruction does not cohere particularly well with prominent solar minima and large volcanic eruptions over the past 650 years. The MXD-based temperature reconstruction from the CNP represents a new benchmark for high-resolution, pre-instrumental climate variability in the southwestern Mediterranean region, though users of this record are reminded of the large uncertainties, particularly before 1600 CE.