Taxonomic and functional diversity of North American vegetation during the last interglacial–glacial cycle

Abstract

We synthesized pre-last glacial maximum pollen records to reconstruct North American pollen diversity since ca. 130 ka. Using taxonomic diversity (a measure of the number and abundance of taxa) and functional diversity (a measure of the number and abundance of different phenotypes) we identified temporal and spatial diversity trends for six North American bioregions: Arctic, Intermountain West, Mexico, Pacific Northwest, Southeast, and Yucatán. Reconstructed taxonomic temporal and spatial trends vary among bioregions, with regional diversity patterns captured in the functional metric, suggesting shifts in species composition coincide with shifts in ecosystem function. However, significant shifts in taxonomic pollen diversity differed in frequency, magnitude, and timing from their functional counterparts. Variations in both regional taxonomic and functional diversity response to global and regional temperature trends were evident, suggesting temperature alone does not fully explain changes in species composition. Regional richness estimates exhibited higher stability relative to the weighted diversity estimates indicating low levels of species turnover through Late Quaternary warming–cooling phases. Shifts in regional diversity did not predictably respond to stadial and interstadial transitions. Instead, North American patterns of plant diversity over the last ca. 130 ka differ geographically, likely responding to regional rather than global climate change.

Introduction

Studying past climatic variations to understand the processes that determine plant community composition and resilience to large-scale compositional shifts is central to paleoecology (Tausch et al., 1993; Willis et al., 2010; Birks, 2019). North American climate–vegetation relationships from the last glacial maximum (LGM) to the present are well explored (e.g., Bartlein et al., 1998; Williams et al., 2004). Longer timescale relationships, such as the conditions during the last interglacial–glacial cycle, are less known but do provide an opportunity to examine multiple examples of vegetation responses to climate shifts, such as those that occur between Marine Oxygen Isotope stages (MIS).

The last interglacial–glacial cycle (MIS 5–2) contains several warm (interstadial) and cold (stadial) periods useful for examining vegetation responses to rapid, global-climate shifts. Climate driven responses of vegetation have been explored using European pollen records (Guiot et al., 1993; Cheddadi et al., 2005, Lucchi, 2008; Jiménez-Moreno et al., 2010b), and regionally in North America (Glover et al., 2021). Compiling pre-LGM pollen records at continental scales provides both a novel dataset for paleoecological study and for reconstructing vegetation dynamics independent of human influence through multiple climate scenarios between stadial and interstadial periods (Tzedakis et al., 2009). European records suggest MIS 5e was >4.3℃ warmer than present with higher plant community turnover in Central Europe compared to the north and south (Wohlfarth, 2013; Felde et al., 2020; Wilcox et al., 2020). Europe transitioned between warm, forested periods during MIS 5 and cool, open vegetation from MIS 4 through MIS 2 (Helmens, 2014). In North America few syntheses of last interglacial–glacial vegetation changes exist, despite available data (Glover et al., 2021) (Table 1). Around 50 North American pollen records are reported to pre-date the LGM (Williams et al., 2018). Of these, ∼20 have sampling-resolution and age-model constraints that allow for a pre-LGM analysis of continental-scale, temporal and spatial patterns of diversity change. These data are distributed across several North American ecoregions (Fig. 1). Grouping individual sites into broader ecoregions (Bailey, 2014) allows for an ecosystem-centric examination of vegetation dynamics for the past 130 ka.

Figure 1. Map of the 23 sites used in this study. Site locations are shown as circles with the site names. Sites are colored based on the seven regional groupings determined by the DCA (Supplemental Figure 2). The seven DCA regions encompass nine level-I ecoregions. The represented ecoregions are shown in color and listed in the legend. Several of the seven DCA regions represent multiple level ecoregions. Due to the lack of sites in the northeast, the Eastern Temperate Forest is divided into level-II classifications.

Table 1. Table of sites used in this study and their associated publications. Age Range calculated as the difference between the youngest date and the oldest. Region includes the DCA-based grouping made in this study.

Plant communities vary along environmental gradients at different spatial and temporal scales (ter Braak and Prentice, 2004). These often appear as elevational and latitudinal gradients that result from underlying temperature and moisture gradients, particularly in mountainous regions (Gentry, 1988; von Humboldt and Bonpland, 2010; Lamanna et al., 2014; Felde et al., 2016). As environmental gradients shift, the spatial extent of plant communities shift in response (Comes and Kadereit, 1998; Thuiller et al., 2005). Plant communities can be characterized through estimates of taxonomic diversity, incorporating measures of richness (number of taxa) and evenness (distribution of taxa abundances) (Schuler and Gillespie, 2000; Howard and Lee, 2003; Sax and Gaines, 2003). Studies show high plant diversity maintains ecosystem functions and services including high net primary productivity (NPP) (Tilman et al., 1997; Hector et al., 1999; Isbell et al., 2011). High diversity is supported by biotic and abiotic factors that allow for independent regional responses to climate (Zak et al., 2003; Benton, 2009; López-Angulo et al., 2020). Understanding how plant communities respond to shifts in climate in addition to the conditions that support high diversity can help predict the effect of current warming on ecosystems as they adapt (Stohlgren et al., 2000; Suggitt et al., 2019).

