The age-depth model (Figure 2) predicts an age of 14,025 yr BP (13,390–14,190 yr BP) for the base of the 4.7 m core. Modeled sedimentation rates are ∼0.03 cm/yr for much of the core, but exceed 0.15 cm/yr in the upper 33 cm and between 332 and 357 cm. Mean deposition time is 30 yr/cm and ranges between 2.2 and 51.5 year/cm. The core from Anderson Fen begins in gray clay with a minor sand fraction and 3–6% LOI (Supplementary Figure 2). Basal clays are sharply overlain by organic lake sediments (454–411 cm) containing Potamogeton and Nuphar polysepala seeds with LOI that increases from 20 to 50%. The remainder of the core consists of various types of peat with LOI generally above 90% (Supplementary Figure 2).

In total, 55 different pollen and spore taxa were identified in the Anderson Fen record. Only the most frequent taxa are shown in Figure 3. The mean temporal resolution of the pollen record is 126 cal yr.

Pollen Zone 1A (14,025–13,200 yr BP) corresponds with clay at the base of the core. This pollen zone is dominated by Salix and Pinus contorta, with each accounting for 30–50% of the total sum (Figure 3). Alnus viridis reaches up to ∼10% and low levels (< 5%) of Abies and Shepherdia canadensis also occur. Pollen from herbaceous plants, mostly Cyperaceae and Artemisia, account for up to 10%. Pinus contorta accounts for 80% at the beginning of Zone 1B (13,200–11,360 yr BP), before decreasing to ∼30% as A. viridis increases to 40–60%. Other conifers in Zone 1B include Abies (10–20%), Tsuga mertensiana (up to 10%) and Picea sitchensis (up to 5%). Zone 1C is a 750-year interval (11,360–10,600 yr BP) of high P. contorta percentages (50–75%), with a corresponding decrease in most other taxa, except Pteridium aquilinum, which reaches up to ∼20%. Pinus concentrations are similar to Zone 1B (Supplementary Figure 3), which suggests that the high pine percentages of Zone 1C are a statistical artifact rather than an indication of increasing pine populations.

There is a large turnover of taxa at the beginning of Zone 2 (10,600–8,230 yr BP), with Pseudotsuga menziesii and Tsuga heterophylla accounting for up to 60 and 40%, respectively (Figure 3). Alnus rubra ranges between 10 and 20% and P. aquilinum accounts for up to ∼20%. Increases in Ericaceae, Cyperaceae and Sphagnum at ∼10,400 yr BP correspond with terrestrialization of the lake basin to a Sphagnum-dominated peatland (Lacourse et al., 2019a).

Zone 3A (8,230–5,400 yr BP) is characterized by the decreasing relative abundance of P. menziesii (30–5%) and increasing abundances of T. heterophylla (25–50%) and Abies (10–30%). A large increase in T. heterophylla pollen concentrations suggests it was more abundant on the landscape after 8,000 yr BP (Supplementary Figure 3). Pinus contorta, A. rubra and Pteridium account for up to 10%. Zone 3B (< 5,400 yr BP) is dominated by T. heterophylla (∼40–60%) with 10–20% Abies, ∼10% A. rubra and ∼5% T. mertensiana. Pinus contorta shows a notable increase to ∼20% over the last 2,200 cal yr coincident with minor increases in Pinus monticola and Picea sitchensis. Alnus rubra increases to 56% of the total pollen sum over the last 100 years.

Axes 1 and 2 of the PCA explain 47.5 and 19.9% of the total variation, respectively, that is present in the pollen record from Anderson Fen. PC axis 1 separates taxa abundant in the late Pleistocene such as Pinus contorta, Alnus viridis, Salix and Athyrium, which have positive loadings, from most other taxa, including most notably Tsuga heterophylla (Figure 4). PC axis 2 effectively separates Holocene assemblages: Pseudotsuga menziesii and Pteridium aquilinum, which are most abundant during the early Holocene, have positive loadings, and mid- and late Holocene Tsuga heterophylla-dominated assemblages generally have negative loadings. The ordination of pollen samples (Figure 4) corresponds closely with pollen zones identified independently through optimal splitting (Figure 3). In general, late Pleistocene and late Holocene samples are more tightly clustered, which reflects higher compositional similarity of pollen assemblages at these times relative to the early Holocene, when assemblages are more diverse.

