Understanding Spatiotemporal Variation in Richness and Rate of Within-Site Turnover for Vegetation Communities in Western Eurasia over the Last 4000 Years

: Vegetation communities are intricate networks of co-occurring species. Logistical challenges in collecting primary data means research often utilises short-term data from restricted geographical areas. In this study, we examine spatiotemporal change in richness and turnover of vascular plants and bryophytes over the last 4000 years at 23 sites in western Eurasia using high-resolution palaeoe-cological data. We ﬁnd support for the Latitudinal Diversity Gradient and Altitudinal Diversity Gradient in both the overall vegetation community (arboreal and non-arboreal species) and the shrub and herb sub-community (non-arboreal species only), as well as a signiﬁcant temporal increase in the gradient of both relationships. There was a temporal increase in (alpha) richness; the rate of turnover was high but temporally consistent for the overall vegetation community and high but decreasing over time for the shrub and herb sub-community. The rate of change in turnover was affected by latitude (steeper negative relationship at higher latitudes) and altitude (steeper negative relationship at lower altitudes). The Diversity-Stability Hypothesis was supported: vegetation communities changed from “lower richness, higher turnover” historically to “higher richness, lower turnover” more recently. Causal mechanisms for these complex interlinked biogeographical patterns remain ambiguous, but likely include climate change, non-native introductions, increasing homogenisation of generalist taxa, landscape simpliﬁcation, and anthropogenic disturbance. Further research into drivers of the spatiotemporal patterns revealed here is a research priority, which is especially important in the context of biodiversity decline and climate change.


Introduction
Vegetation communities are intricate networks of co-occurring species that interact with one another and their abiotic environment in complex ways that vary over both time and space.Understanding spatiotemporal patterns is thus integral to the study of community ecology, especially given the pressing contexts of climate change and biodiversity loss [1][2][3].However, such patterns can be challenging to study over ecologically relevant timescales because most primary ecological data are relatively recent [4,5].Several studies have demonstrated the utility of high-resolution palaeoecological data for analysing community dynamics over longer timescales than would have been possible using primary ecological data [6][7][8][9][10][11].In addition to such approaches being helpful in understanding temporal change in ecological communities at specific sites, use of "recent past" palaeoecological data from multiple sites has considerable potential for exploring spatiotemporal patterns in ecological parameters [12,13].This is a powerful approach since it enables consideration of spatiotemporal patterns in community ecology across timespans that are meaningful for natural shift within the taxa under consideration [14][15][16].This also allows patterns in richness and within-site turnover to be studied throughout the transition from relatively low human environmental impacts to the more intensive pressures of the present day [17,18].
Richness is influenced by many different geographical and environmental variables.The Latitudinal Diversity Gradient (LDG) suggests richness is typically highest at the equator and decreases towards the poles, while the Altitudinal Diversity Gradient (ADG) suggests that richness is highest at low/mid altitudes and decreases towards higher altitudes.Evidence to support the LDG occurs in most taxonomic groups but is not universal (see meta-analyses [19,20]).Similarly, numerous studies have confirmed the presence of an ADG, although the exact pattern and specific elevations involved differ depending on location and the taxa studied (review by [21]).Temporal patterns in richness have not received as much attention as spatial patterns [22].Analysis of temporal patterns has, until recently, focused on assessing increases in (alpha) richness due to a net positive in the balance between introductions of non-native species and extirpation of native species (e.g., [23]).Two highly cited studies using primary survey data from constant-effort biodiversity monitoring programmes across the world have, contrary to expectations, concluded that alpha richness is increasing almost universally at local scales [24,25].However, this research has been criticised (e.g., [26,27]) for being short-term and data-poor (typically 13-20 years) and using data across time periods that differ in duration, actual years analysed, and relative position to external stressors such as climatic change.Expanding the time period under consideration using palaeoecological data, Berglund et al. [6] and Colombaroli and Tinner [18] have shown similar patterns of richness increase over time.However, while these studies have a key strength in their duration, their spatial resolution is limited.Multi-site palaeoecological data at inter-continental scale, covering a standardised time period at high resolution, provide a unique opportunity to understand broadscale patterns in community ecology.A dataset that is extensive both spatially and temporally also enables the LDG and ADG to be empirically tested, thereby enabling consideration of potential interactions between spatiotemporal factors.Such interactions are understudied but intriguing results from Silvertown [28] and Giesecke et al. [29] have showed that the strength of the LDG across Europe increased over the last 13,000 years when studied in 1000-year and 3000-year increments, respectively.
Within-site turnover is inherently temporal as it quantifies change in the species that are present in the community at different time periods.However, the rate of within-site turnover could vary over time because of factors that vary with time such as climate change [30][31][32][33][34], direct human influence [35]; and non-native introductions [23,36].Moreover, temporal patterns in the rate of within-site turnover might be mediated by spatial variation.For example, Korhonen et al. [37] demonstrated that interannual turnover in aquatic ecosystems was greater at higher latitudes than it was at lower latitudes.Advancing understanding of such patterns is increasingly important given the pressing contexts of human-accelerated climate change and other anthropogenic impacts.Most previous research has involved short-term studies of primary data, which is problematic as measured turnover can covary with study duration [37].Study duration has previously been extended using innovative space-for-time substitution [30], but a recent study by Mottl et al. [38] has demonstrated that using palaeoecological data to analyse rate of change has untapped potential.
Finally, as well as there being spatiotemporal patterns in richness and the rate of within-site turnover individually, richness and turnover can themselves inter-relate both overall and in ways that covary with time and/or space.The Diversity-Stability Hypothesis (DSH) states that higher richness is associated with higher ecological resilience and thus lower turnover.The DSH was originally postulated by Darwin [39] before being extended by MacArthur [40] and McNaughton [41].The relationship is usually considered to be the result of functional linkages and interactions (especially competition and trophic connections) that are greater when richness is higher.This then decreases invasion and extirpation rates and increases resilience to disturbance, which results in stable ecological states.There is considerable experimental evidence for the DSH, at least over short timescales (e.g., [42,43]).However, the pattern is contentious and complex and an increasing number of studies are finding no significant relationship or a significant opposite trend (see review by Ives and Carpenter [3]).Quantifying patterns at inter-continental scales over longer timescales would be beneficial.
In this study, we use pollen records as a proxy for vegetation over the last 4000 years at 23 sites in western Eurasia.Specifically, we obtain robustly dated high-resolution palaeoecological data on vascular plants and bryophytes to quantify spatiotemporal change in richness and the rate of within-site turnover in the vegetation community.Our aims are threefold: (1) to analyse spatiotemporal change in richness in western Eurasia to test the latitudinal and altitudinal diversity gradients, quantify temporal patterns, and establish whether there are significant interactions between temporal and spatial variables; (2) to quantify temporal change in the rate of within-site turnover to establish whether, on average, turnover is accelerating, slowing, or static given the increasing human pressures on natural systems, and whether temporal patterns are mediated by spatial variables; and (3) to relate the rate of within-site turnover to richness to test the diversity-stability hypothesis on a dataset that is extensive both spatially and temporally.In this way, we seek to improve understanding long term processes in ecology by providing a hitherto understudied historical dimension to the LDG, ADG, and DSH.

