Climate reconstructions based on GDGT and pollen surface datasets from Mongolia and Baikal area: calibrations and applicability to extremely cold–dry environments over the Late Holocene

Our understanding of climate and vegetation changes throughout the Holocene is hampered by representativeness in sedimentary archives. Potential biases such as production and preservation of the markers are identified by comparing these proxies with modern environments. It is important to conduct multi-proxy studies and robust calibrations on each terrestrial biome. These calibrations use large databases dominated by forest samples. Therefore, including data from steppe and desert–steppe sites becomes necessary to better calibrate arid environments. The Mongolian Plateau, ranging from the Baikal area to the Gobi desert, is especially characterized by low annual precipitation and continental annual air temperature. The characterization of the climate system of this area is crucial for the understanding of Holocene monsoon oscillations. This study focuses on the calibration of proxy–climate relationships for pollen and glycerol dialkyl glycerol tetraethers (GDGTs) by comparing large Eurasian calibrations with a set of 49 new surface samples (moss polster, soil and mud from temporary dry ponds). These calibrations are then cross-validated by an independent dataset of top-core samples and applied to four Late Holocene paleosequences (two brGDGT and two pollen records) surrounding the Mongolian Plateau: in the Altai mountains, the Baikal area and the Qaidam basin, to test the accuracy of local and global calibrations. We show that (1) preserved pollen assemblages are clearly imprinted on the extremities of the ecosystem range but mitigated and unclear on the ecotones; (2) for both proxies, inferred relationships depend on the geographical range covered by the calibration database as well as on the nature of samples; (3) even if local calibrations suffer from reduced amplitude of climatic parameters due to local homogeneity, they better reflect actual climate than the global ones by reducing the limits for saturation impact; (4) a bias in climatic reconstructions is induced by the over-parameterization of the models by the addition of artificial correlation; and (5) paleoclimate values reconstructed here are consistent with Mongolia–China Late Holocene climate trends and validate the application of local calibrations for both pollen and GDGTs (closest fit to actual values and realistic paleoclimate amplitude). We encourage the application of this surface calibration method to reconstruct paleoclimate and especially consolidate our understanding of the Holocene climate and environment variations in arid central Asia.

1 Introduction

Since the understanding of the interactions between the paleoclimate proxies, such as pollen or biomarker abundances, and general circulation model outputs became a major issue in future climate change modeling, resolving the issue of climate proxy calibration is crucial (Braconnot et al., 2012). Current climate changes in extremely cold environments (Masson-Delmotte, 2018), such as Mongolia and Siberia (Fig. 1), are amplified compared with other places around the world (Tian et al., 2014), and the drivers of the current degradation of Mongolian environments in diversity and biomass production still need to be understood. From a climatic point of view, Mongolia is at a junction between the westerlies which are driven by the North Atlantic Oscillation (NAO) and the East Asian Summer Monsoon which is linked to the El Niño–Southern Oscillation (ENSO) and the Inter-tropical Convergence Zone (ITCZ; An et al., 2008). The Mongolian Plateau is a hinge area: the high altitude of the Altai range to the west and the Sayan range to the northwest of the country partially block both the westerlies arriving from the northern Atlantic ocean through the Baikal area and the East Asian Summer Monsoon (EASM, Fig. 2; Chen et al., 2009). The Mongolian system is thus driven by a mix of the distant drag of these two main climatic cells. The understanding of the complex interaction of these cells is necessary and paleoclimate studies are needed to better document this region (Braconnot et al., 2012).