Redundancy analysis revealed that the paleoenvironmental variables explain a significant amount of variation (43.6%, p < 0.0001) in the pollen data. Temperature and summer insolation are correlated strongly with axes 1 and 2, whereas precipitation shows a stronger correlation with axis 2 (Figure 4). Late Pleistocene samples dominated by P. contorta, A. viridis and Salix are associated primarily with low temperatures. Early Holocene Pseudotsuga-dominated assemblages are associated with relatively high temperature and summer insolation, and relatively low precipitation. Late Holocene samples dominated by T. heterophylla and Abies are associated primarily with high precipitation. These relationships suggest that PC axis 1 represents a gradient in temperature, whereas PC axis 2 represents variation in precipitation and summer insolation, which are negatively correlated (r = −0.62, p < 0.01).

The multivariate time series plots reveal clear differences in temporal changes in woody plant taxa along the elevational transect of sites (Figure 5 and Supplementary Figure 4). At low and middle elevations, forest composition changes from being dominated by P. contorta in the late Pleistocene, to P. menziesii in the early Holocene, and to T. heterophylla and A. rubra in the mid-late Holocene. At high elevations, plant communities are dominated by P. contorta, Abies and Alnus viridis before ∼8,000 yr BP, followed by T. mertensiana and Abies in the mid-late Holocene. The high-elevation record from Burman Pond provides macrofossil evidence that the Abies time series is composed of both A. lasiocarpa and A. amabilis (Figure 5E; Mazzucchi, 2010). There is also a noticeable, though subtle increase in P. contorta across all elevations over the last 2,000 years (Figure 5A). The boxplots associated with each time series shows that some taxa have been largely confined to low and middle elevations (e.g., P. menziesii, A. rubra) or to high elevations (T. mertensiana) for the duration of the records (Figure 5 and Supplementary Figure 4). Other taxa such as P. contorta, T. heterophylla and Abies have occurred at a wider range of elevations.

RLQ and fourth-corner analyses were used to assess the roles of paleoenvironmental change and species traits in directing changes in forest composition as recorded in the pollen percentages. RLQ axes 1 (94.2) and 2 (5.5%) explained almost all the covariance (i.e., cross matrix) between the paleoenvironmental variables and species traits (Figure 6). Permutation testing indicates a significant relationship between paleoenvironmental variables and species traits (p < 0.0001), although the correlation with the pollen matrix for both axes is low (0.204 and 0.122). The variance of the paleoenvironmental scores is captured effectively on the first two RLQ axes (99%), as is the variance in species traits (79%). Paleotemperature has the largest contribution to total inertia (43%), followed by insolation (34%), which provides a proxy for seasonality in temperature and therefore growing season length. Among species traits, shade tolerance, maturity (i.e., minimum seed-bearing age), relative growth rate and waterlogging tolerance contribute most to total inertia (14, 14, 12, and 12%, respectively). We focus on RLQ axis 1 as it explains 94.2% of the co-variance of the three tables. RLQ axis 1 is correlated most strongly with paleotemperature (Figure 6A) and traits such as relative growth rate, waterlogging tolerance and frost tolerance in the positive direction and traits such as shade tolerance, maturity and height in the negative direction (Figure 6B). This axis separates “cool” late-glacial pollen assemblages with abundant P. contorta and A. viridis with positive loadings from “warm” late Holocene assemblages with abundant T. heterophylla and T. mertensiana with negative loadings (Figures 6C,D). Early Holocene assemblages with abundant P. menziesii plot along RLQ axis 2 in association with warm temperatures, high seasonality (represented by insolation), and low precipitation, and traits such as adaptation to fire and large seeds.

Fourth-corner analysis identified significant correlations between almost all pair-wise comparisons of paleoenvironmental variables and species traits (Table 2). However, all environment-trait correlations were weak i.e., r < 0.2. These weak but significant correlations highlight that interspecific differences in a suite of traits are important in regulating long-term forest community assembly. The strongest positive correlations occurred between paleotemperature and typical adult height, shade tolerance, and lifespan. The strongest negative correlations occurred between temperature and relative growth rate and frost tolerance, and between summer insolation and maturity. Temperature was the only environmental variable with a significant correlation to seed mass.

