Introduction

Invasive plant species are expanding in riparian areas, which leads to catastrophic impacts on ecosystem services such as reduced habitat quality, decreased stream water flow, and lower vegetation cover and diversity (Rundel et al. 2014). Managing invasive species spread in these areas is challenging because mechanical treatments are expensive and may be destructive towards native species, few chemical treatments are available for use in semi-aquatic ecosystems, and heavy anthropogenic and animal use of these regions continually spread invasive plant species across shorelines (Dana et al. 2017). To understand the effects of plant invaders on shoreline plant communities, plant community composition should be quantified through time, i.e., succession, especially following management inputs aimed at reducing invasive species abundance (DeMeester and Richter 2010). This is particularly important for native plant species following mechanical removal because native plant community assembly after disturbance events may follow altered trajectories when an invasive species is present and native plant successional processes may not be sufficient for community recovery (de Abreu and Durigan 2011). Monitoring may be the most efficient and cost-effective way of understanding succession in high-disturbance ecosystems.

Permanent plots can be used to generate hypotheses regarding mechanisms for succession, identify the impacts of invasive species, and determine how plant communities reassemble after invasive species control (Bhaskar et al. 2014). However, succession generally requires long time frames to observe ecological changes (Buzzard et al. 2016) and identifying plant community change requires that the underlying ecological variation be defined and encompassed in replicate regions throughout the study area (Ulrich et al. 2016). Substituting historical space (successional stages) for time has been a common sampling approach to circumvent long-term monitoring (Prach et al. 2014) that can provide rich information in short time frames, especially in the face of major ecosystem changes like climate change (Peters et al. 2012). However, space for time sampling includes the assumption that spatial and temporal variation are similar (Pickett 1989) and thus, will typically average spatial data across the landscape which misses transient dynamics, like climate effects on vegetation production (Pickett et al. 2009). The weaknesses of the averaging approach combined with high likelihood of experimental failure leads us to consider long-term monitoring as the best approach for understanding succession in high-disturbance environments such as along the undammed Yellowstone River.

Monitoring functional groups and species of interest through time are common methods used to track succession, especially in areas affected by biological plant invaders (Byun et al. 2013). Accounting for the dominant functional groups allows researchers to forecast plant community response to management type and intensity (Hazelton et al. 2014). When functional group monitoring is combined with monitoring a species of interest, like a plant invader, researchers may accurately forecast plant community responses to environmental changes and management inputs when funding is limited and uncertain (Luzuriaga et al. 2015).

Russian olive (Elaeagnus angustifolia L.) is a nitrogen-fixer with cascading, ecosystem-wide impacts, including altering beaver population dynamics, watershed food webs, and nitrogen cycling (Collette and Pither 2015b). For example, Lesica and Miles (1999) found that beavers avoided foraging on Russian olive trees as compared to neighboring plains cottonwood trees along the Marias River in northern Montana. Furthermore, Russian olive has high seed dispersal, seed longevity, and resistance to herbivory (Edwards et al. 2014). These properties have allowed Russian olive to effectively spread and colonize many riparian areas, especially throughout the northern Great Plains (Espeland et al. 2014).

Forecasting plant community change in areas affected by invasive species, like Russian olive, requires understanding ecosystem dynamics within invaded and non-invaded ecosystems (Madurapperuma et al. 2013). Our objective was to set up a long-term monitoring project to understand successional transitions in a Russian olive-invaded river to determine post-removal management techniques and prevent reinvasion. Because our end goal is to determine probabilities of successional transitions, we wanted our monitoring plots to be replicated within successional entities, and for these entities to be real, functioning ecological units (Girvan and Newman 2002). To identify successional entities, we explore here the utility of a data-driven approach regarding the plant cover and soil properties that likely differentiate successional states affected by invasive tree species, like Russian olive, along a stretch of the Yellowstone River.

Materials and methods

Study site