Figure 1 (a) Topographic map of Mongolia (from ASTER data) with the location of surface samples and weather stations considered in the present study; (b) mean annual precipitation; (c) mean annual air temperature; (d) focus on the samples surrounding Taatsiin Tsagaan Lake, Gobi desert; (e) focus on the samples along a valley in the Khentii range; (f) localization of Khangai surface samples; (g) focus on the Baikal Lake transect following the Angara valley. The Mongolian GIS data are issued from the ASTER dataset (https://terra.nasa.gov/about/terra-instruments/aster, last access: January 2018), the meteorological dataset from WorldClim2 and infrastructures from public dataset (ALAGaC) (https://marine.rutgers.edu/~cfree/gis-data/mongolia-gis-data/, last access: January 2018).

Figure 2 Eurasian map of all the pollen surface samples included in the database. The color code refers to the biome pollen inferred for each site. The biomes are WAMX, warm mixed forest; WAST, warm steppe; TEDE, temperate deciduous forest; XERO, xerophytic shrubland; COMX, cool mixed forest; HODE, hot desert; CLMX, cold mixed forest; PION, pioneer forest; TAIG, taiga forest; COST, cold steppe; COCO, cold conifer forest; TUND, tundra; ANTH, anthropic environment; CLDE, cold deciduous forest; CODE, cold desert. The thickest points underline the COST samples selected for this study to operate the transfer function method among all COST sites (shown by a lozenge on map). The arrows indicate the main climatic system driving the Mongolian climate: in orange the Westerlies arriving from the North Atlantic ocean and in blue the Indian Summer Monsoon (ISM) and the East Asian Summer Monsoon (EASM). The dashed line represents the EASM limit following Chen et al. (2010), Q. Li et al. (2018) and Haoran and Weihong (2007) for the northernmost part of the boundary. The Mongolian Plateau (MP) is highlighted by the brown shaded area (Windley and Allen, 1993; Sha et al., 2015). The light blue crossed circles localize the five cores used as test benches for the calibrations (top-core and paleosequences).

Lake sediment archives are commonly used to infer past variations of these climate and environmental systems associated with vegetation and human land use (Lehmkuhl et al., 2011; Felauer et al., 2012; Wang and Feng, 2013). Among the proxies available, pollen and geochemical biomarkers are used as past temperature indicators (ter Braak and Juggins, 1993; Weijers et al., 2007b), and the combination of these proxies helps to polish lake sediment shift interpretations (Atahan et al., 2015; Watson et al., 2018; Martin et al., 2020; Kaufman et al., 2020). For decades the pollen signal has been used to display shifts in vegetation composition and structure (Bennett and Willis, 2002) and has allowed quantitative reconstructions of climate parameters such as precipitation regime and temperatures (Birks et al., 2010; Ohlwein and Wahl, 2012; Wen et al., 2013; Cao et al., 2014; Marsicek et al., 2018). Since vegetation structure and pollen production are mainly influenced by climatic parameters (Zheng et al., 2008) in the absence of human influences, the paleopollen signal is very often interpreted as a response to the climate variations through time (Kröpelin et al., 2008; Wagner et al., 2019). Even if human activities influence pollen rain as well (Hjelle, 1997; Hellman et al., 2009a), these empirical observations of the pollen–climate relation lead to the development of semi-quantitative (Ma et al., 2008) and quantitative calibrations (Brewer et al., 2008; Salonen et al., 2019) of the signal. Different methods have been developed to reconstruct past climates (Chevalier et al., 2020): probability density functions, assemblage approaches, transfer functions (TFs) and methods based on vegetation models (Guiot et al., 2000; Birks et al., 2010; Bartlein et al., 2011; Ohlwein and Wahl, 2012). More precisely, these methods are as follows: the inverse modeling method (IM; Guiot et al., 2000), the weighted averaging partial least squares regression (WAPLS; ter Braak and Juggins, 1993; Ter Braak et al., 1993), the artificial neural networks (ANNs; Peyron et al., 1998), the modern analogue technique (MAT; Overpeck et al., 1985; Guiot, 1990; Jackson and Williams, 2004), the response surface technique (RST; Bartlein et al., 1986), probability density functions (PDFs; Kühl et al., 2002; Chevalier, 2019), modified mutual climate range method (MMCRM; Klotz et al., 2003, 2004), Bayesian hierarchical models (BHMs; Ohlwein and Wahl, 2012), the boosted regression trees (BRTs; Salonen et al., 2014), etc. For northern Europe and despite some problems and pitfalls, Seppä et al. (2004) demonstrated that pollen-inferred climate reconstructions are generally consistent with other independent climatic reconstructions. This study encourages us to lead multi-proxy studies to refine climate reconstruction understanding, and especially in tricky and dry contexts such as the Mongolian Plateau (Rudaya et al., 2009).

Among new promising proxies and from the three last decades, biomarkers such as the glycerol dialkyl glycerol tetraethers (GDGTs) have provided new perspectives on continental temperature reconstructions (Naafs et al., 2017a, b). Among the GDGTs, we will focus on two major groups: the isoprenoid-GDGTs (isoGDGTs; Hopmans et al., 2000) and the branched-GDGTs (brGDGTs; Sinninghe Damsté et al., 2000; Weijers et al., 2007a, b; Dearing Crampton-Flood et al., 2020). BrGDGT assemblages reflect archaeal activity in rivers (De Jonge et al., 2014b), soil (De Jonge et al., 2014a) or lake water column (Dang et al., 2018). The GDGT input origin could be traced using the BIT index (branched and isoprenoid tetraether index, Hopmans et al., 2004; Pearson et al., 2011) and the ${\mathrm{III}}_{\mathrm{a}}/{\mathrm{II}}_{\mathrm{a}}$ ratio (Xiao et al., 2016; Martin et al., 2019a; Cao et al., 2014). BrGDGT environmental drivers are linked to climate parameters (Weijers et al., 2007b), soil typology and vegetation cover (Davtian et al., 2016), which in turn imply land cover and land use. Accurate determinations of the relationships between brGDGT assemblages and climate still need some improvements (Naafs et al., 2018; Wang et al., 2019, 2020) and especially on local to regional scales and in extreme environments. It has been shown empirically (Weijers et al., 2004; Huguet et al., 2013) on cultures of pure strains (Salvador-Castell et al., 2019) as well as on meso- and microcosm experiments (Chen et al., 2018; Martínez-Sosa et al., 2020) that organisms adjust their membrane plasticity by the degree of methylation and cyclization of the compounds. Moreover, some studies have focused on variations in the archaeal community structure (Xie et al., 2015), the archaeal group responses to environmental changes (Knappy et al., 2011) and the GDGT occurrences in different archaeal communities (Liu et al., 2012b) to determine the potential effects of community structure on GDGT relative abundances. To evaluate the provenance and the climatic information brGDGTs bear, several indexes have been proposed in the literature (Table S1 in the Supplement). To monitor these changes, the cyclization ratio of branched tetraethers (CBTs) and methylation of branched tetraether (MBT) indexes linked to environmental factors such as climate and soil parameters have been proposed (Weijers et al., 2007b; Huguet et al., 2013). In particular, with regard to the methylation degree, the ratio of 5-, 6- (De Jonge et al., 2013) and 7-methyl isomers (Ding et al., 2016) responds to environment forcing: the 5-methyl brGDGTs mathematically correlate mainly with temperature (R²=0.76; Naafs et al., 2017a), while 6- (R²=0.69) and 7-methyl brGDGTs (R²=0.44) seem to correlates with moisture and pH (Yang et al., 2015; Ding et al., 2016). More specific indexes have been proposed by De Jonge et al. (2014a) to limit the multi-correlation systems with the withdrawal of 5-methyl compounds such as ${\text{MBT}}_{\text{5Me}}^{\prime }$, which is independent of the pH, and CBT_5Me which is more representative of the soil pH than the former version of the index (index formula in Table S1). The statistical relevance of these indexes is a major issue in brGDGT calibration (Dearing Crampton-Flood et al., 2020). Some regional indexes for soil temperature such as Index₁ (De Jonge et al., 2014a) and Index₂ for Chinese soils (Wang et al., 2016) have been explored too, in the context of a strong local calibration demand (Ding et al., 2015; Yang et al., 2015). IsoGDGTs were first attributed to lake water column production (Schouten et al., 2012), but they were also described in significant but lower proportions in soils (Coffinet et al., 2014). The ratio of isoGDGTs to brGDGTs ($R_{\mathrm{i}/\mathrm{b}}$) has been proposed as a reliable aridity proxy (Yang et al., 2014; Xie et al., 2012). It has been shown that a linear relation exists between these GDGT indexes and some climatic features at large regional scales (in the wide Chinese biome range, from tropical forest to central arid plateau, for instance; Yang et al., 2014; Lei et al., 2016).

Since multi-proxy studies become more and more accurate in both temperature and precipitation reconstruction, local to regional calibrations have been proposed for dry areas such as the arid central Asian (ACA) area: pollen semi-quantitative climate reconstruction (Ma et al., 2008), pollen transfer functions (Herzschuh et al., 2003, 2004; Cao et al., 2014; Zheng et al., 2014) and brGDGT regression models (Sun et al., 2011; Yang et al., 2014; Ding et al., 2015; Wang et al., 2016; Thomas et al., 2017). Even if all of these studies focus on areas surrounding the EASM line (Fig. 2; Chen et al., 2010; Q. Li et al., 2018), the understanding of the climate cells interaction and oscillation over time is still lacunary, and especially on the ACA upper edge. In this context, our study took place in the northernmost part of this climatic system (Haoran and Weihong, 2007). Moreover, we propose the first multi-proxy calibration exercise in the ACA area based on pollen and brGDGT fractional abundances retrieved from modern samples (soil, moss litter, pond mud) in semiarid to temperate conditions. The aim of this study is to take advantage of new, modern surface sample datasets in Baikal area and Mongolia to propose an adapted calibration of pollen and archaeal biomarker proxies for cold and dry environments. For that purpose, local calibrations are compared with global calibrations to infer the influence of calibration scale and proxy types on derived climatic parameters. Our approach is summarized in the following steps:

1. collection of a new set of modern surface samples for Mongolia with homogeneous characterization of their bioclimate environment followed by pollen and GDGT pattern characterization;

2. evaluation of the match between actual bioclimate environments and the associated pollen rain and biomarker assemblages based on mathematical criterion without eco-physiological considerations;

3. creation of local Mongolian Plateau (MP) climate calibrations for pollen and GDGTs and comparison of local and global calibrations in the Mongolian case study;

4. a posteriori validation of the inferred relationships between proxies and ecological likelihood based on the currently developed evidences of brGDGT and pollen rain ecological significance;

5. discussion of the implications of the calibration mismatches in terms of climatic reconstructions in arid and cold environments;

6. testing of the new calibrations (pollen and brGDGTs) through their application on four surrounding Late Holocene records: two pollen records, Dulikha bog (Baikal area; Bezrukova et al., 2005; Binney, 2017) and Lake D3L6 (Altai; Unkelbach et al., 2019), and two brGDGT records, NRX (Altai; Rao et al., 2020) and XRD (Qaidam basin; Sun et al., 2019).

2 Mongolian and Baikal study area

2.1 Coring, sampling area and sample types

The study area lies from 52^∘29^′ to 43^∘34^′ N in latitude and from 101^∘00^′ to 107^∘06^′ E in longitude (Fig. 1a). The sample sites (n=49) are listed in Table S2 with a description of the sample type, the applied analyses, the coordinates and the associated ecosystem. For each site, the Garmin eTreX10 was used to record GPS coordinates to 5 m accuracy. The surface samples were collected throughout Mongolia in 2016 following four transects (n=29): in the Khentii mountain range (MMNT1 and MMNT2; Fig. 1e), in the Orkhon valley (MMNT3), and in the Gobi desert and the Gobi–Altai range (MMNT4; Fig. 1d). During the same field trip, a fifth transect has been done in the Sayan range along the Angara valley, Russia (MRUT1, n=12; Fig. 1g). A Khangai mountains field trip from spring 2009 enlarged this set of data with a sixth transect of surface samples (MMNT5, n=6, in Fig. 1f) and two lake coring samples from Ayrag Nuur (MMNT5C12) and Shargyl Nuur (MMNT5C11), both in Fig. 1f. Both of these top-cores were added to the surface pollen database, while only the MMNT5C12 core has been used as cross-value to check the accuracy of the brGDGT climate models. This core has been clipped into 62 samples of which the top-core has been replicated 6 times (samples MMNT5C12-1 to MMNT5C12-6). Into the MMNT5 transect, mud from two temporary dry ponds has been sampled. These surface muds are referred in following figures as mud. Into the other transects and depending on aridity and vegetation at each site, a soil or a moss polster was sampled. In figures, soil refers to the first 3–5 cm of the ground in dry ecosystems, while moss is a mix between soil, litter and a bryophyte (or Cyperaceae) layer in wetter environments. Moss acts as a pollen trap recording a 3- to 5-year mean pollen signal (Räsänen et al., 2004). In drier areas, the soil surface samples have the same function, in spite of a lower pollen conservation and over-representation of some taxa (Lebreton et al., 2010). In parallel with the GDGT analysis and following the calibration approaches presented in De Jonge et al. (2014a), Davtian et al. (2016) and Naafs et al. (2017a, 2018), mud from temporary ponds and soil samples as well as the soil part of moss litter were also used for actual GDGT analysis. To summarize, this study is based on 49 sites, 48 samples in the pollen dataset, 44 in the brGDGT dataset and 6 cross-validation samples to test the brGDGT models. In terms of sample types, the dataset consists of 30 mosses, 15 soils, 2 pond muds and 2 top-cores.

To test the reliability of our modern calibrations, we have finally selected four paleosequences within or close to the MP used as test benches of the calibrations. For the pollen analysis, the cores of D3L6 from Unkelbach et al. (2019) located in the Mongolian Altai range and the Dulikha bog (Fig. 1; Bezrukova et al., 2005; Binney, 2017) are compared to the Xiangride section (XRD) used for the brGDGT sequence from Sun et al. (2019), sampled in the Chinese Qaidam Basin and the NRX peat bog (Chinese Altai, Fig. 2; Rao et al., 2020). These two cores have recorded the paleoenvironmental changes of the Late Holocene period.

2.2 Vegetation and biomes

The central part of the Mongolian Plateau (MP) is characterized by a dry and cold flat mosaic of steppes and deserts with a 1220 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ median elevation (Fig. 1a; Wesche et al., 2016) and is intersected in its northern part by the Sayan and Khentii ranges, in its southern part by the Gobi–Altai and Qilian Shan ranges aligned along a NW–SE direction and in the west by the Altai range (Windley and Allen, 1993; Sha et al., 2015). A wet and cold highland in the Khangai ranges culminates at 4000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, and a flatter and wetter Mongolian area, the Darkhad basin, is located in the north, close to the Russian border on the edge of the southern Siberian Sayan range. In the northernmost part of the MP, the Baikal lake area is characterized by a basin at a lower altitude (around 600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, Fig. 1g; Demske et al., 2005).

The distribution of vegetation and biomes follows a latitudinal belt organization: in the north the boreal forest presents a mosaic of light taiga dominated by Pinus sylvestris mixed with riparian forest dominated by birches (Betula spp.), alders (Alnus spp.) and willows (Salix spp.; Demske et al., 2005). On the MP, the light taiga dominated by larches (Larix sibirica) and a small amount of birches is mixed with dark taiga composed of Siberian pines (Pinus sibirica) and spruces (Picea obovata; Dulamsuren et al., 2005; Schlütz et al., 2008). The Mongolian taiga is constrained to a region spanning from the Darkhad Basin to the Khentii range (Fig. 1a). On the north face of the Khangai piedmont, the vegetation is dominated by a mosaic of forest–steppe ecosystems: the steppe is dominated by the Artemisia spp. associated with Poaceae, Amaranthaceae, Liliaceae, Fabaceae and Apiaceae (Dulamsuren et al., 2005). On these open lands there are some patches of taiga forest, following roughly the broadside and the northern face of the crest leading on to the grasslands in the valley (Dulamsuren et al., 2005). The two last vegetation layers in Mongolia through the elevation gradient are an alpine meadow dominated by Cyperaceae and Poaceae with a huge floristic biodiversity and an alpine shrubland with pioneer vegetation on the summits (Klinge et al., 2018). On the southern slope of the range, the ecotone between the steppe and the desert vegetation extends hundreds of kilometers from the northern part of the Gobi desert (with water supplied by the Gobi lake area in between) to the Gobi–Altai range in the south (Klinge and Sauer, 2019). In the southernmost part of the country, the warm and dry climate conditions favor desert vegetation dominated by Amaranthaceae, Nitrariaceae and Zygophyllaceae. The vegetation cover is lower than 25 % and is mainly composed of short herbs, succulent plants and a few crawling shrubs.

2.3 Bioclimate systems

In the central steppe–forest biome, the vegetation is marked by an ecotone with short grassland controlled by grazing in the valley and larches on the slopes. The forest is gathered in patches constituting between 10 % and 20 % of the total vegetation cover. There are also some patches of Salix and Betula riparian forests among the sub-alpine meadows on the upper part of the range. This vegetation is characteristic of the northern border of the Palearctic steppe biome (Wesche et al., 2016). This biome is characterized by a range of 800 to 1600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, a mean annual air temperature (MAAT; Fig. 1c) between −2 and 2 ^∘C, and a mean annual precipitation (MAP; Fig. 1b) from 180 to 400 mm yr⁻¹ (Wesche et al., 2016 based on Hijmans et al., 2005). In Mongolia, even if the MAP is very low (${\text{MAP}}_{\text{Mongolia}}\in [\mathrm{50};\mathrm{500}]\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$), the major part of the water available for plants is delivered during late spring and early summer, in contrast to Mediterranean and European steppes (Bone et al., 2015; Wesche et al., 2016). These seasons are the optimal plant growth periods. An unknown amount of precipitation is also brought in winter as snowfall (Rudaya et al., 2020), which is not always measured into the weather station MAP. The main part of the MP MAP occurs during the summer (climate diagrams from Dulamsuren et al., 2005). However, the precipitation origin for Mongolia is still under debate (Piao et al., 2018). Mongolian summer precipitation up to the Baikal area (Shukurov and Mokhov, 2017) seems to be controlled by the East Asian Summer Monsoon (EASM) system instead of the Westerlies' winter precipitation stocked onto the Sayan and Altai range (Fig. 2; An et al., 2008). The alternating Westerlies/EASM domination on the MP climate system appears to fluctuate throughout the Holocene depending on the monsoon strength (Zhang, 2021): the weaker is the monsoon, the further the EASM brings precipitation up to the ACA hyper-continental area. The EASM force may variate in function of the MP snow cover (albedo effect on sun radiance impact; Liu and Yanai, 2002) and/or the Pacific surface temperature (Yang and Lau, 1998). Finally, Piao et al. (2018) insist on the importance of the locally evaporated water recycling within the Mongolian MAP amount.

3 Methods

3.1 Pollen analysis, modern pollen datasets and transfer functions

Different chemical processes were performed on the samples: bryophytic parts of the moss samples were deflocculated by potassium hydroxide (KOH) and filtered by 250 and 10 µm sieves to eliminate the vegetation pieces and the clay particles. Then, acetolysis was performed to destroy biological cells and highlight the pollen grains. For the soil and pond mud samples, 2 steps of HCl and HF acid attacks were added to the previous protocol to remove all the carbonate and silicate components. All the residuals were finally concentrated in glycerol and mounted between slide and lamella. The pollen counts were carried out with a Leica DM1000 LED microscope on a 40× magnification lens. The total pollen count size was determined by the asymptotic behavior of the rarefaction curve. This diagram was plotted during the pollen count using PolSais 2.0, software developed in Python 2.7 for this study. The rarefaction curve was fitted with a logarithmic regression analysis. The counter was suspended whenever the regression curve reached a flatter step (Birks et al., 1992). A threshold for the local derivation at $\mathrm{d}x/\mathrm{d}y=\mathrm{0.05}$ was set. For each sample, the total pollen count is usually around $n\in [\mathrm{350};\mathrm{500}]$ grains for steppe or forest and $n\in [\mathrm{250};\mathrm{300}]$ for drier environments such as desert and desert–steppe.

Among all of the pollen-inferred climate methods, the MAT and the WAPLS were applied in this study on four different modern pollen datasets, and on the D3L6 and Dulikha fossil pollen sequence to test the accuracy of these calibrations (Unkelbach et al., 2019, Figs. 1a and 2). The MAT consists of the selection of a limited number of analogue surface pollen assemblages with their associated climatic values (Guiot, 1990); while the WAPLS uses a weighted average correlation method on a limited number of principal components connecting the surface pollen fractional abundance to the associated climate parameters (ter Braak and Juggins, 1993; Ter Braak et al., 1993). The first dataset, called the New Mongolian–Siberian Database (NMSDB), is composed of pollen surface samples analyzed in this study (N=49; Figs. 2 and 3). The second one is the same subset aggregated to the larger Eurasian Pollen Dataset (EAPDB) compiled by Peyron et al. (2013, 2017). From this dataset of 3191 pollen sample sites, a pollen–plant-functional-type method was applied to determine the biome for each sample according to the actual pollen rain (Fig. 2; Prentice et al., 1996; Peyron et al., 1998). Then, only the cold steppe (COST) dominant samples were extracted from the main dataset and aggregated with the NMSDB to produce the COSTDB (N=430 sites, shown by a lozenge in Fig. 2). Finally, a scale-intermediate dataset of samples located within the Mongolian border merged with the new Mongolian dataset is presented as MDB (N=151 sites). The relation between each taxon and climate parameter was checked and then the MAT and WAPLS methods were applied with the Rioja package from the R environment (Juggins and Juggins, 2019).

Figure 3 Simplified surface pollen diagram, bioclimatically sorted, of the Siberian–Mongolian transect. The pollen taxa are expressed in % TP. The ecosystem units were determined with a CONISS analysis. The left-hand colored dots represent the ecosystem for each sample from light-taiga–riparian forest (deep blue), light/dark taiga–birch sub-taiga, steppe–forest, alpine meadow, steppe, steppe–desert and desert (deep red). The color scale is presented in Fig. 5. The MAP and MAAT are extracted from Fick and Hijmans (2017).

3.2 GIS bioclimatic data

Because Mongolia and Siberia have relatively few weather stations (Fig. 1a), climate parameters were extracted with R from the interpolated climatic database WorldClim2 (Fick and Hijmans, 2017). We used mean annual precipitation (MAP; Fig. 1b) and mean annual air temperature (MAAT; Fig. 1c), as well as temperatures and precipitation for spring, summer and winter (T_spr, P_spr, T_sum, P_sum, T_win and P_win); mean temperature of the coldest month (MTCO); and the mean temperature of the warmest month (MTWA) in this study to characterize the actual climate. Because the Mongolian Plateau is poor in weather stations, the WorldClim2 database suffers from interpolation errors. The surface sites presenting inconsistent climate parameters (MAP < 0 or MAP < season precipitation) were removed from the global database. The elevation data and the topographic map originate from the ASTER imagery (Fig. 1a). The biome type for each site derives from the LandCover database (Olson et al., 2001), classification and field trip observations.

3.3 GDGT analysis and calibrations

For consistency with the sampling process and the modeling methodologies developed for pollen analysis, soil parts of the moss polsters, soil samples and pond mud were treated for GDGT analysis. After freeze drying, about 0.6 g of material was sub-sampled. The total lipid extract (TLE) was microwave extracted (MARS 6 CEM) with dichloromethane (DCM) $/$ MeOH (3:1) and filtered on empty SPE cartridges. The extraction step was processed twice. Following Huguet et al. (2006), C₄₆ GDGTs (GDGTs with two glycerol head groups linked by C₂₀ alkyl chain and two C₁₀ alkyl chains) were added as internal standard for GDGT quantification. Then, apolar and polar fractions were separated on an alumina SPE cartridge using hexane $/$ DCM (1:1) and DCM $/$ MeOH (1:1), respectively. Analyses were performed in hexane $/$ iso-propanol (99.8:0.2) by high-performance liquid chromatography with atmospheric pressure chemical ionization mass spectrometry (HPLC-APCI-MS, Agilent 1200) in the laboratory of LGLTPE-ENS de Lyon, Lyon, following Hopmans et al. (2016) and Davtian et al. (2018).

Each compound was identified and manually integrated according to its $m/z$ and relative retention time following the integration descriptions from Liu et al. (2012a), and De Jonge et al. (2014a) for 5- and 6-methyl brGDGTs and Ding et al. (2016) for 7-methyl brGDGTs (the peak chromatogram integration is displayed in Fig. S1 in the Supplement). Statistical treatments on isoGDGT (Fig. 4a) and brGDGT (Fig. 4b) abundances were treated following two methods presented in Deng et al. (2016), Wang et al. (2016) and Yang et al. (2019): compounds were gathered by chemical structures such as cycles (CBT) or methyl groups (MBT; De Jonge et al., 2014a). brGDGTs were expressed as fractional abundance [x_i] (Fig. 4b; Sinninghe Damsté, 2016), as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }f[x{]}_i={\displaystyle \frac{n_i}{\sum _{j=\mathrm{1}}^{N_{\text{brGDGT}}}{x}_j}}.\end{array}$$

