Identifying Global‐Scale Patterns of Vegetation Change During the Last Deglaciation From Paleoclimate Networks

During the last deglaciation (∼19–11 ka before present), the global mean temperature increased by 3–8 K. The concurrent hydroclimate and land cover changes are not well constrained. Here, we use a pollen database to quantify global‐scale vegetation changes during this transitional period at orbital (∼104 years) and millennial timescales (∼103 years). We focus on the proportion of tree and shrub pollen, the arboreal pollen (AP) fraction. Temporal similarities over long distances are identified by a paleoclimate network approach. At the orbital scale, we find coherent AP variations in the low and mid‐latitudes which we attribute to the global climate forcing. While AP fractions predominantly increased through the deglaciation, we identify regions where AP fractions decreased. For millennial timescales, we do not observe spatially coherent similarity structures. We compare our results with networks computed from three deglacial climate simulations with the CCSM3, HadCM3, and LOVECLIM models. Networks based on simulated precipitation patterns reproduce the characteristics of the AP network. Sensitivity experiments with statistical emulators indicate that, indeed, precipitation variations explain the diagnosed patterns of vegetation change better than temperature and CO2 variations. Our findings support previous interpretations of deglacial forest evolution in the mid‐latitudes being the result of atmospheric circulation changes. The network analysis identifies differences in the vegetation‐climate‐CO2 relationship simulated by CCSM3 and HadCM3. We conclude that network analyses are a promising tool to benchmark transient climate simulations with dynamical vegetation changes. This may result in stronger constraints of future hydroclimate and land cover changes.

