A quantitative paleoclimatic reconstruction of the non-analogue environment of oxygen isotope stage 3: new data from small mammal records of southwestern Germany

Ensuring comparability between results is a key goal of all paleoecological reconstructions. Quantitative estimates of meteorological variables, as opposed to relative qualitative descriptions, provide the opportunity to compare local paleoenvironmental records against global estimates and incrementally build regional paleoclimatic records. The Bioclimatic Method provides quantitative and qualitative estimates of past landscape composition and climate along with measures of statistical accuracy by applying linear discriminant functions analysis and transfer functions to faunal taxonomic abundance data. By applying this method to the rodent data from Geißenklösterle and Hohle Fels, two Paleolithic cave sites located in the Ach Valley of southwestern Germany, we classify the regional vegetation according to Walters’ zonobiome model. We also present new estimates of meteorological variables including mean annual temperature, mean annual precipitation, and vegetative activity period of the Ach Valley for the period spanning ~ 60,000 to 35,000 cal BP. The results suggest the Ach Valley contained a non-analogous landscape of arctic tundra and temperate deciduous woodland with occasional arid steppe expansion. Meteorological estimates suggest the climate was significantly colder during the Middle and Upper Paleolithic than today, with higher annual precipitation and dramatic temperature shifts between seasons. These results fit well with climatic reconstructions from Switzerland and the Netherlands based on a variety of proxies. They also provide further evidence of a localized climatic response within southwestern Germany to the stadial-interstadial shifts preceding the Heinrich 4 event. Finally, these results reinforce our previous claims that climatic volatility was not a driving force in the loss of Neanderthal groups throughout the Swabian Jura during OIS 3.