Reconstructions of pollen-based, paleodiversity provide scenarios to study shifts in plant community composition (Birks, 2019). Modern pollen-derived taxonomic diversity estimates have effectively captured the structure and composition of vegetation on landscapes at local scales (van der Knaap, 2009; Matthias et al., 2015; Birks et al., 2016; Felde et al., 2016) and, while caution must be taken in interpreting pollen-derived diversity estimates (Weng et al., 2006; Goring et al., 2013; Väli et al., 2022), diversity estimates can address questions of diversity change from the trends produced. By reconstructing plant paleodiversity trends across North America, we can interpret both ecosystem diversity and dynamics as they responded to past climate changes.

Paleodiversity has been reconstructed via several common indices, including species richness, Shannon’s index, and Simpson’s index (Morris et al., 2014). While these methods are still commonly used, other methods such as Hill numbers are more versatile via their ability to vary between different orders (Jost 2006, 2007). Applying diversity indices in a paleo context can be complicated by several taxonomic-scale-related problems (Cleal et al., 2021), such as varying taxonomic resolution (Birks and Birks, 2000) and less information for taxonomic classifications (Godfray et al., 2004). While much of this is simply the nature of working within the fossil record, incorporating plant functional traits can provide a means to reconstruct ecosystem dynamics where taxonomy alone cannot (Eronen et al., 2010; Barnosky et al., 2017; Chacón-Labella et al., 2023).

In this study, we used pollen-derived taxonomic and functional diversity estimates to quantify North American plant community dynamics from the last interglacial period (MIS 5e) to the present (MIS 1). We hypothesize that global temperature variability and latitudinal temperature gradient compression drove shifts in taxonomic and functional diversity. We examine whether continental and regional shifts in pollen-derived taxonomic diversity are mirrored in functional trait-based diversity for six geographic regions of North America based on 23 fossil pollen records.

Regional Setting

Modern North American ecoregion floristics are represented by high species richness in the tropics of southern Mexico, southeastern United States, and Pacific coast (Ricketts et al., 1999). We can use sub-biome ecoregions (Bailey, 2014) to characterize subcontinental scale regions with similar floristics to scale up from individual sites to broader, floristic heterogeneity at the continental scale. Richness is lowest in the sagebrush steppe and grasslands of the continental interior and Arctic tundra (Cavender-Bares et al., 2020). While plant and pollen diversity are not directly analogous, regional correlations between the two have been established in northern Europe (Reitalu et al., 2019).

To verify whether North American ecoregions reasonably correspond to the pollen record, we grouped 23 North American pollen records, dated to at least 40 cal ka BP, into six regions based on a detrended correspondence analysis (DCA) of pollen type abundances (Supplemental Figs. 1 and 2). These DCA-based groupings align with level-I ecoregions of North and Central America with 9 of 15 represented (Fig. 1) (CEC, 1997; Griffith et al., 1998). DCA clusters grouped similar pollen records based on floristic history rather than modern composition while still being recognizable within modern North American floristics. DCA clusters aligned with ecoregion determinations (CEC, 1997; Griffith et al., 1998), which allowed us to use ecoregion changes from present as an interpretive framework (Fig. 1).

Methods

The overall workflow is outlined in Fig. 2 with each step pertaining to one of three categories: data compilation, in which the datasets used were established, data preparation, in which data were transformed and filtered, and data analysis, in which data were used to generate the results presented.

Figure 2. Flowchart of the workflow to generate taxonomic and functional diversity trends for the six regions in this study using publicly available datasets. Steps are colored based on three categories: data collection (gray), where data for the study were accessed; data preparation (orange), where data were filtered and transformed; and data analysis (blue), where the data presented in this study were created.

Records selection

Published pollen data within North America with dated materials to at least 25 cal ka BP and extrapolated to at least 40 cal ka BP were compiled from the Neotoma Paleoecology Database (Williams et al., 2018). Of the 28 records, three were removed for low resolution (fewer than 10 pollen samples or with data spanning less than 10 ka) and six were removed based on constraints in extrapolated age-depth models. Accepted sites also had extrapolations that followed reasonable sedimentation rates based on their 20–30 ka dated intervals. Four newer records not archived in Neotoma at the time of our study were included for a total of 23 records (Table 1). Raw pollen counts were available for all records used in this study.