Data Acquisition
The 2020 sites listed on the European Pollen Database [44] at the time of the study were screened to identify all sites where high-resolution palynological data over the last 4000 years were available electronically.Specifically, sites selected for inclusion needed to have: (1) been used within a published output in a peer-reviewed journal or peer-examined doctoral thesis; (2) at least two radiocarbon dates within the last 4000 years forming a robust age-depth model of sedimentation to allow precise time-horizons to be identified by the original author(s); (3) high-resolution sampling at ≤5 cm intervals; and (4) data from at least one time period every 100 years over the 4000-year study period (i.e., 0-99 years ago; 100-199 years ago; 200-299 years ago and so on).This last criterion was ultimately relaxed to allow missing data in up to four time periods (10% of the total) to increase the number of eligible sites.This protocol gave 23 sites from 15 countries (latitude: 63 • 53 0 to 37 • 39 55 North; longitude: 7 • 59 35 W to 46 • 37 55 East; altitude: 1 m to 2410 m; Figure 1).Of these sites, eight had data for all time periods while 15 had missing data up to, but not exceeding, 10% of the total.

Data Manipulation
Palynological data (vascular plants and bryophytes) were extracted for the time periods of interest at each of the 23 sites using the taxonomic level given by the original author(s).As the taxonomic level to which pollen and spores could be identified varied for different taxa based on their morphological distinctiveness, the overall dataset contained family-level, genus-level and species-level data (hence we use "taxonomic" and "taxa" rather than "species" throughout this research).It is important to note that there were no systematic differences in data over time (i.e., taxonomic resolution was temporally consistent for each site) such that subsequent analysis of temporal patterns was not an artefact of changing taxonomic resolution.Moreover, all sites contained taxa identified to species-level, genus-level and family-level, such that there was no systematic bias in taxonomic resolution between sites that would confound spatial analysis.Two biodiversity metrics were calculated.Taxonomic richness calculated for each site for each time period as a simple count of taxa recorded.Within-site turnover was calculated on presence/absence data for successive sequential time periods to give a series of up to 39 values for each site describing temporal change in the vegetation community at that site (i.e., 0-99 years compared to 100-199 years, 100-199 years compared to 200-299 years etc); it should be noted that while 39 values were present for sites were there were no missing data, some sites had fewer values where data for specific time periods were missing.It is recognised that using presence/absence data to quantify turnover is less sensitive than measures of turnover that quantify proportional change in taxa.However, this was a deliberate decision to ensure that analysis of change in vegetation communities was due to species invasion and extirpation rather than compositional change in vegetation communities due to change in relative abundance, since this can be substantially skewed by relatively minor changes in taxa that dominate the palynological dataset [45].Thus, although a presence/absence matrix is less information rich, it is ecologically and statistically more robust.All withinsite turnover calculations were based on Jaccard's Coefficient of Community Similarity (CC j ): this approach has been used previously to assess changes in community ecology in both ecological and palaeoecological contexts (e.g., [34,36,46,47]).Here, we calculated the inverse of CC j (sometimes called Jaccard Distance) using Equation (1), where c was the number of taxa common to two successive time periods and S was the sum of the total number of taxa found in both time periods.This gave a measure of vegetation community dissimilarity, running between 0 and 1 where higher numbers were indicative of greater dissimilarity, and thus greater within-site turnover.Sites included in this study.All palaeoecological data were collated from the European Pollen Database [44]; data are accessible as per the data accessibility statement at the end of this paper.