This study occurred in seven floodplains (sites) bordering the Yellowstone River near Miles City, Montana, USA, all of which are either dominated by or susceptible to Russian olive invasion (Fig. 1). Study sites ranged from approximately 0.3–2.5 km apart and encompass 7 river kilometers. Three sites border the north side of the river; East Lower Blakely Creek (Site 1), Lower Cottonwood (Site 2), and the N. Forest Service Pasture (Site 3), while four sites border the south side of the river; Boat Ramp (Site 4), Woods Site (Site 5), West Yellowstone River (Site 6), and East Yellowstone River (Site 7). At all sites, four 0.5 ha plots that were approximately 0.5 km apart and 17 m away from the river bank were established with a permanent stake. At five sites (Site 2, 3, 4, 5, and 7) two plots at each site were in areas 50% or greater cover of Russian olive cover, while two plots were established in areas with less than 50% Russian olive cover. At the two remaining sites, all plots were established in areas with less than 50% Russian olive cover. Prior to analyzing these data, we excluded four of the original plots, (Site 3- Plot 2 and Plot 4, and Site 7-Plot 2 and Plot 3) because Russian olive had recently been mechanically removed at these plots.

Fig. 1
figure 1

Sampling stakes (locations) of the study site. Codes and names correspond as follows: Site 4-Boat Ramp, Site 1-East Lower Blakely Creek, Site 7-East Yellowstone River, Site 2-Lower Cottonwood, Site 3- N. Forest Service Pasture, Site 5-Woods Site. Numbers alongside each site denote plot number

Climate is typical of mid-continental zones, with temperatures that average a high of 14.7 °C and a low of 0.8 °C (W.R.C.C. 2017). Precipitation averages 315 mm and primarily occurs in spring with a flush of precipitation in autumn. Snowfall averages 740 mm in this region and generally occurs between October and May. Across and among sites, soil series vary. The dominant soil series mapped as the Glendive complex, Coarse-loamy, mixed, superactive, calcareous, frigid Aridic Ustifluvents, and dominated 38% of our study plots. Less dominant soil series include the Havre soil series, Fine-loamy, mixed, superactive, calcareous, frigid Aridic Ustifluvents, which mapped at 24% of plots; 24% of plots were also mapped as the Harlake soil series, Fine, smectitic, calcareous, frigid Aridic Ustifluvents; while the Hanly soil series, Sandy, mixed, frigid Aridic Ustifluvents, mapped in 13% of plots. These soil series are all found in 0–6% slopes and differ in the dominant soil texture and percent of calcium carbonates.

While woody vegetation among sites is dominated by the invasive tree, Russian olive, and the native tree, plains cottonwood (Populus deltoides W. Bartram ex Marshall ssp. monilifera (Aiton) Eckenwalder), there was also a mixture of trees and shrubs including native green ash (Fraxinus pennsylvanica Marshall), western snowberry (Symphoricarpos occidentalis Hook.), Woods’ rose (Rosa woodsii Lindl.), and chokecherry (Prunus virginiana L.) species, and the invasive tree Tamarisk (Tamarix ramosissima Ledeb.). Graminoids at this site included native perennial grasses western wheatgrass (Pascopyrum smithii (Rydb.) Á. Löve) and prairie sandreed (Calamovilfa longifolia (Hook.) Scribn.); the invasive perennial grasses smooth brome (Bromus inermis Leyss.), and crested wheatgrass (Agropyron cristatum (L.) Gaertn.); invasive annual grasses cheatgrass (Bromus tectorum L.) and Japanese brome (Bromus arvensis L.); and the invasive forbs spotted knapweed (Centaurea stoebe L.), and Canada thistle (Cirsium arvense (L.) Scop.).

Study design