Data conversions

Pollen taxa identified in each record were filtered based on the level of specificity of the identification. Only terrestrial pollen data were used, excluding spores, aquatic pollen types, and non-pollen palynomorphs to minimize site specific variability (e.g., wetland changes). Pollen taxa with less than 1% abundance in an individual record were excluded. Filtering identified pollen types reduced overall pollen-based, taxonomic richness, but aided in reducing researcher specific biases in pollen identification. The resulting data set contained 142 total pollen types, with 39 family-, 96 genus-, and 7 species-level pollen identifications (Supplemental Table 1).

We interpolated each record to address issues of radiometric dating. Raw pollen counts from each record were converted into relative proportions of the total terrestrial pollen and interpolated using a moving window weighted means interpolation to 1000-year time steps between 0 and 150 ka. Weighted mean interpolations of pollen data were ignored if temporal gaps within the original datasets exceeded ± 5 ka of the centroid of the targeted, 1000-year time step, resulting in discontinuities in our records. Samples within the ± 5 ka window were weighted using a normal curve based on their proximity to the centroid.

DCA and region determination

Similarities of the ecological and taxonomic clustering of the 23 pollen records were verified using DCA, an ordination method appropriate for grouping proportional datasets distributed along broad environmental gradients (Holland, 2008; Oksanen et al., 2013). DCA was conducted on all samples from all records ensuring that complete histories were considered simultaneously. Other ordination approaches for determining pollen–environment relationships (e.g., canonical correspondence analysis [CCA]) did not improve or enhance our analysis.

Hill numbers

Hill numbers are a family of interrelated diversity estimates that vary by the weights of rare and common taxa and are expressed as ‘effective number of species’ (Chao et al., 2014a, b). Hill numbers have proven to be a versatile approach to representing multiple diversity estimates simultaneously (Ohlmann et al., 2019). While commonly used to represent taxonomic diversity, a framework exists for calculating functional diversity with Hill numbers, although it has yet to be adopted in a paleoecological context (Chiu and Chao, 2014). With taxonomic Hill numbers, the ‘species’ that diversity is calculated from is the relative proportion of the taxa. In functional Hill numbers, the ‘species’ is the functional distance between two taxa rather than the proportions (Chao et al., 2014a). Simply, the outputs of taxonomic Hill numbers are in units of ‘effective number of species,’ indicating a value of x is equivalent to x number of equally abundant species within a community. Functional diversity is interpreted as x number of equally abundant and functionally distinct species.

Varying both taxonomic and functional Hill numbers by the variable q adjusts the weights between rare and common taxa or functional distances in an assemblage. Lower q values produce higher diversity estimates through the inclusion of both rare and common taxa/distances while higher q values produce lower estimates by filtering out rare taxa/distances. We calculated three taxonomic and three functional Hill numbers for each sample within a site to model diversity over time within the sites. The taxonomic Hill number of order q (^qD) follow formula Equation 1a (Chao et al., 2014b):

(1a)\begin{equation}{}^{q}D = {\left( {\mathop \sum \limits_{i = 1}^S {p_i}^q} \right)^{\frac{1}{{1 - q}}}}\end{equation}

Where S equals the total number of species in the assemblage, p_i equals the proportion of species i. To represent functional diversity, Chiu and Chao (2014) modified taxonomic Hill numbers (^qD) to include the mean functional distance between any two randomly sampled individuals in the assemblage using Rao’s quadratic entropy (Q) as defined by Equation 2 (Rao, 1982):

(2)\begin{equation}Q = \mathop \sum \limits_{i = 1}^S \mathop \sum \limits_{j = 1}^S {d_{ij}}{p_i}{p_j}\end{equation}

Where S and p_i remain unchanged from Equation 1a, p_j equals the proportions of species j, and d_ij equals the functional distance between species i and j, calculated using Gower distances. Distance measures are calculated between all species but vary in precision based on the number of trait observations. The resulting functional Hill numbers (^qD(Q)) are described by Equation 3a (Chiu and Chao, 2014):

(3a)\begin{equation}{}^{q}D\left( Q \right) = {\left[ {\mathop \sum \limits_{i = 1}^S \mathop \sum \limits_{j = 1}^S \frac{{{d_{ij}}}}{Q}{{\left( {{p_i}{p_j}} \right)}^q}} \right]^{\frac{1}{{2\left( {1 - q} \right)}}}}\end{equation}

Where all variables are the same as defined in Equations 1a and 2. When q = 1 in either metric, the equations are undefined and Equations 1b and 3b, which are equivalent to the limit as q approaches 1, must be used for taxonomic and functional diversity, respectively:

(1b)\begin{equation}{}^{1}D = \,\mathop {\lim }\limits_{q \to 1} {}^{q}D = {\text{exp}}\left( { - \mathop \sum \limits_{i = 1}^S {p_i}\log {p_i}} \right)\end{equation}

(3b)\begin{equation}{}^{1}D\left( Q \right) = \mathop {\lim }\limits_{q \to 1} {}^{q}D\left( Q \right) = {\text{exp}}\left[ { - \frac{1}{2}\mathop \sum \limits_{i = 1}^S \mathop \sum \limits_{j = 1}^S \frac{{{d_{ij}}}}{Q}{p_i}{p_j}\log \left( {{p_i}{p_j}} \right)} \right]\end{equation}

The first three orders of taxonomic diversity collapse down to three common indices: Richness (⁰D), Shannon’s index (¹D), and Simpson’s index (²D). These indices weight all species equally, all species by their relative abundance, and all species by the square of their relative abundance, respectively. For functional diversity, ⁰D(Q) and ²D(Q) correspond to functional attribute diversity (FAD) and Gini–Simpson index. ¹D(Q) does not correspond to an existing index, but for naming consistency we refer to it as functional Shannon’s index, although it is important to reiterate the two are not directly equivalent (Chiu and Chao, 2014). ⁰D(Q) can be thought of as the sum of all functional distances, ¹D(Q) as the sum of all functional distances weighted by their relative abundance, and ²D(Q) as the sum of all functional distances weighted by the square of their relative abundance.

Seven traits were selected for functional diversity analysis: leaf mass per area, leaf lifespan, max plant height, max plant longevity, seed mass, whole plant growth form, and whole plant woodiness. We selected these traits due to their abundant number of observations and ability to represent most plant ecological strategies through the leaf-height-seed scheme (Westoby, 1998). Some traits are auto-correlated (e.g., leaf lifespan and leaf mass per area), affecting d_ij and Q, which in turn alter the magnitude of functional diversity values, but not the overall pattern, which is driven by p_i and p_j.

Trait data were compiled following Brussel and Brewer (2021), using the Botanical Information and Ecology Network (BIEN) database and filtered to include only North and Central America taxa (Maitner et al., 2018). For species-level pollen-types, trait values for all occurrences were averaged and reported. For genus- or family-level pollen types, we averaged the values of all species in those groupings and used the averages as the trait value for the taxon. Families and genera containing species with multiple growth forms (herbaceous and woody taxa) were treated as a different category from woody or herbaceous taxa when calculating functional diversity. To quantify the uncertainty associated with using trait averages, we performed a Monte Carlo simulation. Trait values were resampled to create simulated distributions of 5000 values from which we calculated the mean. This was repeated 10,000 times, creating a distribution of 10,000 means per trait, from which we sampled 100 means and calculated the associated functional diversity estimates.

Distances are calculated from the trait values of each species compared pairwise with all other species in the assemblage. This method is effective in cases when all taxa within an assemblage are identified to the species level and there is little variation in trait values (Pérez-Harguindeguy et al., 2016). Complications arise when taxa within an assemblage are identified at genus or family levels and there is higher variation in trait values. Averaged trait values ignore the natural variability that exists within a genus or family (Roscher et al., 2018; Sotomayor et al., 2022). However, our goal is to use traits to group pollen (plants) independently of their underlying taxonomy rather than use traits for their functional interpretations. The resulting functional diversity estimate loses the ability to draw specific environmental interpretations but provides comparator diversity datasets to analyze along with taxonomic diversity data.

Temporal analysis

Regional diversity trends were interpolated using a generalized additive model (GAM) to smooth the individual trends of each site within a region. The GAM model was calculated using the ‘gam’ function within the mgcv package in R Studio (Wood, 2001). Age was used as the predictor and diversity as the response with a Gaussian family, a normal/identity link function, and a restricted maximum likelihood (REML) smoothing method. Records within each region vary in length, affecting the GAM smoother and decreasing individual site influence towards the present. We selected a GAM to allow for the identification of significant shifts in the diversity time series based on the methods outlined in Simpson (2018). Times of significant change were highlighted on the splines. Site diversity values were normalized around the mean to enable comparisons of shift magnitude between regions and to highlight differences between the temporal patterns of the taxonomic and functional metrics. Without normalization, questions of the validity of diversity estimate values could prevent thorough interrogations of the trends being presented. Raw diversity trends are shown in Fig 6 and Supplemental Figs 3, 4, and 5. GAM splines are plotted alongside regional temperature records (δ¹⁸O) and atmospheric CO₂.

Spatial analysis