The pollen record from Anderson Fen (Figure 3) is similar to records from both lower and higher elevations on Vancouver Island, both in the immediate vicinity (Figure 5) and further afield (Brown and Hebda, 2003; Lacourse, 2005; Lacourse and Davies, 2015), revealing an ecotonal sequence of changing forest compositions as one would expect. Forests between 14,000 and 10,600 yr BP were dominated by P. contorta with local presence confirmed by needles at 11,920 and 10,670 yr BP. Abundant shrubs during this pine-dominated interval include A. viridis, and Salix and the insect-pollinated Shepherdia canadensis before 13,000 yr BP. These early forests occurred under cool climatic conditions (Figure 4). Several other conifers (Abies, T. mertensiana, P. sitchensis) increase in abundance during the particularly cool Younger Dryas chronozone (12,900–11,700 yr BP). Tsuga mertensiana is a well-established marker of the Younger Dryas cooling on the Northeast Pacific Coast (Mathewes, 1993; Lacourse, 2009; Gavin and Brubaker, 2015).

Forest composition near Anderson Fen changes rapidly at ∼10,600 yr BP with the establishment of P. menziesii and T. heterophylla forests (Figure 3). The transition from pine-dominated forests is also accompanied by increases in Cupressaceae, although it is not clear whether this is T. plicata and/or Chamaecyparis nootkatensis, due to the inability to differentiate their pollen. Acer macrophyllum, a low to mid-elevation species that is regularly found in association with P. menziesii on southeastern Vancouver Island, is usually poorly represented in fossil pollen records because it is mostly insect-pollinated, but it is also locally present during the transition (11,100–10,300 yr BP). The transition to P. menziesii forests was facilitated by increased temperature and seasonality and decreased precipitation (Figure 4). More frequent wildfires in the early Holocene would have also favored P. menziesii (Gavin and Brubaker, 2015). Over the remainder of the Holocene as precipitation increased, forests became increasingly dominated by T. heterophylla and Abies with the abundance of P. menziesii decreasing especially after 5,400 yr BP, although it continued to be present, as even low levels of its pollen indicate local presence (Heusser, 1985). The last 2,200 years is marked by increases in P. contorta and Pinus monticola, both of which are present in low abundance at Anderson Fen today. The large increase in Alnus rubra pollen percentages over the last 100 years (Figures 3, 7) primarily reflects disturbance in the region due to widespread logging.

There are few other mid-elevation paleoecological records in the region to compare to Anderson Fen. The record from Rhamnus Lake (Allen, 1995), which sits at 320 m elevation approximately 160 km southeast, is notably similar to Anderson Fen. Forests there transition from being dominated by Pinus 14,000–12,000 yr BP, to P. menziesii in the early Holocene, and then co-dominated by P. menziesii and T. heterophylla forests by 7,000 yr BP. Mid-elevation pollen records from sites on the adjacent mainland of British Columbia (Mathewes, 1973) show a similar sequence of changing forest composition, although P. menziesii is less common and T. plicata plays an important secondary role in mid-late Holocene forests, reflecting the overall wetter climate of that area. Further south on the Olympic Peninsula in Washington, United States, pollen and plant macrofossil records from Yahoo Lake at 710 m (Gavin et al., 2013) are also similar, reflecting the overarching control that macroscale climate plays in driving changes in forest composition. Differences of note relative to coastal BC are the prominent roles of P. sitchensis in late-glacial forests and A. rubra in early Holocene forests as well as earlier expansion of T. heterophylla by a few thousand years.

Comparison of changes in forest composition along the elevation gradient shows clear sorting of woody plant species through time and by elevation (Figure 5 and Supplementary Figure 4). It is also clear that temporal changes in climatic conditions interact significantly with life-history and stress tolerance traits to direct long-term patterns in forest community assembly (Figures 6, 7). Species with similar functional traits tend to be most abundant at similar times, whereas species that are functionally dissimilar are separated in time (Lavorel and Garnier, 2002; Swenson et al., 2011). Below, we focus on woody plant taxa that are well-represented by pollen records on central Vancouver Island to describe elevational sorting of species through time and the role species traits in changing forest composition.

