DOI: 10.5061/dryad.r7sqv9sh8
title: Plant community compositional stability over 40 years in a Fraser River Estuary tidal freshwater marsh

description: 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.
description: 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 (&gt;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<sup>2</sup> 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<sup>2</sup> 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 <em>Carex lyngbyei</em> or
 <em>Juncus</em> 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 &lt; 25% cover, cover class 2 represents
 25–50% cover, cover class 3 represents 50-75% cover, and cover class 4
 represents &gt; 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
 <em>Agrostis</em> species, we assumed <em>Agrostis
 alba</em> L. identified in 1979 and 1999 was synonymous with
 <em>Agrostis stolonifera</em> L. in 2019. All species and
 their synonymous nomenclature from prior data collection years are
 available in Supplemental Table S7.
 We elected to
 classify <em>Phalaris arundinacea</em> 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 <em>P. arundinacea</em> 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
 (<em>Rubus armeniacus</em> 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.
description: 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.
description: # Plant community compositional stability over 40 years in a Fraser River
 Estuary tidal freshwater marsh | Name | Description | | ---- | -----------
 | | d1979\_formatted.csv | Plant community data collected in 1979 | |
 d1979\_formatted- METADATA.pdf | Metadata for plant community data
 collected in 1979 | | d1999\_formatted.csv | Plant community data
 collected in 1999 | | d1999\_formatted- METADATA.pdf | Metadata for plant
 community data collected in 1999 | | d2019\_formatted.csv | Plant
 community data collected in 2019 | | d2019\_formatted- METADATA.pdf |
 Metadata for plant community data collected in 2019 | |
 d2019\_position.xlsx | Location data for plant sampling conducted in 2019
 | | d2019\_position-METADATA.pdf | Metadata for location data for plant
 sampling conducted in 2019 | | CommunityStability.Rmd | All code used to
 conduct analyses found in manuscript and publication |