Data Manipulation
Palynological data (vascular plants and bryophytes) were extracted for the time periods of interest at each of the 23 sites using the taxonomic level given by the original author(s).As the taxonomic level to which pollen and spores could be identified varied for different taxa based on their morphological distinctiveness, the overall dataset contained family-level, genus-level and species-level data (hence we use "taxonomic" and "taxa" rather than "species" throughout this research).It is important to note that there were no systematic differences in data over time (i.e., taxonomic resolution was temporally consistent for each site) such that subsequent analysis of temporal patterns was not an artefact of changing taxonomic resolution.Moreover, all sites contained taxa identified to species-level, genus-level and family-level, such that there was no systematic bias in taxonomic resolution between sites that would confound spatial analysis.Two biodiversity metrics were calculated.Taxonomic richness calculated for each site for each Equation (1): To allow better understanding of change over time, taxonomic richness and within-site turnover were calculated for the overall vegetation community and also, separately, for the shrub and herb (hereafter SAH) sub-community.This allowed separate consideration of change in non-arboreal taxa (mainly shorter-lived r-strategists) compared to the overall vegetation community, which is strongly influenced by arboreal taxa (longer-lived,

K-strategists).
There is inevitable variability in the number of pollen grains counted by different researchers, which can confound between-site comparisons as the number of species found (taxonomic richness) is fundamentally dependent on sample effort (number of grains counted) [48,49].There are "gold standard" methods to address these limitations, involving rarefraction or random sub-sampling based on the lowest count.However, such approaches are time consuming, result in the loss of a large amount of data, and, most importantly, are not always necessary as long as the lowest sample effort is above a statistically derived minimum threshold to be considered robust and spatiotemporally comparable.Research in the early 1990s [50,51] suggested that ≥ 300 pollen grains should be counted for an accurate vegetation profile.More recently, statistical modelling by Djamali and Cilleros [52] has shown that a count of 230 grains yields 95% reliability in the pollen profile (i.e., using the concept of α = 0.05, counting above 230 grains has a very limited effect on the vegetation profile).All sites included in this study had a high pollen count (number of grains counted per time period per site: average = 2175; minimum = 336), Given this, and because taxonomic richness and within-site turnover metrics used presence/absence data rather than abundance data, variation in total number of pollen grains counted per time period was not considered likely to bias analysis and thus rarefraction or random sub-sampling was not deemed necessary.