Figure 4 Fractional abundances of (a) isoGDGTs and (b) brGDGTs for moss polsters (green), soil surface samples (orange) mud from temporary dry ponds (blue) and the full sequence of the Lake MMNT5C12 as paleo-brGDGTs comparison (grey). The punctuation marks ^′ and ^′′ refer to 6- and 7-methyl, respectively.

To infer temperatures from brGDGT abundances, two types of model were applied: linear relationships between temperature and MBT–CBT indexes, and multiple regression (mr) models between one climate parameter and a proportion of multiple brGDGT fractional abundances. For the simple linear regression model, a correlation matrix between climate parameters and indexes was calculated using the corrplot Rcran library. For mr models, we developed in the R environment a stepwise selection model (SSM; Yang et al., 2014) to determine the best-fitting model connecting climate parameters with brGDGT fractional abundances. Then we gathered some of the climate–GDGT linear relations established in previous studies (De Jonge et al., 2014a; Naafs et al., 2017a, b, 2018; Sinninghe Damsté, 2016; Yang et al., 2014, 2019) focusing on a single climatic parameter, MAAT (Table S1). These models were clustered into three categories, by sample type (mosses, soils or pond muds), geographical area (regional or worldwide scale) and the statistical model (MBT-CBT based on multiple regression models). According to the type of environment from which the samples originated, there was peat, soil and lake-inferred modeling. All these models were applied to the Baikal area–Mongolian surface samples, compared with the actual MAAT value at each site and applied to the brGDGT XRD section (Fig. 2; Sun et al., 2019) and the NRX bog (Figs. 1a and 2; Rao et al., 2020).