For each region, we calculated the average of non-transformed taxonomic and functional diversity of each MIS. We then interpolated average diversity using empirical Bayesian kriging, chosen for its ability to quantify standard error and assumption of spatial correlation between inputs (Krivoruchko, 2012). Resulting maps use a bivariate color scale in which high values in both taxonomic and functional diversity are displayed in blue and high values in either taxonomic or functional diversity are displayed in black or orange, respectively. Results for each region are mapped to indicate the distribution of and change in diversity between MISs.

Results

Regionalization of records

Pollen records were grouped into six North American ecoregions: the Arctic (4 sites), Pacific Northwest (3 sites), Intermountain West (10 sites), the Southeast (2 sites), Mexico (3 sites), and Yucatán (1 site). The Yucatán site was an outlier in the DCA (Supplemental Figs. 1 and 2). The floristic composition of the Yucatán record is primarily tropical taxa, distinguishing it from all other records.

Temporal diversity trends

Monte Carlo simulations show minimal variation in the magnitude and no variation in the trend of functional diversity, indicating the estimates do not vary significantly with uncertainty in trait averages. We found trends for richness (⁰D) and FAD (⁰D(Q)) to be distinct from the trends for Shannon’s (¹D) and functional Shannon’s indices (¹D(Q)) and Simpson’s (²D) and Gini–Simpson’s indices (²D(Q)), while the trends of the latter two are similar. Richness and FAD are distinct from the other two orders because of its exclusion of abundance (Fig. 3). In most cases, trends for Simpson’s and Gini Simpson’s indices have steeper increases and decreases, appearing as Shannon’s and functional Shannon’s indices trends with higher amplitude (Figs. 4 and 5). In all three orders of q, diversity trends differ by region. Although continental-scale patterns of diversity change are not apparent, regional diversity show examples of long-term trends and the influence of stadial and interstadial periods.

Figure 3. Normalized regional plant diversity of order q = 0 time series with climate. (A) Oxygen isotope record for the Arctic (Cronin et al., 2019) with arctic (ARCT) regional plant diversity. (B) Oxygen isotope record the Eastern Pacific Ocean (Herbert et al., 2001) with Pacific Northwest (NW), Intermountain West (IMW), and Mexico (MEX) regional plant diversity. (C) Oxygen isotope record for the Caribbean (Schmidt et al., 2004) with Southeast (SE) and Yucatán (YUC) regional plant diversity. Global atmospheric CO₂ concentration from the EPICA Dome ice core (Bereiter et al., 2015) is plotted in all climate plots in blue. Taxonomic richness (⁰D) is shown in black. Functional attribute diversity (FAD; ⁰D(Q)) is shown in orange. Individual site trends are plotted as thin lines. Gray vertical bars represent stadial (cool) MIS and white bars represent interstadial (warm) MIS (Lisiecki and Raymo, 2005). The dotted line at y = 0 indicates the mean of each site following normalization. Bold lines represent GAM splines displaying the interpolated diversity of order q = 0 trends for each region with the 95% confidence interval shaded. Thickened segments of the spline indicate sections of significant slope.

Figure 4. Normalized regional plant diversity of order q = 1 time series with climate following the format of Fig 3. Shannon’s index (¹D) is shown in black. Functional Shannon’s index (¹D(Q)) is shown in orange. ARCT = Arctic; NW = Pacific Northwest; IMW = Intermountain West; MEX = Mexico; SE = Southeast; YUC = Yucatán.

Figure 5. Normalized regional plant diversity of order q = 2 time series with climate following the format of Fig 3. Simpson’s index (²D) is shown in black. Gini–Simpson’s index (²D(Q)) is shown in orange. ARCT = Arctic; NW = Pacific Northwest; IMW = Intermountain West; MEX = Mexico; SE = Southeast; YUC = Yucatán.

Richness (⁰D) and FAD (⁰D(Q)) trends follow each other well with significant shifts in one commonly reflected in the other. Both metrics remains stable around the mean in all regions except Mexico. Richness trends in Mexico show relative increases, with below average richness at the beginning of the record and average richness by the end. All other regions fluctuate around the mean but do not exhibit long-term shifts in magnitude. Richness and FAD vary inconsistently from global temperature as represented by atmospheric CO₂ with regions responding independently and within region responses being unpredictable during stadial and interstadial conditions. Regional climate records align more closely with richness and FAD however there remains are large degree of variance in response to shifts in regional climate.

Shannon’s (¹D) and functional Shannon’s indices (¹D(Q)) trends are similar showing less stability around the mean than richness and FAD. Multiple regions fluctuate around the mean (e.g., Arctic, Pacific Northwest, Intermountain West, and Mexico), while others decrease (e.g., Southeast) or experience large shifts in diversity (e.g., Yucatán). Shannon’s and functional Shannon’s indices follow similar patterns to richness and FAD when compared with the regional climate records and atmospheric CO₂. However, shifts in regional climate correspond more consistently with shifts in the Shannon’s and functional Shannon’s indices than with richness and FAD.