Data Analysis
Univariate relationships between taxonomic richness, latitude and altitude were visualised and analysed using simple bivariate regressions that were undertaken on raw data.Then, to analyse the combined effects of space and time a Generalised Linear Model (GLM) was performed.This GLM used taxonomic richness as the dependent variable with latitude, altitude and time period entered as continuous factors: the interactions between each spatial variable and time were also added (i.e., latitude*time and altitude*time).This analytical framework allowed for the relative importance of spatiotemporal variables to be ascertained as latitude, altitude and time were entered into a single model; the inclusion of the interactions terms was necessary to allow quantification of temporal change in the LDG and ADG, respectively, which was one of our focal aims.The model used a Poisson distribution with a log link function, which was appropriate for the count data being analysed and was the optimal model based on comparison of Akaike's Information Criterion scores relative to other statistically valid model distribution/function combinations such as negative binomial models [53,54].A second GLM was created for taxonomic richness for the SAH sub-community only.
Temporal relationships for within-site turnover were plotted graphically and analysed using simple bivariate regressions that were undertaken on raw data.Then, to determine statistically whether the rate of within-site turnover was accelerating, slowing, or static given the increasing human pressures on natural systems, and the possible mediating effects of latitude and altitude on any temporal patterns, two GLMs were performed.These models used within-site turnover data for the overall vegetation community (model 1) and within-site turnover data for the SAH sub-community (model 2) as the dependent variable modelled against time, time*latitude, and time*altitude.Scale models were created because, although turnover was technically proportional, each proportion was the result of a complex calculation (Jaccard Distance) rather than arising out of binomial events (such that a binomial distribution was not appropriate); both dependent variables were also normally distributed.A gamma distribution with a log link function was used as this provided a better fit (considerably lower Akaike's Information Criterion scores) relative to a standard linear scale model for both overall vegetation community and the SAH sub-community.
Finally, to relate within-site turnover to richness, and thus test the Diversity-Stability Hypothesis, we calculated the relationship between the average turnover scores for each time period and regressed these values against taxonomic richness for the preceding time period.This provided an initial understanding of the relationship between turnover and richness in western Eurasia over the last 4000 years.To explore this further, Principal Components Analysis (PCA) was undertaken to describe the relationship between turnover (y) and richness (x) so that this relationship could itself be related to time (z) to quantify temporal mediating effects.All analyses were conducted in SPSS version 27 (USA: IBM).

Spatiotemporal Change in Taxonomic Richness: Univariate Patterns
There was a significant negative relationship between taxonomic richness and latitude as predicted by the LDG for the overall vegetation community (F 1,866 = 92.405;p < 0.001; R 2 = 0.139) and the SAH sub-community (F 1,865 = 90.430;p < 0.001; R 2 = 0.149).The gradient of the relationship was slightly steeper for the overall vegetation community compared to the SAH sub-community (Figure 2a).There was also a significant negative relationship between taxonomic richness and altitude as predicted by the ADG for the overall vegetation community (F 1,866 = 60.978.405;p < 0.001; R 2 = 0.084) and the SAH subcommunity (F 1,865 = 84.081;p < 0.001; R 2 = 0.121).The gradients of the relationships were approximately equal for the overall vegetation community and the SAH sub-community (Figure 2b), such that the lines of best fit were parallel.
Hypothesis, we calculated the relationship between the average turnover scores for each time period and regressed these values against taxonomic richness for the preceding time period.This provided an initial understanding of the relationship between turnover and richness in western Eurasia over the last 4,000 years.To explore this further, Principal Components Analysis (PCA) was undertaken to describe the relationship between turnover (y) and richness (x) so that this relationship could itself be related to time (z) to quantify temporal mediating effects.All analyses were conducted in SPSS version 27 (USA: IBM).

Spatiotemporal Change in Taxonomic Richness: Univariate Patterns
There was a significant negative relationship between taxonomic richness and latitude as predicted by the LDG for the overall vegetation community (F1,866 = 92.405;p < 0.001; R 2 = 0.139) and the SAH sub-community (F1,865 = 90.430;p < 0.001; R 2 = 0.149).The gradient of the relationship was slightly steeper for the overall vegetation community compared to the SAH sub-community (Figure 2a).There was also a significant negative relationship between taxonomic richness and altitude as predicted by the ADG for the overall vegetation community (F1,866 = 60.978.405;p < 0.001; R 2 = 0.084) and the SAH subcommunity (F1,865 = 84.081;p < 0.001; R 2 = 0.121).The gradients of the relationships were approximately equal for the overall vegetation community and the SAH sub-community (Figure 2b), such that the lines of best fit were parallel.There was a significant positive relationship between taxonomic richness and time (measured as years before present in 100-year intervals), such that richness was significantly higher in more recent times than it was historically.This trend held for both overall vegetation community (F 1,866 = 30.323;p < 0.001; R 2 = 0.622) and the SAH sub-community (F 1,865 = 37.471; p < 0.001; R 2 = 0.745).The gradient of the line of best fit, which was indicative of the magnitude of change, was almost identical: 0.0028 higher richness per 100 years (11.2 taxa over the entire study period) versus 0.0029 higher richness per 100 years (11.6 taxa over the entire study period), respectively, such that the lines of best fit were parallel (Figure 3).As regards site-specific trends (not shown graphically), there was a significant positive relationship between taxonomic richness and time within the overall vegetation community (i.e., higher richness in more recent time periods) at 17 of 23 sites (74%); the remaining sites were split evenly between weak positive trends (three sites; 13%) and weak negative trends (three sites; 13%), all of which were non-significant.For the SAH sub-community, there were significant positive relationships at 12 of 23 sites (52%), with the remaining 11 sites (48%) being non-significant.
maining sites were split evenly between weak positive trends (three sites; 13%) and weak negative trends (three sites; 13%), all of which were non-significant.For the SAH subcommunity, there were significant positive relationships at 12 of 23 sites (52%), with the remaining 11 sites (48%) being non-significant.

Spatiotemporal Change in Taxonomic Richness: Multivariate Modelling
A multivariate analysis for taxonomic richness that considered all spatiotemporal variables and the interactions between them revealed complex patterns (Table 1).As expected given the univariate modelling, time was significant as was latitude, both in the expected directions.Once latitude and time were factored in, there was support for the ADG for the overall vegetation community but not the SAH sub-community, suggesting that this pattern was driven by changes in arboreal taxa rather than non-arboreal taxa.Interestingly, there was also considerable evidence of spatiotemporal interactions,

Spatiotemporal Change in Taxonomic Richness: Multivariate Modelling
A multivariate analysis for taxonomic richness that considered all spatiotemporal variables and the interactions between them revealed complex patterns (Table 1).As expected given the univariate modelling, time was significant as was latitude, both in the expected directions.Once latitude and time were factored in, there was support for the ADG for the overall vegetation community but not the SAH sub-community, suggesting that this pattern was driven by changes in arboreal taxa rather than non-arboreal taxa.Interestingly, there was also considerable evidence of spatiotemporal interactions, whereby time changed the strength of spatial relationships: both the LDG and ADG were steeper in more recent times than they had been historically.

Spatiotemporal Change in Rate of Within-Site Turnover
There was no relationship between the rate of within-site turnover in the overall vegetation community relative to time (F 1,813 = 1.127; p < 0.001; R 2 = 0.059), with a line of best fit that was almost flat (Figure 4).However, there was a statistically significant temporal trend in the rate of within-site turnover within the SAH sub-community (F 1,813 = 9.477; p < 0.001; R 2 = 0.283), such that the rate of turnover was significantly lower in more recent times than it was historically.The gradient of the line of best fit for the SAH sub-community was −0.0012 (i.e., within-site turnover changed by 0.12% for every 100 years or ~5% over the entire study period; Figure 4).Undertaking analysis on the rate of within-site turnover temporally for each of the 23 sites individually (not shown graphically) resulted in a mixed picture, with direction and magnitude of change varying considerably between sites.Gradients ranged from −0.0103 to 0.0044 (mean = −0.00004)for the overall vegetation community and from −0.0126 to 0.0051 (mean = −0.0012)for the SAH sub-community).Thus, while the average rate of turnover was significantly lower in more recent times for the SAH sub-community (Figure 4), two sites (8%) had a significant opposite trend whereby rate of turnover was higher in more recent times, and eight sites (35%) had no significant pattern.This suggested that spatial variables might be mediating temporal change.This was confirmed via multivariate modelling that showed significant interactions for time*latitude and time*altitude (Table 2); time itself was only significant for the SAH sub-community as indicated by the baseline pattern shown in Figure 4.
Table 2. Generalised Linear Models of within-site turnover in relation to time and mediating effects of latitude and altitude.Chi square values are likelihood ratio for overall models and Wald for fac- Undertaking analysis on the rate of within-site turnover temporally for each of the 23 sites individually (not shown graphically) resulted in a mixed picture, with direction and magnitude of change varying considerably between sites.Gradients ranged from −0.0103 to 0.0044 (mean = −0.00004)for the overall vegetation community and from −0.0126 to 0.0051 (mean = −0.0012)for the SAH sub-community).Thus, while the average rate of turnover was significantly lower in more recent times for the SAH sub-community (Figure 4), two sites (8%) had a significant opposite trend whereby rate of turnover was higher in more recent times, and eight sites (35%) had no significant pattern.This suggested that spatial variables might be mediating temporal change.This was confirmed via multivariate modelling that showed significant interactions for time*latitude and time*altitude (Table 2); time itself was only significant for the SAH sub-community as indicated by the baseline pattern shown in Figure 4. Negative in both models: general negative effect of time on rate of turnover was greater for sites at lower altitudes (steeper negative relationship) and weaker for sites at higher altitudes (shallower negative relationship).At very high altitudes, the relationship became positive in some cases.