Pinus contorta is most abundant at low to mid-elevations, especially during the relatively cold interval 14,000–10,000 yr BP (Figures 5–7). This is consistent with numerous studies along the northwest coast of North America that have documented the expansion of P. contorta populations following deglaciation with its presence on the North Pacific coast predating the records from central Vancouver Island (e.g., Heusser, 1960; Lacourse, 2009; Gavin and Brubaker, 2015). The life history and stress tolerance traits of P. contorta allowed this species to be among the first conifers to establish widespread presence and continue to facilitate its establishment on disturbed, rocky and waterlogged sites today. Relative to other trees, it has a “fast” life history (Charnov, 1993; Franco and Silvertown, 1996): early sexual maturity, high relative growth rate, low shade tolerance, and short life span (Supplementary Table 2). It has wide ecological amplitude and high tolerance to waterlogging, drought and frost, relative to other trees in this region, and as a result, it often occurs in habitats that are unfavorable to potential competitors. For example, P. contorta is often present in a stunted form in ombrotrophic bogs on the North Pacific coast, where high water tables preclude the establishment of other trees. There is an increase in P. contorta over the last 2,000 years across elevations on central Vancouver Island (Figure 5A) as well as at several other paleoecological sites in the wider region (e.g., Heusser, 1960; Cwynar, 1987; Lacourse, 2009; Lacourse and Davies, 2015). This region-wide increase in P. contorta is likely linked to the late-Holocene expansion of peatlands that support shrubby P. contorta var. contorta. With future changes in climate, P. contorta is likely to continue to expand on disturbed, xeric and waterlogged sites and to benefit from increased wildfire, although at regional scales its overall abundance may change little in the decades ahead due to losses in some areas and gains in others (Hamann and Wang, 2006; MacKenzie and Mahony, 2021).

Pseudotsuga menziesii is also most abundant at low to mid-elevations on central Vancouver Island, reaching its maximum abundance after P. contorta, between 11,000 and 8,000 yr BP. Previous studies have documented the northward migration of P. menziesii from the south (e.g., Heusser, 1960; Lacourse, 2009; Gavin and Brubaker, 2015). Pseudotsuga arrived first at low elevations and migrated upslope in the early Holocene as temperature and seasonality increased and precipitation decreased (Figures 5B, 7). More frequent fire in the early Holocene would have favored fire-resistant P. menziesii over most other conifers (Brown and Hebda, 2002; Gavin and Brubaker, 2015; Halofsky et al., 2018). Pseudotsuga then decreased in abundance at mid-elevations to become most abundant at low elevations, as it is today especially to the south in drier, rainshadow forests on southern Vancouver Island and in Washington state. In combination, these studies demonstrate both northward and upslope migration of P. menziesii in response to early Holocene warming, followed by contraction as temperatures cooled after 6,000 yr BP. In the RLQ analysis, P. menziesii plots away from all other species in association with traits such as adaptation to fire, high seed mass and dispersal capacity, large height, and long lifespans as well as early Holocene climatic conditions of relatively high temperatures, low precipitation, and lengthened growing seasons as reflect by high insolation seasonality (Figures 6, 7). This combination of traits gives P. menziesii a competitive advantage in dispersal, establishment and persistence (Weiher et al., 1999; Leishman, 2001). Given responses of P. menziesii to past climate change, northward and upslope range expansion is likely given projected future warming and increases in wildfires. Bioclimate-envelope and process-based models similarly predict that changes in climate and increased wildfires will expand coastal P. menziesii northward and upslope with a net increase in abundance and area occupied (Hamann and Wang, 2006; Gray and Hamann, 2013; Halofsky et al., 2018; MacKenzie and Mahony, 2021).

There are two common species of hemlock in temperate rainforests on the North Pacific coast and their long-term paleoecological histories on central Vancouver Island are similar. Both Tsuga species are most abundant after ∼8,000 yr BP, with increasing abundance since then at most sites; however, T. heterophylla has occurred primarily at low and mid-elevations and T. mertensiana is almost exclusively confined to high elevation forests (Figures 5C,F). The expansion of hemlocks over the mid-late Holocene was facilitated by higher precipitation over this interval (Figure 7). This similar paleoecological history is reflected in the RLQ analysis with both hemlocks plotting close to each other in the ordination space (Figure 6C), close to mean annual precipitation (Figure 6A) and with similar life history and stress tolerance traits (Figure 6B). For example, both hemlocks reach sexual maturity at about 20 year, and have the same or similar seed mass, dispersal capacity, relative growth rate, lifespan, and tolerance to shade, drought and waterlogging (Supplementary Table 2). Relative to P. contorta, both species are characterized by a “slow” life history. However, T. mertensiana has substantially higher frost tolerance, which accounts for its abundance at high elevations where temperatures are cooler, snow cover persists longer, and growing seasons are shorter. This also helps explain its presence in low and mid-elevation forests during the cool late-glacial period (Figures 3, 7; Mathewes, 1993; Gavin and Brubaker, 2015). Climate models predict increases in year-round temperatures and winter precipitation for much of Pacific Canada (Wang et al., 2016) that are likely to facilitate expansion of T. heterophylla northward and to higher elevations in response (Hamann and Wang, 2006; Gray and Hamann, 2013; Halofsky et al., 2018; MacKenzie and Mahony, 2021). For T. mertensiana, increases in temperature are likely to override any benefit of increased available moisture and potentially lead to contraction of its overall range. Indeed, Hamann and Wang (2006) and MacKenzie and Mahony (2021) note that T. mertensiana forests are particularly vulnerable to ongoing changes in climate. Northward and upslope shifts for the species are expected to be accompanied by substantial decreases in total area occupied, especially on Vancouver Island where T. mertensiana forests could be eliminated altogether within the next 50 years (Hamann and Wang, 2006; MacKenzie and Mahony, 2021).

