Global biome patterns of the Middle and Late Pleistocene

Our primary aim was to assess the hypothesis that distinctive features of the patterns of vegetation change during successive Quaternary glacial–interglacial cycles reflect climatic differences arising from forcing differences. We addressed this hypothesis using 207 half‐degree resolution global biome pattern simulations, for time slices between 800 and 2 ka, made using the LPJ‐GUESS dynamic global vegetation model. Simulations were driven using ice‐core atmospheric CO2 concentrations, Earth's obliquity, and outputs from a pre‐industrial and 206 palaeoclimate experiments; four additional simulations were driven using projected future CO2 concentrations. Climate experiments were run using HadCM3. Using a rule‐based approach, above‐ground biomass and leaf area index of LPJ‐GUESS plant functional types were used to infer each grid cell's biome. The hypothesis is supported by the palaeobiome simulations. To enable comparisons with the climatic forcing, multivariate analyses were performed of global vegetation pattern dissimilarities between simulations. Results showed generally similar responses to glacial–interglacial climatic variations during each cycle, although no two interglacials or glacials had identical biome patterns. Atmospheric CO2 concentration was the strongest driver of the dissimilarity patterns. Dissimilarities relative to the time slice with the lowest atmospheric CO2 concentration show the log‐linear relationship to atmospheric CO2 concentration expected of an index of ecocarbon sensitivity. For each simulation, extent and total above‐ground biomass of each biome were calculated globally and for three longitudinal segments corresponding to the major continental regions. Mean and minimum past extents of forest biomes, notably Temperate Summergreen Forest, in the three major continental regions strongly parallel relative tree diversities, hence supporting the hypothesis that past biome extents played an important role in determining present diversity. Albeit that they reflect the climatic consequences only of the faster Earth system components, simulated potential future biome patterns are unlike any during the past 800 ky, and likely will continue to change markedly for millennia if projected CO2 concentrations are realised.

ditional simulations were driven using projected future CO 2 concentrations. Climate experiments were run using HadCM3. Using a rule-based approach, above-ground biomass and leaf area index of LPJ-GUESS plant functional types were used to infer each grid cell's biome. The hypothesis is supported by the palaeobiome simulations. To enable comparisons with the climatic forcing, multivariate analyses were performed of global vegetation pattern dissimilarities between simulations. Results showed generally similar responses to glacial-interglacial climatic variations during each cycle, although no two interglacials or glacials had identical biome patterns. Atmospheric CO 2 concentration was the strongest driver of the dissimilarity patterns. Dissimilarities relative to the time slice with the lowest atmospheric CO 2 concentration show the log-linear relationship to atmospheric CO 2 concentration expected of an index of ecocarbon sensitivity. For each simulation, extent and total above-ground biomass of each biome were calculated globally and for three longitudinal segments corresponding to the major continental regions. Mean and minimum past extents of forest biomes, notably Temperate Summergreen Forest, in the three major continental regions strongly parallel relative tree diversities, hence supporting the hypothesis that past biome extents played an important role in determining present diversity. Albeit that they reflect the climatic consequences only of the faster Earth system components, simulated potential future biome patterns are unlike any during the past 800 ky, and likely will continue to change markedly for millennia if projected CO 2 concentrations are realised.