Relationship between Richness and Rate of Within-Site Turnover: Is the Diversity-Stability Hypothesis Supported?
There was a significant negative relationship between taxonomic richness and withinsite turnover for both the overall vegetation community (Figure 5a) and the SAH subcommunity (Figure 5b).This meant that turnover was lower when richness was higher, exactly as predicted by the DSH.The relationship was steeper and stronger for the SAH sub-community, with a gradient of −0.0044 (or 0.44% less turnover per additional taxon in the community) and an R 2 of 0.279, compared to the overall vegetation community (gradient = −0.0018;R 2 = 0.101).
Based on visual analysis of Figure 5, the driver for the significant relationship between turnover and richness (and thus the driver for the support found in this study for the DSH) appeared to be time-or, possibly more accurately, other variables that change temporally (see Discussion).The position of the datapoints along the line of best fit strongly linked to time, with older samples being "higher turnover, lower richness" and more recent samples being "lower turnover, higher richness".To explore this further, the diagonal position of the datapoints was quantified by using Principal Components Analysis to create a new axis, PC1, which approximated the original line of best fit by passing through the centroid of the data cloud on the longest axis (i.e., such that the PC1 score for each datapoint was broadly akin to the position of that datapoint on the lines of best fit shown in Figure 5a,b).PC1 was statistically negatively related to time expressed as years before present for both overall vegetation community (F 1,37 = 25.045,p < 0.001, R 2 = 0.404) and the SAH sub-community (F 1,37 = 64.658,p < 0.001, R 2 = 0.636): this meant that the relationship between richness (x) and turnover (y) was causally related to time (z).This result is consistent with both the temporal pattern for richness (higher in more recent time periods; Figure 3) and turnover (lower in more recent time periods; Figure 4).became positive in some cases.