3.4 Statistical analyses

GDGTs and pollen data were analyzed with a principal component analysis (PCA) to determine the axes explaining the variance within the samples. The biotic values (pollen and GDGTs) were also compared to abiotic parameters (climate, elevation, location and soil features) by the way of a redundancy analysis (RDA). The regression models were run with the p value < 0.05 (model relevance), the R² (correlation level between the variables), the root mean square deviation (RMSE, error on parameter reconstruction) and Akaike's information criterion (AIC, effect of over-parameterization on multiple regression models;  Arnold, 2010; Symonds and Moussalli, 2011). A cross-validation test was performed for all the brGDGT calibrations (from this study and from the literature) using an independent set of six lacustrine samples from the lake MMNT5C12 top-core. Statistical analyses were performed with the Rcran project, using the ade4 package (Dray and Dufour, 2007) for multivariate analysis. All the plots were made with the ggplot2 package (Villanueva and Chen, 2019) or the Rioja package (Juggins and Juggins, 2019) for the stratigraphic plot and the pollen clustering using the CONISS analysis method (Grimm, 1987).

4 Results

4.1 Pollen, climate and ecosystems relations

4.1.1 Modern pollen rain and vegetation representation

The pollen rain (Fig. 3) is dominated by six main pollen taxa: Pinus sylvestris, Betula spp., Artemisia spp., Poaceae, Cyperaceae and Amaranthaceae. The pollen diagram, sorted by bioclimate from the wet and relatively warm Baikal area on the upper part to the dry–warm Gobi desert on the bottom, presents a net arboreal pollen (AP) decrease from 85 % to 5 %. A total of 34.26 % of the variance is explained by PC₁ extending from positive values associated with non-arboreal pollen (NAP; Amaranthaceae, Poaceae and Artemisia spp.) to negative values associated with AP (Pinus undet., Betula spp. Picea obovata and Larix sibirica, in Fig. 5c). This trend shows the transition between ecosystems, marked by the seven main CONISS clusters (Fig. 3) and PC₁ and PC₂ variations (Fig. 5c). Below are the over-representative main taxa for each of the Siberian–Mongolian transect ecosystems:

1. light taiga–riparian forest dominated by Pinus sylvestris (>70 %), Pinus sibirica and very low NAP (<5 %);

2. mixed light/dark taiga–birches sub-taiga with an assemblage of Larix sibirica, Picea obovata, Pinus sylvestris and P. sibirica;

3. forest–steppe ecotone same AP assemblages that the light/dark taiga ecosystem with 20 % of Artemisia spp., plus occurrence of Poaceae, Cyperaceae, Thalictrum spp. and Convolvulus spp.;

4. steppe still dominated by Artemisia spp. (30 %) and rising Poaceae (25 %) and Brassicaceae;

5. alpine meadow overpowered by Cyperaceae up to 50 %, Poaceae, Brassicaceae, Amaranthaceae and Convolvulus spp.;

6. steppe–desert ecotone highlighting by the transition between Amaranthaceae–Caryophyllaceae community and Poaceae–Artemisia spp. assemblages;

7. desert dominated by Amaranthaceae (from 25 % to 65 %) and by rare pollen-type Caryophyllaceae, Thalictrum spp., Nitraria spp. and Tribulus spp.

Figure 5 Multivariate statistics for the proxies clustered by ecosystems: (a) principal components analysis (PCA) and (b) redundancy analysis (RDA) for brGDGT fractional abundances; (c) PCA and (d) RDA for pollen fractional abundances. The variance percentage explained is displayed on the axis label; the size of the dataset (n) and the RDA linear regression (R²) are inserted in each plot area.

4.1.2 Pollen–climate interaction

The pollen rain trends follow similar variations to bioclimate parameters in MAP, MAAT and elevation (Fig. 3). The highest AP values are correlated to large MAP (up to 500 mm yr⁻¹) and relatively high MAAT (around 1 ^∘C), in the low-elevation Baikal area. Then the transition between AP and NAP dominance is marked by decreases in both MAAT (−1 ^∘C) and MAP (300 mm yr⁻¹) connected to the high-elevation Khangai range. Finally, the dominance of NAP in the Gobi desert area is linked to very arid conditions (MAP < 150 mm yr⁻¹) and relatively warm MAAT (up to 4 ^∘C). The correlation between the taxa themselves and climate parameters is R²=0.38 (RDA; Fig. 5d). The rise in MAAT is associated with that of Amaranthaceae, Poaceae, Sedum-type and Caryophyllaceae percentages. On the contrary, the decrease in MAAT is associated with a rise in the AP and Cyperaceae, Artemisia spp. and Brassicaceae percentages. MAP, closely related to RDA₁, rises with AP and decreases with NAP (Fig. 5d). Finally, the elevation gradient favors Artemisia spp. and Cyperaceae for NAP and Salix spp. and Larix sibirica for AP (Fig. 5d).

4.1.3 Pollen-inferred climate reconstructions: MAT and WAPLS results

To reconstruct climate parameters from pollen data, MAT and WAPLS methods were applied on the four scales, modern pollen datasets and the 10 climate parameters (Table 1). All these methods can be run with $n\in [\mathrm{1};\mathrm{10}]$ parameters: the number of analogues for MAT and the number of components for WAPLS. The best transfer functions among all of them were selected by the following approach: in a first step, for each climate parameter the methods fitting with the higher R² and the lower RMSE were selected. Then, in case the highest R² and the lowest RMSE were not applied for the same number of analogues or components, we chose the method presenting the lower number of parameters. Despite the small number of parameters relative to the number of observations, the method fits well (Table 1, Arnold, 2010). MAT method gives better R² in bigger DB than in smaller ones. Fitting increases with the diversity and the size of DB, since MAT is looking for the closest value between climate and pollen abundance. By contrast, WAPLS fits better on the local scale and especially with a smaller number of sites. In this case, the pull of data is largest and the variance is largest (ter Braak and Juggins, 1993). WAPLS also leads to better values of RMSE than R², in contrast to MAT. For temperature, pollen fits better with T_spr or MTWA in Mongolia. Temperatures of the warmest months indeed control both vegetation extension and pollen production (Ge et al., 2017; Li et al., 2011) and especially in very cold areas such as Mongolia. For precipitation, the significant season is the one associated with the summer monsoon system in Mongolia (Wesche et al., 2016). Almost all the Mongolian precipitation falls during the spring and the summer (Wang et al., 2010), and the amount of precipitation controls, among other parameters, the openness of the landscape in Mongolia (Klinge and Sauer, 2019). To simplify the confrontation of the diverse models, the MAAT and MAP were isolated from the rest of the climate parameters. Even if these two climate parameters are not the best-fitting pollen methods, they are the easiest to interpret and are comparable with the GDGT regression models commonly based on MAAT and MAP. In any case, these models are mitigated by the spatial autocorrelation affecting any models made on ecological database (Legendre, 1993) and especially the MAT method more than the WAPLS (Telford and Birks, 2005, 2011).