The last deglaciation, approximately from 19 to 11 ka before present (BP), is the transition from the Last Glacial Maximum (LGM, ∼21 ka BP) to the current warm period, the Holocene (starting at 11.7 ka BP). It is the most recent period with large natural changes in radiative forcing ( Figure 1a) and thus suitable to benchmark the ability of Earth System Models to simulate the response to such changes. The deglaciation is externally forced by a spatial and seasonal redistribution of solar radiation due to the varying orbital configuration of the Earth (Berger, 1978). This orbital forcing is amplified by responses of the Earth system, in particular, a rise in the carbon dioxide (CO 2 ) concentration from roughly 185 ppm to around 280 ppm (Köhler et al., 2017), the melting of large ice sheets, and changes in the distribution of land, ice, and ocean (Ivanovic et al., 2016).
Resulting environmental changes include an increase in temperature Snyder, 2016), a predominant increase in precipitation amounts McGee, 2020), and an expansion of forests (Binney et al., 2017;Prentice et al., 2000;Tian et al., 2018). The non-linear warming trend of the deglaciation, which started earlier in the Southern than in the Northern Hemisphere , was modulated by abrupt events that occur on millennial timescales. Notable events were the Bølling-Allerød-age warming (14.7-12.9 ka BP) and the Younger Dryas-stadial cooling (12.9-11.7 ka BP) in the North Atlantic region, and the Antarctic cold reversal which proceeded out-of-phase with the North Atlantic events . In this study, we distinguish between the trends, which characterize the glacial-interglacial transition, and the collection of events that deviate from the trends and last a few millennia at most. We call the former orbital-scale variations, because the (slower) glacial-interglacial transition occurs on the order of 10 4 years, and the latter millennial-scale variations, as the individual (faster) events are on the order of 10 3 years.
While temperature increased by 3-8 K in the global mean and almost everywhere on the globe since the LGM (Cleator et al., 2020;MARGO Project Members, 2009;Shakun et al., 2012;Snyder, 2016), reconstructed as well as simulated spatio-temporal patterns of precipitation change are more heterogeneous Cleator et al., 2020;Kageyama et al., 2021;McGee, 2020). Globally averaged, precipitation and temperature are consistently positively correlated in simulations with typically 2% of precipitation increase per Kelvin of temperature rise (Allen & Ingram, 2002;Held & Soden, 2006; G. Li et al., 2013;Rehfeld et al., 2020). For  Berger, 1978), CO 2 concentration (Köhler et al., 2017), and simulated decadal-average anomalies with respect to 6 ka BP of global mean surface air temperature (ΔGMST) and global mean annual precipitation (ΔGMPR) from TraCE, HadCM3, and LOVECLIM. (b) Geographical location and shape-coded median inter-sample range (ISR) of ACER AP records and EPICA Dome C (EDC, EPICA Community Members, 2004) and NGRIP (NGRIP members, 2004) ice cores. The explained variance of the ACER AP signal is color-coded for records that are used in our analysis. The background layer shows the Last Glacial Maximum (21 ka BP) land-sea-ice-mask from the ICE-5G ice sheet history (Peltier, 2004), interpolated to the TraCE horizontal grid. Modern coastlines are displayed in black. Exemplary time series of ACER AP, and TraCE tree and shrub fraction (TSF), annual precipitation (PR), and mean annual surface air temperature (TS) box-aggregated to the time axis of ACER records are displayed as insets (see Section 2.2 for processing of model output). Time series labels correspond to the respective site ID of the ACER database. Displayed records are Toushe Basin (4), Kohuora (15), Lake Masoko (21), Monticchio (37), Potato Lake (39), and ODP1234 (69).
Vegetation is an interactive component of the Earth system, which responds to atmospheric forcing but also serves as boundary condition for the atmosphere, leading to land-atmosphere feedbacks (e.g., Crucifix et al., 2005;Davies-Barnard et al., 2017;Kubatzki & Claussen, 1998;O'Ishi & Abe-Ouchi, 2013). In this study, we focus on the coupling from the atmosphere to the vegetation and how vegetation composition responds to climate changes. Major climatic factors which limit vegetation growth are moisture availability, growing season length and temperature, winter temperature, radiation, and CO 2 concentration (e.g., Crucifix et al., 2005;Harrison et al., 2010;Seddon et al., 2016). Globally aggregated, changes of these factors led to an expansion of forests during the last deglaciation (Prentice et al., 2000(Prentice et al., , 2011. However, the existence of areas with decreasing tree cover has been noted on a local level (e.g., Heusser, Heusser, Mix, et al., 2006;Jiménez-Moreno et al., 2010). Modeling studies show that the importance of the limiting factors vary in space and time, depending on the local climate state (e.g., Claussen et al., 2013;Crucifix et al., 2005;O'Ishi & Abe-Ouchi, 2013;Woillez et al., 2011).
The listed bioclimatic factors have a direct mechanistic effect on plant growth. In turn, they are influenced by other climate variables. For example, moisture availability depends on precipitation and evaporation, with increased precipitation leading to a higher moisture availability and increased evaporation to a lower moisture availability. As evaporation is higher for higher temperatures, increased temperatures can have a positive or negative influence on plant growth depending on the local environment. Statistical analyses of limiting factors are additionally hampered by the correlation between different bioclimatic variables in space and time. The explanatory value of the indirect limiting factors mean annual temperature and precipitation has been shown by observational as well as modeling studies (e.g., Brovkin et al., 1997;O'Ishi & Abe-Ouchi, 2013).
Past vegetation changes can be diagnosed from pollen samples that are extracted from sediment cores or peat bogs (Bradley, 2015). In turn, past climatic conditions can be inferred by determining the limiting factors that govern the reconstructed vegetation variations (e.g., Birks et al., 2010;Chevalier et al., 2020). The interpretation of pollen assemblages is complicated by the complex high-dimensional raw signals as typically tens-to-hundreds of different taxa are contained in a pollen assemblage (Bradley, 2015). Therefore, the raw signal is often reduced to a low-dimensional or univariate signal using ecological principles and statistical algorithms (e.g., Harrison et al., 2010;Legendre & Legendre, 2012;F. Li et al., 2010;Prentice et al., 1996). In this study, we use the arboreal pollen fraction (AP), which is the fraction of tree-and shrub-pollen in an assemblage. AP can be extracted from numerous pollen records across the globe. It provides an estimate of the tree cover at a given time and location, but this estimate is biased by various factors. Pollen samples average signals over space and time but catchment areas and temporal averaging scales vary strongly between records, especially when comparing between marine and terrestrial sediment cores. Pollen productivity varies between taxa, with anemophilous species typically producing orders of magnitude more pollen than zoophilous species (Bradley, 2015). Additionally, dispersal efficiency depends on the dispersal agent of the specific taxa. Differences in the geometric characteristics lead to different pollen transport ranges for different taxa. These factors make a global-scale comparison of the AP fraction challenging as reported metadata is often insufficient to explicitly model record-specific characteristics and because correction factors are missing for most taxa and regions. Comparing relative AP changes is more promising than comparing absolute differences between locations. This is because relative changes only rely on the assumption that the confining factors do not change significantly over time while permitting changes in space. Quantitative approaches based on pollen assemblage characteristics are prone to missing important biological components because they overlook contributions from rare taxa. Therefore, statistical approaches such as this study need to be accompanied by global-scale analyses based on expert elicitation (e.g., Nolan et al., 2018), and ideally, the two complement each other.
In this study, we characterize long-distance similarity structures in AP records covering the last deglaciation. We extract the records from the ACER pollen and charcoal database . They are distributed over all continents except for Antarctica ( Figure 1b). We quantify the temporal similarity of the AP signals with a pairwise correlation measure for irregular time series (Rehfeld et al., 2011). To summarize the spatial structure of the temporal similarities, we use paleoclimate networks, which are a useful tool to identify spatio-temporal patterns in proxy record compilations (Bühler et al., 2021;Konecky et al., 2014;Oster & Kelley, 2016;Rehfeld et al., 2013). We assign proxy record locations as nodes of the network. Network links represent significant correlations of the proxy fluctuations. Network measures allow to quantify the overall similarity, cross-regional similarities, and intra-regional homogeneity of proxy signals (Rehfeld et al., 2013). Here, we focus on long-distance similarities which are more likely to result from coherent changes in the governing climatic factors. We complement the network analysis with a principal component analysis (PCA), which is a more established technique of multivariate signal analysis.
We compare the AP networks with networks computed from simulations of the last deglaciation generated with three Earth System Models of varying complexity ( Figure 1): CCSM3 (Liu et al., 2009), HadCM3 (Armstrong et al., 2019), and LOVECLIM (Menviel et al., 2011). We exploit similarities between AP and surrogate networks to assess how well different climatic factors can explain vegetation variations diagnosed from the AP network. The surrogate networks additionally help to analyze the vegetation-climate-CO 2 relationship in the deglacial simulations. Transient simulations of past climate require new evaluation techniques. Weitzel et al. (2019) called for approaches to compare spatio-temporal patterns beyond the classical point-to-point comparison. In this sense, our work is a step towards process-based evaluations of transient simulations against proxy record compilations. In particular, our results do not rely on a specific realization of internal variability in the simulations.
In the following, we first introduce the pollen database and climate simulations employed in our analysis (Section 2). We describe the statistical methods in Section 3. The results based on pollen, simulated climate and surrogate vegetation networks are given in Section 4. We discuss limitations and coherency of our results with earlier studies in Section 5 and conclude with Section 6.

Pollen Records
The ACER pollen and charcoal database  contains 93 high-resolution pollen records which cover parts of the last Glacial cycle (130-0 ka BP). The data set contains ∼20,300 pollen samples with ∼2,400 distinct taxa. We first extract the records which contain samples in the period 22-6 ka BP. To quantitatively analyze millennial-to-orbital-scale features, a sub-millennial resolution is required . Therefore, we select only those records which have a median inter-sample range of less than 500 yrs within the period 22-6 ka BP. This excludes very low-resolution records. 63 records in the ACER database with ∼6,100 samples fulfill this criterion and are therefore chosen for our analyses (see Table S1 in Supporting Information S1).
The locations of the 63 selected records spread from 67°N to 43°S (Figure 1b). The records are distributed over all continents except for Antarctica with the highest spatial density in Europe, North America, and Eastern Asia. Most records are located in the low and mid-latitudes. Large gaps exist in continental Asia, North Africa, and in areas covered by ice sheets during the LGM such that a study of the extension of boreal forests during the deglaciation is not possible. The records stem from 42 terrestrial and 21 marine sediment cores. The median inter-sample range averaged over this set during the period 22-6 ka BP is 228 years. In preparing the ACER database, Sánchez  globally harmonized plant species, families, and genera compared to the original publications to facilitate comparative analysis of large-scale vegetation characteristics. We harmonize reported taxa counts and percentages by computing taxa percentages for all samples.
The ACER database provides original and harmonized chronologies. The original age models employ various methods and radiometric calibration curves. Therefore, Sánchez  constructed harmonized chronologies for 86 of the 93 records with the software package CLAM (Blaauw, 2010); see the Supporting Information of Sánchez  for details. They harmonized radiometric ages to the IntCal13 calibration curves (Reimer et al., 2013), while adopting event ages to the GICC05 (Wolff et al., 2010) and AICC2012 (Veres et al., 2013) chronologies. Here, we use the harmonized chronologies whenever they are provided and original ages for the remaining samples. An update of the chronologies to IntCal20 (Reimer et al., 2020) is desirable for future use of the ACER database, but is outside the scope of this study.

Climate Simulation Output
We use output from published simulations with three Earth System Models of various complexity for comparison with the pollen records: The TraCE transient deglaciation simulation with CCSM3, a series of concatenated time slice simulations with HadCM3, and a transient deglaciation simulation with LOVECLIM.
The TraCE simulation is a transient simulation of the last 22 ka BP with CCSM3 (He, 2011;Liu et al., 2009). CCSM3 is a global, coupled ocean-atmosphere-sea ice-land surface model (Collins et al., 2006). The atmospheric component CAM3 was run with T31 (∼ 3.75 • ) horizontal grid-spacing (Yeager et al., 2006). Vegetation changes were modeled dynamically with CLM-DGVM (Levis et al., 2004). The simulation started from an equilibrium LGM simulation (Otto-Bliesner et al., 2006). Orbital parameters and greenhouse gas concentrations followed Berger (1978) and Joos and Spahni (2008). Ice sheets and topography were updated every 500-1,000 years according to the ICE-5G ice sheet history (Peltier, 2004). Ocean freshwater influx was adjusted manually to reproduce reconstructions of sea-level rise and ice sheet melting (He, 2011).
Second, we use time slice (snapshot) simulations with the coupled climate model HadCM3 (Armstrong et al., 2019;Davies-Barnard et al., 2017;Singarayer & Valdes, 2010). The atmospheric component HadAM3 has a horizontal grid spacing of 2.5° × 3.75° . Vegetation variations were simulated with the dynamic vegetation model TRIFFID (Cox, 2001;Cox et al., 1998). Time slice simulations were produced every 1 Kyr from 22 ka BP to 6 ka BP with insolation, greenhouse gas concentrations, topography, and ice masks adapted as described in Armstrong et al. (2019). We create a pseudo-transient simulation by concatenating the time slice simulations.
The third simulation stems from the intermediate complexity model LOVECLIM (Menviel et al., 2011). The atmospheric component ECBilt has a T21 horizontal grid spacing (∼ 5.6 • ) and three vertical layers of which only the lowest performs moisture storage and transport Goosse et al., 2007;Opsteegh et al., 1998). The transient simulation was started from an LGM equilibrium run (Menviel et al., 2008;Timm & Timmermann, 2007) and continued in transient mode from 18 to 6.2 ka BP. Insolation changes followed Berger (1978), CO 2 concentrations followed Monnin et al. (2001), and the ice sheet topography was updated every 100 years following Peltier (1994). Freshwater pulses were induced in the North Atlantic and Southern oceans to reproduce millennial-scale events (Menviel et al., 2011).
We extract mean annual surface air temperature and annual precipitation amount from all simulations for our analyses. We use mean annual surface air temperature because seasonal temperatures are highly correlated in deglaciation simulations (see Figure S32 in Supporting Information S1). Therefore, effects specific to either growing season or winter temperature are difficult to isolate. We use precipitation instead of moisture availability to separate the temperature and precipitation effects on moisture availability. Because temperature changes tend to be better constrained for the deglaciation, quantifying the importance of the less understood precipitation changes is beneficial for future model-data comparison efforts. More discussion on the choice of variables is provided in the Supplement (Text S18 in Supporting Information S1). For TraCE and HadCM3, output from the dynamical vegetation models is publicly available. As a summary measure of the simulated vegetation evolution, we compute surrogate tree and shrub fractions (TSF) by combining all simulated plant functional types, which contain trees or shrubs. TSF serves as an additional variable for the comparison of TraCE and HadCM3 against the pollen records.