Relationship between Richness and Rate of Within-Site Turnover: Is the Diversity-Stability Hypothesis Supported?
There was a significant negative relationship between taxonomic richness and within-site turnover for both the overall vegetation community (Figure 5a) and the SAH sub-community (Figure 5b).This meant that turnover was lower when richness was higher, exactly as predicted by the DSH.The relationship was steeper and stronger for the SAH sub-community, with a gradient of −0.0044 (or 0.44% less turnover per additional taxon in the community) and an R 2 of 0.279, compared to the overall vegetation community (gradient = −0.0018;R 2 = 0.101).Based on visual analysis of Figure 5, the driver for the significant relationship between turnover and richness (and thus the driver for the support found in this study for the DSH) appeared to be time-or, possibly more accurately, other variables that change temporally (see Discussion).The position of the datapoints along the line of best fit strongly linked to time, with older samples being "higher turnover, lower richness" and more recent samples being "lower turnover, higher richness".To explore this further, the diagonal position of the datapoints was quantified by using Principal Components Analysis to create a new axis, PC1, which approximated the original line of best fit by passing through the centroid of the data cloud on the longest axis (i.e., such that the PC1 score for

Discussion
Analysis of paleoecologically derived data from 23 sites in western Eurasia and spanning a 4000-year period support some of the major spatial relationships (LDG, ADG) and interrelationships (DSH) within community ecology.More importantly, they also demonstrate the integral importance of temporal processes, either directly (taxonomic richness and rate of within-site turnover) or as mediating effects on LDG, ADG and DSH relationships to create the complex spatiotemporal patterns summarised in Figure 6.
Our data support the LDG and ADG, with taxonomic richness being greater at lower latitudes and altitudes.Such patterns are common although not universal, as previously shown by Hillebrand [19], Kinlock et al. [20] and Fischer et al. [21] for primary ecological data and Reitalu et al. [55] for palaeoecological data.Where these patterns occur, they can be driven by greater loss of taxa at higher latitudes and altitudes and/or greater gain of taxa at lower latitudes and altitudes [2].However, our analysis also highlights that time interacts with these relationships, as they have both become steeper over time.The casual mechanisms for this are unclear.In one of the very few previous analyses of temporal change in the LDG, Silvertown [28] found the latitude-richness gradient became steeper over time for woody taxa across Europe over the last interglacial period when analysed in 1000-year increments.Silvertown ascribed natural spatiotemporal change in climate to be the causal mechanism and concluded that the LDG steepening was due to flora readjusting after the end of the last ice age.This was tangentially supported by research that demonstrated the lack of any LDG immediately after prehistoric mass extinction events (MEEs), which were followed by strengthening gradients over time [2,56].This includes patterns for North American mammals following the Cretaceous/Palaeocene MEE [57] and marine taxa following the Permian/Triassic MEE [58].However, as our current study spans 4000 years of relatively stable climate (albeit with human-accelerated climate change effects becoming apparent over the last couple of hundred years), climate change cannot fully explain the observed patterns and combinations of natural and anthropogenic processes are likely to be responsible as discussed below.
ods; Figure 3) and turnover (lower in more recent time periods; Figure 4).