Table 1 Statistical results of the MAT and WAPLS methods applied to four surface pollen datasets and 10 climate parameters^a. The lower number of parameters leading to the best performance is highlighted in bold.

^a Only the five better-fitting regression models for each climate parameter are shown. ^b The climate parameters correspond to mean annual air temperature (MAAT), mean temperature of the warmest (MTWA) and the coldest (MTCO) months, spring temperature (T_spr), mean annual precipitation (MAP), and precipitation for summer (P_sum) and spring (P_spr).
^c Corresponding to the number of parameters used in the model inferring the best R² and RMSE. ^d Number of parameters, R², RMSE of the finally selected model.

4.2 GDGT–climate calibration

4.2.1 GDGT variance in the dataset

In the MMNT5C12 sediments, isoGDGTs are dominated by GDGT-0 and crenarchaeol (74.6 % and 9.8 % in relative abundances, respectively, in Fig. 4a, grey boxplots). These compounds, mainly lake-produced (Schouten et al., 2012), are also present in the moss samples (32.7 % and 31.3 %, green boxplots) and in soils (57.4 % and 26.7 %, orange boxplots). The variations of their fractional abundance in our soil and moss polster dataset are discrete and poorly linked to climate parameters (Fig. 4a). IsoGDGT patterns in lake sediments do not really diverge from soil samples which can lead to postulation that the in situ production of isoGDGTs in shallow and temporary lakes like MMNT5C12 is reduced (Fig. 4a). At least, it may show that the archaeal community both in lakes and in soils is dominated by methanogenic Euryarchaeota more than Thaumarchaeota (Zheng et al., 2015; Y. Li et al., 2018; Besseling et al., 2018). Then, it appears (Fig. 4a) that the isoGDGTs produced in soils are dominated by crenarchaeol in accordance with studies on high alkalinity of the soil (Y. Li et al., 2018) linked to the impact of aridity (Zheng et al., 2015). However, no relationship exists between the crenarchaeol concentration and MAP (R²=0.14, p value > 0.005). The putative regio-isomer response to MAP (Buckles et al., 2016) is not evidenced in NMSDB.

The average [brGDGT]_tot concentrations differ depending on the sample type:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }[{\text{brGDGT}}_{\text{tot}}{]}_{\text{sed}}=\mathrm{25.8}\pm \mathrm{8.1}{\text{ ng\hspace{0.17em}g}}_{\text{sed}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }[{\text{brGDGT}}_{\text{tot}}{]}_{\text{moss}}=\mathrm{23.2}\pm \mathrm{26.8}{\text{ ng\hspace{0.17em}g}}_{\text{moss}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }[{\text{brGDGT}}_{\text{tot}}{]}_{\text{soil}}=\mathrm{0.3}\pm \mathrm{0.14}{\text{ ng\hspace{0.17em}g}}_{\text{soil}}^{-\mathrm{1}},\\ \text{(2)}& {\displaystyle }& {\displaystyle }[{\text{brGDGT}}_{\text{tot}}{]}_{\text{all}}=\mathrm{16.7}\pm \mathrm{23.6}{\text{ ng\hspace{0.17em}g}}_{\text{sample}}^{-\mathrm{1}}.\end{array}$$

brGDGT fractional abundances are consistent with each sample type: the predominant compounds are the I_a, ${\text{II}}_{\mathrm{a}}^{\prime }$, II_a and III_a (Fig. 4b). These compounds mostly explain the total variance (Fig. 5a). Particularly, the PC₁ represents 22.8 % of the total variance and distinguishes two contrasted poles: the 5-methyl group (mostly with PC${}_{\mathrm{1}}>-\mathrm{0.3}$) associated with steppe–forest and forest sites and the 6- and 7-methyl group on the far negative PC₁ values associated with steppe and desert sites. Even if the 7-methyl brGDGTs appear to have weak significance in the brGDGT variance explanation (Fig. 5a), the surface sample 7-methyl average fractional abundance around 4.6 % is following the normal order of magnitude (4.3 % in Cameroon lakes and 6.2 % for Chinese lakes; Ding et al., 2016).

The sediment samples from the lake MMNT5C12, used for past sequence comparison, are more homogeneous than the surface samples, especially when compared with the moss polsters that present a wide variability (Fig. 4b). On this figure it appears that, globally, soil samples are more relevant analogues to sediments than moss polsters (especially the [III${}_{\mathrm{a}}^{\prime }$], [II_a] and [I_a] fractions in Fig. 4). This variability shows an influence of the sample type on brGDGT responses. On the other hand, sample type also supports first-order climate and environment information, since soil and moss polsters originate mainly from steppe to desert environments and forest/alpine meadows, respectively. About the pond mud samples, the BIT and ${\text{III}}_{\mathrm{a}}/{\text{II}}_{\mathrm{a}}$ indexes (Fig. S2) show that a coherent soil origin is leading the brGDGT input instead of a lacustrine one (Pearson et al., 2011; Martin et al., 2019a; Cao et al., 2020).

4.2.2 Climate influence on brGDGT indexes

The brGDGT/climate RDA shows that the brGDGT variance is dominated by the MAP as the first component (Fig. 5b: RDA₁=10.01 %). The negative values show higher precipitation and uncyclized 5-Methyl GDGTs, such as I_a, II_a and III_a, while the lower MAP match with 6- or 7-Methyl GDGTs, such as ${\text{III}}_{\mathrm{a}}^{\prime }$, ${\text{II}}_{\mathrm{a}}^{\prime }$ and ${\text{II}}_{\mathrm{a}}^{\prime \prime }$ in accordance with De Jonge et al. (2014a). The RDA₂ is slightly more connected to MAAT as opposed to elevation, also clustering the methylated and cyclized GDGTs to the higher MAAT. As in the pollen-climate response, the elevation is a second driving factor not to be neglected. The correlation between relative abundance of methylated and cyclized brGDGTs with climate parameters was not strong (Weijers et al., 2004; Huguet et al., 2013). All the MBT, MBT', ${\text{MBT}}_{\text{5Me}}^{\prime }$ and CBT, CBT', CBT_5Me relations with climate parameters were tested and showed a very low correlation with $R^{\mathrm{2}}\in [\mathrm{0.1};\mathrm{0.35}]$ (Fig. S3b). Once the MBT (Fig. S4a) and the ${\text{MBT}}_{\text{5Me}}^{\prime }$ indexes (Fig. 6) are compared with the world database (Yang et al., 2014; Naafs, 2017; Dearing Crampton-Flood et al., 2019) it appears that the NMSDB set is consistent with known values instead of a wide sample dispersion.

Figure 6 ${\text{MBT}}_{\text{5Me}}^{\prime }$–MAAT relation comparison between the NMSDB surface samples and the world peat and soil database (Dearing Crampton-Flood et al., 2019).

4.2.3 Multi-regression models

The stepwise selection model (SSM) for climate–brGDGT modeling was applied only on the 5- and 6-methyl, because 7-methyl brGDGTs show weak significance in the variance explanation (PCA; Fig. 5a). To guarantee the homogeneity of the calibration, the SSM has been applied on the total surface dataset excluding the two pond mud samples (even if their GDGT input seems to be validated by the BIT and III${}_{\mathrm{a}}/$II_a indexes in Fig. S2). The N_SSM different combinations of the 15 brGDGT compounds result in $N_{\text{SSM}}={\mathrm{2}}^{\mathrm{15}}=\mathrm{32}\mathrm{768}$ models possible for each climate parameter. Even the models including some minor compounds ([br]_i<5 %) have been considered since, in the NMSDB, brGDGT fractional abundances are more fairly distributed than in the global database, in which few compounds overlap the majority of the compounds (Fig. S4b). Indeed, the cumulative fractional abundance curve (Fig. S4c) is slightly faster in reaching the asymptote line for the world peat (blue curve) and the world soil (in brown) than in the NMSDB. The world peat database needs only three brGDGTs to explain more than 85 % of the fractional abundance against more than 10 compounds in the NDMSDB soils. Then, the better-fitting equations (with low RMSE and AIC and high R²) were selected for each number of parameters (number of brGDGTs issued in the linear regression) for both MAAT and MAP. Within the 15 models (one model for each parameter addition), the 9 more contrasted ones were selected for discussion (Table S3). The models with the best statistical results were comprised of between 5 to 12 parameters and present a $R^{\mathrm{2}}\in [\mathrm{0.62};\mathrm{0.76}]$, a RMSE around 1.1 ^∘C or 68 mm yr⁻¹, and an ${\text{AIC}}_{\text{MAAT}}\in [\mathrm{149};\mathrm{156}]$ or ${\text{AIC}}_{\text{MAP}}\in [\mathrm{503.8};\mathrm{503.1}]$. The earlier a parameter is used in the mr models, the greater is its influence. For both MAAT_mr and MAP_mr models, III_a, ${\text{III}}_{\mathrm{a}}^{\prime }$, III_b and ${\text{III}}_{\mathrm{b}}^{\prime }$ are the most relevant compounds for the climate reconstruction (Table 2) which is consistent with the PCA and RDA observations displayed (Fig. 5a and b). Both the MAAT_mr models infer a positive contribution of ${\text{III}}_{\mathrm{a}}^{\prime }$ and a negative contribution of III_a, which confirms these models as eco-physiologically consistent with the RDA results. Moreover, except for ${\text{II}}_{\mathrm{b}}^{\prime }$, all the compounds positively correlate with MAAT and negatively with MAP, in accordance with the MAP–MAAT anti-correlation. The ΔT values closest to 0 reveal the best-fitting model on each point (Fig. 7, left panels). Then, the box plot (Fig. 7, right panels) summarizes the best-fitting model at a regional scale.