Simpson’s (²D) and Gini–Simpson’s indices (²D(Q)) trends again follow each other at a similar level to Shannon’s and functional Shannon’s indices. However, the two metrics mismatch to a higher degree in Mexico and the Arctic. Within these two regions, Gini–Simpson’s index (²D(Q)) still has the same general pattern as Simpson’s index (²D), but there are differences in magnitude and timing of the shifts within the patterns. Simpson’s and Gini–Simpson’s indices show nearly identical patterns of diversity shifts as Shannon’s and functional Shannon’s indices, differing only in magnitude, mirroring the atmospheric CO₂ and regional climate patterns previously described.

Geographic diversity distribution

Taxonomic (^qD) and functional diversity (^qD(Q)) align in terms of the geographic distribution of high and low regional diversity (Fig. 6). However, anticipated continental gradients from low- to high-diversity regions are not consistent between the two metrics (e.g., Lamanna et al., 2014). Variation in the spatial extent of low taxonomic or functional diversity was not consistent continent-wide. Overall, eastern paleoecological records indicate higher taxonomic pollen diversity compared to western North America. North–south patterns show the Yucatán Peninsula exhibiting the highest diversity in all three orders of q. Prior to MIS 2, the Pacific Northwest shows similarly high diversity to the Yucatán followed by the Southeast. These patterns of taxonomic pollen diversity are generally captured by the functional metrics with variance in the magnitude of functional diversity relative to taxonomic diversity.

Figure 6. Bivariate maps of the geographic distribution of plant diversity in North America. Taxonomic diversity is indicated by the gray shading while functional diversity is orange shading. High values are indicated by dark shading and low values by light shading. Blue shading indicates overlap between taxonomic and functional diversity. Site locations are indicated by points colored by interpolation input values. Mismatches between site color and interpolated diversity are due to a combination of interpolation, including diversity of adjacent regions and the rounding of values. Ecoregions lacking data are overlain in white. Color schemes for both taxonomic and functional diversity are geometrically scaled (color intervals are based off a geometric series). Taxonomic richness (⁰D) is scaled from |minimum value – 17.9|17.9 – 24.6|24.6 – maximum value|. Functional attribute diversity (⁰D(Q)) is scaled from |minimum value – 25.2|25.2 – 34.0|34.0 – maximum value|. Shannon’s index (¹D) is scaled from |minimum value – 4.8|4.8 – 6.6|6.6 – maximum value|. Functional Shannon’s index (¹D(Q)) is scaled from |minimum value – 7.3|7.3 – 9.2|9.2 – maximum value|. Simpson’s index (²D) is scaled from |minimum value – 3.3|3.3 – 4.7|4.7 – maximum value|. Gini–Simpson’s index (²D(Q)) is scaled from |minimum value – 5.0|5.0 – 6.0|6.0 – maximum value|. The extent of the Laurentide Ice Sheet is shown in light blue (Batchelor et al., 2019). Continent coastlines based on the reconstructed sea-level trends from Cutler et al. (2003).

Discussion

North American plant taxonomic and functional diversity underwent several asynchronous changes in the last ca. 130 ka. We reject our hypothesis that global temperature variability and latitudinal temperature gradient compression drove shifts in taxonomic and functional diversity. Other global mechanisms, such as CO₂ variability, also do not appear to be drivers (Figs. 3–5). Changes in moisture flux into North America could explain regional diversity changes through the last interglacial–glacial cycle, but independent region-to-continental precipitation reconstructions are not available for comparison for this time period. Punctuated patterns of cooling and warming over the last 130 ka did not consistently correspond to changes in taxonomic or functional diversity. When significant diversity changes occurred, we note misalignments in the temporal range between taxonomic and functional trends. With few exceptions, the trends for taxonomic pollen and functional diversity are similar, suggesting that these diversity calculations may be interchangeable.

Richness versus diversity

The trends for richness (⁰D) and FAD (⁰D(Q)) in each region are generally more stable than the trends for both Shannon’s (¹D) and functional Shannon’s indices (¹D(Q)) and Simpson’s (²D) and Gini–Simpson’s indices (²D(Q)). Many regions experience periods of significant change in richness and FAD, but significant changes are less frequent and smaller compared to diversity orders of q = 1 and q = 2, likely reflecting low regional species turnover over time. The pattern of significant change is consistent with species introductions and extirpations that require larger environmental shifts to change plant composition (Dakos et al., 2019). The assumption of an individualistic response to climate in paleoecological records is critical and implicitly includes other abiotic and biotic interactions, such as soil or competition conditions (Bartlein et al., 1986; Davis and Shaw, 2001). While difficult to quantify, it is possible that changes in richness and FAD require larger magnitude environment changes than are required to shift diversity orders of q = 1 and q = 2, particularly on millennial time scales if extirpation and immigration rates are equal. This is specifically true with pollen data where ‘extirpation’ of a pollen type must overcome local presence and potential residual long-distance transport from extra-local populations.