Discussion
Analysis of paleoecologically derived data from 23 sites in western Eurasia and spanning a 4000-year period support some of the major spatial relationships (LDG, ADG) and interrelationships (DSH) within community ecology.More importantly, they also demonstrate the integral importance of temporal processes, either directly (taxonomic richness and rate of within-site turnover) or as mediating effects on LDG, ADG and DSH relationships to create the complex spatiotemporal patterns summarised in Figure 6.Our data support the LDG and ADG, with taxonomic richness being greater at lower latitudes and altitudes.Such patterns are common although not universal, as previously shown by Hillebrand [19], Kinlock et al. [20] and Fischer et al. [21] for primary ecological data and Reitalu et al. [55] for palaeoecological data.Where these patterns occur, they can be driven by greater loss of taxa at higher latitudes and altitudes and/or greater gain of taxa at lower latitudes and altitudes [2].However, our analysis also highlights that time interacts with these relationships, as they have both become steeper over time.The casual mechanisms for this are unclear.In one of the very few previous analyses of temporal change in the LDG, Silvertown [28] found the latitude-richness gradient became steeper over time for woody taxa across Europe over the last interglacial period when analysed in 1000-year increments.Silvertown ascribed natural spatiotemporal change in climate to be the causal mechanism and concluded that the LDG steepening was due to flora There are complex interactions between space and time in relation to richness [29].Here, not only does time mediate the spatial relationships in taxonomic richness, it has an important link to richness itself.The finding that alpha richness has increased over the recent past towards the present day supports the contentious previous studies using primary ecological survey data [24,25], but using much more robust datasets (longer time, higher resolution, more data rich, and consistent time periods analysed at all sites).This accords with longer-term palaeoecological studies [6,18,29,59,60] that have shown similar patterns of richness increase across the Holocene (the last 11,500 years).However, it should be noted that all of these studies using palaeoecological data-and indeed our analysis-are at a temporal resolution that will not allow detection of very recent (>100 years) taxonomic richness declines.Moreover, this baseline temporal increase should emphatically not be taken as evidence against the biodiversity crisis since site-specific (alpha richness) trends do not to preclude catastrophic loss of biodiversity over a wider scale (gamma richness) [26,27,61].Trends could also be transitory as communities often gain new taxa faster than they lose existing taxa [62].Moreover, it is possible that a component of the temporal increase in richness is non-desirable biotic homogenisation [23,29,63,64], as well as being complicated by extirpation lag [62].It is also important to highlight that temporal alpha richness increases are not universal: richness did not increase over time at 26% of sites (overall vegetation community) and 48% of sites (SAH sub-community).Some non-significant patterns tended towards the negative, especially at higher latitudes and/or altitudes).Because taxa differ in preservation potential, it is impossible to exclude the possibility that some of the apparent increase in richness over time could be an artefact rather than a genuine pattern.However, all palynological reconstructions are fundamentally based on the assumption that any bias in preservation is minimal, such that vegetation communities can be meaningfully reconstructed over much longer time periods than considered here [65].
The factors that covary with time, and thus might be responsible for the temporal patterns in taxonomic richness (and the mediating effects of time on the LDG and ADG), are unclear.Thus, while intriguing spatiotemporal patterns have been revealed by this paper, the causal mechanisms are more ambiguous.Speculative possibilities include increases in non-native plants [23] (although if this was the sole explanation an exponential pattern would be expected, with the rate of change increasing towards more recent times).A second possibility is change in landscapes caused by human activity is creating habitat mosaics that increase richness, as found in relation to more open landscapes being created in the mid-Holocene [18,59].A third possibility is an increase in the number of generalist taxa within communities over time, either with or without loss of more specialist taxa [63,64]: such a pattern could result from increasing landscape simplification [66] or anthropogenic disturbance [67].Finally, it is likely that this pattern is at least partly driven by humanaccelerated climatic change altering fundamental and realised ranges, especially the loss of taxa at high latitudes/altitudes and the poleward/upward expansion of other taxa into Europe at middle latitudes and altitudes) [68,69].The potential for climate change to be an important driver of the steepening LDG and ADG is tangentially supported by Rohde [70] and Erwin [71].The former suggested the causal mechanism for the LDG/ADG to be spatial differences in temperature that influence evolutionary speed (warmer temperatures at lower latitudes and altitudes = faster evolutionary rate = higher diversity), while the latter postulated a temporal relationship between temperature and richness.Taken together, it is possible that recent rising temperatures could be acting to steepen the LDG and ADG All of the above processes would involve arrival of new taxa at study sites over time and/or slower rates of loss (i.e., giving a net positive balance).For all of these speculative drivers, it might be expected that temporal change in richness would be non-linear (rather than linear across the 4000-year study period) as the drivers themselves are non-linear.It is possible, however, that biotic responses are possibly complicated by lag times meaning that, at any one time, a site might be experiencing extinction debt or immigration credit as per Jackon and Sax [62], or that multiple drivers are occurring together in ways that produce more complex temporal relationships.Although in our study rate of within-site turnover is decreasing temporally for the SAH sub-community, rates still remain at ~40% per hundred years.This means that temporal decrease in within-site turnover is not at odds with processes that involve turnover being responsible for richness increases over time-and indeed it could even be argued that a slowing of turnover would be expected if biotic homogenisation (which, by definition, cannot be infinite) is starting to plateau.
Given that higher richness is thought to drive higher stability (and thus lower turnover) [40][41][42][43], and that our data support this, it is not surprising that within-turnover has decreased over time while richness has increased over the same time period.Moreover, the main driver for the DSH in our data appears to be time.While natural mechanisms (functional diversity, competition and trophic connections) could be important as they decrease invasion and extirpation risk, human activity might also be a contributory factor, with increased management of ecosystems potentially partly supressing natural turnover.This accords with evidence from Vogel et al. [72], who demonstrated that a positive relationship between richness and stability in grassland communities only occurred in ecosystems that were highly managed, not those that were dominated by natural processes and interactions.It should be noted that there is also a potential artefact whereby turnover estimates might be higher when richness is lower as change in one or two taxa has a higher relative impact.
Throughout this study, and indeed all studies using palaeoecological data, there is an implicit assumption that vegetation diversity (richness and community composition) is equivalent to palynological diversity.Where this has been empirically tested, agreement between vegetation community and pollen community varies.For example, in tropical environments of Africa and South America there is no significant agreement between vegetation community and pollen community [73], whereas where this has been tested in northern Europe there are significant positive relationships between these two measures of vegetation [55].Although it recognised that pollen records will never be fully representative of vegetation communities, they are one of the best data sources to evaluate inter-continental changes in community ecology over periods longer than of those covered by primary ecological data.It should also be noted that there are inherent differences in the morphological distinctiveness of different pollen and spores, which affects the level of taxonomic resolution that is possible during microscopic identification, and thus all sites in this study contained taxa identified to species-level, genus-level and family-level.There was no systematic spatial or temporal bias in taxonomic resolution that would unduly bias results.However, "taxonomic richness" calculated in this way is inevitably lower than species richness and thus the spatiotemporal changes identified in the paper could be underestimated relative to reality.Overall, our findings support and extend previous work, including widely-recognised biogeographical relationships such as the LDG and ADG as well as more recent (and contentious) work on temporal increases in alpha richness, using a robust dataset that is extensive both spatially and temporally.The casual mechanisms, however, remain ambiguous and although this is indicative of the inherent complexity of ecological systems, we suggest that more research into the drivers of these patterns is undertaken as a matter of priority.This is especially important given the contexts of widespread biodiversity decline at national and global scales and strengthening extrinsic changes including those to climate, non-native introductions, increasing homogenisation of generalist taxa, landscape simplification and anthropogenic disturbance.