Table 2 Statistical values and equations of the best brGDGT MAAT_mr and MAP_mr models.

Figure 7 Validation of brGDGT-climate models on the study sites: reconstructed values for literature MAAT (a), NMSDB mr–MAAT (b) and MAP (c). Models are tested on the NMSDB sites (1) and the box-plot statistics (2) are provided. Sites are clustered in four groups: cross-value on the six first samples of the independent core MMNT5C12, Arkhangai; moss polsters from Mongolian steppe–forest; Gobi steppe–desert soil samples; and moss polsters from Baikal area. Values are plotted in anomaly.

5 Discussion

5.1 Issues in modeling Mongolian extreme bioclimate

5.1.1 Appraisal modeling in arid environments

According to Dirghangi et al. (2013) and Menges et al. (2014), the commonly used brGDGT indexes (MBT and CBT) are not relevant for arid areas with $\text{MAP}<\mathrm{500}\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ because the relation between low soil water content and soil brGDGT preservation and conservation interferes in the brGDGT/climate parameters (Dang et al., 2016). The MAAT models based on MBT and MBT' indexes provide colder reconstructions (Fig. 7c2) as shown by De Jonge et al. (2014a), because arid soils favor 6-methyl brGDGTs (by pH raising due to the low weathering effect of the weak precipitation; Dregne, 1983; Haynes and Swift, 1989) and drive the MBT to decrease towards zero. This explains the colder MAAT_Ding and ${\text{MAAT}}_{{\text{MBT}}_{\text{DJ}}^{\prime }}$ reconstructed values compared to the modern ones. Moreover, the main issue in climate modeling in Mongolia is the narrow relationship between MAAT and MAP. Because of the climatic gradient from dry deserts in the southern latitudes to wet taiga forests in the northern ones, MAAT and MAP maps are strongly anti-correlated (Figs. 1b, c and S5). If this correlation is not statistically determined on the range of the global database (R²=0.35, p<0.005), the impact is significant on the range of the Mongolian sites (R²=0.91, p<0.005). This correlation could be a bias resulting from the interpolation method of the WorldClim2 database. In fact, there are very few weather stations (Fig. 1a; Fick and Hijmans, 2017), and their distribution on the large MP is interrupted by mountain ranges. According to Fick and Hijmans (2017) the interpolation model used in the ACA area(which includes our study area) gives a strong correlation (R²=0.99) and a little error ($\text{RMSE}=\mathrm{1.3}{}^{\circ }\mathrm{C}$) for MAAT and R²=0.89 and $\text{RMSE}=\mathrm{23}\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$ for MAP. Whenever the Baikal area–Mongolian calibrations are used for paleoclimatic reconstructions, the RMSE of the climate parameters has to be added to the RMSE model. Moreover, the relevance of the interpolation models suffers from the transition threshold made in Mongolia between the EASM and the Eurasian Westerlies (Fig. 2; An et al., 2008) and reinforced by the topographic break (Fig. 1a). Because the mr–GDGT models have been compiled with the group of Baikal sites which are out of the MAAT–MAP strong auto-correlation range (Fig. S5), the reliability of the independence of the MAAT and MAP models seems to be guaranteed.

The topographic fence in Mongolia also affects the pollen and brGDGT distributions by itself, as seen in both RDA analyses (Fig. 5b and d) where elevation appears to be a major ecophysiological parameter. Elevation affects vegetation and pollen rain not just because of its influence on local MAAT and MAP but also because it drives other ecophysiological parameters such as O₂ concentration, wind intensity, slopes and creeping soils, snow cover, and exposure (Stevens and Fox, 1991; Hilbig, 1995; Klinge et al., 2018). Elevation as one of the main brGDGT drivers could also be explained by the archaeal community responses to pH, moisture and soil compound variations along the altitude gradient (Laldinthar and Dkhar, 2015; Shen et al., 2013; Wang et al., 2015) and the vegetation shifts (Lin et al., 2015; Davtian et al., 2016; Liang et al., 2019). In any case, a better understanding of the archaeal community's response to ecophysiological parameter variations will considerably improve the brGDGT calibration process (Xie et al., 2015; Dang et al., 2016; De Jonge et al., 2019).

5.1.2 Particularity of the southern-Siberian–Mongolian climate system

Both GDGT and pollen calibrations show that the precipitation calibrations are more reliable than temperature ones (Tables 1 and 2, Figs. 3, 8 and 7), reflecting that the southern-Siberian–Mongolian system seems to be mainly controlled by precipitation. This dominance of precipitation could be due to seasonality. Even if the brGDGT production is considered to be mainly linked to annual temperature means (Weijers et al., 2007a, b; Peterse et al., 2012), the high-pressure Mongolian climate system (Zheng et al., 2004; An et al., 2008) favors a strong seasonal contrast: almost all the precipitation and the positive temperature values happen during the summer (Wesche et al., 2016). Consequently, for the NMSDB pollen transfer functions, the seasonal parameters such as MTWA, T_sum and P_sum better describe the pollen variability than MAAT and MAP climate parameters (better R² and RMSE in Table 1). The opposite is found for EAPDB and COSTDB models and calibrations made on large-scale databases. The Mongolian permafrost persists for half of the year in the northern part of the country (Sharkhuu, 2003) and acts on vegetation cover and pollen production (Klinge et al., 2018). Furthermore, the effects of frozen soils on soil bacterial communities and GDGT production are thought to be important (Kusch et al., 2019) since the archaeal community seems to be shifting with abrupt temperature modifications (De Jonge et al., 2019). This seasonality leads to a quasi-equivalence between MAP and P_Sum (if P_win≈0 then MAP≈P_sum) while MAAT is torn apart by the large T_Sum−T_win contrast (because the MAAT is an average value and not a sum as for MAP). The mathematical consequences of the seasonality on these two climate parameters are not the same. Finally, the MAP appears to be the most reliable climate parameter for southern-Siberian–Mongolian climate studies according to the NMSDB sites (with $\text{MAAT}<\mathrm{5}{}^{\circ }\mathrm{C}$). Even if the brGDGTs seem to respond to summer temperature (Wang et al., 2016; Kusch et al., 2019), the summer mr models are not significantly improving the calibration compared to the MAAT_mr ones. For instance, the best $T_{{\text{sum}}_{\text{mr}}}$ is selected by its AIC, and $T_{{\text{sum}}_{\text{mr6}}}$ is inferred using 6 brGDGTs' fractional abundance (R²=0.63 and $\text{RMSE}=\mathrm{1.53}{}^{\circ }\mathrm{C}$). This lack of seasonality effect, expected in such cold areas, is consistent with temperate Chinese sites (brGDGT reconstructions; Lei et al., 2016).

Figure 8 Statistical values plotted against the number of parameters of the different mr–GDGT models: the R², the RMSE normalized on the highest RMSE value and the AIC also normalized. (a) Selection of the two best multi-regression models for each number of parameters; (b) combination of the R² (b) and the AIC (c) values for of all the mr models. The blue dots are for the MAP models, orange dots for the MAAT one.

5.1.3 Extreme bioclimatic condition modeling lead to a better global climate understanding

To reduce the signal $/$ noise ratio, a wider diversity of sample sites should be added as initial inputs in the models. This raises the question of the availability of reliable samples in desert areas. The soil samples in the steppe to desert biomes are often very dry, and these over-oxic soil conditions are the worst for both pollen preservation (Li et al., 2005; Xu et al., 2009) and GDGT production (Dang et al., 2016). brGDGT concentrations in moss polsters and temporary dry pond muds are thus higher than in soils in our database (Eq. 2 and Fig. 4). The explanation of the signal difference between the three types of samples could also originate from the in situ production of brGDGTs inside the moss predominant over the wind-derived particles brought to the moss net. As well as this, it seems that the pool of moss polster is associated with a similar trend to the worldwide peat samples from Naafs et al. (2017b) (Figs. 6 and S4a). Moreover, in the steppe or desert context of poor availability in archive sites, the edge clay samples or top-cores of shallow and temporary lakes could be a solution for paleosequence studies. The two pond mud samples of the NMSDB are included within the soil–moss trend for all models (Figs. 5 and S4 and S3 in the Supplement). Even if the brGDGT production and concentrations are different in soils than in lakes due to lake in situ production (Tierney and Russell, 2009; Buckles et al., 2014), this effect is a function of the lake depth (Colcord et al., 2015) and consequently negligible for shallow lakes, and almost absent for lake edge samples as shown by Birks, H. J. B. (2012) for Lake Masoko in Tanzania.