The two deciduous nitrogen-fixing species, A. viridis and A. rubra, are also generally sorted by elevation through time. Alnus viridis occurs mostly at mid- to high elevations, reaching its highest abundance during the cool late-glacial period in association with P. contorta, and then moving upslope to occupy higher elevations almost exclusively for much of the last 10,000 years (Figures 5D, 7). In contrast, A. rubra has been an almost exclusively low-elevation species that increases in abundance over the Holocene (Supplementary Figure 4), along with species such as T. heterophylla and increases in available moisture (Figure 7). Both species of alder are characterized by a relatively “fast” life history and have similar shade, drought and waterlogging tolerances (Figure 5 and Supplementary Table 2). However, as with the two species of hemlock, they differ most notably in their tolerance of frost with the higher frost tolerance of A. viridis explaining its abundance at higher elevations. Based on their paleoecological histories and stress tolerance traits, A. viridis is likely to decrease in frequency on the North Pacific coast, whereas A. rubra is likely to increase. Bioclimate-envelope models support these inferences (Hamann and Wang, 2006; MacKenzie and Mahony, 2021). Increased disturbance from wildfire and human activities such as logging is likely to continue to benefit Alnus rubra in particular, as it readily invades open sites with bare soil.

Abies is most abundant at middle and high elevations but is more variable through time than other taxa (Figure 5). This variability is undoubtedly linked to the current inability to identify Abies pollen to species. On the North Pacific coast, Abies pollen derives from up to three different species, making it difficult to untangle their paleoecological histories. Abies amabilis is most common today at middle elevations, A. lasiocarpa is found in present-day high elevation forests, and A. grandis is present at low-middle elevations but is generally infrequent. Many paleoecological records in the region, including the Anderson Fen record (Figure 3), show an increase in Abies pollen in the cool late-glacial period between ∼13,000 and 12,000 yr BP, coincident with an increase in T. mertensiana. The late-glacial increase in Abies likely represents expansion of cool-adapted A. lasiocarpa (Heinrichs et al., 2002). In contrast, it is likely that expansion of A. amabilis populations occurred later, especially at mid-elevations (Gavin and Brubaker, 2015). The record from Burman Pond (Mazzucchi, 2010) provides macrofossil evidence for this species turnover at high elevations between 12,000 and 11,000 yr BP (Figure 5E). Changing pollen abundances at Anderson Fen (Figures 3, 5E) and Grant’s Bog (Lacourse et al., 2019b) suggest a similar turnover of species at lower elevations with Abies pollen in the late-glacial period likely deriving from A. lasiocarpa and Abies pollen after ∼8,000 yr BP likely deriving from A. amabilis. Abies has low loadings on both axes of the RLQ analysis (Figure 6C), which may be to some extent an artifact of being composed of two different species that are most abundant at different times. Overall, Abies is associated most strongly with high precipitation (Figure 6A), reflecting the overall preference of both Abies species for relatively high available moisture. Abies amabilis and A. lasiocarpa have similar life history and stress tolerance traits, although A. lasiocarpa has notably smaller seeds and higher frost tolerance (Supplementary Table 2). The latter explains its success in the late-glacial period and at high elevations, as is the case with T. mertensiana. Additional Abies macrofossil records or methods to differentiate Abies pollen to species are needed to resolve species dynamics of true firs throughout northwestern North America.