Introduction
The last 20 years of archaeological research have seen a proliferation of efforts at paleoenvironmental reconstruction and climate modelling. Much of this work has been focused on the Last Glacial period in an attempt to identify terrestrial records reflecting the drastic climatic instability documented in Arctic ice cores during this time (North Greenland Ice Core Project members Rasmussen et al. 2014). The varying success of these studies has illuminated the stochastic nature of local ecological response to continental climatic change across Europe. Various factors may affect vegetative composition at a regional level, including the presence of geographic barriers to migration (Müller et al. 2003) and the presence of biomass suppressing consumers such as megaherbivores, humans and fire (Bond 2005). As such, robust terrestrial records of local floral and faunal community structure during OIS 3 (~ 60,000-25,000 BP) are required to tease out the role that regional characteristics played in mitigating the effect of continental climatic events. Such records are also needed to assess the continued applicability of environmentally deterministic models of Neanderthal extinction and early modern human expansion throughout Central Europe (Discamps and Royer 2017;Rhodes et al. 2019). Terrestrial paleoenvironmental records are based on a variety of proxy data sources including stable isotopes, pollen spectra, mammalian community structure, lake sediments, and speleothems. These proxies provide qualitative and/or quantitative climatic data at different scales (i.e. site, regional, and/or continental) covering different time frames (i.e. isolated events, cumulative decades, palimpsests of varying or unknown duration). Unfortunately, many of these records are limited to dichotomous relative categories when describing past climates (warm vs. cool, wet vs. dry). The application of quantitative statistical methods from other fields within the natural sciences has the potential to increase the resolution and precision of our paleoclimatic data, as well as its comparability across space and time (Lyman 2016).
Paleoecological transfer functions have been used for decades to measure past atmospheric and oceanic conditions from fossil planktonic foraminifera, variation in tree ring widths, and pollen assemblages derived from lake sediments (Sachs et al. 1977). The commonality between these applications lies in the use of calibrated functions derived from multivariate regression analysis of modern multi-species (or multi-tree ring width) distributional data. These transfer functions are applied to older samples to estimate environmental parameters (Sachs et al. 1977). Transfer functions were perhaps most famously used by the CLIMAP group (CLIMAP 1976(CLIMAP , 1981 in their study of ocean surface temperatures reaching back to the Last Glacial Maximum. The Quantitative Bioclimatic Model (QnBM) (Hernández-Fernández and Peláez-Campomanes 2005) uses transfer functions based on modern distributions of groups of terrestrial fauna to estimate past vegetative biome presence and various meteorological phenomena such as mean annual temperature (MAT), mean annual precipitation (MAP), and vegetative activity period (VAP). It builds upon the Qualitative Bioclimatic Model (QlBM) which uses weighted averages and linear discriminant function analysis to determine the dominant vegetative biome(s) of a locality based on the associated mammal fauna (Hernández-Fernández 2001). The QlBM reconstructs past landscape vegetation and climate following Heinrich Walters' international climatic zonobiomes model (Walters 1970). In this way, the QlBM is most similar to the Habitat Weighting Method (Evans et al. 1981;Andrews 1996) but utilizes well-known, comparable, and standardized terminology (Table 3). Furthermore, cross-validation studies using modern data from transitional ecotones (Hernández-Fernández and Peláez-Campomanes 2005) suggests that the QnBM can produce accurate estimates past non-analogue environments.
Despite the potential strength of this method in reconstructing past climatic variables, it has been sparsely applied to contexts outside of Spain. However, since 2005, a number of small mammal researchers have used the Bioclimatic Method to estimate the local paleoecology of Paleolithic sites from Portugal, France, Italy, Germany, Belgium, and Poland (Fernández et al. 2007;Gąsiorowski et al. 2014;Socha 2014;López-García 2015, 2016, 2017a, 2017b, 2018a, 2018bPinero et al. 2016;Rey-Rodríguez et al. 2016;Discamps and Royer 2017;Berto et al. 2018Berto et al. , 2019Fagoaga et al. 2018;Wong et al. 2017;Royer 2020;Wong et al. 2020a). In this paper, we apply both the Qualitative Bioclimatic Model (QlBM) and the Quantitative Bioclimatic Model (QnBM) to recently complied rodent data from the sites of Geißenklösterle and Hohle Fels in the Ach valley of the Swabian Jura (Rhodes et al. 2018Rhodes 2019).
The Swabian Jura is an important region with regards to the Middle to Upper Paleolithic transition of Central Europe due to its numerous well-stratified archaeologically rich cave sites. Previous paleoecological reconstructions suggest this valley experienced a gradual cooling of the climate and an associated extension of arctic tundra environments beginning in the late Middle Paleolithic (MP) and continuing well into the Upper Paleolithic (UP) (Miller 2015;Riehl et al. 2015;Rhodes et al. 2018Rhodes et al. , 2019. Furthermore, they suggest that the valley may have had a non-analogous landscape pattern in the past (Rhodes, 2019). As such, this study provides a rigorous test of the Bioclimatic Models ability to work with complex records from Paleolithic sites. This is the first study to derive quantitative measures of meteorological phenomena from the Middle Paleolithic (MP) and Upper Paleolithic (UP) in the Ach Valley, extending the precision and comparability of the already robust paleoecological record of this region (Niven 2006;Krönneck 2012;Kitagawa 2014;Miller 2015;Riehl et al. 2015;Wong et al. 2017Wong et al. , 2020aRhodes et al. 2018Rhodes et al. , 2019Baumann et al. 2020;Starkovich et al. in press).

Site description and materials
The Swabian Jura is a karstic plateau of Jurassic limestone in southwestern Germany bounded by the Neckar Valley to the north and the Danube Valley to the south. It contains three river valleys of archaeological importance, the Ach Valley, the Lone Valley, and the Lauchert Valley (Fig. 1). Hohle Fels Cave sits at 534 m a.s.l. and faces northwest overlooking the Ach Valley. The cave comprises a 29-m-long corridor (where modern excavation is focused) and a 23-m-long grand hall, both of which contain sediments accumulated primarily through a now closed chimney at the back of the cave. Archaeological investigations at this site have been ongoing since 1977, first under the direction of Joachim Hahn and subsequently by Nicholas Conard, both of the Universität Tübingen. Combined, these excavations have produced a 5-m-thick stratigraphic record spanning the Middle Paleolithic through to the Holocene which is divided into geological horizons (GHs) designated by Arabic numerals. This study uses previously published data on the rodent taxonomic composition of the MP and UP periods from two 1 m 2 excavation units : Table 2) which is detailed in Table 1. The excavation and recovery methods are published in detail elsewhere (Conard et al. 2013;Rhodes 2019) and include water screening of bulk sediment and sorting of the bone material by trained volunteers.
Geißenklösterle Cave is also located in the Ach Valley and was first excavated by Eberhard Wagner in 1973. Campaigns directed by Joachim Hahn from 1974 to 1991 and Nicholas Conard from 2000 to 2002 produced the material used in this study and previous studies identifying the small mammal taxonomic composition and taphonomic history (Rhodes et al. 2018;Ziegler 2019). The MP and UP deposits at Geißenklösterle are interstratified by a nearly culturally sterile geogenic horizon (GH 17) and document occupation of the cave by Neanderthal and anatomically modern human groups spanning the period from ~ 48,000 to 33,000 cal BP. The rodent remains which come from three 1 m 2 excavation units were collected using the same recovery methods applied at Hohle Fels and the previously published taxonomic composition of the material (Rhodes et al. 2018: Table 2) is detailed in Table 2. Along with Sirgenstein, Brillenhöhle, Vogelherd, Holenstein-Stadle, and Bockstein caves, these sites have produced some of the earliest and most impressive evidence of UP symbolic behaviour including works of portable figurative art and musical instruments (Conard 2009;Conard et al. 2009). The long record of human occupation, well stratified deposits, and large sample sizes of small mammal remains available from both sites are ideal for highresolution diachronic recreations of the paleoenvironment of the Ach Valley.

Qualitative bioclimatic model (QlBM)
The QlBM is an ordination method used to estimate the probable zonobiome (herein referred to as biome)  characterizing any locality as defined within Walter's (1970) typology. It is based on the assumption that faunal community structure is shaped by both climatic variables and the associated local vegetative properties (Hernández Fernández 2001a, b). Heinrich Walter's zonobiome classification is well suited to this method as it incorporates both climate and vegetation type (Breckle and Walter 2002). Cross-validation of the model using different mammal groups (all mammals, macromammals, micromammals, and micromammals and bats) against world-wide modern localities (n = 50) with known biome zones has shown that micromammal groups (rodents excluding insectivores and bats) predict the biome with 97% accuracy, second only to complete mammal assemblages (Hernández Fernández and Peláez-Campomanes 2003). We applied the QlBM to the most recent taphonomically assessed rodent records from Hohle Fels and Geißenklösterle (Rhodes et al. 2018Rhodes, 2019), as complete faunal data was not yet available. We were able to reproduce Hernández Fernández and Peláez-Campomanes' (2003) cross validation results and confirm that of the 50 modern locality spectra used only three were miscategorized. In these three cases the second highest likely biome categorization matched the actual modern biome. As such, we feel confident assigning a 94% accuracy to the results of the QlBM biome predictions for the archaeological spectra, however the method may produce results which are actually closer to 97% accurate, as suggested by Hernández Fernández and Peláez-Campomanes (2003). Furthermore, by using only the rodent records we are able to apply the QlBM to each (GH) individually and derive a diachronic record of climatic change through time. We assigned each rodent species a value for each biome based on its modern geographic distribution. A species is considered to occupy a biome if it encompasses 15% or more of the species modern geographic range (Hernández-Fernández 2001). If modern populations of a species do not occupy a particular biome the value is 0; if the species does occupy the biome then the value is 1/n with n being the total number of biomes in which the species is present. Hernández-Fernández refers to this as the climatic restriction index (CRI) and a detailed description of the process used to determine the biomes inhabited by a species can be found in Hernández Fernández (2001a, b). The CRI is analogous to the weightings of the Taxonomic Habitat Index (Evans et al. 1981) and the Habitat Weighting Method (Andrews 1996), but differs in that a species score is evenly distributed across all biomes, even those in which the species is only rarely sighted (Hernández Fernàndez 2001). The CRI values for the species identified from Hohle Fels and Geißenklösterle are included in Supplementary Material Table 1 and 3, respectively. Our CRI values follow earlier studies (Hernández Fernández 2001a, b) with the exception of Lemmus lemmus, which is known to inhabit the same biome as Dicrostonyx torquatus/gulielmi (van Kolfschoten 2014; Rhodes et al. 2018). Specimens attributed to the family level only (i.e. Gliridae and Muridae) were excluded from our analysis as the species within these familes inhabit very divergent climatic zones including savanna, grassland, tundra, desert, and most forest biomes. Without more precise genus or species level identifications, these specimens are unlikely to provide any meaningful data on past climate or landscape composition.
The CRI values are then used to calculate the total representation of each biome within a locality, known as the bioclimatic component (BC) using the following formula: where i is the biome and S is the total number of species within that locality (Hernández Fernàndez 2001). In this way, a total weighted average of the representation of each biome is calculated for each locality or, as in the current study, for each geological horizon. The BC values for each horizon are included in Supplementary Material Table 2 and 4.
We input the bioclimatic components of each GH and aggregated spectra for the complete Middle Paleolithic, Aurignacian, and cultural hiatus periods (a total of 24 spectra) in SPSS version 24 along with the spectra of 50 modern localities compiled by Hernández Fernández and Peláez-Campomanes (2003). The biome categorization of the 50 modern localities was used as the independent variable to teach the linear discriminant functions analysis to classify the highest likely and second highest likely biome categorization for the archaeological spectra.

Quantitative bioclimatic model (QnBM)
Although not described as such by Hernández Fernández and Peláez-Campomanes' (2003), the QnBM utilizes transfer functions to extrapolate meteorological variables from the bioclimatic spectra of a locality. These transfer functions are based on the relationship between modern climate variables and terrestrial faunas (Hernández-Fernández and Peláez-Campomanes 2005) and have been derived by applying multivariate linear regression to a world-wide sample (n = 50) of mammal faunas from known climates with associated modern meteorological data. Hernández Fernández and Peláez-Campomanes (2005) then cross-validated this method using a second group of modern faunas and applied the transfer functions to a Pleistocene assemblage, the results of which were found to be similar to the palynological record from the same site (Hernández-Fernández and Peláez-Campomanes 2005). This process also produced measures of standard error and coefficients of determination for each meteorological variable estimated. Like in the QlBM, complete mammalian faunas and assemblages of rodents (without insectivores and bats) produced the most accurate results (Hernández- Fernández and Peláez-Campomanes 2005). We used transfer functions for rodent assemblages (Hernández Fernández and Peláez-Campomanes 2003) to derive estimates of past mean annual temperature (MAT), mean annual precipitation (MAP), the mean temperature of the warmest month (MTW), the mean temperature of the coldest month (MTC), and an estimate of the vegetative activity period (VAP; broadly defined as the number of months in which temperatures exceed 7 °C). A basic linear regression formula was applied to the QlBM bioclimatic spectra (with roman numerals indicating different biomes following Table 3): [predicted = B+(IV*aIV)+(VI*aVI)+(VII*aVII)+(VII I*aVIII)+(IX*aIX)] Using the slope (aIV, aVI, aVII, aVIII, and aIX) and intercept (b) values provided by Hernandez-Fernández and Peláez-Campomanes (2005) for each meteorological variable (also included in Supplementary Materials Table 2 and 4). Hernandez- Fernández and Peláez-Campomanes (2005) found that MAT, MTC, and VAP were strongly correlated with rodent fauna spectra (r 2 = > 0.930) while MAP and MTW have weaker coefficients of determination (r 2 = 0.746). These, along with the standard error of the estimates and other relevant variables are reported in Supplementary Materials Table 2 and 4. Our estimates of MAT, MTW, and MTC are reported degrees Celsius, MAP in millimeters, and for VAP as number of months.

Qualitative bioclimatic model results
Overall, the bioclimatic spectra from both Hohle Fels and Geißenklösterle are broadly similar (Tables 4 and 5).
Three biomes dominate all horizons-typical temperate with deciduous forests (VI), cold-temperate boreal climate (VIII), and arctic tundra (IX). In most horizons, arctic tundra dominates the distribution. This is particularly true at Hohle Fels, where only the oldest Middle Paleolithic horizon, GH 12 (> 60,000 BP), has temperate deciduous forest spectra values exceeding those of biome IX. Our linear discriminate function analysis also suggests that GH 12 and the Middle Paleolithic as an aggregated assemblage at Hohle Fels would be more correctly classified as biome VI-typical temperate with deciduous forests (Table 4). This is likely driven by the presence of the pine vole (Microtus subterraneous) within these assemblages. Otherwise, arctic tundra is the most likely biome represented by the small mammal assemblages from Hohle Fels according to the discriminate function analysis, with temperate and deciduous forests being the second most likely.
This pattern is also seen at Geißenklösterle, with arctic tundra landscapes dominating the bioclimatic spectra and results of the discriminate functions analysis of most assemblages. Geological horizons 11, 19 and 22 are exceptions (Table 5). In the Aurignacian (GH 11) the exceptionally low sample size within this horizon (NISP = 20) is likely skewing the analysis towards a deciduous forest signal, partially due to an underrepresentation of lemming species relative to other periods. In GHs 19 and 22, however, it would appear that the presence of the pine vole (Microtus subterraneous) and the loss of the tundra vole (Microtus oeconomus) respectively, are driving a more forested signal. Although not strong enough to change the discriminate functions results, the presence of the ground squirrel (Spermophilus superciliosus) is clearly reflected in the high values of biome VII-arid temperate steppe landscape within GH 18, 22 and the aggregated Middle Paleolithic assemblage. At Hohle Fels, similarly high values for biome VII in GH 8, 12, and the aggregated MP assemblage are also driven by the presence of the ground squirrel. As in previous studies (Rhodes et al. 2018, this species is a key indicator of the presence of an otherwise rarely denoted biome -the arid steppe.

Quantitative bioclimatic model results
Our reconstructed meteorological results also reflect the dominance of the arctic signal throughout the period studied. At Hohle Fels, the reconstructed MATs range from 7.6 °C to 11.1 °C colder than the Baden-Württemberg historical average (Table 6). At GK our MATs are equally as disparate, ranging from 8.2 °C to 11.7 °C colder than the historical average ( Table 7). The MAP results indicate lower precipitation amounts throughout both sequences, with only GH 19 at Geißenklösterle and GH 7 at Hohle Fels exceeding the historical average by modest amounts. However, the coefficient of determination for the MAP regression model is 0.746, which suggests that a substantial portion of the variation in this measure is due to factors other than the bioclimatic spectra values. Hernandez Fernández and Peláez-Campomanes (2005) suggest that the fact that precipitation influences mammal faunas indirectly (i.e. through vegetation) may explain this low r 2 score. This measure may also be affected by soil water retention and other sedimentary characteristics.
Our mean average warmest month (MTW) temperatures indicates there may have been a gradual increase in summer temperatures from the MP through the early UP of ~ 3 °C at both Hohle Fels and Geißenklösterle. Overall, this meteorological variable appears to mirror shifts in MAT, with periods of increased average cooling (i.e. GH 7, 7a/7aa, 10, and 11 at Hohle Fels and GH 13, 19, 20, and 21 at Geißenklösterle) also experiencing increased summer temperatures (Figs. 2 and 3, respectively). The MTC results suggest that winter temperatures were more stable through time in the Ach Valley. At Hohle Fels, after a drop from an average of -10 °C to -14 °C in the MP (GH 12 to 11), winter averages remain steady with only minimal changes in the earliest Aurignacian period (Fig. 1). At Geißenklösterle two shifts in MTC have been identified, the first between GH 20 and GH 19 in the MP when winter temperatures rose by just over 3 °C and the second between GH 13 and GH 12 in the Aurignacian when a period of gradual cooling is interrupted by a 2 °C increase in average temperature during the coldest period of the year (Fig. 2). Interestingly, the warmer winter signal in GH 19 is accompanied by estimates of an increase Table 4 Bioclimatic Spectra for Hohle Fels by geological horizon. The yellow bars illustrate the % or loading that each climatic zone holds out of the total spectra (100). Linear Discriminate Function 1 st and 2 nd biome classification predictions are included Table 5 Bioclimatic spectra for Geißenklösterle by geological horizon. The yellow bars illustrate the % or loading that each climatic zone holds out of the total spectra (100). Linear Discriminate Function 1st and 2nd biome classification predictions are included in annual precipitation of almost 200 mm. Both of these meteorological shifts appear to have been short-lived.
The final meteorological measure reproduced is the number of months of vegetative activity or productivity (VAP). This variable measures the number of spring and early summer growing months and, conversely, the length of the vegetative dormancy period. As the length of these periods is directly related to temperature, this measure can also be described as the number of months during which temperatures exceed 7 °C (Gavilán and Fernández González 1997). Hernandez Fernández and Peláez-Campomanes (2005) report a strong determination coefficient (r 2 = 0.955) for this meteorological variable and rodent faunal taxonomic composition. At Hohle Fels, this period ranges from 2.4 to 3.7 months, while at Geißenklösterle VAP ranges from 2.3 to 3.6 months. This suggests that in both records the period of vegetative dormancy may have lasted between 8 and 10 months of the annual year.

Discussion
Based on the results of the QlBM and QnBM, it is clear that the Paleolithic Ach Valley landscape can be characterized by more than a single modern biome at any given time. At Geißenklösterle, both the MP and UP were largely arctic tundra landscapes, however a signal for temperate deciduous forests is also present (primarily as the 2nd most predicted group by the discriminant functions analysis). Twice in the MP record arid steppe environments (biome VII) encroach on this pattern. Cold-temperate boreal forests (biome VIII) appear around Geißenklösterle in the early MP and late early UP, suggesting significantly colder temperatures than during previous periods. The paleoecological record at Hohle Fels exhibits greater stasis, with arctic tundra environments extending from the MP through to the UP. However, the dominance of deciduous forests (biome VI) as the 2 nd predicted grouping and the boreal forest (biome VIII) signal in the aggregated MP assemblage suggest the presence of a boreo-nemoral ecotone (Breckle and Walter 2002) nearby.
Consistency in our paleoecological records from the region is clear, and the QlBM results described herein broadly match the results from other proxy material studies. The dominance of arctic tundra landscapes is documented in other small mammal (Rhodes et al. 2018Ziegler 2019) and botanical (Riehl et al. 2015) signals, while the presence of lemmings and reindeer in both the Hohle Fels and Geißenklösterle assemblages supports the presence of arctic climates in the valley during this time (Starkovich et al. in press;Rhodes 2019). Our study also evidences the  (Beck et al. 2016;Viñas et al. 2016) and both boreal and deciduous forest environments are also suggested by the presence of shrews (Rhodes et al. 2018. The periodic encroachment of arid steppe vegetation during the MP at Geißenklösterle and MP and early UP at Hohle Fels may indicate the presence of a more typical steppe-tundra landscape (Chytrý et al. 2019) akin to the well-known Mammoth Steppe (Guthrie 2001(Guthrie , 2013. However, the lack of a consistent signal for steppe vegetation in the Ach Valley throughout the period studied suggests this landscape may have differed from the classical description of the mammoth Steppe in the past. Therefore, the Ach Valley was likely not a homogeneous arctic tundra climatic zone, but rather a transitional landscape between deciduous and coniferous forests, similar to the boreo-nemoral ecotone VI/VIII (Breckle and Walter 2002), or between boreal forest and treeless tundra [i.e. the forested tundra ecotone VII/IX (Breckle and Walter 2002)] at different points in time. In this way, our results are similar to other applications of the QlBM in France (Discamps and Royer 2017) and Belgium (López-García et al. 2017b) from the same time period. The periodic presence of arid steppes documented in the Geißenklösterle MP record also fits with this interpretation, as arid steppes are known to border boreo-nemoral ecotones in modern regions with continental climates (Breckle and Walter 2002;pg. 414). The cyclical expansions of boreal forest and steppe biomes during interstadials indicated by pollen record of the Bergsee (Duprat-Oualid et al. 2017) and at Le Grand Pile in France (Woillard 1978;De Beaulieu and Reille 1992) provide further support to this interpretation. As such, the Ach Valley landscape during this time fits well with Guthries' revised description of the Mammoth Steppe when he stated "whatever it was it was a mosaic" (2001; pp. 565).
In terms of climate, the Ach Valley was colder and dryer than modern averages during both the MP and UP. There is also a clear pattern of decreasing VAP through time alongside a slight increase (~ 3 o C) in MTW. However, it is important to note that Hernández Fernández and Peláez-Campomanes, (2005) suggest a minimum of 100 specimens Fig. 3 Geißenklösterle meteorological estimates by geological horizon (GH). Standard error shown in grey dashes and the modern MAT indicated by the yellow bar is necessary to produce reliable past climatic estimates, and this threshold has not been met within many of the horizons studied, particularly in the Aurignacian layers, at both Geißenklösterle and Hohle Fels. The degree to which sample size has affected the representation of rare species in both assemblages has been quantified elsewhere (Rhodes et al. 2018 and suggests that between 2 and 3 species are underrepresented or absent from most UP horizons. As such, the Aurignacian results must be interpreted cautiously and their implications on our understanding of the past landscape and climate of the Ach valley considered critically. Compared to modern measures, our MAT estimates for the Ach Valley are similar to modern coastal arctic tundra temperatures from Sweden and the boreal zone of northern Europe (i.e. Archangelsk at 0.4 o C) and Siberia (i.e. Irkutsk at -1.3 o C) (Breckle and Walter 2002). The boreal zonobiome is characterized by long, cold winters and short summers (with associated short VAP periods) and extreme temperature fluctuations throughout the year. This fits with our VAP results which indicate vegetation growth during 2-4 months per year throughout the entire sequence. This short VAP period is also seen in arctic tundra regions (Breckle and Walter 2002). Our estimated precipitation amounts are between 4 and 5.5 times higher that of modern arctic continental values (Breckle and Walter 2002) and are similar to modern precipitation amounts from beech forests in Luxembourg (MAP 739 mm) (Breckle and Walter 2002).
Comparing our results to other European applications of the QnBM we still see that our raw estimates are substantially different. For example, in Belgium, MAT during OIS 3 was found to differ from modern averages by between 5 and 6 °C at Caverne Marie Jeanne (López-García et al. 2017b) and by 10.7 °C at Scladina Cave (López-García et al. 2017a). Similar to our study, these MAT estimates are all colder than modern MATs. In Poland, Socha (2014) reports estimates of MAT ranging between 1 °C and 5.5 °C. Our results are substantially lower, with MAT ranging from 0 °C to − 3.1 °C. They are also substantially colder than results from Mediterranean locals like Fumane (López-García et al. 2015); however, this is to be expected.
Yet, our MTC results do correlate with estimates from other material records specifically those of Niederweningen, in Switzerland, where pollen, plant, macrofossils and beetle remains suggest MTC temperatures between -20 o C and − 5 o C, and at Gossau, where lake peat deposits reconstruct MTC estimates between − 15 o C and − 7 o C (Heiri et al. 2014). Both of these records come from the lowlands of the Swiss alps, so altitudinal influence is possible but not extreme in its effects. Our MTW results, however, far exceed the range of 8-13 o C reported from the same Swiss proxy records. In the Netherlands, sedimentary and paleobotanical records suggest MTC temperature ranging from − 11.5 to -13 o C (Huijzer and Vandenberghe 1998), which are also consistent with our findings. Overall, we get a picture of a much colder and moister environment with short spring-summer periods and long winters with dramatic drops in the daily temperature. Despite being somewhat contradictory, our meteorological estimates do fit what could be expected for a forest tundra ecotone (VIII/IX). The presence of deciduous forests on the landscape during most periods is made possible by the low minimum winter temperatures endured by deciduous forests, which range from between -17 o C and -12 o C.
Looking at our climatic estimates through time, we see a number of minor shifts in the meteorological variables which follow the general pattern described above by the climatic zone proportions. At Geißenklösterle we see a decrease in MAT in the mid-MP (GH 21 to 20) followed by warmer MTC estimates in the late MP (GH 19; Fig. 3). Combined with the increased presence of woodland and steppe environments in GH 19, this may indicate a shift from a stadial event to an interstadial period (GI 11). We also see a decrease in MAT and VAP in GH 13 preceding a return to previous MAT levels in GH 12, which following the radiocarbon chronology for the site (Higham et al. 2012) may signal the height of Heinrich 4 (Rhodes 2019). At Hohle Fels there is a similar signal of two particularly cold events indicated by decreased MAT and MTC in GH 11 and MAT and VAP in GH 7 (Table 6). That these two periods may correspond to the end of the Heinrich 5 and 4 events, respectively, fits well with the most recently calibrated radiocarbon dates from this sequence (Bataille and Conard 2018).
The magnitude of the centennial to millennial scale climatic changes of OIS 3 are well documented at northern latitudes, although the degree to which this is reflected in the Ach Valley climatic record and other terrestrial records is arguable (Stringer et al. 2003;Dalen et al. 2012;Staubwasser et al. 2018). Greenland ice core data suggests Dansgaard-Oeschger (D-O) events are associated with fluctuations of up to 15 o C (Johnsen et al. 1992 in Arppe andKarhu 2010) although central European enamel-derived 18 O estimates tend to record changes closer to 5-7 °C between cold and warm climate states (Arppe and Karhu 2010). Within the Ach Valley, the greatest fluctuation in MAT is 2.2 o C, occurring in the MP (GH 12 to 11) at Hohle Fels. While the standard error of our MAT estimate (3.6 °C) does bring this and other fluctuations into the lower range of temperature changes from terrestrial records attributed to continental D-O shifts, other applications of the QnBM have produced records comprised of similarly low-magnitude temperature fluctuations (López-García et al. 2015, 2017b, 2018a, which may reflect the fact that rodent communities are responding to changes in vegetation and snow cover more directly than shifts in temperature or precipitation. However, the lowerthan-expected amplitude of meteorological change during this period in the Ach Valley may also be due to various characteristics of the rodent assemblages themselves. In particular, the a) mode of accumulation of the material and the site formation history, b) the potentiality of non-analogous community structures in the distant past and c) the occurrence of localized or regionally restricted ecological responses to the continental climatic changes during this time. These complicating factors are discussed in more detail below.
A) Accumulation and/or site formation Many factors act upon fossil assemblages during accumulation and fossilization which can modify the identifiability of the remains and, ultimately, the taxonomic composition of the collection. Rodent assemblages are primarily accumulated in caves by predatory birds and small carnivores (Denys 1985;Andrews 1990;Fernández-Jalvo and Avery 2015) and this process introduces unique situations within which the materials reflection of the original rodent community structure may become skewed. Predator dietary preferences and hunting behaviours may lead to a lack of representation of key species in the fossil assemblage, or selective destruction of identifiable bones and teeth during consumption or shortly thereafter. To account for this, a taphonomic analysis was completed on both the Hohle Fels and Geißenklösterle rodent assemblages to quantify the potential biases introduced to the taxonomic record during material accumulation (Rhodes et al. 2018;. The results of this analysis showed the material was accumulated by various predatory birds including the Eurasian eagle owl, great grey owl, snowy owl, kestrel and by small carnivores, most likely the red and/or arctic fox. These predators are all generalists who hunt in a variety of ecozones. The Eurasian eagle owl has the widest hunting range, which can reach 42.5 km 2 (Cantrell 2004). As such, it is fair to assume that the rodent material comes primarily from within the river valley. Some biases in the taxonomic composition of the assemblages may have been introduced, however, including a lower proportion of murids than is present on the landscape, a lower proportion of burrowing and/or nocturnal species such as the ground squirrel and European mole, and a higher proportion of the water vole than would otherwise be expected (based on the Eurasian eagle owls preference).
However, modelling of the effect of time averaging, predator behaviour and prey population dynamics in small mammal assemblages has shown that ~ 140 years of accumulation is all that is needed to provide a stable signal of relative taxonomic abundance that accurately reflects long-term prey community dynamics (Terry 2008). The effect of sample size on the representation of rare taxa is also a potentially confounding issue (Rhodes et al. 2018 as is the history of the sedimentary deposits themselves. Miller (2015) has documented stratigraphic incon-gruities, such as the occurrence of a warm and wet erosional event at the top of GH 17 at Geißenklösterle, and other evidence of a warm karstic environment throughout the MP record of both sites. These events may have modified the accumulating rodent assemblage to some undefinable degree. A recent study of the MP lithic material from this site suggests that the stratigraphic integrity is sound . Furthermore, the temporal correlation between biome and MAT changes and the timing of hemispheric climatic events suggests that any movement of the small rodent material that may have occurred in the past did not dramatically bias the paleoenvironmental signal of the faunal record. The correlations between the faunal and botanical records also suggests post-depositional movement did not dramatically affect either material record. B) Lack of modern analogues Paleoecological transfer functions cannot predict climatic situations which are not present in the modern data used in the regression model (Sachs et al. 1977). Hernandez Fernández and Peláez-Campomanes (2005) recognize this short-coming in their quantitative model and attempt to address it in two ways. In the first cross-validation of the QnBM the localities chosen (5 per climatic zone = 50 total) exhibited the average meteorological values for those biomes. However, the meteorological values were collected from a 10,000 m 2 range in an effort to include variation. In later work the authors include 16 new localities which are situated within ecotones to increase the power of the model in identifying such environments. The Ach valley rodent assemblages do contain non-analogous species compositions, specifically the co-occurrence of steppe, tundra, and rocky species (i.e. Clethrionomys nivalis, Microtus gregalis, and Lemmus lemmus). This is common in European Paleolithic rodent and large faunal records and may reflect the presence of non-analogous environments in the past (Stewart 2005). Claims that the OIS 3 ecology was non-analogue in nature are far from new (Guthrie 2001(Guthrie , 2013Stewart et al. 2003Stewart et al. , 2019Stewart 2005). Although the likelihood of non-analogous environments may increase with age (Stewart 2008), the QlBM and QnBM model has been successfully applied to wide number of archaeological assemblages ranging in age from OIS 7 to the Holocene (Socha, 2014;López-García et al. 2015, 2016, 2017a, 2017b, 2018a, 2018bPinero et al. 2016;Rey-Rodríguez et al. 2016;Discamps and Royer 2017;Berto et al. 2018Berto et al. , 2019Fagoaga et al. 2018). The results of these studies successfully correlate with paleontological, paleobotanical, and geoarchaeological evidence for past environment, for the most part. While the possibility remains that these different material records are producing simul-taneously biased yet consistent climatic and ecological signals, it appears more parsimonious to us to assume that this method, and those applied to the other proxy material records, are correctly deriving signals of the past paleoenvironment. Furthermore, the non-analogous nature of OIS 3 may have been less pronounced in the past than some paleoecological reconstructions may be leading us to believe. The fact that the QlBM, and similar methods, are measuring species modern realized niches (even on a global scale) rather than potential niches may be inflating this issue. While some suggest that these results are caused by 'mixtures of fossils from different times' (Coope and Angus 1975;Coope 2000in Stewart et al. 2019 this seem unlikely given the amount of consistency between records derived from different proxy material types. However, various processes have undoubtedly skewed the paleoclimatic signature of proxy material records to different extents, including time averaging of deposits, the error ranges inherent in our absolute dating methods, and preservation issues, among others. Still, boreal, deciduous, steppe and tundra ecozones were clearly present across southern Germany during this period, and this may reflect an ecotone not present in modern day environmental variability. C) Regional and/or local response to OIS 3 climate This hypothesis has been put forth based on previous qualitative paleoenvironmental reconstructions of the Ach Valley landscape ) and the fact that geographic features of Central Europe may have shaped the vegetative and faunal response to climatic changes (Guthrie 2001). Genetic studies of forest adapted plants (Kolář et al. 2016) and insects (Drees et al. 2016) suggest that the mountainous regions north of the Swiss Alps may have harbored cryptic temperate forest refugia (Stewart and Lister 2001) during interstadial periods. Genetic studies of the recolonization of Central and Western Europe by the common vole (Microtus arvalis) also suggests the persistence of isolated populations within glacial refugia from central Germany to Spain (Tougard et al. 2008). Royer et al. (2016) found evidence for the presence of refugia populations in both isolated areas of favorable climate and pockets of mosaic environments with generally unfavorable but not lethal climates (i.e. within species potential niche ranges). This scenario would fit with our QlBM results if these refugia persisted in the form of transitional ecotones during stadial events. However, if the Ach Valley were a true forest refugia we would expect the presence of forest adapted murids (Apodemus or Eliomys species) and a more consistent presence of temperate voles (Microtus subterraneous and Myodes glareolus) throughout all assemblages.

Conclusions
The QlBM suggests the Ach valley landscape was a combination of arctic tundra and boreal and deciduous forests with periodic expansion of nearby steppe environments. This is consistent with the meteorological estimates produced by the QnBM, as well. Overall, this study suggests the climate was significantly colder than modern regional temperature ranges during the Paleolithic with much higher annual precipitation rates. The Ach Valley Paleolithic landscape had between 2 and 4 months of active vegetation growth a year as well as dramatic temperature shifts between winter and summer periods. Two periods of climatic change are recognizable in the Hohle Fels records characterized by the presence of steppe environments and decreases in arctic tundra landscapes. These periods may correlate to interstadial periods predating the Heinrich 4 event (Bataille and Conard 2018). At Geißenklösterle, two periods of steppe expansion are documented in the MP record, with related decreases in the presence of arctic tundra zones. These periods are tentatively correlated to interstadial events directly before and after Heinrich 5 (Higham et al. 2012). In the Aurignacian, two periods of increased deciduous forest and decreased arctic tundra (without steppe expansion) are also suggestive of warmer continental events straddling what is likely a local response to Heinrich 4.
These results present a paleoecological baseline for future studies examining the role of OIS 3 climate on various aspects of human behaviour during this period, including shifts in subsistence strategies, settlement dynamics, and material culture expression. The presence of a nonanalogue ecotonal environment in the Ach Valley and the evidence of oscillations between cold tundra environments and woodland habitats suggests the local inhabitants must have adapted substantially to landscape variability during both the MP and UP. Finally, the lack of any clear signal of climatic amelioration or declining temperatures through time lends further support to previous claims that climatic volatility and/or degradation did not play a large role in the abandonment of the region by Neanderthal populations.
Funding Open Access funding enabled and organized by Projekt DEAL. This research was partially funded by a DAAD Long-Term Research Grant awarded to S.E. Rhodes.
Data availability All data is provided as supplementary materials attached to this submission.

Conflict of interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.