Methods
In this section, we introduce the methods used for data processing, constructing paleoclimate networks, and quantifying the similarity of networks. An overview of the main steps of the data processing and analysis workflow is given in Figure 2.

Extracting Arboreal Pollen Fractions From Pollen Records
To analyze spatio-temporal similarities in the pollen records, we extract AP fractions from all records by summing up the fractions of tree and shrub taxa contained within each sample. The AP fraction is a univariate and ecologically meaningful signal, which is interpretable on a global scale. It is a continuous variable, that is more straightforward to employ for quantitative comparisons compared to categorical signals such as the identification of dominant biomes. To semi-automatize the computation of AP fractions, we harmonize the taxa on a coarse level (mostly genus), remove very rare taxa, and finally select all arboreal taxa. Over all selected records within 6-22 ka BP, this approach leaves only 0.6% of pollen counts as unclassified (see Supplement for details of the procedure). The resulting time series are in good agreement with original publications and our analysis methods are robust to rare misclassifications, which can occur in our automated approach (not shown). Figure S1 in Supporting Information S1 shows the resulting AP fraction time series, which we call the AP signal in the following.
We compute the explained variance of the time series of AP fraction by performing a redundancy analysis (Legendre & Legendre, 2012) for each record (see Supporting Information S1 for details). The univariate AP signal explains on average 44% of the variance in the pollen records, with particularly high values in Europe and North America ( Figure 1b). Lower values cluster in East Asia and Africa. This shows that the AP signal is a good representative of the deglacial vegetation evolution in most regions. We apply a probit transform to all AP time series, which reduces potential edge effects in the similarity estimates caused by AP fractions near zero and one. Our results are not sensitive to this transformation (see Supporting Information S1). As the AP signal mixes information from multiple timescales, we additionally investigate similarity structures of millennial-scale AP variations. To this end, we apply a Gaussian band-pass filter suitable for irregularly sampled time series (Rehfeld et al., 2011, which preserves frequencies from 2 to 8 Kyr. The filter removes orbital-scale signals and potential non-climatic high-frequency fluctuations. This facilitates a focused analysis of millennial-scale variations. We refer to the resulting time series as "filtered AP signals" in the following.

Extracting Surrogate Time Series From Model Simulations
To compare climate simulation output and pollen records, we create surrogate records for each of the real-world pollen records. Therefore, we first extract the time series of all selected variables (temperature and precipitation for TraCE, HadCM3, and LOVECLIM, TSF for TraCE and HadCM3) at the nearest grid box of the pollen records. For marine sediment cores, we choose the closest terrestrial grid box. This approach guarantees consistency between simulated precipitation, temperature, and TSF. Then, we block-sample the time series to the resolution of the respective pollen record to imitate the limited resolution, irregular sampling, and integrated nature of the AP signals. This means that we cut each time series into non-overlapping slices determined by the midpoints between sample ages, average the signal contained in each slice, and assign this average to the date of the corresponding sample. For the analysis of millennial-scale signals in the simulations, we apply a Gaussian band-pass filter to the surrogate TraCE and LOVECLIM time series similar to the AP signals. We do not include HadCM3 in the millennial-scale analysis. The individual "snapshot" equilibrium simulations cannot show realistic millennial-scale events when concatenated (see Section 2.2). We therefore only assess the orbital trend here.
To study the effects of precipitation, temperature, and CO 2 on TSF more directly, we additionally train a statistical emulator to the simulation output of TraCE and one to HadCM3 output. Similar to Wei et al. (2020), we use generalized additive models (Wood, 2017), which are non-parametric, non-linear statistical models, in which TSF is the response variable and mean annual temperature, annual precipitation, and CO 2 concentration are predictors. The emulators explain the major vegetation changes over the deglaciation as simulated by TraCE and HadCM3. We simulate sets of surrogate TSF time series using transient climate forcing by all three variables (EM ALL), and only by precipitation (EM PR), by temperature (EM TS), and by CO 2 (EM CO2). More information on the emulators are given in Text S5 in Supporting Information S1.

Construction of Paleoclimate Networks
Paleoclimate networks are a tool to analyze spatio-temporal similarities between proxy records (Rehfeld et al., 2013. They consist of nodes, which in our case are the locations of the pollen records, and links between them. Unweighted links are assigned based on classifying the value and sign of the pairwise similarity of time series. We gather the set of links in an adjacency matrix . It is used to compute network measures that characterize the structure of the network . We examine the similarity of the time series described in Sections 3.1 and 3.2 to create networks for the unfiltered as well as the millennial-scale AP signals, and each of the surrogate signals. We estimate the similarity of two time series using Gaussian kernel correlation (Rehfeld et al., 2011). This is a weighted moving-window technique that is directly applicable to irregular time series. Samples are weighted proportional to the time difference from the window mean. The aggregation of all windows results in the final correlation estimate. Following Rehfeld et al. (2011), we choose the kernel width h = Δt/4, with Δt being the larger of the mean inter-sample ranges of the compared time series. Gaussian kernel correlation was shown to be more accurate than similarity estimators based on a priori interpolation and requires fewer samples than non-linear measures such as Mutual Information . The estimator returns NA if the overlap of the two time series is insufficient to accurately estimate the similarity. Here, we choose the minimal overlap to be 10 samples in both time series. This removes 5 of the 63 previously selected records because they do not overlap sufficiently with other records to result in a valid similarity estimate. Gaussian kernel correlation compares centered and standardized time series. Thus, it measures the similarity of relative fluctuations from the mean of the respective record. This means that it assesses the synchronicity between periods of relatively high and low AP concentrations compared to the mean in the respective time series.
To test the significance of the estimated correlations, we compare the values against a statistical null model, which consists of regionally homogeneous and temporally auto-correlated vegetation variations. We choose the null model parameters such that the temporal decorrelation length roughly replicates the time needed by vegetation to adapt after perturbations (e.g., Reick et al., 2013). We select the spatial decorrelation length in line with estimates from multi-millennial climate simulations without strong globally synchronous temperature changes (Reschke et al., 2019) (see Supporting Information S1 for a detailed description of the null model and significance test). In accordance with the aims of this study, our null model enables the identification of long-distance coherent vegetation variations. The signals described in Sections 3.1 and 3.2 can be positively or negatively correlated while the null model only exhibits positive correlations. Therefore, we consider a correlation estimate for a pair of time series as significant if the absolute correlation ranks above the 95% percentile of the null model. We call a correlation estimate significantly positive if the correlation is significant and above zero, significantly negative if the correlation is significant and below zero, and insignificant otherwise. We categorize these correlation estimates in an adjacency matrix to create the paleoclimate network. Adjacency matrix entries are set to −1 if the records are significantly negative (negative links), 0 if they are insignificant, and +1 if they are significantly positive (positive links). We further investigate sub-networks, which consist of a subset of nodes (e.g., all records in South America) and the links between these nodes, as well as links between two sub-networks (e.g., all links between South American and European records).
We visualize paleoclimate networks in two ways: as a complete link-by-link visualization, which shows the full geographical structure of the network, and in an aggregated view focusing on inter-continental similarity patterns. The network visualizations are based on the tidygraph (Pedersen, 2019b) and ggraph (Pedersen, 2019a) packages within the statistical programming language R (R Core Team, 2020). The link-by-link rendering builds upon hierarchical link bundling (e.g., Holten, 2006). For the aggregated view on inter-continental links, we count positive and negative site-to-site links between pairs of continents and represent links as an undirected graph. Definitions of aggregated regions, locations of continent-level nodes, and strength of link bundling can be adapted manually. Technical details are provided in Text S8 in Supporting Information S1 and in the code stack provided in support of our study.

Network Measures
To quantitatively study and compare the network characteristics, we use three summary measures: link density, node degrees, and cross-link fractions (CLFs). These have been tested for paleoclimate data in  and characterize different aspects of spatio-temporal dynamics. The link density (LD) is the ratio of present links (represented by non-zero entries in the adjacency matrix) over the number of links in a network where all nodes are connected with each other (fully connected network) Here, A is the adjacency matrix of the network and N V is the number of nodes. The link density takes values between zero and one, with higher values corresponding to more significant links. The node degree (NDG) of a node V i measures its connectivity. It is given by the sum of all links involving V i Node degrees can take values between 0 and N V − 1. Values are high for nodes, which evolve synchronously with many other time series, and low for those, which are significantly correlated with only a few other records. CLFs quantify the aggregated similarity between two sub-networks, for example, between records in North America, which form one sub-network, and records in Africa, which form a second sub-network. To compute CLFs, we consider the nodes V 1 in the first sub-network and V 2 in the second sub-network and compare the number of links between those sub-networks with the number of links in a fully connected network CLFs take values between 0 and 1. The higher the value, the more similar the time series signals are between two sub-networks.
We calculate all three network measures separately for positive and negative links in the adjacency matrix A, and for all links combined. Thanks to the linearity of the network measures, the total of any of the measures is the sum of the respective positive and negative portions. For interpretation of the network measures, note the integrated nature of link densities and CLFs. The link density is proportional to the summation over all node degrees, whereas a CLF is proportional to the sum of node degrees when only links between two sub-networks are considered.
To quantify the spatial alignment of networks with comparable network measures, we compute Cohen's kappa coefficients κ. κ quantifies the agreement of adjacency matrix entries compared to randomized networks with the same number of positive and negative links. κ takes values between −1 and +1. Positive values signify that more entries than expected between random networks coincide and negative values mean that less values than in the random networks coincide.

Principal Component Analysis
PCA is an established method for analyzing correlated multivariate data (Abdi & Williams, 2010). Its goal is to represent the information contained in a data set in a basis of orthogonal variables, the principal components (PCs). The PCs identify statistically independent modes of variability, which explain the highest fractions of variance of the initial multivariate signal (Vidal et al., 2016). We compute PCAs for the ACER AP time series and for TraCE temperature, precipitation, and TSF surrogate time series. To compute the PCA, we interpolate the records to a common time axis. We standardize all time series which avoids over-representation of records with large variability magnitudes. We explain technical details in Supporting Information S1.
Unlike the network analysis, PCA requires interpolation of records to a common time axis which can lead to statistical artifacts. In addition, all records need to cover the full-time period of interest in the PCA. Due to these disadvantages, our PCAs include less records than the network analyses (35 instead of 58 in the networks analysis). The network analysis makes it simpler to account for the non-uniform spatial coverage of records as similarity of variations is quantified for pairs of records and not for all records simultaneously. Therefore, different spatial scales and domains can be separated more easily. An advantage of the PCA is that the common signal of a spatial pattern is directly given by the principal component which has the form of a time series.

Results
We aim at identifying spatio-temporal similarity structures in the vegetation evolution during the last deglaciation and attributing them to governing climatic factors. In the following, we first present results from the AP network. Then, we compare these results with surrogate temperature, precipitation, and vegetation networks from the three Earth System Models. We further compare AP and model networks with networks based on the emulated TSF response in TraCE and HadCM3 to temperature, precipitation, and CO 2 forcing. Finally, we study pollen and surrogate networks derived from signals which are filtered to extract millennial-scale fluctuations.

Spatio-Temporal Similarity Structures in Deglacial Arboreal Pollen Records
The spatio-temporal similarity structures in the AP records are visualized by the ACER AP network layout (Figures 3a and 3b, see Supporting Information S1 for description of the visualization) and summarized in the corresponding network measures (Figures 4 and 5a). The network features numerous inter-continental links (Figure 3b). Here, an inter-continental link is defined as a correlation between time series from different continents which is significant compared to the null model. The total link density of 0.31 quantifies the high fraction of links that are significantly stronger correlated than in the null model ( Figure 4). The node degrees, which measure the connectivity of a location, and CLFs, which quantify inter-continental similarity structures, further show that records across all combinations of continents are linked to each other (Figure 5a). Total node degrees averaged over continents are between 15.1 and 24.9 for all continents but Asia (8.4). On average, records outside of Asia are linked to 34.8% of all other nodes (note that the maximum node degree is 57, that is, the number of records minus the node itself, see Figure 3a). Total CLFs are between 0.31 and 0.60 for all continent combinations not involving Asia. The lower number of links connected to Asia might be due to the comparably little explained variance of the AP signal in the Asian pollen records (see Figure 1b and Section 5.1).
Notably, the sign of the links, corresponding to significantly positive or negative correlation estimates, varies with 54.9% positive and 45.1% negative links (Figure 4a). This shows that regionally heterogeneous AP variability patterns exist. All CLFs and average node degrees feature a substantial portion of positive and negative links (Figure 5a), which indicates that regions with opposing variability patterns exist on most continents. The existence of such regions is supported by the PCA (Figure 5b). The first principal component explains 48% of the variance in the AP records and consists of a monotonic trend starting around 17.5 ka BP. However, the direction of this dominant mode varies regionally with 54.3% of the records having a positive and 45.7% a negative loading (as described in Section 3.5, only 35 records are included in the PCA). For both signs, records with high loadings exist. A positive loading implies an increasing AP trend and a negative loading implies a decreasing AP trend over the deglaciation.
The network analysis and PCA complement each other. Both show a coherent deglacial AP evolution across the low and mid-latitudes as quantified by the high total link densities and explained variance of the first principal component. We attribute this coherence to the dependence of vegetation on the deglacial climate evolution. However, a majority of records with increasing AP trends are accompanied by regions with negative AP trends on most continents. These anti-correlated AP trends manifest in the substantial fraction of negative inter-continental links and the opposite sign of first principal component loadings. In particular, mid-latitude western edges of continents feature anti-correlated behavior compared to the predominant AP evolution with sharp geographic transitions along the west coasts of South and North America (Figure 5b).

Comparison of Arboreal Pollen and Surrogate Networks
To attribute the AP similarity structures to climatic factors, we compare the ACER AP network with the surrogate networks. Additionally, this comparison provides insights into the degree of realism of the climate-vegetation-CO 2 relationships of the dynamical vegetation models incorporated in CCSM3 and HadCM3. Figures 3c-3e show the inter-continental links of the TraCE TSF, precipitation, and temperature networks. The HadCM3 and LOVECLIM networks are visualized in Supporting Information S1. Link densities of all surrogate networks are given in Figures 4a and 4b. Figure 5 shows CLFs and average node degrees of TraCE TSF, precipitation, and temperature, and we provide analogous figures for the other networks in Supporting Information S1. The continent-averaged node degrees are predominantly higher in networks with higher link densities (Figure 5a for TraCE, Figure S21 for HadCM3, Figure S22 for LOVECLIM in Supporting Information S1). This shows that the additional links are distributed relatively uniform in space. The CLFs, which tend to be highest in temperature networks across most continent combinations (Figures 5a, S9, and S10 in Supporting Information S1), support this interpretation. For several individual CLFs, the ratio of positive to negative links is less balanced in the TraCE precipitation than in the ACER AP network (Figure 5a). All 15 CLFs of the ACER AP network have a negative link portion of at least 1/3, but only 12 do so for TraCE precipitation, 12 for HadCM3 precipitation, and 8 for LOVECLIM precipitation. This suggests spatially more clustered regions of anti-correlated signals in the precipitation networks compared to the AP network. We note that the three precipitation networks are spatially not well-aligned as the κ-coefficients between each of them are below 0.1 (not shown). Due to these inter-model differences, an exact spatial alignment of simulated and observed networks cannot be expected and we have to rely on statistics of large-scale properties such as the network measures.   2). This could contribute to the apparent differences in the TSF variation patterns (see Text S19 in Supporting Information S1, for further discussion).
To further disentangle the temperature and precipitation effects in TraCE, we compute κ-coefficients to compare the similarity of the adjacency matrices of simulated TraCE TSF and emulated TraCE TSF with precipitation or temperature forcing only. As we assume perfect spatial alignment for different variables within one simulation, κ-coefficients provide good evidence whether temperature or precipitation is a better predictor of TSF in TraCE. We find that the agreement of the emulated TSF network with temperature forcing with the simulated TSF network is almost indistinguishable from the agreement random networks with the same number of positive and negative links would have (κ = 0.04). In contrast, TraCE TSF has a significant spatial agreement with emulated TraCE TSF with precipitation forcing (κ = 0.19, see Figure S25 in Supporting Information S1). The κ-coefficient increases to κ = 0.35 when insignificant entries are removed from the adjacency matrices ( Figure S25 in Supporting Information S1). Thus, precipitation explains the TraCE TSF variations better than temperature and CO 2 , as it matches the network characteristics more closely than the emulated TSF with CO 2 forcing and because its spatial agreement with TSF is superior to the agreement with emulated TSF forced by temperature. A substantial part of the variability in the surrogate TraCE TSF time series is not explained by the emulated joint response to temperature, precipitation, and CO 2 (κ = 0.17). This could be due to internal variability that is not modeled by the statistical emulator or confounding variables.
The similar characteristics of the precipitation and ACER AP networks make precipitation a candidate to explain the diagnosed AP variability patterns through a monotonic relationship. In contrast, temperature network measures deviate strongly from the ACER AP network such that the large-scale AP characteristics cannot be explained through a monotonic relationship with temperature. However, a more complex relationship due to the negative effects of temperature on moisture availability cannot be studied by this direct comparison. For this, we turn to the TSF networks. With the high percentage of positive links, the HadCM3 TSF network differs strongly from the ACER AP network. The characteristics of the TraCE TSF network are more similar to the ACER AP network.
The dominant mode of variability in the TraCE TSF PCA is a monotonic trend over the deglaciation as it is in the ACER AP network, although the locations of decreasing TSF and AP trends do not always coincide (Figure 5b). Assuming that the AP patterns are representative for tree cover changes (see Section 5.1 for a discussion of confining factors), TraCE TSF reproduces the diagnosed deglacial vegetation changes more realistically than HadCM3 TSF. Combined with the results of the emulated TSF networks presented above, this indicates that precipitation changes explain the diagnosed patterns of AP changes better than changes in temperature and CO 2 .

Millennial-Scale Similarity Structures in Arboreal Pollen and Surrogate Networks
As we expect that the similarity structures in the networks presented in Sections 4.1 and 4.2 are predominantly the result of orbital-scale variations, we additionally study networks with signals filtered to extract millennial timescales. We note that the link densities of all millennial-scale networks except for LOVECLIM precipitation are much lower than in the respective unfiltered networks (Figure 4). This shows that the spatial coherence of millennial-scale fluctuations during the deglaciation was much weaker than the overall coherence. Meanwhile, the fraction of positive to negative links is on the same order as in the respective unfiltered network in all millennial-scale networks. The link density of the millennial-scale ACER AP network is around the level (LD ≤ 0.05) that would be expected by chance if the statistical null model was true. Thus, we cannot confidently diagnose the existence of coherent vegetation structures from this network.
Comparing the millennial-scale TraCE networks, we find the most spatial coherence for surrogate temperatures, followed by precipitation and TSF (Figure 4). However, this hierarchy is absent in LOVECLIM, where the coherence of signals is similar for temperature and precipitation. The inter-continental coherence, as quantified by the CLFs, are mostly unstructured for ACER AP, TraCE TSF, and TraCE precipitation ( Figure S20 in Supporting Information S1). A clearer pattern can be identified for TraCE temperature. Links are predominantly positive between Northern Hemispheric records but mostly negative between Oceania and the Northern Hemisphere.
The δ 18 O oxygen isotopic records from the Greenland ice core NGRIP (NGRIP members, 2004) and the Antarctic ice core EPICA Dome C (EPICA Community Members, 2004) feature distinct high-resolution signals of millennial-scale variations in both hemispheres during the deglaciation (Figure 6). We test for synchronous millennial-scale fluctuations at different nodes with either of these records to further characterize the network layouts. To this end, we enhance the millennial-scale networks by adding the δ 18 O records from NGRIP and EPICA Dome C as external nodes and their significant correlations with AP records as links. Similarly, we append the nodes and links resulting from simulated temperature and precipitation variations at the ice core locations to the respective surrogate networks. For ACER AP (Figure 6a), we cannot identify homogeneous patterns, which additionally underlines the notion that the links could appear by random chance. For TraCE precipitation (Figure 6b), mostly positive links between Greenland and Europe as well as Alaska are found. The links from Antarctica to Europe, Greenland, and Alaska are negative. Additionally, negative links between Greenland and the west coast of North America are found. For TraCE temperature (Figure 6c), the Northern Hemisphere is mostly positively correlated to Greenland, whereas Oceania and extra-tropical South America are positively correlated to Antarctica. For tropical locations, the structures are more heterogeneous.
In contrast to the strong spatial coherence of AP signals on orbital timescales (Section 4.1), the absence of coherent patterns in the millennial-scale ACER AP network hints at predominant local processes on millennial timescales. However, other archive-specific factors could also play a role as discussed in Section 5.3. In contrast, surrogate temperature and, less pronounced, precipitation networks exhibit spatially homogeneous variations of extra-tropical locations with a distinct hemispheric signature.

Discussion
Our paleoclimate network approach reveals a coherent deglacial AP evolution in the low and mid-latitudes, but with regionally anti-correlated trends. In this section, we first evaluate the impact of record quality on these findings. Since the comparison against surrogate networks suggests that the identified structures are better explained by precipitation than temperature and CO 2 variations, we subsequently compare our findings to previous literature on the deglacial vegetation and climate evolution. In contrast to the analysis of the unfiltered signals, our focused study of millennial-scale AP fluctuations does not feature spatially homogeneous structures. We discuss potential reasons for the absence of such structures in the final part of this section.

Influence of Record Quality and Site Type on the Network Measures
The network layouts and measures depend on several factors: the spatio-temporal similarity of the governing climate and vegetation variations, the representativeness of the AP signal for the vegetation evolution at a given location, site type (marine vs. terrestrial provenance of records), pollination mechanism, pollen productivity and transport, taphonomic processes, the period covered by a record, the temporal resolution of records, the precision of an archive's age-depth model, and the spatial distribution of proxy locations. To understand the implications of our results for the spatio-temporal vegetation evolution and the governing climatic processes, we need to evaluate the importance of these influencing factors. As link density and CLFs derive from summations of node degrees (see Section 3.4), we focus on the relationship between the mentioned influencing factors and node degrees.
The impact of some influencing factors can be quantified by comparing the node degrees in the ACER AP network with the respective factors in the records (Figure 7a). We find that the period covered by a record (R 2 = 0.213) and the explained variance of the AP signal (R 2 = 0.110) have a minor impact on the total node degree. Here, the coefficient of determination R 2 quantifies how well node degrees can be predicted from an influencing factor. The impacts of temporal resolution (median inter-sample range), chronological precision (mean uncertainty of the CLAM age models where available), and spatial layout (mean distance from other records) are negligible (Figure 7a). None of the tested measures has a substantial impact on the fraction of positive compared to negative node degrees ( Figure 7a). Based on these results, we conclude that the network measures are mainly a result of the spatio-temporal coherence of the underlying climate and vegetation variations. Only the comparably low number of links connected to Asia might originate from the lower explained variance of AP (Figure 1b). The records, which cover shorter periods, are not clustered in a particular region ( Figure S26 in Supporting Information S1). Thus, this factor is unlikely to bias CLFs and continent-averaged node degrees substantially.
Our methodology does not explicitly account for differences in the spatial and temporal averaging scales between records, as they occur particularly between marine and terrestrial sites. To study these effect, we analyze two AP sub-networks containing only marine or only terrestrial records. The absolute link densities of marine and terrestrial networks are very similar (marine: LD = 0.32, terrestrial: LD = 0.31), which means that the degree of large-scale similarity between the records in the two sub-networks is almost equal (Figure 7b). This suggests that the imprint of the large-scale climate changes on the AP fraction is stronger on the long timescales, which we consider, than the effects from varying temporal and spatial averaging scales of the records. The percentage of negative links is larger in the marine (53.3%) than in the terrestrial network (35.5%). As the spatial distribution of the records in the two sub-networks is different with most marine records being located along mid-latitude western edges of continents, we compare the differences against randomized networks and find that none of the positive or negative link densities in the sub-networks is outside the 5th or 95th percentile of the randomized networks (see Supporting Information S1 for details). While we cannot preclude that the spatio-temporal averaging scale of a record has a substantial influence on the fraction of negative links in a network, these significance tests indicate that other factors such as the spatial distribution of the records have a stronger effect than the classification into marine and terrestrial records.
Differences in pollen productivity and dispersal efficiency between taxa, particularly due to varying pollination mechanisms, and taphonomic processes could additionally influence the ACER AP network. However, we observe a decrease of link density from the TraCE temperature and precipitation networks to the TraCE TSF network, which features realistic network characteristics and is not influenced by pollen-specific processes. Therefore, the complex non-climatic processes that influence vegetation could be sufficient to explain smaller link densities for AP compared to temperature or precipitation networks. Accounting for varying pollen productivity between taxa and spatio-temporal averaging scales of records is desirable for future research to confirm that the diagnosed AP change patterns are indeed reflecting tree cover variations. This step requires the development of pollen productivity estimates for more taxa and continents. ). The statistics are (columns from left to right) median inter-sample range (ISR), explained variance of the AP signal, mean uncertainty of CLAM age models (where available), length of the covered time period, and mean distance to other records. The coefficient of determination R 2 of each quantity is given in the top right corners. All statistics are computed for the interval 22-6 ka BP. For link densities, black dots and bars indicate median, 5% and 95% percentiles of randomized networks (see Supporting Information S1). The link density of the combined network as in Figure 4 (ALL) is shown for comparison.

