Long-term data sets documenting temporal changes in vegetation communities are uncommon, yet imperative for understanding trends and triggering potential conservation management interventions. For example, decreasing species diversity and increasing non-native species abundance may be indicative of decreasing community stability. We explored long-term plant community change over a 40-year period through the contribution of data collected in 2019 to two historical datasets collected in 1979 and 1999 to evaluate decadal changes in plant community biodiversity in a tidal freshwater marsh in the Fraser River Estuary in British Columbia, Canada. We found that plant assemblages were characterized by similar indicator species, but most other indicator species changed, and that overall α-diversity decreased while β-diversity increased. Further, we found evidence for plant assemblage homogenization through the increased abundance of invasive species such as yellow flag iris (Iris pseudacorus), and reed canary grass (Phalaris arundinacea). These observations may inform concepts of habitat stability in the absence of direct anthropogenic disturbance and corroborate globally observed trends of native species loss and non-native species encroachment. Our results indicate that within the Fraser River Estuary, active threat management may be necessary in areas of conservation concern in order to prevent further native species biodiversity loss.

Sampling design & harmonization between observations

Our main goal was to sample the vegetation in a representative way to allow comparison with the datasets collected in 1979 (Bradfield & Porter, 1982) and 1999 (Denoth & Myers, 2007). This publication will reference dates the data were collected, rather than publication dates of the corresponding studies.

In the original 1979 study, eight transects (Q-X) were laid out in a subjective fashion to cross through the main features of vegetation diversity at Ladner Marsh (Fig. 1 in Bradfield & Porter, 1982). All transects spanned a similar elevation range across the marsh platform, with the three main plant assemblages (Sedge, Fescue, Bogbean) separated by apparent changes in hydrological conditions along transects.

In the 1999 study, Bradfield & Porter’s (1982) Fig.1 was used to visually approximate the locations of transects to repeat the vegetation survey (Denoth & Myers, 2007). In this study (2019 survey), transect locations were determined by overlaying Bradfield & Porter’s (1982) Fig. 1 on a georeferenced basemap, aligning prominent features such as tidal channel tributary junctions, marking GPS locations in Avenza Maps (Avenza Systems Inc., Ontario, Canada, v. 3.2), and finding these points in the field (Fig.1C). Difficulties arising from the inexact relocations of transects in the 1999 and 2019 surveys, and aggressive shrub encroachment along transect Q, resulted in an incomplete resampling of all eight transects from the original 1979 survey (further explained below). To evaluate the potential for differences in transect relocation to affect trends observed in the data, or to evaluate marsh-wide spatial trends in plant composition, we calculated the percentage of plots clustered in each assemblage group for each transect.

All three studies used a semi-systematic approach for determining locations of 1x1 m quadrats along transects. In the 1979 study, quadrats were mainly located at 10 m intervals along transects although this varied in places from 2–20 m depending on local changes in the vegetation (Bradfield & Porter, 1982). In the 1999 study, an attempt was made to follow the quadrat spacing shown in Bradfield & Porter’s (1982) Fig. 3 regardless of perceived vegetation changes along transects. For the 2019 study, quadrat placement was guided by visual assessment of vegetation patchiness along transects. If patches dominated (>50 % cover) by one or two species (not necessarily the three assemblage identifiers) continued more than 10 m of transect length, or if no dominant species was evident, we sampled every 10 m of transect length (Fig. 2D). No patches were so small that the 1 m² plot was less than 1 m from the boundary of the next patch. Such fine-scale variations in decisions for quadrat placement among the three studies were considered inconsequential for the broader scale assessments of assemblage changes over time. 

Plot-scale sampling

Species were recorded if their most basal stem originated within the 1 m² quadrat, and cover within the plot was considered for all above-ground vegetation that occurred within the quadrat boundary; vegetation overhanging the quadrat from individuals whose basal stems originated outside the quadrat boundary was not considered. In the instance where the basal stem was inside the plot, but aerial vegetation extended beyond the boundary of the quadrat, we only considered vegetation cover for portions of the plant within the boundary of the quadrat. We treated each ramet of rhizomatous species such as Carex lyngbyei or Juncus sp. as individuals, rather than attempting to delineate extent of each continuous rhizome of the genetically distinct individual. For these species, whenever the quadrat fell on top of an individual ramet, the ramet was considered in the plot if more than 50% of the leaves emerging from the ramet were immediately under or inside the quadrat boundary. Aerial plot cover was estimated by modified Braun-Blanquet cover classes used by Bradfield & Porter (1982) and Denoth & Myers (2007), where cover class 0 = 0% cover (absent), cover class 1 represents < 25% cover, cover class 2 represents 25–50% cover, cover class 3 represents 50-75% cover, and cover class 4 represents > 75% cover. Owing to consultation with one of the co-authors (Gary Bradfield) in the 1999 and 2019 studies, differences in between-observer cover estimation were considered minimal.

Vegetation identification