Plots were established in 2014 and annually sampled for soil and vegetation in July 2014, 2015, and 2016. For this long-term monitoring study, we used a quick and cheap approach to ensure annual sampling regardless of varying budget cycles. To quantify understory species cover and density (plants m−2), we used the Daubenmire method (Daubenmire 1966) to sample three 20 cm × 50 cm Daubenmire frames that were randomly chosen using a random number generator to identify values ranging from 1 to 100 and pacing these values north, south, east, or west from a centrally located stake. In addition, 3–3 × 20 × 10 cm soil core samples were collected from the soil surface at each Daubenmire sampling location, mixed, and sifted through a coarse 2 mm screen to exclude detritus and roots, and sent to A&L Eastern Laboratories to identify the soil properties of the site. Soil properties included soil type (% sand, silt, or clay), organic matter (OM) (% and lb ac−1), soil nutrients (ppm and  % base saturation of phosphorus (P), potassium (K), calcium (Ca), and magnesium (Mg), pH, and cation exchange capacity (CEC) (meq 100 g−1) (©A&L Eastern Laboratories, Richmond, VA). Soil mineral nitrogen, total soil nitrogen, and total soil carbon abundance were also measured, but these data had little to no variation across sampling sites and thus were not included in our analysis.

In November 2014, July 2015, and July 2016 we returned to each site and quantified the dominant woody type and cover by visually estimating the cover of each tree and shrub species at each plot into four categories, 0–25%, 25–50%, 50–75%, or 75–100% cover. We also used the canopy-gap intercept method (Herrick et al. 2005) in November 2014 to quantify the frequency of each tree and shrub species at each plot using three randomly located 50-m transects throughout each plot area and recording the species type covering the transect line with a minimum gap size of 20 cm.

Data analyses

To quantify the spatial variation within and among sites, we used linear and non-parametric statistical learning techniques in JMP, such as principal components analysis (PCA), multivariate correlations, clustering techniques, discriminant analysis, CART-classification and regression trees, recursive partitioning rank variable tests, regression, and ANOVA models (Version 13; SAS Institute Inc., Cary, NC).

In an initial analysis, we used principal components model to identify major axes of plant and soil variation among and within sites. These analyses indicated that plant and soil variables had high variance both among and within sampling sites, (i.e., 85% variation among sites and 57% variation within sites). Because of the high among and within site variation, we chose to classify plots into groups (i.e., plant communities) based on over and under story plant cover using recursive partitioning. Recursive partitioning (RP) is a non-parametric classification method (CART-Classification and Regression Trees) that uses a decision tree algorithm which partitions/splits the input response data by a binary split at each step into sub-classes based on error reduction, resulting in terminal nodes of the tree representing clusters of observations with similar attributes, i.e., plant communities (Rouget et al. 2001). We used understory and woody plant cover data to classify plots into classification groups (or plant communities) with similar vegetation composition. The determination of the final number of nodes was based on observations per node and in this dataset. The outcome of our partition was five Leaves that were each uniquely classified by a tertiary species’ cover. A Wilks’ Lambda test was used to test model significance for this partition and treatments were considered significant when α was ≤0.05.

To evaluate soil variability among the 5 plant community types, we used a discriminant analysis. A log transformation was performed on calcium and magnesium (ppm) to improve model fit, while all other data were normally distributed. The RP analyses coupled with the discriminant analysis were used to identify plant-soil relationship patterns along the Yellowstone River.

Basal plant cover was determined using a mixed model ANOVA in JMP (Version 13; SAS Institute Inc., Cary, NC). This model tested the year, Leaf, and year by Leaf interactions on basal plant cover data collected throughout this study. Species were pooled by functional group into one of five categories; basal woody, invasive annual grass, invasive perennial grass, native perennial grass, or forb cover and we also took the sum of all plant cover for a total plant cover category. Total plant cover was normally distributed, but all functional groups required a square-root transformation to meet the standards of normality prior to analysis. Only three of the functional groups; invasive perennial grass, native perennial grass, and forb cover, and total plant cover had significant interactions from the Leaf, year, and Leaf by year treatments (P ≤ 0.10). Tukey’s HSD was then used to determine differences in Leaf and year treatments and treatments were considered significant when α was ≤0.10.

To understand trends in understory species distributions, we produced a correlation matrix with the ten most frequently sampled understory species. Species in this matrix included the exotic species crested wheatgrass, Japanese brome, smooth brome, Kentucky bluegrass (Poa pratensis L.), spotted medic (Medicago arabica (L.) Huds.), sweet clover (Medicago officinalis (L.) Lam.), and Russian olive seedlings; and the native species prairie sandreed, western wheatgrass, and western snowberry. Because all species were not always present in all sampling plots, our matrix had a low sample size (N = 18) and thus, we note all correlations that were significant at the P ≤ 0.15 level (Table 1; P ≤ 0.15).

Table 1 Correlation table of ten most frequent species found at sampling plots
Table 2 Four letter symbols, scientific names, and common names of the ten most frequent species growing at sampling stakes

Results

Five plant communities were selected from the RP (Fig. 2), because more than five partitions did not greatly improve the percentage of variance explained in the model (final R 2 = 0.57; Fig. 2). Each split was based on the third most dominant tree/shrub species (indicator species) cover at each plot (woody cover 3) and plant communities were uniquely identified by indicator species, represented in italics (Fig. 2). The first partition, which split shrubs (sagebrush, snowberry) from tree (Russian olive, juniper, tamarisk) explained 26% of the variation among all plots (Fig. 2). The second split separated the tree types; juniper and Russian olive from plots with no third level woody plant cover (open canopy) or were dominated by Tamarisk (plant community 3; 6 plots from 3 sites); this second partition explained an additional 13% resulting in a total of 39% of the variation explained. The third partition split the shrub node into Sagebrush (plant community type 1; 2 plots 2 sites) and Snowberry (plant community type 4; 4 plots 2 sites), these three splits (plant community types 1–4) accounted for 51% of the total variation of the whole model. The final split separated Russian olive (plant community type 2; 2 plots from 2 sites) from Juniper (plant community type 5; 4 plots from 3 sites), these four splits created 5 vegetation leaves and explained 57% of the variation in vegetation among sampling sites based on the tertiary species.

Fig. 2
figure 2

Tree diagram of LEAF partitions at reference plots. LEAVES represent sampling stakes (locations) that were sorted by plant cover (understory and woody) using recursive partitioning. All partitions were based on their tertiary woody plant cover and are represented by Sagebrush, Russian olive, Tamarisk/open canopy, Snowberry, or Juniper as tertiary species. Text above partitions identifies the percent variation explained by the partition. Text inside boxes describes the LEAF number (bold), LEAF name (italics), representative vegetation of the LEAF, and the representative sampling site (S) and plot (P) within that site. Note, site 1 and site 6 only have one plot location and, thus, are only represented by site (S) in this figure

Fig. 3
figure 3

Discriminant graph representing the classifications between LEAVES (tertiary dominant cover type- big circles) and soil properties (lines extending from center indicate strength of the variable—rays) at reference plots. LEAVES represent the dominant species present in each of the LEAVES; Tamarisk or. (Tamarisk; LEAF 1), silver sagebrush (Sagebrush; LEAF 2), western snowberry (Snowberry; LEAF 3), Russian olive (Russian olive; LEAF 4), or western juniper (Juniper; LEAF 5). Soil properties include percent sand, clay, OM, and canopy cover, soil pH, ppm of calcium (Ca), magnesium (Mg), phosphorus (P), and potassium (K), and the cation exchange capacity (CEC). Symbols represent the sampling stakes (locations) and their designated LEAF

Discriminant analyses of the five plant community types by soil variables (%sand, %clay, %organic matter, calcium, magnesium, potassium, phosphorus, cation exchange capacity, pH, and total canopy cover) produced a model with 26% misclassified plots (Fig. 3; P < 0.0001). The misclassification results from overlapping of Snowberry and Russian olive plots in the soil parameter space. Two canonical axes were satisfactory for this model and significantly explained 83% of the plot variability. The first canonical discriminant function is an axis dominated by % sand, CEC, and Ca and explained 58% of the total variation in the model. Along this axis, Juniper and Tamarisk are differentiated from Sagebrush. The second canonical axis associated with soil P and pH and explained an additional 24% of the total variation. Sagebrush and Juniper plant community types have similar values along axis 2; yet, Sagebrush had high total canopy cover, % clay, and P and low levels of Mg and K, while Juniper had high CEC and low levels of Ca and OM. Russian olive and Snowberry plant communities had high levels of Ca, pH, and OM and low CEC and % sand. Alternatively, the Tamarisk plant community type had high Mg and K, low canopy cover and low % clay.

Fig. 4
figure 4

Basal plant cover for each designated plant community. Different letters refer to differences in plant communities across years at α ≤ 0.10

Basal plant cover differed by Leaf and year in total plant cover and in three of the five sampled functional groups; invasive perennial grass, native perennial grass, and forb cover (Fig. 4; P ≤ 0.10). 2016 had about 20% higher plant cover than 2014 or 2015, which was likely due to higher precipitation during the 2016 growing season compared to those of 2014 and 2015 (Supplementary Fig. 1). However, where total plant cover was lowest in 2014 and 2015, i.e., the Russian olive and Juniper plant communities, it was highest in 2016 (Fig. 4a; P = 0.0930). Like total plant cover, native perennial grasses had the highest cover in 2016 and were highest in the Tamarisk and Russian olive plant communities compared to the Sagebrush, Snowberry, or Juniper plant communities (Fig. 4b; P = 0.0099). Across all years, Sagebrush plant communities had the highest cover of invasive perennial grasses and in 2016 invasive perennial grass cover was also higher in the Snowberry and Juniper plant communities compared to the Tamarisk and Russian olive plant communities (Fig. 4c; P = 0.0006). While forb cover did not greatly vary among sites or years; forb cover was higher in 2014 compared to 2015 in the Snowberry plant community (Fig. 4d; P = 0.0006).

When comparing plant species associations within plots, we found that most correlations were negative, and six of the seven significant correlations were between an invasive and native plant species (Table 1). A strong negative correlation of r = −0.58 occurred between the introduced grass, crested wheatgrass (AGCR), and the native, western wheatgrass (PASM) (P = 0.012). Crested wheatgrass (AGCR) also had a negative association with the native shrub, western snowberry (SYOC) (r = −0.35, P = 0.148), while western wheatgrass (PASM) was negatively correlated with the exotic forb, spotted medic (MEAR) (r = −0.39, P = 0.110). Alternatively, we found a positive association among the invasive grass Japanese brome (BRJA) and the native grass prairie sandreed (CALO) (r = 0.36, P = 0.145), a positive association among the exotic forb spotted medic (MEAR) and the native shrub western snowberry (SYOC) (r = 0.36, P = 0.137), and a positive association between the invasive Kentucky bluegrass (POPR) and the native western wheatgrass (PASM) (r = 0.44, P = 0.064). When two invasive plant species were correlated, it was between Russian olive (ELAN) and Kentucky bluegrass (POPR) and was negative (r = −0.38, P = 0.115).

Discussion

For this study, we strategically chose four plots in seven floodplains to represent areas a gradient of Russian olive plant cover and thus, characterize ecosystem properties in invaded compared to non-invaded Russian olive plots. However, our initial analyses revealed that the variation within and among sites was high, suggesting that a finer ecological scale was driving variation across these sites. While our original choice of plots based on dominant tree cover was conceptually well-supported by other studies (see Gonzalez-Rodriguez et al. 2011; Montoya et al. 2012); our statistical learning techniques revealed unexpected underlying ecological patterns involving less dominant plant species and their soil relationships (Figs. 2, 3). For example, our finding that plant community partitions were based upon the tertiary woody plant species cover suggests that ecosystem properties below the scale of biological invasions likely predict successional pathways after control activities, such as in Blank et al. (2015) and Hirsch-Schantz et al. (2014). In the case of Russian olive-invaded areas, our data show that succession will likely be driven by the interactions of soil properties with less dominant plant species.

Plant community types

In this study, we identify five plant community types, which were distinguished by tertiary woody plant cover and associated soil properties (Figs. 2, 3). Identifying plant community types by their tertiary woody plant species may be more practical in riparian regions invaded by Russian olive because the two primary woody plant species, cottonwood and Russian olive, tend to be generalists that grow throughout these riparian zones (Uresk 2015), while lesser dense species, i.e., tertiary woody species, likely correlate with the ecosystem properties occurring within these sites (Cushman et al. 2010). For example, we found that plant community types tended to differ by soil type; such as the tree species Juniper and Tamarisk plant communities that correlated with sand, and the shrub communities Sagebrush and Snowberry plant communities that correlated with clay. Alternatively, Russian olive plant communities initially partitioned with tree species (Fig. 2), yet had overlapping soil preferences with Snowberry plant communities (Fig. 3). This difference is likely because Leaf partitions, which identified distinct plant communities, were driven by the aboveground functional similarity of tertiary woody species, yet these tertiary woody species appear to have preferred soil growing conditions, which differ among species, regardless of these species aboveground functional similarity (Grime 2002). For example, it has been well documented that both western snowberry and Russian olive prefer heavier loam and clay textured soils (Bai et al. 2009; Fensham et al. 2015; Gaddis and Sher 2012), while tree species, like juniper and tamarisk generally prefer sandy soils (Nagler et al. 2011). Furthermore, Russian olive communities appear to be associated with high pH soils and that the probability of presence increased with increasing topsoil pH is consistent with Collette and Pither (2015a) who found that topsoil pH was an important variable in predicting Russian olive’s presence in the landscape. Soil pH is generally higher in clay and loamy compared to sandy soils because of negatively charged clay ions attracted to positively charged hydrogen ions, thus increasing soil pH (Fensham et al. 2015). Thus, Russian olive likely persists in areas with higher pH or moderate to high clay content.

Long lived species, like cottonwood trees, tend to indicate late successional plant communities (Klaar et al. 2015), while plant communities dominated by grasses and forbs are representative of early successional plant communities (Edmondson et al. 2014). In this study, plant communities with the highest cover of cottonwood trees, i.e., late successional plant communities, included those in Russian olive, Sagebrush, and Juniper leaves (Supplementary Table 1). Alternatively, where cottonwood is limited, such as in the Tamarisk or Snowberry plant communities, it seems likely that these regions may be at an earlier successional stage. Along riverbanks, plant communities are frequently disturbed by flooding events and commonly identified as early successional plant communities (Egger et al. 2015) and tamarisk is a fast-establishing species that generally prefers growing in frequently disturbed regions (Glenn and Nagler 2005). However, our finding that total basal plant cover did not differ among plant communities (Fig. 4a) and that the Tamarisk plant community had higher cover of native perennial grasses than Sagebrush, Snowberry, and Juniper plant communities in 2016 (Fig. 4b), seems to indicate that tamarisk invasion could be greater in regions dominated by native perennial grasses compared to invasive grasses. Furthermore, the higher cover in invasive perennial grasses in Sagebrush plant communities (Fig. 4c) could simply indicate that these woody species can easily establish and persist among seeded crested wheatgrass. In a recent study throughout the Great Basin, Williams et al. (2017) found a strong correlation between sagebrush growth and seeded invasive perennial grass (i.e., crested wheatgrass) plant communities. Consequently, while these plant communities may be in different successional states, it will likely take more years of continued monitoring to fully support that conclusion.

Managing plant community types

Correlations among commonly found plant species revealed that native and invasive species were generally negatively correlated (Table 1) and the strongest (negative) correlation was between the invasive crested wheatgrass and native western wheatgrass species. These species are both cool-season grasses and are arguably the most common native and invasive grasses in upland northern mixed-grass prairie ecosystems (Sanderson et al. 2015; Toledo et al. 2014). Crested wheatgrass effectively outcompetes western wheatgrass because of its faster growth rates, only requiring 82 days to produce a leaf compared to 98 days for western wheatgrass (Frank and Ries 1990). Furthermore, western wheatgrass has better water use efficiency than crested wheatgrass, which may allow western wheatgrass to persist in upland mixed-grass prairie ecosystems despite crested wheatgrass interference (Hamerlynck et al. 2016). This finding seems to suggest that species distributions are determined by plant traits and interactions with other species (Dawson et al. 2017) and that many invasive species likely fill the same ecological niche as natives (Yannelli et al. 2017). Furthermore, while our study did identify other invasive and native species functioning like crested and western wheatgrass, such as Japanese brome and prairie sandreed or western wheatgrass and Kentucky bluegrass (Table 1), the small sample size and P value restricted our abilities to extract meaningful and robust conclusions for these relationships.

Conclusions

Understanding plant community dynamics in diverse and dynamic ecosystems, like riparian areas, requires long-term monitoring to generate robust datasets across plant community types (Laine et al. 2015). While establishing plots based on species dominance is well supported (Gonzalez-Rodriguez et al. 2011; Montoya et al. 2012), this practice did not accurately represent our sites because of the high variation both within and among sampling sites. Instead, we identified five plant community types, all distinguished by their tertiary woody plant cover. Tertiary woody plant cover is likely a better predictor of plant community dynamics in these regions because dominant species tend to be widely distributed and locally abundant throughout these riparian areas, while underlying species appear to have preferred growing conditions determined by soil properties that may be spatially linked to successional processes. Quantifying underlying ecosystem variation prior to undertaking a long-term monitoring study should allow researchers to better understand initial differences between and forecast successional dynamics of at risk regions, like riparian areas invaded by Russian olive.