| INTRODUC TI ON
The Middle and Late Pleistocene, the last 780 ky (Gibbard et al., 2005), is characterised by alternating glacial and interglacial stages, with glacial terminations (rapid transitions from glacial to interglacial conditions) occurring at ca. 100 ky intervals. Hays et al. (1976), using stable oxygen isotopes from marine microfossils, demonstrated that these climatic fluctuations were driven by variations in the Earth's orbit (Milankovitch cycles). Given that the orbital forcing of climate during each glacial-interglacial cycle was unique, so too were the resulting climatic conditions as recorded by various palaeoclimatic proxies ( Figure S1). Palaeovegetation records spanning two or more glacial-interglacial cycles, although sparsely distributed globally, similarly evidence vegetation differences, especially between successive interglacials, at any one location (e.g. Miyoshi et al., 1999;Tzedakis, 1994). In Europe, which has the highest density of such records, these also evidence the unique temporal development of vegetation across the continent during different interglacial stages. For more than three decades, the generally accepted hypothesis has been that the unique vegetation characteristics of each glacial-interglacial cycle resulted from climatic differences, and that these, in turn, resulted from differences in orbital forcing and their effects on, for example, ice-sheet extent and albedo, as well as on atmospheric concentrations of naturally occurring greenhouse gases (carbon dioxide, methane and nitrous oxide). As Watts (1988) wrote, 'each well-studied interglacial period reveals a coherent and protracted development of vegetation in response to climatic change' (p. 160), the diversity of expression of interglacials being 'most probably a response to the interrelation of the three main orbital parameters which may result in distinctive forcing patterns and seasonality for each separate interglacial ' (p. 161).
Despite its paradigm status, there are several reasons why this hypothesis, that the differing spatiotemporal vegetation patterns of different interglacials reflect differing climatic patterns as a consequence of novel orbital forcing, has not to-date been systematically tested. The principal problem is obtaining independent continental palaeoclimatic records for successive glacial-interglacial cycles with which objectively to compare palaeovegetation records. Although pollen records can be used to infer palaeoclimate (e.g. Allen & Huntley, 2009;Newnham et al., 2017), this leads to circularity of argument if the inferred palaeoclimate is subsequently compared with the palaeovegetation evidence. An alternative approach, applied to the period since the last glacial maximum (LGM), is to simulate palaeoclimate using a general circulation model (GCM) and compare the simulated climate with the palaeovegetation evidence.
Comparison may be achieved either by driving a vegetation model using the simulated climate and comparing the simulated vegetation with that recorded by pollen data (see e.g. Bigelow et al., 2003;Kaplan et al., 2003), or by quantitatively reconstructing palaeoclimate from palaeovegetation data and comparing the results with the GCM-simulated climate (e.g. Harrison et al., 2014). Whilst both approaches have been applied successfully to a number of regions, and to the period since the LGM, their utility is limited in regions where, or during times when, palaeovegetation records are sparse. Both are also difficult to apply to intervals beyond ca. 50 ka, the effective limit of 14 C (radiocarbon) dating, without a potential for circularity if palaeovegetation records are dated by 'tuning' to marine oxygen isotope (MOI) records of ice volume and/or to orbital forcing (see e.g. Torres et al., 2013).
Effective application of either data-model comparison approach also requires a large number of palaeovegetation data points.
Palaeovegetation records spanning two or more glacial-interglacial cycles at a single locality, however, are globally extremely sparse in occurrence. Although a larger number of records, each representing a single interglacial stage, is available, independently dating such records is even more difficult than is dating records spanning two or more glacial-interglacial cycles. This is well illustrated by the historical recognition in Europe of only two Middle/Late Pleistocene interglacials, referred to as the Eemian (last) and Holsteinian (penultimate). For many decades, all records, most spanning only a single interglacial or part of an interglacial, were assigned to whichever of the Eemian or Holsteinian their palaeovegetation sequence showed greater similarity (see e.g. Hall, 1980). When records were obtained spanning two or more glacial-interglacial cycles, however, these provided definitive evidence of two additional glacial-interglacial cycles interposed between the Holsteinian interglacial, as generally recognised, and the last interglacial (LIG), conventionally labelled the Eemian and correlated with MOI Sub-stage 5e. Key records in this context include that from the Ioanina basin, north-west Greece (Tzedakis, 1994), extending continuously back to MOI Stage 12, the overlapping records from the Velay, Massif Central, France, each spanning two or three interglacial stages (Reille et al., 2000), and independently dated records from various sites in England (Penkman et al., 2011;Schreve, 2001Schreve, , 2019. The Holsteinian thus corresponds to MOI Stage 11, and a number of interglacial records from Europe formerly assigned to either the Holsteinian or the Eemian are now recognised as representing MOI Stages 9 or 7. Despite these advances, independently assigning an age to a record from Europe representing a single interglacial stage remains difficult (Tzedakis et al., 2001). Elsewhere in the world it is often virtually impossible, unless, for example, tephra layers are present in the sediment sequence containing the interglacial deposits and can be dated independently (e.g. Westgate & Pearce, 2017), or the record can be correlated unambiguously with a long continuous palaeovegetation record (e.g. Kershaw et al., 2007;Litwin et al., 1997;Singh et al., 1985;Torres et al., 2013), which is rarely achievable.
Notwithstanding these difficulties, our initial objective was to assess whether and how global vegetation patterns, and their temporal development, differed between successive Middle and Late Pleistocene glacial-interglacial cycles, and whether any differences related to differences in climatic boundary conditions, and hence climatic patterns, as the paradigm hypothesis proposes. Given that the sparse palaeovegetation data preclude a temporally and spatially comprehensive synthesis of global palaeovegetation patterns, even for the most recent glacial-interglacial cycle, and given also the dating and other difficulties described above, we adopted a modelling approach. The LPJ-GUESS dynamic global vegetation model (DGVM; Smith et al., 2001Smith et al., , 2014 was used to simulate global vegetation characteristics for a series of time slices for which GCM simulations of palaeoclimate were available. LPJ-GUESS outputs were then used to infer global biome patterns. Simulations were made for the 'present' and for 206 past time slices between 800 and 2 ka, thus spanning the nine glacial-interglacial cycles commencing with MOI Stage 20. Ten intervals with interglacial characteristics are represented 17,15e,15a,13,11,9e,7e, 5e and 1), as are nine intervals with global ice volume indicating substantial glaciation, albeit not always as extensive as at the LGM (MOI Stages 2,6,8,10,12,14,16,18 and 20). In addition, substantial periods of intermediate ice volume, including a number of interstadials (e.g. MOI sub-stages 3, 5a, 5c and 9a), are represented. Mapping and analysis of the simulation results enabled assessment of the nature of global biome patterns, and their relationship to climatic boundary conditions, since 800 ka.
Given evidence from the simulations of the changing extent and location of biomes since 800 ka, we also explored the extent to which past variations in extent of biomes on the principal continental regions may account for present diversity patterns of key organism groups in those biomes, as argued by Huntley (1993) for northern hemisphere temperate forests. By performing equivalent simulations of potential global biome patterns for projected future climatic conditions and atmospheric carbon dioxide concentration ([CO 2 ] atm ), we also explored the extent to which future patterns may differ from any since 800 ka.

| Climate experiments
An internally consistent series of 206 palaeoclimate experiments, plus a 'pre-industrial' experiment (using mid-twentieth century ice-sheet and orbital boundary conditions, but pre-industrial [ca. 1800 CE] atmospheric concentrations of CO 2 , CH 4 and N 2 O, see 0 ka experiment, Table S1), as well as four experiments representing potential future conditions in 2050 and 2100 under the RCP4.5 and RCP8.5 pathways (van Vuuren et al., 2011), were performed using the HadCM3 fully coupled atmosphere-ocean general circulation model (AOGCM). Design of the experiments followed Singarayer and Valdes (2010) with respect to the boundary conditions and forcings applied (Table S1), but with two important exceptions. First, rather than using a static pre-industrial land surface for all experiments, the AOGCM was coupled to Triffid (Cox, 2001), a simple DGVM with five plant functional types (PFTs), enabling vegetation structure and associated climatically important land-surface properties to interact dynamically with the atmosphere. Second, ice sheets were specified following de Boer et al. (2013 who modelled the history of the four major global ice sheets (Antarctic, Eurasian, Greenland and North American) over the past 5 My, providing simulated extents and thicknesses for each at 2 ky intervals (de Boer, 2014). Simulations were run for a minimum of 500 years to obtain a climate in quasiequilibrium with the imposed boundary conditions. Simulations did not include millennial scale forcing (e.g. meltwater pulses) and hence do not include Heinrich and Dansgaard-Oeschger events. The precise model configuration is HadCM3B-M2.1aD, as defined by Valdes et al. (2017) who showed that this configuration of the model outperforms many higher-fidelity CMIP5 models. More details of model setup can be found in Davies-Barnard et al. (2017).
For each palaeoclimate and potential future experiment, anomalies relative to the pre-industrial experiment were calculated for monthly mean temperatures (°C), monthly precipitation totals (mm) and monthly mean cloudiness (% potential bright sunshine).
Anomalies were calculated using the climatology of the last 100 years of the experiment in each case and interpolated to the centres of the cells of the 0.5 × 0.5° grid used with LPJ-GUESS. Following Huntley et al. (2013), the potential 1901-1930 annual time series of monthly mean temperature, monthly precipitation and monthly mean cloudiness for 0.5 × 0.5° grid cells on shelf areas exposed by lowered glacial sea levels were estimated using thin-plate spline surfaces fitted to data for the 0.5 × 0.5° land grid cells in the CRU TS 3.0 dataset (Harris et al., 2014). Following the approach of Miller et al. (2008), anomalies derived from the palaeoclimate and potential future experiments were applied to the 1901-1930 annual time series of climatic data for each CRU TS 3.0 land grid cell and each shelf grid cell, monthly mean temperature time series first being detrended to remove the warming trend over those three decades. Atmospheric carbon dioxide concentration for each past time slice was determined according to Antarctic ice-core data (Loulergue et al., 2008;Lüthi et al., 2008;Petit et al., 1999) using the EDC3 chronology (Parrenin et al., 2007). Obliquity was specified following Laskar et al. (2004). Boundary conditions for potential future simulations were specified as in Huntley et al. (2021). Soil textures were specified for all time slices using present-day values (Haxeltine & Prentice, 1996), shelf cells being specified as having a mediumcoarse soil texture.

| LPJ-GUESS simulations
LPJ-GUESS simulations were performed as described by Allen et al. (2020) using the carbon-only version of the model (Smith et al., 2014). In all, 20 PFTs were used, comprising 11 tree, 7 shrub and 2 grass PFTs, attributes of which were defined by Allen et al. (2020). In total, 25 replicate 0.1 ha stands, a number more than sufficient to achieve a relatively stable 'equilibrium' PFT composition (Smith et al., 2001), were simulated for each 0.5 × 0.5° grid cell. The model was spun-up from bare ground for 500 simulated years, a sufficient period for the simulated vegetation to reach approximate equilibrium with the climate and [CO 2 ] atm (Hickler et al., 2012), and then run for a further 90 simulated years. The were available, as well as for the 'present' and four RCP scenarios.
For the 'present' simulation, 1901-1930 climatic data were used without applying anomalies but with detrended monthly temperatures, atmospheric carbon dioxide concentration was 280 ppmv and obliquity specified at its contemporary value. Two variables computed by LPJ-GUESS for each PFT were used to infer biome type for each grid cell at each time slice: cmass-above-ground vascular plant biomass (referred to hereafter as carbon mass); and lai-leaf area index (LAI). Grid cell total carbon mass for each PFT, and overall, was computed using estimated ice-free land area of the grid cell at that time (see Appendix S1).

| Biome inference
The biome occupying each 0.5° grid cell during each time slice was inferred from simulated LAI and carbon mass values for the 20 PFTs using the rule-based procedure described by Allen et al. (2020). In all, 21 biomes were recognised (Table S2), the small minority of grid cells not matching the criteria for any biome being categorised as Unclassified. Global biome distributions were mapped for all time slices using ArcMAp 10.5®, ice-sheet and ocean masks being overlain onto each map, and the maps compiled as a slideshow (see Supporting Information, Slideshow). The following measures were calculated for each time slice, globally and for three longitudinal segments (135° W-30° W-Americas; 30° W-65° E-Europe-Africa; and 65° E-135° W-Asia-Beringia): extent of each biome, both as total area (10 3 km 2 ) and as percentage of ice-free land area; number and percentage of grid cells with ice-free land assigned to each biome; total carbon mass (kg) for each biome; total ice-free land area; number of grid cells with icefree land; and overall total carbon mass. For each biome, global total carbon mass of each PFT and percentage of the biome's total carbon mass contributed by each PFT were also calculated for each time slice.

| Multivariate analyses
To quantify differences between biome patterns for pairs of time slices, a measure of dissimilarity was first calculated for each grid cell. This was computed as the sum of two components, the first reflecting the difference in bioclimatic zone, according to the simulated vegetation, and the second the difference in vegetation structure. Each biome was assigned to one, or in two cases two, of six bioclimatic zones, and to one, or in one case two, structural types (Table S2). Dissimilarity between adjacent bioclimatic zones was defined as 20 units, dissimilarities between non-adjacent zones being appropriate multiples of 20 units (Table S3). For structural types, dissimilarities used a non-linear scaling designed to emphasise the greater magnitude of structural difference between forest and non-forest biomes, and also to assign a small dissimilarity even when two biomes' structural types were the same, thus avoiding zero dissimilarity between biomes of the same structural type in the same bioclimatic zone (Table S4). Because the Desert biome encompasses both polar and warm temperate/sub-tropical deserts, when computing the matrix of dissimilarities between biomes (Table S5), dissimilarities between Desert and other biomes were calculated after first assigning the Desert biome either to the Arctic or to the Warm temperate/Sub-tropical bioclimatic zone, whichever was closer to the bioclimatic zone of the biome with which it was being compared. Semi-desert biome dissimilarities were calculated as the mean of two values taking its structural type as Grassland and as Shrubland, this biome being intermediate between these structural types. Similarly, Temperate Broadleaved Evergreen Forest biome dissimilarities were calculated as the mean of two values taking its bioclimatic zone as Temperate and as Warm temperate/Subtropical, this biome spanning these bioclimatic zones.
Using the dissimilarity matrix between biomes, grid-cell dissimilarities between the biomes occupying each grid cell for each of a pair of time slices were assigned for all grid cells with non-zero ice-free land area for both time slices being compared, excluding grid cells for which the biome was Unclassified. The global weighted mean, weights being the ice-free land areas of the grid cells, of these grid-cell dissimilarity values was then used as the overall dissimilarity between global biome patterns for the two time slices. Dissimilarity values were computed in this way for all possible time-slice pairs and the resulting dissimilarity matrix used as input to the multivariate analyses.
To enable comparison with the dissimilarities in biome patterns, dissimilarities in climatic forcing between time slices were computed as the standardised Euclidean Distance using three forcing factors: northern hemisphere summer insolation, using the value for the month June 21st to July 20th at 65° N, computed following Laskar et al. (2004); global ice volume, as reflected by the LR04 Global Pliocene-Pleistocene Benthic δ 18 O stack (Lisiecki & Raymo, 2005); and atmospheric carbon dioxide concentration, derived from Antarctic ice-core measurements (Bereiter et al., 2015;Lüthi et al., 2008, data Hammer et al., 2001) was used to perform a multiple linear regression with these three variables as independent variables and biome dissimilarity as the dependent variable.
The biome pattern dissimilarity matrix was also subjected to multivariate analysis, both ordination and classification methods being applied. Ordination was performed using both principal coordinates analysis (PCoord), as implemented in PAST, and non-metric multidimensional scaling (NMDS), using the 'metaMDS' function in the R-package 'vegan' Version 2.5-7 (R Core Team, 2021) The dataset was classified using K-means cluster analysis, a non-hierarchical method, as implemented in pASt.
Similarly, there is a good overall match between our inferred global biome pattern at the LGM and biomes inferred from palynological data (e.g. Ciais et al., 2012;Prentice et al., 2011). Furthermore, LPJ-GUESS, in the variant that we use without a simulated nitrogen cycle, is generally somewhat less sensitive to changes in [CO 2 ] atm than many similar models, similarly somewhat under-estimating vegetation sensitivity compared to Free Air CO 2 Enrichment (FACE) experiment results (see Terrer et al., 2019, figure 3a). This low sensitivity to [CO 2 ] atm , despite the standard Farquhar-based photosynthesis module, results from a strong co-limitation of photosynthesis by plant-root-available water if plants are water-stressed . Given this evidence, our maps can be viewed as reasonably accurate representations of likely global biome patterns, even though we acknowledge that the simulated plant-physiological responses to [CO 2 ] atm changes are somewhat uncertain.

| Interglacial global biome patterns
Overall global biome patterns for the 10 time slices falling closest to ice-volume minima comparable to that of the Holocene, and thus representing interglacials ( Figure S3), are broadly similar to the present potential biome pattern ( Figure S2). However, closer inspection reveals that no two are identical, with many regional differences in biome patterns; four examples serve to illustrate the range of these variations ( Figure 1). Notable amongst these are varying relative extents of Boreal Evergreen Needle-leaved Forest and Boreal Parkland in both Eurasia and North America. Boreal Evergreen Needleleaved Forest is most extensive during MOI Stages/Sub-stages 1 (8.9 Mm 2 ), 15e (10.3 Mm 2 ) and 19 (9.5 Mm 2 ), and least extensive during 13 (5.1 Mm 2 ), 15a (5.8 Mm 2 ) and 17 (5.2 Mm 2 ), whereas Boreal Parkland is most extensive during 5e (12.2 Mm 2 ), 11 (11.7 Mm 2 ), 15a (13.3 Mm 2 ) and 17 (11.3 Mm 2 ) and least extensive during 1 (7.1 Mm 2 ), 15e (4.6 Mm 2 ) and 19 (5.4 Mm 2 ). Tropical Evergreen Forest extent in the Europe-Africa longitudinal segment ranges between 3.4 Mm 2 during MOI Stage 1 and 2.5 Mm 2 during MOI Sub-stage 15a, the latter value being only 69% of the former. In the Asia-Beringia segment, that includes Australasia, the extent of Desert varies very considerably between interglacials, reaching 5.3 Mm 2 during MOI Sub-stage 15e but only 1. 15a.
An exhaustive listing of further regional differences is not appropriate here, but their overall magnitudes are summarised by biome pattern dissimilarities between the 10 interglacials and the present (Table 1). These identify MOI Stage 17 and Sub-stage 15a as those having global biome patterns most different from that of the present and MOI Sub-stage 9e as that with a global biome pattern most similar to the present.
Examination of the complete series of global biome maps (Supporting Information, Slideshow) reveals further regional differences between interglacials with respect to simulated temporal patterns of biome changes. Such differences in temporal patterns between interglacials are consistent with evidence from those few locations providing records of two or more interglacial palaeovegetation sequences (e.g. Allen et al., 2002;Allen & Huntley, 2009;Miyoshi et al., 1999;Tzedakis, 1994).

| Glacial global biome patterns
Once again, overall global biome patterns are broadly similar for the time slices closest to the ice volume maxima of the last nine glacial stages ( Figure S4). The 668 ka time slice, although falling well before the MOI Stage 16 ice volume maximum, shows a similar biome pattern and is included in Figure 2 and Figure  Steppe are, respectively, 6.8 and 6.6 Mm 2 , 4.7 and 4.6 Mm 2 and 1.8 and 1.8 Mm 2 , whereas for glacial time slices their mean and median extents are, respectively, 8.7 and 9·0 Mm 2 , 5.5 and 5.5 Mm 2 and 4.0 and 4.2 Mm 2 . That Steppe shows the greatest contrast in simulated extent is also consistent with the palaeovegetation record, on the basis of which large parts of the temperate latitudes are inferred to have been occupied by Steppe during glacial times (see e.g. Prentice et al., 2011).
Furthermore, the simulated LGM global biome pattern (MOI Stage 2, 18 ka, Figure 2d) has good visual matches to the maps of 'megabiomes' inferred from palynological data by Ciais et al. (2012 Figure S2) and of biomes inferred from the same data by Prentice et al. (2011).
Dissimilarities in biome patterns between glacial stages and the present, and between earlier glacial stages and MOI Stage 2 (Table 2), identify MOI Stage 14 as the least dissimilar from the present and MOI Stage 2 as the most dissimilar, apart from 668 ka BP which has the greatest dissimilarity from the present of all past time slices simulated, also being the time slice with the lowest [CO 2 ] atm . F I G U R E 1 Simulated selected interglacial biome patterns. Inferred biomes are mapped for four selected marine oxygen isotope stages/ sub-stages that approach or reach global ice volume minima similar to that of the Holocene, and that thus can be considered interglacials. Time slices shown are those closest in age to the ice-volume minimum of each marine oxygen isotope stage/sub-stage. The four marine oxygen isotope stages/sub-stages shown are selected to illustrate the range of variation in interglacial biome patterns. (Maps drawn using the Robinson Projection.) As with the interglacials, examination of the complete series of global biome maps (Supporting Information, Slideshow) reveals further regional differences between glacials with respect to simulated temporal patterns of biome changes. Such differences in temporal patterns are again consistent with palaeovegetation evidence from those rare locations spanning two or more glacial stages.    Laskar et al., 2004); global ice volume, as reflected by a global benthic δ 18 O stack (Lisiecki & Raymo, 2005); and atmospheric CO 2 concentration (Bereiter et al., 2015), both dissimilarity value series being plotted against age for the 206 past time slices.  Points are coloured to indicate to which of the six clusters, identified using K-means cluster analysis, they belong: cluster 1-blue; cluster 2-cyan; cluster 3-yellow; cluster 4-magenta; cluster 5-orange; cluster 6-pink.

| Overall history of global biome patterns 800 ka to present
relatively high global ice volumes ( Figure 5), typical of stadial, but not full glacial, conditions.

| Global biome extents and above-ground vascular plant biomass 800 ka to present
In addition to changes in location and regional extents of major biomes (Figures 1 and 2 Although we estimate total ice-free land area to have varied by <±3% since 800 ka (Figure 6), and Tropical Evergreen Forest, accounting for 38% of simulated total global carbon mass under present conditions, is simulated to have ranged only between 79% and 107% of its present simulated extent, simulated total global carbon mass ranges between 54% and 102% of its present simulated value Globally, simulated carbon mass per unit area, however, ranges between 53% and 102% of its present value, with a range of 73-103%

| Past biome patterns and present biological diversity patterns
Previous studies have reported relationships between contemporary species richness and past biome constancy (Colville et al., 2020;Huntley et al., 2016Huntley et al., , 2021 and have identified a relationship between tree species-richness in the temperate forests of North America, Europe and eastern Asia and glacial forest biome extents in those regions (Huntley, 1993). Our simulations quantify biome extent variations since 800 ka (Table 5; Table S7), allowing exploration of the degree to which contemporary species richness is more generally related to past biome extent variation. Some obvious parallels emerge between relative species richness and past global biome extent variation. Tree species richness (Keppel et al., 2021), for example, is F I G U R E 5 K-means cluster membership in relation to global ice volume. Time series of K-means cluster membership of time slices, each time slice being represented by a dot, coloured according to cluster membership as in Figure 4, and plotted on a vertical scale reflecting cluster number, the time slices being connected in temporal order by the thin black line. To provide context, the global ice volume, as reflected by the LR04 benthic δ 18 O stack (Lisiecki & Raymo, 2005), is shown by the red line.   Four progressively darker shades of orange indicate values in the ranges: >105% -150%; >150% -250%; >250% -500%;
a Zero extent simulated for present prevents calculation of past relative values. and decreases for Tropical Evergreen Forest, Temperate Shrubland,

Temperate Needle-leaved Evergreen Forest and Boreal Evergreen
Needle-leaved Forest (Figure 6; Figure S7) contribute to the large dissimilarities. Although there is also an increase in unclassified grid cells, this has very limited impact because even at their maximum, in the 2100 projection under RCP8.5, they account for only 1.8% of global ice-free land area.

| Global biome patterns and the paradigm hypothesis
Our results support the hypothesis that ascribes the differing spa-

| Global biome extents and above-ground vascular plant biomass
The direct physiological effects upon plants (Urban, 2003), and hence vegetation, of changes in [CO 2 ] atm have received considerable attention in the context of the accelerating increase in [CO 2 ] atm over the past two centuries, and especially since the mid-twentieth century, resulting from anthropogenic activities, especially burning of fossil fuels (Keeling & Graven, 2021). These effects have also received attention from Quaternary palaeoecologists, although it is the impacts of lower [CO 2 ] atm during glacial stages that has been their principal focus. Lower [CO 2 ] atm favours C4 species relative to C3 species, leading to increased extents of C4 grasslands in at least some tropical habitats during glacial stages (Urban et al., 2015). It also favours non-woody relative to woody taxa, leading to a shift in their competitive balance, especially in semi-arid or seasonally-arid regions where wildfire is a characteristic feature of the ecosystems (Bond et al., 2003;Bond & Midgley, 2012). In addition, these effects, especially the relationship between [CO 2 ] atm and plant water use efficiency, result in biases in palaeoclimate reconstructions made from palaeovegetation data if they are not taken into account (Prentice et al., 2017). LPJ-GUESS has been shown to reproduce well both FACE experiment results (e.g. Hickler et al., 2008;Zaehle et al., 2014) and estimated vegetation productivity changes 1901-2015 (Kastner et al., 2022). However, without nitrogen limitation, the version used here may overestimate the positive response of photosynthesis to increased [CO 2 ] atm . Nonetheless, photosynthesis increases almost linearly with leaf-internal CO 2 pressure over the glacial-interglacial range of [CO 2 ] atm (Gerhart & Ward, 2010;Kgope et al., 2010), whereas the CO 2 effects show increasing saturation under present and future [CO 2 ] atm Long et al., 2004). Thus, substantial plant-physiological CO 2 effects on global biome patterns are likely over the past 800 ky, especially where C4 grasses compete with C3 trees.
Our results are consistent with palaeovegetation evidence of shifting balances between C4 and C3 grasses and between woody and non-woody plants as [CO 2 ] atm varied during the mid-and late-   (Table S2)  LGM and Holocene global vegetation patterns and carbon mass.
Our simulations, spanning the nine glacial-interglacial cycles since 800 ka, considerably extend and reinforce this evidence.

| Past biome extents and contemporary biodiversity
Species origination through evolutionary diversification and species loss through extinction are two principal, and opposing, mechanisms contributing to determining species richness in any given region, biome or taxonomic group. Higher diversities of species will accumulate where origination rates are high and/or extinction rates are low. Varying mid-and late-Pleistocene global climate and biome patterns likely had substantial but regionally varying impacts upon both rates. Fragmentation of species' distributions during intervals when climatic conditions suitable for their persistence were available only in restricted and isolated geographical areas could potentially lead F I G U R E 8 Non-metric multidimensional scaling analysis of global biome pattern dissimilarities, including projected future biome patterns. Plot of non-metric multidimensional scaling axis 1 (horizontal) versus axis 2 (vertical) for an analysis including the potential future biome patterns for 2050 and 2100 under the RCP4.5 and RCP8.5 concentration pathways, as well as all past time slices. Each point represents a time slice or potential future and is labelled with its age (ka) or concentration pathway and year (CE). Points representing past time slices are coloured to indicate to which of the six clusters, identified using K-means cluster analysis, they belong: cluster 1-blue; cluster 2-cyan; cluster 3-yellow; cluster 4-magenta; cluster 5-orange; cluster 6-pink. Points representing potential future biome patterns are coloured red for 2050 projections and purple for 2100 projections. RCP, representative concentration pathway.
to evolutionary divergence of the isolated populations, and hence increase species origination rate (Kadereit & Abbott, 2021). Less obviously, intervals offering suitable climatic conditions more extensively potentially allowed previously isolated but related species to exchange genetic material. Especially in the case of plants, new species can arise from such hybridisation, once again increasing species origination rate (Kadereit & Abbott, 2021).
Species extinctions can result from various mechanisms related to changing environmental conditions (Huntley, 1998). In the present context, two mechanisms are perhaps most relevant. First, species unable to track the shifting geographical location of suitable conditions are very likely to face a severely increased extinction risk.
Furthermore, this risk increases the more rapidly the environment, and hence the geographical location of suitable conditions, changes.
Extinction risk increase is also likely to be greater for species with more specialised habitat requirements, such as those limited to a single biome (Huntley et al., 2016). Long-term survival of such species will be greatest, and hence species richness tend to be higher, where the required habitat or biome is able to persist, notwithstanding changing environmental conditions (Colville et al., 2020;Huntley et al., 2016Huntley et al., , 2021. The second mechanism relates to both the widely observed species-area relationship and the increased extinction risk faced by a species if its geographical range, and hence its population, becomes restricted. Thus, according to the species-area relationship, species richness of a biome in a continental region where its mean mid-and late-Pleistocene extent was less than its present extent is expected to be less than predicted by that present extent. Species richness of the same biome in a continental region where its mean mid-and late-Pleistocene extent was similar to, or greater than, its present extent is expected to conform to, or exceed, that predicted from its present extent. In addition, extinction risk will have been higher for species associated with a biome that, in a given continental region, had a substantially reduced extent during at least one interval of the mid-and late-Pleistocene, and less for species associated with a biome that has not suffered any such severe reduction(s) in extent. Such contrasting species extinction risks likely resulted in differing rates of species loss. Hence, lower than expected species richness is predicted for biomes in continental regions where they have suffered one or more extreme reduction(s) in extent, compared to the same biome in continental regions where such extreme reduction in extent has not occurred.
Our results support the contributions to determining present species richness of the species-area relationship and mean past biome extent in different continental regions, as well as of the differential species extinction risk in biomes experiencing differing degrees of past extent reduction in different continental regions.
These mechanisms are independent of species' relative capacities to shift their ranges in response to climatic changes, albeit that changes in a biome's extent in a continental region will often be associated with changes in its geographical location. This is especially the case in higher northern latitudes where biome constancy is generally low, as is species richness, contributing to the lack of any relationship between biome constancy and species richness in less diverse regions (Huntley et al., 2021). In these latitudes, a biome's past extent in a given continental region apparently has had a more marked influence on its present species richness than has its lack of areas of constancy. Nonetheless, our results also show past extent is important in regions where biomes have had substantial constancy.
Tropical Evergreen Forest, for example, has large areas of high constancy (Huntley et al., 2021), but also contrasting degrees of past reduction in extent, and parallel patterns in present tree diversity, in its three principal continental regions (Keppel et al., 2021).
They found the LGM simulation had the greatest dissimilarity from the present simulation, followed by the LIG and mid-Holocene in the opposite direction to that of the latter time slices (Figure 8).
One potential reason for this contrast is the differing [CO 2 ] atm values used in the simulations being compared. Whereas the present and doubled [CO 2 ] atm simulation pairs used by Sykes et al. (1999) used 300 and 600 ppmv (Manabe & Wetherald, 1987), and 326 and 652 ppmv (Schlesinger & Zhao, 1989) Long et al., 2004) will contribute to lesser simulated responses of global biomes to these future increases in [CO 2 ] atm than to equivalent past increases. Nonetheless, the large magnitude deviations from the regression line for future biome patterns suggest that this results principally from future climate projections taking into account only relatively faster Earth system components, and excluding slow components, in particular ice-sheet dynamics, that are likely to take millennia to reach a new equilibrium, and in doing so will lead to substantial further climatic and global biome pattern changes. Simulated potential future biome patterns, that are without any analogue during at least the past 800 ky, would have substantial negative implications for global biodiversity. The likelihood of even greater longer-term changes raises the possibility of very severe losses of global biodiversity because temperate and boreal biomes would likely be markedly reduced in extent, if not lost completely.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The following files can be accessed at https://doi.org/10.5061/dryad.