The soils of the Gobi desert also have a high salinity level which is also a parameter of control on brGDGT fractional abundances (Zang et al., 2018). This taphonomic bias (also climatically induced) could explain part of the histogram variance of Fig. 4 related to the sample type as well as the shift of the soil cluster from the regression line in the cross-value plot of brGDGT MBT'/CBT models in Fig. S3 in the Supplement. Even if the impact of salinity on sporopollenin is not well understood, salt properties may affect pollen conservation in soils (Reddy and Goss, 1971; Gul and Ahmad, 2006).

Finally, the saturation effect of the proxies when they reach the limits of their range of appliance is also to be taken into consideration. Since both pollen and brGDGT signals are analyzed in fractional abundance (i.e., percentage of the total count of concentration), these proxies evolve in a [0;1] space. The saturation effect appears when extreme climatic conditions are reached (Naafs et al., 2017a, b). For instance, in a tropical context, temperature values are too high to be linearly linked to fractional abundances (Pérez-Angel et al., 2019). Considering pollen–climate relationships, the inferior limit of pollen percentage is critical: for the majority of pollen types, whenever MAAT or MAP reaches a very high or low threshold, the pollen fractional abundance approaches zero (Fig. 9). These limit areas need to be closely investigated, which legitimizes the local calibration methods.

Figure 9 Relationships between the eight major pollen taxa (% TP) and MAP (mm yr⁻¹, upper part of the facet plot) and MAAT (^∘C, lower part). The black line is the linear fitting for all samples (EAPDB), the orange for all the samples from steppe biome (COSTDB) and the blue only for the NMSDB samples presented in this article.

5.2 Statistical tools for best model selection

5.2.1 Over-parameterization and best models selection

Among the possible methods, statistical values help to select the most reliable paleoclimate reconstruction. However, the correlation (R²) and errors (RMSE) are not informative enough to discriminate between methods and to point to the most suitable ones for paleoclimate modeling. This is especially true for the multi-parameter methods (such as brGDGT multi-regression models and pollen transfer functions). Indeed, the more input parameters in the method, the more accurate it is (Tables 1 and S3 in the Supplement and Fig. 8a and b). All the regression models improve with parameter additions, and in particular the less well fitting methods improve exponentially (lower limit of the R² area; Fig. 8b). The best R² models for each parameter number (Fig. 8a) correspond to the upper limit of the R² area (Fig. 8b). This figure shows that the R² vs. parameter number trend follows a logistic regression (both for MAAT_mr and MAP_mr models). However, and especially for MAAT_mr regression models, this logistic curve becomes asymptotic early, similar to the RMSE decrease. The over-parametrization of the models has proven to produce artifacts in ecological modeling (Arnold, 2010; Symonds and Moussalli, 2011). The issue is thus to identify the threshold in the parameter numbers selected. We used Akaike's information criterion (AIC) to determine the better model without over-parameterization for brGDGT regression models: the lower the AIC, the better the model (Tables 2 and S3). The trend of AIC vs. the parameter number is, however, more complex (Fig. 8c). For MAAT_mr, the regression model becomes more accurate from one to five parameters rapidly but then slowly decreases. The AIC curve takes an asymmetrical hollow shape around five parameters with a steeper slope on the left side (Fig. 8a). The AIC values for MAAT_mr6 and MAAT_mr7 are almost identical (Fig. 8a). The MAP_{mr 6, 7, 8} have almost equivalent AIC values, while the AIC curve shapes differ for the other MAP_mr models (asymmetrical hollow shape around five with a steeper slope on the left side; Fig. 8a). To summarize, the most universal models are MAAT_mr5 and MAP_mr7 (Table 2), but the closed models are also valuable in some local contexts, and especially in similar dry–cold regions. We need to determine the cross-values of these models to select the appropriate ones for the southern-Siberian–Mongolian context.

5.2.2 Assessment of the calibration feedback

The cross-values of the nine best MAAT_mr regression models (Fig. 7a1 and a2) and the best MAP_mr regression models (Fig. 7b1 and b2) were tested. The MAAT reconstructions provide different responses to the three main bioclimate areas (parcel A1): if they properly estimate temperatures in the Baikal area, they overestimate and underestimate them for the center of the northern Mongolian mountains and the Gobi desert, respectively. For precipitation (parcel B1), all the MAP_mr calculated with local to regional databases also misrepresent the extreme values: precipitation values are too high and too low for the Gobi desert and the Baikal area, respectively. To conclude, the wider the dataset extension, the more alleviated the extreme values.

For both MAAT_mr and MAP_mr models, the 95 % interval shrinks with parameter addition, but the mean values do not necessarily get closer to the measured value of the climate parameter (the dashed line in Fig. 7a2 and b2). Therefore, if the tests on the AIC point toward the MAAT_mr4 and the MAP_mr6 regression models, the back-cross plots suggest the MAAT_mr5, MAAT_mr6, MAP_mr5 and MAP_mr7 regression models (Table S3 in the Supplement, colored in blue, and Table 2) provide the best estimates for climate reconstruction in lacustrine archives (ΔMAP=0 and best-fitting temperature for the mean value of all samples; Fig. 7b2 and b1).

5.2.3 Global vs. local calibration

Whatever proxy is used, when reconstructing temperatures and precipitation from past records in a given location, there is the issue of basing reconstructions on calibrations based on local or global datasets (among others; Tian et al., 2014; Cao et al., 2014; Ghosh et al., 2017; Dearing Crampton-Flood et al., 2019). We tested both approaches on our datasets with a cross-value run on the NMSDB-independent set of MMNT5C12 core samples. The global brGDGT–climate calibration artificially reaches higher R² than local ones due to the larger range of values of the involved climate parameters. Since the world soil database in Naafs et al. (2017b, a) covers a wide temperature range ($\text{MAAT}\in [-\mathrm{5};\mathrm{30}]$), counter to the NMSDB ($\text{MAAT}\in [\mathrm{0};\mathrm{5}]$), then its signal $/$ noise ratio gets lower (Fig. 6). Despite the relatively lower R² of 0.66 scored by the MAAT_mr5 compared with world calibrations (Pearson et al., 2011; De Jonge et al., 2014a; Naafs et al., 2017a, b), the boxplots for all MAAT_mr calculated from the NMSDB are mostly centered on the MAAT average value with the shortest variance spreading for all the sites (Fig. 7c1 and c2). These local calibrations fit best with the MAAT_Ding from Ding et al. (2015), which is also a local calibration made on the Tibet–Qinghai plateau database. The global databases made on worldwide sites (De Jonge et al., 2014b, a; Naafs et al., 2017a, b) provide MAAT_model>MAAT_real and large standard deviation (SD). These global calibrations also attenuate the extreme MAAT values: the very cold MP sites are reconstructed with temperatures up to +5 to +10 ^∘C warmer, while the warm Gobi desert sites are down by up to −3 to −5 ^∘C. On the other hand, local calibrations performed on subtropical to tropical Chinese transects (Yang et al., 2014; Thomas et al., 2017) have smaller SD but largely overestimate MAAT values due to the warmer conditions of the initial database sites. In brief, the lake core sediment samples match the best to the modern MAAT and MAP value with the brGDGTs mr models, which invites us to consider that these local brGDGT calibrations present a robust way to approach past climate.

Similarly, for pollen transfer functions, the geographic range of the surface samples on which the calibration relies is a relevant parameter to take into account for the reliability of the paleoclimate reconstructions. The choice of the maximum value of this geographic range has been discussed previously for vegetation modeling, for example, the Relevant Source Area of Pollen (RSAP; Prentice, 1985; Hellman et al., 2009a, b; Bunting and Hjelle, 2010). For MAT and WAPLS regression models, the same issue holds true. The responses of the eight over-represented taxa to climate parameters are different in the three geographic ranges (NMSDB, COSTDB and EAPDB). The linear tendency allows for checking the main trends between taxa distribution and climate parameters, despite the weak linear regressions (p value > 0.005 and R²<0.4, in Fig. 9). For the majority of these taxa, the trend is the same, independent of the database size (Larix spp. and Cyperaceae percentages increasing with weaker MAAT, or Amaranthaceae and Pinus sylvestris percentages increasing with higher MAAT). However, due to the shift between pollen types and their associated vegetation (i.e., Poaceae-pollen signal similar for a wide diversity of Poaceae communities with very contrasted ecophysiological features), trends are controlled in some peculiar cases by the geographical clipping of the DB. Thus, Poaceae have a positive response to MAP on the global scale but not inside the Mongolian area. The human influence on pollen rain is also dependent on the biogeographical context; thus, Artemisia spp. are not considered to be influenced as much by human activity in the Asian steppe environment (Liu et al., 2006) than in the European one (Brun, 2011).

Concerning transfer functions, WAPLS performs better for the local database than for the COST and EAP databases (Table 1). On these subsets, the WAPLS RMSE and R-square values are even higher than for the MAT approach. The major difficulty resides in the reconstructions of precipitation. Even if the RMSE and R² values are higher for all models of MAP than MAAT, the influence of precipitation on vegetation cover is not well understood. In Mongolia it is clear that the precipitation controls the treeline in mountainous areas (Klinge and Sauer, 2019) and the global openness in the steppe–forest ecotone (Wesche et al., 2016) as well as human land use (Tian et al., 2014), but the risk of autocorrelation between MAAT and MAP signals is important, even if the RMSE and R² values are higher for MAP regression models than for MAAT ones (Telford and Birks, 2009; Cao et al., 2014). Tangibly, for the two proxies, even if the global calibrations can operate on our study area, the local calibrations reach higher accuracy.