For all sampling years, observation of vascular plant species was conducted in early summer when species are identifiable by sexual reproductive traits, but before senescence (approx. June – July). In all datasets, most plants were identified to species according to Hitchcock & Cronquist (1973), although a few were identified at higher taxonomic levels due to insufficient identifying characteristics (n = 6 to genus, n = 2 to Family; see Table S7). To account for changes in nomenclature revision over time, all datasets were harmonized to use the most recently accepted species name as reported in the PLANTS Database of the United States Department of Agriculture, Natural Resources Conservation Science [USDA NRCS]. For example, in the instance of Agrostis species, we assumed Agrostis alba L. identified in 1979 and 1999 was synonymous with Agrostis stolonifera L. in 2019. All species and their synonymous nomenclature from prior data collection years are available in Supplemental Table S7.

We elected to classify Phalaris arundinacea as non-native to align our treatment of the species with the designation provided by the British Columbia Ministry of Environment Species & Ecosystems Explorer (B.C. Conservation Data Center, 2023), which is the authoritative source for species conservation information for the province. While molecular analysis has confirmed P. arundinacea was native to North America prior to European colonization (Anderson et al., 2021), regional pollen studies have demonstrated some evidence for its absence in wetlands around the Salish Sea (Townsend & Hebda, 2013). Perhaps most important to consider is that hybridization of native with introduced varieties has resulted in aggressive invasive attributes, resulting in this species being of high management concern in Salish Sea Ecosystems (Sinks et al., 2021).

Differences between datasets

In 1999 and 2019, some plots were omitted due to access or relocation issues (Table S1). Most notably, transect “Q” (n = 7 plots) was omitted in 1999 and 2019 due to inaccessibility. In 1979, this transect was placed within approx. 100 m of Ferry Rd, which forms the eastern boundary of the marsh by an approx. 2-m elevated grade to keep the road above high tide elevations. In this portion of the marsh, riparian forest with an understory of non-native Himalayan blackberry (Rubus armeniacus Focke) grew so densely that by 2019, access to the transect would have required significant and costly vegetation removal. The encroachment of blackberry and riparian thicket was also a challenge preventing surveyors in 1999 from accessing this area. Thus, data from transect Q are not included in the present analyses.

The 1999 survey approximately located all plots from transects R-X from the original 1979 survey; however, the number of plots along these transects differed in 2019. This is partially due to the surveyors in 1999 seeking to exactly relocate original plot locations, while in 2019, our objective was to place the plots according to visual perceptions of shifts in dominant species. Besides the plots omitted by not sampling transect Q, we noted a total of 20 fewer plots surveyed in 2019 (Table S1). This is most likely due to our methods in 2019 placing plots to characterize patches dominated by distinct species, resulting in a number of plots being contingent on vegetation composition rather than spatial accuracy. Additionally, we acknowledge that spatial inaccuracy of transect relocation would result in different total transect lengths, and thus a different number of plots to be sampled along the transect. We also speculate there may have been some bank erosion resulting in wider channel mouths where some transects originated or ended, resulting in shorter transects overall. Visual comparison of satellite imagery suggests that erosion would have been minor, but not absent. To reconcile these differences, we excluded 1–4 plots per transect from the 1979 and 1999 datasets that had the least potential for spatial proximity to plots sampled in 2019 in order to compare an equal number of plots between sampling years along a similar length of transect (Table S1). 

Analyses

All analyses were performed in R v. 4.2.1 (R Core Team, 2022). We performed cluster analysis on species composition and abundance at the plot scale for each dataset. Similar to Bradfield & Porter (1982), we used Euclidean distance as the measure of plot dissimilarity (“stats,” R Core Team). Following Legendre & Legendre (2012), we also performed cluster analysis using Bray-Curtis dissimilarity to compare with Euclidean distance and found no meaningful difference in results from the two distance measures (Fig. S1).

For each dataset, we identified three assemblages, defined by the 3-group cut-level for each dendrogram, to facilitate direct comparisons of changes in vegetation properties over time. Species indicator analysis was used to determine which species abundance characterized each assemblage (“indicspecies,”R package De Cáceres & Jansen, 2016). Indicator Value (IndVal) association indices between species and clustered assemblages were calculated using an abundance-based point biserial correlation coefficient (multipatt func = “r.g”), and significance of associations was tested by permutational analysis (Dufrêne & Legendre, 1997). We also performed indicator analysis on the Bray-Curtis clusters to confirm that significant indicator species were comparable between the two distance measures (Table S2).  All species’ mean cover abundance is summarized in Table S6.

Community diversity calculations for each year of observation followed Whittaker (1975), with α-diversity calculated as the mean number of species per plot within an observation year and assemblage, and β-diversity calculated as the total number of species within the assemblage divided by α-diversity. These calculations were also performed on all data recorded for each observation year to generate community-wide measures of diversity. To address inconsistent numbers of plots grouped into assemblages each year, diversity metrics were bootstrapped 10 times using the minimum number of plots observed in an assemblage each year (n = 18) (Table S3). 

Community turnover for each assemblage was measured using the “codyn” R package (Hallett et al., 2016). Total species turnover (total magnitude of change), species gained (appearances), and species lost (disappearances) were calculated as a percent change for each assemblage between 1979–1999, and 1999–2019. Total turnover was calculated as a ratio of the absolute value of species gained and lost to the total number of species observed in both timepoints.

All analyses were performed in R Studio v. 4.2.1. Please see ReadMe, then metadata files for each data file. R code is annotated with notes about the analysis process and often links to resources used to produce analyses.