Site Selection and Study Area
Representative field sites were chosen from a pool of over 30 available sites based on location and size of invasive species populations, common native plant assemblages, site layout, and accessibility. Only non-tidal wetland mitigation sites displaying dominant patches of target organisms were used. Field sites were assigned age classes consistent with DeBerry and Perry (2012) due to the documented importance of age for plant community structure on wetland mitigation sites. On multi-user sites, distinct “phases” or areas constructed during separate time periods were treated independently so long as they fell into different age classes (this was a common condition on study sites operating as mitigation banks). The age classes were determined from site records on the number of complete growing seasons after site construction and included: 1-2 years old; 3-5 years old; 6-10 years old; 11-15 years old; and >15 years old.
Among the sites screened during site reconnaissance, 23 met suitability criteria and were selected for the study. Site ages ranged from 1 to 23 years post-construction and were evenly distributed across the Piedmont (11 sites) and Coastal Plain (12 sites) in Virginia (Figure 2). Most sites were either mitigation banks or in-lieu fee sites (i.e., sites established under the Virginia Aquatic Resources Trust Fund).
Research Design and Data Collection
Sampling transects consisted of five identical 4m2 (2m x 2m) vegetation plots, randomly assigned to an area that captured the gradient from completely invaded (i.e., the invasive species was considered dominant, or comprising at least 20% of the overall relative dominance of the community) to uninvaded (i.e., the invasive species was absent or not comprising more than 5% relative dominance).
Transect Configuration and Plot Locations
The randomization procedure for transect/plot layout involved identifying the center of an invasive species population within a given site and establishing a 4m2 grid with 9 vertices (Figure 3). Using a random numbers generator, a random number between 1 and 9 was selected, and its location on the grid was defined as the center of the vegetation plot for the most invaded site (Plot A). From that point, the direction of the transect was initially determined by defining an arc through which all possible transects could be delineated that would lead toward an uninvaded section of the site with similar environmental conditions. The length of this arc was taken as the domain for another random numbers draw, this time with the value representing the compass bearing from the center of Plot A to the edge of the invasive species population. At the edge of the population, another 4m2 grid was established and another random vertex was drawn, this one representing the center of Plot C. From this point, a straight line was defined from the center of Plot A to the center of Plot C and then an equivalent distance beyond the edge of the invasive population to delineate the remainder of the sampling transect. The center of Plot B (“second most invaded”) was then defined at half the distance between Plots A and C. The center of Plot E (uninvaded) was established at the far end of the transect, and the center of Plot D (“second least invaded”) was established at half the distance between C and E. This procedure resulted in five plots along the invasion gradient from most invaded (Plot A) to edge of invasion (Plot C) to uninvaded (Plot E) (Figure 3). Transect length varied among sites but typically ranged from 50 to 100 meters.
Soil Sampling
The center of each plot was GPS-located in the ESRI-based Collector application for iPad, then a soil sample was taken to a depth of 10cm using a 6cm-diameter soil corer. Soil samples were textured on-site using field methods (Ritchey et al. 2015), then shipped to the Virginia Tech Soil Testing Lab where soil chemical variables were measured with Mehlich extractions for P, K, Ca, Mg, CEC, Zn, Mn, Cu, Fe, and B, and Elementar high-temperature combustion for total values of C and N. Finally, an automated pH analyzer was used to measure pH values of wet samples at a 1:1 soil:water ratio (Maguire and Heckendorn 2019).
Canopy Cover
Canopy cover was measured by taking a skyward, hand-leveled photograph from the center of each plot using a 180-degree hemispheric lens adapter for iPad. Photographs were taken from 1 meter above the ground in Arthraxon and Microstegium plots, and from 2 meters above the ground in Typha plots. These photograph heights allowed us to capture canopy cover skyward of the target organisms while avoiding any potential effects of self-shading. Photographs were analyzed using ImageJ (Rueden et al. 2017) and the package Hemispherical 2.0 (Beckschäfer 2015) to obtain a ratio of open sky to canopy cover.
Vegetation Sampling
We quantified vegetation abundance using cover estimates for all species within each of four 1m2 subplots nested in the 4m2 plots. Cover estimates were based on a modified Daubenmire cover class scale with midpoints used for analysis (Mueller-Dombois and Ellenberg 1974). The cover classes, with midpoints in parentheses (rounded to the nearest whole integer), included: 0-1% (1%); 1-5% (3%); 5-25% (15%); 25-50% (38%); 50-75% (63%); 75-95% (85%); and 95-100% (98%). Cover classes were recorded for each species and then averaged across the four 1m2 subplots. Identifications of all vascular plants were either obtained onsite or samples were gathered and preserved for later verification. Intact collections were deposited at the William & Mary Herbarium (WILLI) following confirmation of identity by a senior botanist. Nomenclature follows Weakley et al. (2020).
Hydrology
Following transcription of the vegetation data, prevalence index (PI) values were calculated for use as a proxy of relative wetness (hydrology) between wetland sites (Atkinson et al. 1993; Tiner 2017). PI values are calculated from the wetland indicator status values for all species recorded within a plot. Wetland indicator status values are numbers assigned to wetland indicator status codes in accordance with the National Wetland Plantlist (Lichvar et al. 2016). The values include: 1=obligate wetland species (OBL); 2=facultative wetland species (FACW); 3=facultative species (FAC); 4=facultative upland species (FACU); and, 5=obligate upland species (UPL). Each species’ indicator status value is multiplied by the relative abundance of that species within the plot then summed to produce a weighted average index between 1 and 5. Plots closer to 1 are considered to have wetter conditions and plots closer to 5 are drier (Tiner 2017).
Statistical Analysis
Data analysis was completed using R version 4.0.3 (R Core Team 2020) including the packages vegan, Hmisc, and BiodiversityR (Kindt and Coe 2005; Borcard et al. 2018; Harrell et al. 2020; Oksanen et al. 2020). The datasets for each invader were analyzed separately due to expected variation in their relative tolerances for environmental stressors and discrepancies among growth requirements (Zedler and Kercher 2004; Swearingen et al. 2010). The correlation between relative abundance of each invader and variables in the environmental matrix was calculated using the nonparametric Spearman rank-order correlation test. The Spearman test was chosen due to its robustness to deviations from normality, as well as its ability to detect both linear and monotonic relationships, without appreciable loss of statistical power in comparison with parametric tests (Legendre and Legendre 2012).
Canonical Correspondence Analysis (CCA) (ter Braak 1986) was used to evaluate the overall community response to environmental variation along the invasion gradient. Prior to CCA analysis, rare species were removed from the abundance matrix of each dataset due to the outsized influence that rare species have on the Χ2 distance used in CCA (Legendre and Gallagher 2001; Peck 2016). Rare species reduction followed the Borcard method, which uses a stepwise approach based on the correspondence analysis (CA) component of CCA to evaluate the effect of progressive species removals (Legendre and Legendre 2012). Final CCA models were chosen with a combination of forward and backward model selection using the ordistep() function of vegan, which eliminates environmental variables based on significance of permutation tests in combination with the Akaike Information Criterion (AIC) (Borcard et al. 2018), in addition to variance inflation factor (VIF) to identify and remove highly correlated variables (Legendre and Legendre 2012). This procedure results in a parsimonious model when all environmental variables retained in the model are statistically significant and the adjusted R2 for the final model doesn’t exceed the adjusted R2 for the global model (global model = all environmental parameters included) (McCune and Grace 2002; Borcard et al. 2018). All permutation tests of significance were set at 1000 iterations, and statistical analyses were evaluated at α = 0.05.