Diversity patterns of the entire community (order of q = 1) and the dominant taxa (order of q = 2) are similar, indicating diversity shifts are occurring evenly throughout different floristic communities. Compared to richness (⁰D) and FAD (⁰D(Q)), Shannon’s (¹D), functional Shannon’s (¹D(Q)), Simpson’s (²D), and Gini–Simpson’s (²D(Q)) indices all experience higher degrees of diversity fluctuations (Figs. 3–5). These fluctuations indicate regional community composition varies while large species turnover is uncommon as indicated by the richness stability. Shifts in richness indicate extirpation and immigration as species ranges change. Our results show that the composition of floristic communities across the continent have been pliable in response to regional climate change.

Taxonomic versus functional diversity

Taxonomic and functional diversity trends are similar over the last ca. 130 ka, indicating paired analysis is not necessary if taxonomic or functional data are available. The two metrics differ in two key ways: when significant shifts occur in the time-series and the amplitude of trends. The use of Hill numbers allows us to be confident that similarities and differences between the metrics are a result of data differences and not an artifact of the selected indices. Differences between the two metrics could be due to shifts in the dominant plant functional types without corresponding shifts to the actual taxa present in these ecosystems (Doxa et al., 2020). For example, it’s possible that with expansion of the Laurentide Ice Sheet, the southward displacement of conifers shifted the functional diversity in many regions that were previously dominant in broad-leaved trees (Williams et al., 2004).

Both taxonomic (^qD) and functional (^qD(Q)) metrics follow similar temporal and spatial patterns (Figs. 3–6), with both metrics capturing similar patterns of plant diversity across the continent. The similar overall patterns point to taxonomic and functional diversity being inherently linked. Ecosystems are high or low in both metrics with mismatches uncommon (Fig. 6). Similarities between the two metrics show reconstructions of high taxonomic diversity can be considered as productive ecosystems without the need for functional diversity reconstructions. Some communities may experience a large turnover in the species pool with a limited effect on the functional diversity (e.g., Arctic q = 2 – MIS 3) (Villéger et al., 2010; Ibsen et al., 2020). In the context of climate change, shifts in ecosystem services as characterized by functional diversity are linked to shifts in taxonomic diversity but triggered differently. Functional diversity experienced fewer significant diversity shifts than taxonomic diversity, highlighting the possibility of increased resilience to environmental change. Functional diversity shifts occur both prior to and following species composition shifts, suggesting uncertainty in predicting responses to future climate change. As with taxonomic diversity, regional history plays the most important role in determining the degree to which functional diversity, and by extension ecosystem services, may be affected by climate change.

Regional climate

We expected diversity to increase with distance from the continental ice sheets based on compression of latitudinal temperature gradients. Changes in expected diversity should reflect admixtures of regional floras as they responded to shifts in temperature based on individualistic response of taxa to climate (Bartlein et al., 1986; Minckley et al., 2019). However, our results show regional responses to global temperature changes are difficult to predict due to mismatches between regional and global climate variability. Significant shifts happened since the last interglacial, but not proximal to MIS boundaries (Figs. 3–5). Age uncertainties and sampling resolution of older pollen records may create enough error in these records to blur direct mechanistic links between MIS transitions and regional vegetation changes. Previous studies have shown vegetation composition to reliably track climate at millennial scales (Harrison and Sanchez Goñi, 2010; Jiménez-Moreno et al., 2010a), but our results suggest major floral transitions represented by significant changes to pollen diversity lack a direct association with global climate shifts. This discrepancy is likely due to our incorporation of diversity metrics (^qD and ^qD(Q)) rather than focusing on dominant plant taxa. However, describing how correlated the two are to one another will be important for future studies. Regionally, changes in diversity are not uniform in magnitude or trend. In many cases, shifts observed in one region are not paralleled in the trends of another, even when geographically adjacent.

Understanding the role regional climate plays in influencing plant diversity has important implications for how we predict future diversity trends. The regional variability of diversity shifts is one of the difficulties revealed by our study. Plant diversity metrics (^qD and ^qD(Q)) for the Intermountain West indicate a high degree of stability over the last ca. 130 ka. Individual sites within the region experienced transitions between more open sagebrush steppe and closed spruce–fir forests in line with changes from warm–dry to cool–wet climates (Beiswenger, 1991; Herring and Gavin, 2015). Despite changes to its’ floral composition, the overall diversity of the Intermountain West remains relatively low compared to the rest of the continent, indicating that observed compositional shifts are relatively minor compared to more palynologically diverse regions (Fig. 6).