Figure 1 .
Figure 1.Sites included in this study.All palaeoecological data were collated from the European Pollen Database [44]; data are accessible as per the data accessibility statement at the end of this paper.

Figure 1 .
Figure 1.Sites included in this study.All palaeoecological data were collated from the European Pollen Database [44]; data are accessible as per the data accessibility statement at the end of this paper.

Figure 2 .
Figure 2. Schematic of the relationship between taxonomic richness in the overall vegetation community (solid line) and the shrub and herb sub-community (dashed line) at 23 sites in western Eurasia for: (a) latitude; and (b) altitude.The equations were calculated using regression analysis on raw data (n = 867 for overall vegetation community and 866 for shrub and herb sub-community) and values show the change in richness per 1 • increase in latitude and per 100 m increase in altitude.

Figure 3 .
Figure 3. Schematic of the relationship between taxonomic richness in the overall vegetation community (solid line) and the shrub and herb sub-community (dashed line) at 23 sites in western Eurasia in relation to time.The equations were calculated using regression analysis on raw data (n = 867 for overall vegetation community and 866 for shrub and herb sub-community) and values show the change in richness per 100 years.

Figure 3 .
Figure 3. Schematic of the relationship between taxonomic richness in the overall vegetation community (solid line) and the shrub and herb sub-community (dashed line) at 23 sites in western Eurasia in relation to time.The equations were calculated using regression analysis on raw data (n = 867 for overall vegetation community and 866 for shrub and herb sub-community) and values show the change in richness per 100 years.

Figure 4 .
Figure 4. Schematic of the relationship between the rate of within-site turnover (higher = greater) in the overall vegetation community (solid line) and the shrub and herb sub-community (dashed line) at 23 sites in western Eurasia in relation to time.The equations were calculated using regression analysis on raw data (n = 813 in both cases) and values show the change in richness per 100 years.

Figure 4 .
Figure 4. Schematic of the relationship between the rate of within-site turnover (higher = greater) in the overall vegetation community (solid line) and the shrub and herb sub-community (dashed line) at 23 sites in western Eurasia in relation to time.The equations were calculated using regression analysis on raw data (n = 813 in both cases) and values show the change in richness per 100 years.

Figure 5 .
Figure 5. Relationship between taxonomic turnover (higher = greater) between 39 sequential 100year time periods and taxonomic richness in the earlier of those time periods for sites in western Eurasia for: (a) overall vegetation community; and (b) shrub and herb sub-community.In both cases, data points for sequential 100-year periods falling between present day and 1,000 years ago are shaded black, those 1000-1999 years ago are shared dark grey, those 2000-2999 are shaded light grey and those 3000-3999 are open circles.Note that scales differ between graph panels.

Figure 5 .
Figure 5. Relationship between taxonomic turnover (higher = greater) between 39 sequential 100-year time periods and taxonomic richness in the earlier of those time periods for sites in western Eurasia for: (a) overall vegetation community; and (b) shrub and herb sub-community.In both cases, data points for sequential 100-year periods falling between present day and 1000 years ago are shaded black, those 1000-1999 years ago are shared dark grey, those 2000-2999 are shaded light grey and those 3000-3999 are open circles.Note that scales differ between graph panels.

Figure 6 .
Figure 6.Summary of patterns between taxonomic richness and within-site turnover relative to spatiotemporal variables (time, latitude, altitude).Solid lines denote direct effects, dashed lines denote indirect effects where one factor mediates relationships between other variables.LDG = Latitudinal Diversity Gradient; ADG = Altitudinal Diversity Gradient; DSH = Diversity Stability Hypothesis; SAH sub-community = shrub and herb sub-community.

Figure 6 .
Figure 6.Summary of patterns between taxonomic richness and within-site turnover relative to spatiotemporal variables (time, latitude, altitude).Solid lines denote direct effects, dashed lines denote indirect effects where one factor mediates relationships between other variables.LDG = Latitudinal Diversity Gradient; ADG = Altitudinal Diversity Gradient; DSH = Diversity Stability Hypothesis; SAH sub-community = shrub and herb sub-community.

Table 1 .
Generalised Linear Models of taxonomic richness in relation to spatiotemporal variables.Chi square values are likelihood ratio for the overall models and Wald for factors and interaction terms.

Table 2 .
Generalised Linear Models of within-site turnover in relation to time and mediating effects of latitude and altitude.Chi square values are likelihood ratio for overall models and Wald for factors and interactions.