5.2.4 Test bench of the local calibrations on four paleosequences

To test the reliability of our local calibrations, the pollen transfer function and the brGDGT mr models have been applied on four paleosequences. Because there is still no available core analyzed for both pollen and brGDGTs either in ACA area or in the MP, the Dulikha bog (pollen, Baikal; Bezrukova et al., 2005; Binney, 2017), the D3L6 lake (pollen, Altai; Unkelbach et al., 2019), the NRX bog (Altai, GDGTs; Rao et al., 2020) and the XRD section (Qaidan, GDGTs; Sun et al., 2019) are used. The actual values of the climatic parameters are first compared to the top-core reconstructed climatic parameters (Fig. 10, dashed lines). The amplitude of the variations through time has then to be assessed with regards to the expected regional ranges (Zheng et al., 2004). Finally, reconstructions on known short-term climate events are tested for the last 5000 years. They are namely the Little Ice Age, Warm Medieval Period, Dark Ages Cold Period, Roman Warm Period, 3.5 ka cooling and 4.2 ka event (respectively LIA, WMP, DACP, RWP, 3.5 and 4.2 ka; Zhang et al., 2008; Chen et al., 2015; Aichner et al., 2015).

Figure 10 ACA climate reconstruction for the 5000 yr cal BP. (a) Climate–pollen inferred for the Dulikha peat bog (Bezrukova et al., 2005; Binney, 2017) and (b) Lake D3L6 (Unkelbach et al., 2019) comparing two transfer function methods (WAPLS and MAT) and the four databases (EAPDB, COSTDB, MDB and NMSDB); (c) climate–brGDGT inferred from the XRD section (Sun et al., 2019) and (d) the NRX peat bog (Rao et al., 2020) comparing local (NMSDB) and global calibrations (De Jonge et al., 2014a; Naafs et al., 2017a, b). The climate periods correspond to Little Ice Age (LIA), Warm Medieval Period (WMP), Dark Ages Cold Period (DACP), Roman Warm Period (RWP), 3.5 and 4.2 ka event according to Zhang et al. (2008), Aichner et al. (2015) and Sun et al. (2019). Dashed lines represent the actual surface climate parameters for each cores (from WorldClim2, Fick and Hijmans, 2017).

For the pollen transfer function (Fig. 10a and b), the inferred reconstruction approaches display similar trends during the 5000 years with larger amplitude for MAT than for WAPLS. MAT is also more sensitive to the initial calibration dataset selected than WAPLS (Fig. 10a and b), which is in agreement with previous multi-method studies (Brewer et al., 2008). Indeed, the COSTDB and EAPDB are over-reactive unlike the local calibration, while the WAPLS display the same amplitude for each calibration with a different offset. In particular, the NMSDB consistently displays the values closest to the actual climate parameter values both for MAT and WAPLS. For the D3L6 precipitation, all models seem to drive away from the actual value. This shift highlights the WorldClim2 interpolation issues in the Altai mountains (few weather stations and not accounting for snow melt). For the Dulikha bog precipitation, only the NMSDB has the same trend for WAPLS and MAT. The other datasets are more sensitive to the method used.

In Fig. 10c and d, with regard to brGDGTs, the local calibrations (NMSDB applied to XRD section and NRX peat) provide the closest surface reconstructed temperature values to the actual ($\mathrm{\Delta }\text{MAAT}<\mathrm{2}{}^{\circ }\mathrm{C}$; Fig. 10c and d in the left panels) compared to the global calibration ($\mathrm{\Delta }\text{MAAT}\in \text{[}\mathrm{6};\mathrm{10}\text{]}$ ^∘C; De Jonge et al., 2014a; Naafs et al., 2017a, b). For the XRD section, this is explicable by the high similarity between the type of sediment from the XRD section (dry taphonomy) and a large proportion of the NMSDB surface samples (especially from the Gobi desert). For NRX peat, it is due to the similar high elevated and arid conditions for both NRX and NMSDB. Moreover, the NMSDB calibrations present a more realistic amplitude: 6 ^∘C over 5000 years (similar to pollen-inferred amplitude) as opposed to 10 to 20 ^∘C amplitude for global calibrations in XRD and +4 ^∘C for NMSDB against around +7 ^∘C for global databases in NRX. With regard to the precipitation, the brGDGTs mr models show a decreasing trend throughout the Holocene, particularly well-marked between 1000 and 2000 BP followed by a bounce on the last 1000 years in both brGDGT records (tendency consistent with the WAPLS pollen MAP). All brGDGT calibrations exhibit consistent shifts during the LIA (cold–wet), the WMP (warm–dry) and the 4.2 ka event (cold–wet). These variations are also exaggerated with global calibrations. To conclude, general trends are consistent for all calibration datasets, except for the drying–warming trend inferred by the calibration from De Jonge et al. (2014a).

Even if the pollen-based and brGDGT-based climate reconstructions were not conducted on the same core, the D3L6 and NRX records (Fig. 10b and c) are close to each other within the Altai range. These two signals allow the pollen/brGDGT difference in reconstruction output to be discussed. Except for a few globally calibrated reconstructions (MAP_WAPLS-EAPDB for pollen and MAAT_mr-DJ for brGDGT), the climate inferred follows a similar trend. The main difference between the proxy is the amplitude of the climate shifts: brGDGT-inferred models seem to provide higher-amplitude responses for the Late Holocene trends than for the centennial oscillations.

Figure 10 shows that realistic reconstructed surface values are consistent with literature Holocene trends and validate the application of local calibrations for both pollen and brGDGTs. Furthermore, abrupt oscillations and overall amplitudes of temperature and precipitation variations are realistic, in accordance with regional appraisal (Wu et al., 2020). These results permit us to improve our understanding of the MP and ACA Late Holocene climate variations. Overall, on the 5000-year period the climate appears to follow a drier–warmer trend. More precisely, except for Dulikha bog, the local calibration shows an anti-correlation between temperature and precipitation short-period oscillations for both pollen and brGDGTs: the LIA seems to be colder and wetter than the warm and dry WMP. The same cold–wet behavior is observed for the 4.2 ka event. This conclusion is important and situates the Dulikha, D3L6, NRX and the XRD sequences in the same trend as the majority of the ACA paleosequences (Chen et al., 2010, 2015; Wu et al., 2020), connecting the Altai range and the Qaidam basin to the EASM vs. Westerlies Holocene oscillations. In contrast, the Dulikha bog in the Baikal area is connected with the Siberian Highs and the northern part of the EASM front (Wang and Feng, 2013; Feng et al., 2017).

6 Conclusions

The paleoenvironmental and paleoclimatic signals may present several uncertainties (differential production, preservation, etc.) which can misguide the interpretation of past variations. This study shows how both a multi-proxy approach and an accurate calibration are important in preventing these biases. We propose a new calibration for mean annual precipitation (MAP) and mean annual air temperature (MAAT) from brGDGTs as well as a new pollen surface database available for transfer functions. The correlations between pollen rain and climate on the one hand and brGDGT soil production and climate on the other are visible but are still mitigated by the complex climate system of arid central Asia and the diversity of soils and ecosystems. Precisely, each of our proxies seems to be more narrowly linked to precipitation (MAP) than temperature (MAAT) counter to the majority of calibrations in the literature. This is validated on both modern and past sequences for pollen and brGDGTs. The nature of the samples considered (soil, moss polster and mud from temporary dry ponds) also greatly affected these correlations. The calibration attempt for the extreme bioclimates of the MP is difficult because of the low range of climate values, despite the climate diversity ranging from cold and slightly wet (north) to the arid and warm (south) conditions. Even if global and regional calibrations could be applied in such a setting, local calibrations provide enhanced accuracy and specificity. The MAAT and MAP values do not remarkably spread in the vectorial space, which makes it harder to distinguish the linear correlation against variance noise. Moreover, this range of values is close to the lower saturation limit of the proxies, which makes the accurate local calibration tricky but necessary. The local calibrations also suffer from the reduced size and small geographic extent of the dataset. The vegetation cover, extending from a high-cover taiga forest to bare-soil desert cover, also buffers the climate signal and the GDGT/pollen response. The correlations between climate parameters and GDGT/pollen proportion are therefore lower than they could be at global scale. Nonetheless, and despite the lower correlation of the local calibration, these local approaches appear to be more accurate to fit the actual climate parameters than the global ones: both for pollen transfer functions and brGDGT multiple regression models. It is especially the case during the short Late Holocene period which is not suffering from abrupt ecosystem changes. These positive model results have to be considered in light of over-parameterization limits. Too many parameters in mr–brGDGT models or in pollen MAT or WAPLS transfer function can add artificially to the linear relation between climate and proxies and lead to misinterpretation of paleoclimate records. Akaike's information criterion combined with RMSE and R² values is a fair way to select the best climate model. These local calibrations applied to Dulikha, D3L6, NRX and XRD paleosequences highlighted the temperature and precipitation variation throughout the Late Holocene. The next step will be to test our calibrations on pollen and GDGT records available from the same core. We encourage wider application of this local multi-proxy calibration for a more accurate constraint of these central Asian climatic systems, a crucial improvement to properly model the fluctuations of the monsoon line since the Holocene optimum.