In contrast, the Yucatán, a diversity hotspot for North America, fluctuates more in trend with orders of q = 1 and q = 2 (Figs. 3–5). While the lower diversity of the Intermountain West is partially a result of the biases of the pollen record with family-level pollen types (e.g., Poaceae, Asteraceae, Amaranthaceae) more common in the West than genus-level pollen types (e.g., Betula, Tilia, Acer, etc.) identified in eastern North American pollen records (Qian, 1998; Pinto-Ledezma et al., 2018; Cavender-Bares et al., 2020). With the flora of the Intermountain West adapted to grow in both hot and freezing arid conditions, changes to the temperature of the environment may not have been severe enough to shift the overall communities, but rather altered the latitudes and elevations that different communities could grow (Breshears et al., 2008; Lavin et al., 2013).

Diversity may tend to decrease into stadial periods and increase into interstadials, but this pattern is not universal (e.g., Mexico q = 1 – MIS 4–2). With the range of environmental conditions contained within the seven North American regions, it is difficult to determine external forcings outside of increased uniformity in hemispheric climate, possibly because of maritime effects. As with global temperature and atmospheric CO₂, regional temperature (δ¹⁸O) alone does not explain all observed diversity shifts. Many instances exist where diversity increases with increasing temperature (e.g., Southeast q = 0 – MIS 2–1, Yucatán q = 1 – MIS 2–1, Mexico q = 0 – MIS 4–3) and decrease with decreasing temperature (e.g., Pacific Northwest q = 0 – MIS 3–2, Yucatán q = 1 – MIS 5a–4, Mexico q = 0 – MIS 3–2). Mismatches between both conditions occur, with scenarios of increasing diversity with decreasing temperature (e.g., Arctic q = 0 – MIS 3, Southeast q = 1 – MIS 2, Pacific Northwest q = 0 – MIS 4) and decreasing diversity with increasing temperature (e.g., Arctic q = 0 – MIS 2–1, Southeast q = 1 – MIS 2–1, Intermountain West q = 0 – MIS 4–3). While cooler temperatures do not correlate to lower diversity and vice versa (Scheiner and Rey-Benayas, 1994; Donoghue, 2008), the lack of clear patterns determining whether a decrease in temperature is associated with a diversity increase, decrease, or no change suggests other climatic factors play a major role. The most likely factor is moisture availability as it is an important environmental gradient for global plant niches (Silvertown, 2004). However, highly resolved paleohydrology records older than the last glacial period are not readily available, making explorations of this driver of past diversity for North America not currently possible. Regional temperature records align with plant diversity shifts to a better degree than global temperature, but again, without independent moisture data, continental diversity patterns cannot be fully explained. Increased research focusing on isotopic biomarkers (Seki et al., 2011; Feakins et al., 2016; Strobel et al., 2022) could result in the development of regionalized moisture records and complete the picture of regional climates during the last interglacial–glacial cycle.

Conclusions

The reconstructions of geographic and temporal shifts in plant diversity outline the continental dynamics of vegetation over the last ∼130 ka. We examined changes in floristic diversity of the past 130 ka (MIS 5e–MIS 1) using 23 pollen records from North America via the implementation of Hill numbers, allowing for a simultaneous examination of species richness and FAD based on presence-absence (orders of q = 0), rarity (orders of q = 1), or dominance (orders of q = 2). Both Shannon’s (¹D) and functional Shannon’s indices (¹D(Q)) and Simpson’s (²D) and Gini–Simpson’s indices (²D(Q)) fluctuate to a higher degree than richness (⁰D) and FAD (⁰D(Q)), indicating that the compositions within ecosystems were highly variable across the record, while species pools remained relatively constant.

The incorporation of functional diversity (^qD(Q)) alongside taxonomic diversity (^qD) supported our interpretations of diversity dynamics during the last interglacial–glacial cycle and shows the connectivity between the two metrics. However, desynchronies between taxonomic and functional diversity shifts point to differences in resilience between species composition and ecosystem function. Continued mechanistic investigations into individual region shifts in both taxonomic and functional plant diversity can identify relationships between global climate and regional vegetation to further understand Quaternary plant dynamics.

Our results highlight the relationship between global temperature, regional climate, ecosystem function, and species composition. While shifts in global temperature may correspond with shifts in the regional climate, how plant diversity and the broader ecosystem respond is unpredictable. However, the integration of both taxonomic and functional diversity into our study emphasizes the inseparability between the two. Scenarios exist in which one metric experiences a shift prior to or without the other, yet both maintain similar temporal and spatial patterns to one another throughout the cooling and warming cycles of the last 130 ka.