Implications for the Deglacial Vegetation and Climate Evolution
The significant inter-continental coherence diagnosed from the ACER AP network (Section 4.1) results mainly from synchronous orbital-scale variations. These are likely driven by the concurrent global climate changes. A monotonic trend is the dominant mode of AP variability throughout the deglaciation. A similar trend was also found for temperature and precipitation across all studied regions by Clark et al. (2012). Similar to the global precipitation PCA by Clark et al. (2012), the direction of the trend in the AP network varies between regions, whereas it is consistently increasing in the global temperature PCA in Clark et al. (2012), see Text S20 in Supporting Information S1 for further discussion. However, temperature increases can lead to decreasing AP through increased evaporation and thus decreased moisture availability. Therefore, this difference alone does not preclude an explanatory value of temperature for the deglacial vegetation evolution (see Section 4.2). It should be noted that, except for the Western United States, the regions of negative precipitation trends in the PCA analysis by Clark et al. (2012) do not coincide with the regions of decreasing AP trends in our analysis. We do not find an anti-phased monsoon signal in the AP records between Northern and Southern Hemisphere monsoon systems, as it was diagnosed by Cheng et al. (2012) from speleothems. This is likely due to the spatial distribution of AP records, which are mostly located outside the core monsoon regions. In addition, a monsoon signal would require precipitation to be a limiting factor in the respective regions and several modeling studies found that the ecophysiological CO 2 effect is the main driver of increased tree cover in the South American monsoon region (Claussen et al., 2013;Crucifix et al., 2005;O'Ishi & Abe-Ouchi, 2013).
The increased global mean precipitation and temperature led to a wetting trend in many terrestrial areas (Cleator et al., 2020;Prentice et al., 2000). As a consequence, tree vegetation predominantly increased from the LGM to the early Holocene (Binney et al., 2017;Prentice et al., 2011). As our analysis indicates that precipitation changes are likely a better predictor of AP variations in the low and mid-latitudes than temperature or CO 2 concentration changes (see Section 4.2), we attribute most negative links in the network to the existence of regions where precipitation decreased through the deglaciation. This led to reduced moisture availability. For records in Chile, at the west coast of North America, and in Southeastern Africa, this interpretation is in concordance with previous findings from individual records (e.g., DeBusk, 1998;Heusser, 1998;Heusser, Heusser, & Pisias, 2006;Vincens et al., 2007). As an explanation of the anti-correlated patterns along the American Pacific coast, pole-ward shifts of westerly wind systems during the deglaciation have been proposed (Thompson et al., 1993;Heusser, 1998;Heusser, Heusser, Mix, et al., 2006). The spatial patterns of anti-correlated AP signals, which we find on all continents, suggest that similar atmospheric and oceanic circulation shifts are an important driver of the distribution of positive to negative links in the ACER AP network. Inter-model and model-data differences of the regions, where precipitation decreased over the deglaciation, can be used to investigate potential model biases in future research. Better harmonized simulation protocols as they were proposed by the PMIP4 last deglaciation working group (Ivanovic et al., 2016) would improve the comparability of simulations to preclude that inter-model differences result from varying experimental designs.
The spatial distribution of the ACER records restricts our findings to the low and mid-latitudes. Therefore, our analysis does not cover the northward shift of boreal forests, for which previous studies found strong impacts of the CO 2 and temperature increase (Harrison & Prentice, 2003;O'Ishi & Abe-Ouchi, 2013). Additionally, only a few ACER records are located in regions with tropical rain forests. Therefore, our results do not contradict the dominant impact of ecophysiological CO 2 effects for the expansion of tropical forests (O'Ishi & Abe-Ouchi, 2013; Prentice et al., 2011). Unlike previous modeling studies, our approach is not decomposing the total simulated vegetation changes into responses to different factors. Instead, we identify which of the three studied predictors can explain the diagnosed network structures best. Therefore, our results do not directly contradict the found strong ecophysiological CO 2 impacts on simulated vegetation changes in different Earth system models (Claussen et al., 2013;O'Ishi & Abe-Ouchi, 2013;Woillez et al., 2011), including HadCM3 (Crucifix et al., 2005). For example, the rise of CO 2 concentrations could drive or amplify increasing AP trends in regions with increasing precipitation through the deglaciation. On the other hand, increasing evaporation from increasing temperature can amplify decreasing AP trends in regions with decreasing precipitation. Our conclusions rely on the assumption that the large-scale AP characteristics are representative of the TSF evolution despite the confining factors discussed in Section 5.1. This assumption needs to be tested in future research.

Potential Reasons for Inconsistent Patterns on Millennial Scales
As shown in Section 4.3, millennial-scale coherence of AP and surrogate signals is weaker than the respective unfiltered signals, which are dominated by orbital-scale variations. This is in agreement with the general notion that climate (and vegetation) fluctuations are spatially more coherent with increasing timescales. However, the imprint of events such as the Bølling-Allerød-age warming and Younger Dryas-stadial cooling has been identified in proxy records covering large parts of the globe (e.g., Clark et al., 2012;Renssen et al., 2018), which suggests that millennial-scale fluctuations occurred synchronously on inter-continental scales, albeit not globally homogeneous . As complex non-climatic processes influence the vegetation cover, it is possible that the used proxy records of vegetation composition do not reflect millennial-scale climate fluctuations as prominently.
The surrogate networks support this explanation given that temperature networks feature the spatially most coherent patterns, followed by precipitation and TSF (Figure 4c). Indeed, the link density in the millennial-scale TraCE TSF network is close to the AP network. Millennial-scale vegetation changes could also be more subtle than on orbital scales and therefore not be recorded by the AP signal. However, the AP signal still explains more than 25% of the variance for most band-pass filtered records ( Figure S13 in Supporting Information S1).
Another potential cause for undetected synchronous millennial-scale fluctuations is the limited temporal resolution and dating precision of records. We do not find systematic influences of the record quality on the node degree in the millennial-scale ACER AP network ( Figure S27 in Supporting Information S1). However, this could be because there are too few records in the database with high temporal resolution and small dating uncertainties, since both records need a sufficient quality to identify synchronous fluctuations with Gaussian kernel correlations. Therefore, alternative similarity measures that require a lower record quality might be more suitable to study millennial-scale structures. More high-resolution sediment cores and less uncertain chronologies could potentially resolve these questions. A re-calibration to IntCal20 could reduce age uncertainties and biases as it is substantially more detailed during the deglaciation than IntCal13 (Reimer et al., 2020), see Text S21 in Supporting Information S1 for more details. However, this is only the case if sufficient radiometric dates are available.
Additional multi-proxy studies could help to understand if the non-existing inter-continental coherence of millennial-scale fluctuations diagnosed in this study is a result of the spatial heterogeneity of millennial-scale climate variations or of local non-climatic processes.

Conclusions
We studied the global-scale vegetation evolution during the last deglaciation as recorded by AP fractions in sediment cores. A paleoclimate network approach quantifies spatio-temporal similarity patterns in the record set. We find a strong coherence in the low and mid-latitudes, which likely results from orbital-scale variations that respond to the global climate evolution. Complementing our results with a PCA, we identify two anti-correlated patterns: On the one hand, the predominant temporal evolution is characterized by increasing AP fractions during the deglaciation. On the other hand, distinct areas where AP fractions decreased exist on each continent. The pronounced geographical contrasts between these regions indicate that they are the result of shifts in the atmospheric and oceanic circulation, as previously proposed for the west coasts of South and North America (Heusser, 1998;Heusser, Heusser, Mix, et al., 2006) and diagnosed in model simulations of the Last Glacial Maximum (Braconnot et al., 2012;Kageyama et al., 2021).
Our analysis does not detect spatially coherent millennial-scale structures in the AP records. This could be due to a lack of synchronicity in tree cover variations on millennial timescales, as suggested by results from surrogate networks constructed from the TraCE simulation. Alternatively, the temporal resolution and dating precision could be insufficient to identify similarity structures with our quantitative comparison method. Future work is required to resolve these issues, for example, by testing alternative similarity measures and constructing networks from multi-proxy compilations.
We compare the AP network with surrogate tree and shrub fraction networks from simulations with two Earth System Models. While the TraCE tree and shrub fraction network reproduces the pollen network well, the characteristics of the HadCM3 network deviate. Sensitivity experiments with statistical emulators indicate that these discrepancies are mostly due to differences in response to deglacial CO 2 changes. Our results suggest that our paleoclimate network approach could be a promising tool for benchmarking vegetation models under strongly varying climate conditions. Network and principal component analyses complement each other due to different strengths and weaknesses.
Surrogate precipitation and temperature networks from simulations with three Earth System Models, and sensitivity experiments with statistical emulators of tree and shrub fraction forced by precipitation, temperature, and CO 2 indicate that precipitation variations explain the diagnosed vegetation change patterns in the low and mid-latitudes better than temperature and CO 2 . We find consistent temperature and precipitation network properties between the three simulations, indicating that our results are representative for a large range of Earth System Models. For future research, we recommend to reassess our statistical emulation results using sensitivity experiments with dynamical vegetation models and by repeating the analysis with more transient Earth System Model simulations of the last deglaciation as they become available. Developing methods which account for varying pollen productivity of different taxa and spatio-temporal averaging scales of records is necessary to confirm that the diagnosed patterns of AP variations are indeed reflecting tree cover variations. Our study indicates that network analyses can identify spatio-temporal hydroclimate and vegetation variations on orbital timescales, which provides opportunities for the evaluation of simulations with transient boundary conditions.

Data Availability Statement
R code for data processing and to reproduce all analyses from this study is archived on zenodo (DOI: https://doi. org/10.5281/zenodo.5742817). All data sets used in this study are available in public repositories. The ACER pollen and charcoal database is available on PANGAEA: Sánchez Goñi, Desprat, Daniau, Bassinot, Polanco-Martínez, Harrison, and Yamamoto (2017)