Management Implications
Restoration of sites previously invaded by non-native species is often challenging, as invasive plants easily reestablish, native plant populations may take years to recover, and plant communities may shift without active management. Management of invaded areas typically involves eradicating the invasive plant species and planting native species. However, this is a labor-intensive and expensive process that could require repeated management. In addition, these management techniques only account for aboveground changes following invasion, while research has shown that belowground processes may be transformed by invasion. Legacy effects are the result of invasion-induced changes to belowground biotic and abiotic processes that can last for years and prevent native plant communities from recovering. In this study, we examined the potential for legacy effects in two regions where Linaria vulgaris (yellow toadflax) occurs to determine whether these effects correspond with invasive spread. Managers in Colorado should consider how legacy effects could promote reinvasion and suppress native species. While arbuscular mycorrhizal fungal communities were not significantly different between invaded and uninvaded soil, other fungi or bacteria could contribute to the legacy effects we observed, and invaded sites may benefit from soil inoculation from uninvaded areas. This study sheds light on the changes that L. vulgaris induces to soils in Colorado, although future experiments are required to determine what the best management practices are to restore invaded soils.
Introduction
Invasive species reduce species diversity worldwide because of their ability to dominate natural areas (Mack et al. Reference Mack, Simberloff, Lonsdale, Evans, Clout and Bazzaz2000). It has been theorized that invasive species’ pervasiveness comes from their lack of predators and pathogens, their ability to colonize disturbed sites, or their novel interactions outside of their native environments (Keane and Crawley Reference Keane and Crawley2002; Simberloff Reference Simberloff2010). However, invasion success can be variable, with only a fraction of species becoming invasive and only in certain areas outside their endemic ranges (Mack et al. Reference Mack, Simberloff, Lonsdale, Evans, Clout and Bazzaz2000). Environmental factors, such as climate and nutrient availability, as well as the interactions between invasive and native species can all affect the range and severity of a species’ invasion (Hayes and Barry Reference Hayes and Barry2007; Yang et al. Reference Yang, Carrillo, Jin, Shang, Hovick, Nijjer, Gabler, Li and Siemann2013). Recent research has increased emphasis on the role of soil fungal communities in invasion success (Reinhart and Callaway Reference Reinhart and Callaway2006).
Interactions between native soil fungi and non-native plants are important determinants of whether a species becomes a successful invader. Plant–fungal interactions range from mutualistic to parasitic, often depending on whether plants gain or lose nutrients from the interaction (Johnson et al. Reference Johnson, Graham and Smith1997). For example, interactions with mycorrhizal fungi, which often increase plant nutrient uptake, may be responsible for a plant species’ success in a nutrient-limited environment (Richardson and Simpson Reference Richardson and Simpson2011; Yoshida and Allen Reference Yoshida and Allen2004). In nutrient-rich environments, mycorrhizal fungi may be parasitic, which could inhibit a plant species’ dominance in the community (Bever Reference Bever2002). Invasive plants often interact with soil fungi, including mycorrhizae, within their novel range, and these interactions have proven important in enhancing or inhibiting invasion in several species (Reinhart and Callaway Reference Reinhart and Callaway2006; Yang et al. Reference Yang, Carrillo, Jin, Shang, Hovick, Nijjer, Gabler, Li and Siemann2013).
While plants can have large ranges spanning continents, soil fungal communities are highly variable and depend on their plant host identity and local environmental conditions (Afkhami et al. Reference Afkhami, McIntyre and Strauss2014; Van Der Heijden Reference Van Der Heijden, van der Heijden and Sanders2002). When an invasive plant expands beyond its historic range, it may either escape its coevolved parasitic and pathogenic fungi or enter novel interactions that enhance or suppress its growth, influencing whether the species becomes invasive (Klironomos Reference Klironomos2002; Maron et al. Reference Maron, Klironomos, Waller and Callaway2013; Shah et al. Reference Shah, Reshi and Khasa2009). Several studies have demonstrated that invasive plants coexist with suppressive fungal associations in their native ranges, and beneficial fungal associations within their introduced ranges, suggesting that pathogen escape and/or new associations with mutualistic fungi may favor an invasive plants’ success in certain areas over others (Callaway et al. Reference Callaway, Thelen, Rodriguez and Holben2004; Dostálek et al. Reference Dostálek, Münzbergová, Kladivová and Macel2016; Reinhart et al. Reference Reinhart, Packer, Van der Putte and Clay2003).
Some invasive plants, such as scotch broom [Cytisus scoparius (L.) Link], have adapted to drive shifts in the fungal community around them, which could expand their range by increasing suitable habitat or give them an edge over native species (Davis et al. Reference Davis, Grove, Haubensak and Parker2021). By exuding chemical compounds through their roots or leaf litter, these plants alter the soil community by selecting their preferred mutualists or suppressing those of their competitors (Shah et al. Reference Shah, Reshi and Khasa2009; Zubek et al., Reference Zubek, Majewska, Blaskowski, Stefanowicz, Nobis and Kaputsa2016). Some of the most aggressive invasive plants in North America (e.g., garlic mustard [Alliaria petiolata (M. Bieb.) Cavara & Grande]) exert changes on the soil fungal community around them, indicating that this strategy is key to their success (Lankau Reference Lankau2010). Plant-induced changes to soil fungal communities may also have a lasting effect, which can influence plant community dynamics for many years (Gibbons et al. Reference Gibbons, Lekberg, Mummey, Sangwan, Ramsey and Gilbert2017). These lasting effects, referred to as legacy effects, can drive the future of a plant community by differentially influencing the abundance and dominance of different plant species.
Legacy effects of invasion often include changes to both fungal communities and nutrient dynamics important to plant growth. Plant-induced changes to the fungal community can have indirect effects on nutrient dynamics if plants influence the abundance of fungi important to nutrient uptake. Invasive plants can also directly alter nutrient availability by preempting or assimilating more (or fewer) nutrients than native species or by changing the soil texture, which could influence the rate of leaching from the soil (Corbin and D’Antonio Reference Corbin and D’Antonio2004; Evans et al. Reference Evans, Rimer, Sperry and Belnap2001). For instance, if an invasive plant changes either the nitrogen (N) availability in the soil or the mycorrhizal community, which is often important in facilitating plant N uptake, its invasion could have downstream effects on native plant growth.
Yellow toadflax (Linaria vulgaris Mill., Plantaginaceae), also commonly referred to as butter and eggs, is a forb native to Eurasia that is invasive throughout North America (Sing et al. Reference Sing, De Clerck-Floate, Hansen, Pearce, Randall, Tosevski and Ward2016). It reproduces both sexually from seed and asexually from rhizomes, making it extremely difficult to eradicate (Sing et al. Reference Sing, De Clerck-Floate, Hansen, Pearce, Randall, Tosevski and Ward2016). Linaria vulgaris is widespread throughout the United States but is particularly successful in western states throughout the Rocky Mountain region (Sing et al. Reference Sing, De Clerck-Floate, Hansen, Pearce, Randall, Tosevski and Ward2016). Outside its clonal reproduction, little is known about why L. vulgaris is such a tenacious invader in the western United States. Ambient soil N is historically lower in the Rocky Mountain region, where Linaria aggressively proliferates, relative to soil N in the midwestern United States, where Linaria is nonaggressive (Bowman et al. Reference Bowman, Murgel, Blett and Porter2012; Jelinski and Kucharik Reference Jelinski and Kucharik2009). This suggests that L. vulgaris may alter the soil fungal community or outcompete native plants for N in the Rocky Mountain region to invade more aggressively. Alternatively, L. vulgaris could be suppressed by soil fungi in the midwestern United States, thereby making it an ineffective invader.
This study addresses the biogeography of L. vulgaris invasion by asking why this species dominates western subalpine meadows and not midwestern prairies. Specifically, the objectives are to determine (1) whether L. vulgaris invasion leaves legacy effects that influence native plants and its own success, (2) whether soil microbial communities from Colorado (CO) versus Illinois (IL) influence L. vulgaris success, and (3) whether soil N availability and net uptake rates (NUR) explain the difference in L. vulgaris invasion aggression between CO and IL. To do so, we examined plant biomass accumulation and NUR in L. vulgaris and two native plants grown in previously invaded and uninvaded soil collected from IL and CO, in the presence or absence of N fertilization. Using soil samples collected from L. vulgaris populations at the Rocky Mountain Biological Lab (RMBL) in Gothic, CO (aggressive) and Cook County, IL (nonaggressive), we grew plants in all possible treatment combinations in growth chambers for 120 d. We also used next-generation sequencing to characterize arbuscular mycorrhizal fungal (AMF) communities in previously invaded and uninvaded soils from CO and IL.
Materials and Methods
Study Sites
Fieldwork was conducted from July to August 2020 at RMBL in Gothic, CO (38.9592°N, −106.9898°W) and in the Forest Preserves of Cook County in Cook County, IL (41.8897°N, −87.8057°W). RMBL is a research station located in Gothic, CO, within the Gunnison National Forest. Bulksoil (basel sandy loam) and rhizomes were collected from subalpine meadows (∼2895 m asl) characteristically made up of diverse grasses, forbs, and shrubs (Price and Waser Reference Price and Waser1998). Linaria vulgaris is an aggressive invader at RMBL, and many populations are sprayed with herbicide (glyphosate, Roundup®, Bayer, Leverkusen, DE) to control the spread. Bulk soil and rhizomes were collected from populations that had no history of herbicide use. The Forest Preserves of Cook County (∼178 m asl) are semi-urban natural areas in the Chicagoland area made up of prairies, woodlands, and wetlands (Forest Preserve District of Cook County 2020). Soil (orthents, clayey, nearly level) and rhizomes were collected from L. vulgaris populations without visible herbicide use, but herbicide history was not available at the time of collection.
Study Species
Linaria vulgaris is a forb native to Eurasia that has spread throughout North America since the 1600s (Parker and Gassmann Reference Parker and Gassmann2021; Sing et al. Reference Sing, De Clerck-Floate, Hansen, Pearce, Randall, Tosevski and Ward2016). It is an aggressive invasive in the western United States and Alaska, but it can be found throughout North America (McCartney et al. Reference McCartney, Kumar, Sing and Ward2019). In Colorado, L. vulgaris is listed as a noxious weed list B species, meaning it is managed to slow or stop its spread, whereas in Illinois, it is listed as an exotic species, which is not regulated (Colorado Natural Areas Program 2024; Illinois Department of Natural Resources 2024). Linaria vulgaris can reproduce both from seed and clonally from rhizomes, but the latter is believed to be much more prevalent, as its seed viability is ≤25% (Parker and Gassmann Reference Parker and Gassmann2021). Reproducing clonally allows L. vulgaris to grow in dense patches and subsequently exclude native species. Species in the genus Linaria have both vesicles and arbuscules characteristic of AMF association (Harley and Harley Reference Harley and Harley1987). Linaria vulgaris is an undesirable species for cattle forage, and it invades both rangeland and agriculture lands (Parker and Gassmann Reference Parker and Gassmann2021). While the stem-mining weevil (Mecinus spp.) has shown some success as a biological control, management typically consists of herbicide application (Parker and Gassmann Reference Parker and Gassmann2021; Toševski et al., Reference Toševski, Sing, De Clerck-Floate, McClay, Weaver, Schwarzländer, Krstic, Jovic and Gassman2018).
Field Soil and Rhizome Collection
We identified populations of L. vulgaris from research records in Gothic, CO, and iNaturalist reports in Cook County, IL (https://www.inaturalist.org). In the summer of 2020, soil samples were collected from 10 paired plots of L. vulgaris–invaded and adjacent uninvaded plots in each site. Paired plots were at least 75 m apart and were demarcated using N-S and E-W tape measurements to calculate plot area. In plots greater than 100 m2, we sampled from the centermost 5-m2 area to capture the area with the longest invasion history. For plots smaller than 100 m2, we sampled throughout the plot. To locate uninvaded plots, we selected a random intercardinal direction from the paired invaded plot and established an uninvaded plot 5 to 15 m from the invaded plot and of the same size. To ensure random sampling throughout plots, we used a random number generator to select points on the tape measures to collect each sample. In each site (IL and CO), we collected 160 soil cores (152.4-mm depth, 101.6-mm width; n = 320): 8 samples from each of 10 invaded plots and 10 adjacent uninvaded plots, which were maintained as independent samples. In invaded plots, soil was collected from directly beneath the nearest L. vulgaris individuals, and in uninvaded plots, soil was collected beneath the nearest plants. We also collected 32 rhizomes (each 2.4-cm long) of L vulgaris in each invaded plot, for a total of 640 rhizomes. The soil and rhizomes were used in growth chamber experiments, soil nutrient analyses, and next-generation sequencing at the Chicago Botanic Garden.
Soil Edaphic Factors
Soil samples from each field site were analyzed for texture; plant-available levels of NO3, NH4, P (as orthophosphate); and pH. For soil texture, soil samples were pooled per sampling plot and by invasion (IL invaded, IL uninvaded, CO invaded, CO uninvaded) with 10 samples in each combination (total n = 40) analyzed using the hydrometer method (Bouyoucos Reference Bouyoucos1962). Plant-available nutrients were analyzed for 320 soil samples collected from CO and IL (160 each; 80 invaded and 80 uninvaded). Inorganic N and P were extracted from each soil sample using deionized water (100g soil/L water; pH 6.8) and by vigorous shaking for 30 min. Extracts were then filtered and analyzed colorimetrically on a microplate spectrophotometer (BioTek Instruments, Winooski, VT, USA) using the vanadium (III) reduction method for NO3 (Doane and Horwath Reference Doane and Horwath2003), phenol-hypochlorite method for NH4 (Weatherburn Reference Weatherburn1967), and the malachite green method for PO4 (Baykov et al. Reference Baykov, Evtushenko and Avaeva1988). Each extract was tested for pH using a pH pen (Fisher Scientific, Hampton, NH, USA).
AMF Community Sampling
We randomly selected two soil samples for every invasion history treatment in each plot across both sites (n = 80) for amplicon sequencing. Genomic DNA was extracted from 0.25 g of soil using the DNeasy PowerSoil Kit and following manufacturer’s instructions (MO BIO, Qiagen, Carlsbad, CA, USA). Extracts were cleaned using Cytiva Sera-Mag Select beads before amplification (Global Life Sciences Solutions, Marlborough, MA, USA). A 550-bp section of small-subunit ribosomal DNA (18S rRNA) was amplified using the forward primer NS31 and the reverse primer AML2 (Lee et al. Reference Lee, Lee and Young2008). PCR reactions entailed 2 min at 94 C, 30 cycles of 30 s each at 94 C, 1 min at 59 C, 2 min at 72 C, and 10 min at 72 C. DNA amplification was verified by gel electrophoresis on a 1.5% agarose gel, which was then visualized on a Kodak Electrophoresis Documentation and Analysis System 290 (EDAS 290, Eastman Kodak, Molecular Imaging Systems, Rochester, NY, USA). Amplified DNA was cleaned again using Cytiva XP beads. Sample concentrations were then quantified in Qubit using a dsDNA HS assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). Samples were indexed using a Nextera XT Index Kit (Illumina, San Diego, CA, USA). Equimolar aliquots of each sample were pooled and sequenced at on an Illumina MiSeq sequencer (v. 3, 2 × 300-bp paired-end reads) at the Northwestern University NUseq core. The first sequencing run yielded low-complexity sequences, and as a result, samples were resequenced on the Illumina MiSeq v. 3 using 2 × 150 bp paired-end reads. Fragment size was measured using an Agilent 2100 Bioanalyzer System (Agilent Technologies, Waldbronn, Germany).
Sequence Processing and Analysis
Samples were demultiplexed in Illumina BaseSpace and processed in QIIME2 (Bolyen et al. Reference Bolyen, Rideout, Dillon, Bokulich, Abnet, Al-Ghalith, Alexander, Alm, Arumugam, Asnicar, Bai, Bisanz, Bittinger, Brejnrod and Brislawn2019). DADA2 was used to denoise samples, and reads were truncated at 120 bp where quality dropped below a q value of 20 (Callahan et al. Reference Callahan, McMurdie, Rosen, Han, Johnson and Holmes2016). In addition, only forward sequences were retained for further analyses, as previous studies have shown that forward sequences are sufficient to capture AMF diversity (Davison et al. Reference Davison, Öpik, Zobel, Vasar, Metsis, Moora and Gilbert2012). This processing resulted in a total of 7,092,099 reads. Sequences were included if they were reported at least 10 times and in 3 or more samples. Operational taxonomic units (OTUs) were clustered at 97% (96 OTUs), and chimeras and borderline chimeras were removed using q2-vsearch. Taxonomy was assigned using the q2-feature-classifier (Bokulich et al. Reference Bokulich, Kaehler, Rideout, Dillon, Bolyen, Knight, Huttley and Caporoso2018) naïve Bayes classifier with MaarjAm reference sequences (Öpik et al. Reference Öpik, Vanatoa, Vanatoa, Moora, Davison, Kalwij, Reier and Zobel2010). OTUs with less than 80% alignment to MaarjAm reference sequences were removed. OTUs were then compared against the MaarjAm Database (Öpik et al. Reference Öpik, Vanatoa, Vanatoa, Moora, Davison, Kalwij, Reier and Zobel2010) and GenBank (Benson et al. Reference Benson, Karsch-Mizrachi, Clark, Lipman, Ostell and Sayers2012) to determine AMF sequence identity.
Sequences that could not be identified to AMF family level were aligned in CLUSTAL Omega using default parameters, and the percent identities among sequences were determined. Pairwise sequences or multiple-sequence alignments with >97% similarity were classified as the same OTU. This resulted in 47 unique sequences for phylogenetic analyses. Unique sequences were then concatenated with a reference set of ribosomal small subunit (SSU) sequences (Krüger et al. Reference Krüger, Krüger, Walker, Stockinger and Schüssler2012) for major AMF genera (Scutellospora, Gigaspora, Acaulospora) as well as Glomus species and allied taxa (Archaeospora, Claroideoglomus, Ambispora, Funneliformis, Paraglomus, Rhizophagus). Alignments were trimmed, and a maximum-likelihood phylogenetic tree was inferred for the concatenated dataset using RAxML with default options. Eighteen sequences could be placed within the Glomus clade, while the remaining 29 sequences were consistent with Glomeraceae but showed no specific phylogenetic placement. Clustal Omega and RAxML analyses were undertaken on the CIPRES gateway (v. 3.3) at the San Diego Supercomputer Center (www.phylo.org).
Growth Chamber Experimental Design
The experiment consisted of a complete factorial design (two by two by two) in which we measured the growth and N uptake in Linaria and two native plant species (blackeyed Susan [Rudbeckia hirta L.] and fall witchgrass [Digitaria cognata (Schult.) Pilg.]) in response to combinations of three factors: (1) soil source, either CO or IL; (2) invasion history, uninvaded or L. vulgaris–invaded soil; and (3) two levels of soil N, ambient or N fertilized (N kg ha−1 supplied as NH4NO3). Pots contained either L. vulgaris or Rudbeckia and Digitaria grown together. This combination resulted in eight treatments with 10 replicate pots per treatment and plant species (n = 160). Soil samples were maintained independently in the experiment. To prepare the soil for the experiment, soil from each location and invasion type was gently mixed 1:1 (v/v) with coarse sterile sand (autoclaved). Pots (10.1 cm3) were filled with prepared soil and planted with either L. vulgaris alone (two rhizomes); L. vulgaris with native plants (two rhizomes and three seeds per native species), or native plants only (three seeds per native species). There was poor germination and rhizome sprouting in pots with both L. vulgaris and native seeds, so we excluded them from analyses.
The two native species used in the growth chamber experiment, R. hirta (Asteraceae) and D. cognata (Poaceae), are native to IL and CO and commonly found in disturbed sites. Additionally, as R. hirta is a forb, whereas L. cognata is a grass, they may respond differently to soil legacy effects, so by growing both we hoped to capture broader impacts of L. vulgaris invasion on native species’ success. Seeds for the experiment were donated by the Dixon National Tallgrass Prairie Seed Bank at the Chicago Botanic Garden. While R. hirta and D. cognata seed from IL was readily available, obtaining and cold-stratifying seed native to CO would have delayed planting. Seeds were cool (4 C), dry stratified for 30 d and then surface sterilized in a bleach solution (10% v/v) before planting. Linaria vulgaris rhizomes were washed in deionized water to remove excess soil. Rhizomes collected from CO and IL were planted in their soil of origin across all treatments.
Plants were grown in growth chambers using a 12/12-h photoperiod (day/night), 25/15 C (day/night), and 75% relative humidity for 4 mo. Plants were N fertilized every 30 d at rates of 8.9 kg N ha−1 yr−1 (6.1 kg N ha−1 yr−1 as NO3, 2.8 kg N ha−1 yr−1 as NH4) or were unfertilized (control). This level of N fertilization represents double the average N input from anthropogenic N deposition in Gothic, CO (4.05 kg N ha−1 yr−1, 3:1 NO3:NH4, 20-yr average; NADP 2022; Yoshida and Allen Reference Yoshida and Allen2004) and more than double the average N input from anthropogenic N deposition in Cook County, IL (3.8 kg N ha−1 yr−1, 3:1 NO3:NH4, 22-yr average; BassiriRad et al. Reference BassiriRad, Lussenhop, Sehtiya and Borden2015).
Stable Isotope and Biomass Analysis
At the end of the growing period, 5 ml of 14.37 mM 15NH4-15NO3 solution (>98.9 at %15N; Cambridge Stable Isotopes, Tewksbury, MA, USA) was injected into the root zone of each pot using 100-mm 11-gauge biopsy needles (Fisher Scientific, Waltham, MA, USA) (Dresboll and Thorup-Kristensen Reference Dresboll and Thorup-Kristensen2012; Yoshida and Allen Reference Yoshida and Allen2004), leading to a total of 2.012 mg N added per pot. After 24 h, plants were destructively harvested with each plant separated into root and shoot fractions, dried in biomass ovens (72 h, 80 C), and weighed. After weighing, leaf samples (n = 68) were collected from all treatments in which there was sufficient leaf material for analysis. Each leaf sample was ground to a fine powder using ceramic beads in Fastprep-24 MG with subsamples (1 to 1.5 mg), packed into aluminum capsules and analyzed for 15N atom% and N concentrations in a Costech 4010 ECS (Costech Analytical Technologies, Valencia, CA, USA) coupled to an ion ratio mass spectrometer (IRMS; Delta V Plus, Thermo Scientific) at the Northwestern University Stable Isotope Laboratory.
The 15N data were used to calculate the rate of N uptake per day (as NO3 + NH4). The 15N atom percent excess (APE) was calculated as the difference in atom% 15N between labeled and nonlabeled plants from the literature (typically atom% 15N < 0.1; Torres-Poché et al. Reference Torres-Poché, Mora, Boutton and Morrow2020). The net uptake rate (NUR, µg N g−1 shoot d−1) of N was calculated following McKane et al. (Reference McKane, Johnson, Shaver, Nadelhoffer, Rastetter, Fry, Giblin, Kielland, Kwiatowski, Laundre and Murray2002) as:
$$NUR = {{{{\left[ {{\rm{Plant\;N\;content\;}}\left( {{\rm{\mu g}}/{\rm{g}}} \right) \times \left( {{{{\rm{APE}}}}\over{{100}}} \right)} \right]}}}\over{{\left[ {{\rm{time\;}}\left( {{\rm{days}}} \right) \times \left( {{\rm{atom\% }}15{\rm{N}}/100} \right) \times {\rm{shoot\;mass\;}}\left( {\rm{g}} \right)} \right]}}}$$
For our calculations of N uptake, we assumed that the shoot 15N levels in each plant species were proportional to the amount of 15N label taken up from the soil. While this approach simplifies any physiological differences in plant N allocation, our experiment was conducted when plant growth (native plants) and reproductive demands (Linaria, in flower) were high, so that N was more likely allocated to above- rather than belowground structures. We also determined bulk %N from the same leaf samples used for IRMS to compare the effect of treatments on plant total N content for L. vulgaris and native plants.
Statistical Analysis
Statistical analyses were performed in R v. 3.6.2 (R Core Team 2019). We used linear models and linear mixed-effects models to determine the effects of state, invasion history, and fertilization as well as their interactions on biomass, plant presence, plant NUR, and plant %N (Supplementary Table S1). Model fit was verified using histograms and the qqPlot function in the car package (Fox and Weisberg Reference Fox and Weisberg2019). Statistical significance was designated by a P-value ≤ 0.05. We also performed Tukey post hoc tests to determine which treatments were significantly different. Results were then visualized using the ggplot function in the ggplot2 package (Wickham Reference Wickham2016).
Presence and absence of plants was designated with a 1 or 0, respectively, across all pots by treating all pots with shoot biomass greater than 0 as a 1 and all other pots as a 0. Mean plant presence was used as the response variable in models for each treatment. To compare presence/absence among treatments, we used generalized linear models with binomial distributions across all possible combinations of variables (fertilization and invasion history) within all datasets (site and plants grown).
To compare plant biomass across treatments, we divided root and shoot biomass into separate datasets. Biomass was analyzed separately among plants grown (L. vulgaris or native species) and for each combination of states and invasion histories due to varying sample sizes. We also removed pots without biomass from these analyses, as these were evaluated in presence/absence analyses. We used linear mixed-effects models to compare biomass between site and invasion histories and plot as a random effect using the lmer function in the lme4 package (Bates et al. Reference Bates, Maechler, Bolker and Walker2015). Biomass data for each combination of state and invasion histories were transformed by taking either the log or square to account for nonnormally distributed data. The package emmeans was used to generate Tukey pairwise comparisons (Lenth Reference Lenth2021).
We log transformed plant NUR and used linear models to compare NUR across invasion histories and fertilization treatments. We created linear models to compare plant bulk %N across treatments. We also generated linear mixed-effects models to test differences in soil pH, PO4, NO3, NH4 (log transformed), and percent clay, silt, and sand (untransformed) across sites and invasion histories with plot as a random effect.
We performed permutational multivariate ANOVAs (PERMANOVAs) using the adonis2 function in the vegan package to test whether AMF communities were influenced by invasion history, state, and individual edaphic factors (Oksanen et al. Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn, Minchin, O’hara, Simpson, Solymos, Stevens, Szoecs and Wagner2020). We performed an indicator species analysis using the multipatt function in the indicspecies package (Cáceres et al. Reference Cáceres, Legendre and Moretti2010). We also compared individual OTU and family mean relative abundances across states and invasion histories using Kruskal-Wallis tests. Relative abundances of AMF OTUs and edaphic factors were visualized on nonmetric multidimensional scaling plots using the metaMDS function in the vegan package (Oksanen et al. Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn, Minchin, O’hara, Simpson, Solymos, Stevens, Szoecs and Wagner2020). We used linear mixed-effects models to compare OTU and family richness between states and invasion histories.
Results and Discussion
Soil Chemistry and Texture
Soils in CO were more acidic and contained significantly higher levels of NH4 and silt than those in IL (Table 1; Figure 1). However, there were no significant differences in soil NH4 or NO3 availability between invaded and uninvaded soils. This was unexpected, given that soil NH4 is much greater in IL than CO, meaning that other local abiotic or biotic site characteristics may be more important determinants of invasion than simply N availability (Baron et al. Reference Baron, Reuth, Wolfe, Nydick, Allstott, Minear and Moraska2000; NADP 2022). Due to the scarcity of L. vulgaris populations in IL, our sampling plots were located along trails, roadsides, and train tracks where soil may be enriched by chemical runoff or physically disturbed. Plots in CO were located within a national forest where anthropogenic impact is relatively low but large vertebrate activity is high (e.g., white-tailed deer (Odocoileus virginianus)). This activity may have resulted in higher N inputs.
Table 1. Summary of soil chemistry and texture variables by state of collection (Colorado or Illinois) and invasion history including means ± SE and best-fit model results a
a Significant P-values are denoted in bold.
Figure 1. Nonmetric multidimensional scaling (NMDS) plots demonstrating overall differences in arbuscular mycorrhizal fungal (AMF) communities (A) between CO and IL soils, (B) between Linaria vulgaris–invaded and uninvaded soils in CO; and (C) between invaded and uninvaded soils in IL. Vectors denote significant edaphic factors (P < 0.05).
Aside from differences in N availability, the two study sites showed strong differences in soil chemistry and texture, with higher pH in IL soil and greater silt content in CO soil. These differences could also be influential in L. vulgaris invasion. For example, soils with higher silt content tend to have a strong water-holding capacity, but the high particle density can inhibit plant root growth and thus increase plant water stress. Differences in soil pH may exacerbate nutrient deficiencies and reduce AMF colonization, thereby inhibiting L. vulgaris invasion (Soti et al. Reference Soti, Jayachandran, Koptur and Volin2015). The variables that were not different between sites were PO4, clay, sand, and NO3, which is consistent with a previous survey (McDowell et al. Reference McDowell, Noble, Pletnyakov and Haygarth2023).
More interestingly, soil PO4 and clay were the only tested abiotic variables that differed between invasion treatments, suggesting that the observed legacy effects may be associated with differences in soil abiotic properties. One possibility is that L. vulgaris plants preferentially establish in microsites with high soil P and clay content. Alternatively, soil P may modulate AMF community composition (Singh et al. Reference Singh, Wu, Shao and Zhang2022) and/or plant growth (Richardson and Simpson Reference Richardson and Simpson2011), while clay may act as an absorbing surface for allelochemicals. In addition, soil PO4 may be altered by L. vulgaris root exudates or rhizosphere microbes following invasion. Further, the relationship between soil PO4 and invasion history differed between sites, with greater PO4 in uninvaded than invaded soil from CO and the reverse in IL, which could contribute to the divergent legacy effects we observed. Future experiments should test whether L. vulgaris directly influences soil PO4 or if plants preferentially colonize areas with higher clay and PO4.
AMF Community Composition
Next-generation sequencing resolved 61 OTUs. Of these 61 OTUs, 29 could not be identified to the family level. The remaining 32 OTUs belonged to 7 genera and 6 families. Average family abundances per sample in order from highest to lowest are: Glomeraceae (23 OTUs, 50.0% of total reads), Clarideoglomeraceae (2 OTUs, 20.6% of total reads), Paraglomeraceae (3 OTUs, 12.2% of total reads), Diversisporaceae (1 OTU, 4.5% of total reads), Archaeosporaceae (2 OTUs, 1.8% of total reads), and Ambisporaceae (1 OTU, <1% of total reads). The unidentified OTUs comprised the remaining 10.7% of total reads. Members of the Gigasporaceae were not detected in either site. ANOVAs demonstrated that there was no significant difference in either OTU or family richness between states, invasion histories, or the interaction between state and invasion history (Table 2).
Table 2. Results of ANOVA comparing arbuscular mycorrhizal fungal (AMF) operational taxonomic unit (OTU) and family richness and their interactions with site, invasion, and site and invasion together a
a Significant P-values are denoted in bold.
PERMANOVA revealed that there was a significant difference in AMF community composition between CO and IL (R2 = 0.12, P = 0.001; Figure 1). This was expected, as previous research has shown that AMF communities are localized based on geographic distance, climate, and plant community type (Kivlin et al. Reference Kivlin, Hawkes and Treseder2011). Such differences may enhance or suppress the growth of L. vulgaris and thus may be an important factor in the biogeography of its invasion (Shah et al. Reference Shah, Reshi and Khasa2009). Kruskal-Wallis tests on the mean relative abundance of each AMF family (Figure 2) revealed that Diversisporaceae were significantly more abundant in IL than CO soil (χ2 = 35.2, P <0.001), while Archaeosporaceae were significantly more abundant in CO than IL soils (χ2 = 10.8, P < 0.001). Previous studies have indicated that plant inoculation with members of the Diversisporaceae resulted in greater plant biomass and P uptake, whereas inoculation with Archaeosporaceae provided little benefit compared with control treatments (Säle et al. Reference Säle, Palenzuela, Azcón-Aguilar, Sánchez-Castro, Alves da Silva, Seitz, Sieverding, van der Heijden and Oehl2021). This points to a more mutualistic AMF community in IL than CO, which contradicts the trend in L. vulgaris biomass. Alternatively, IL soils may carry a heavier pathogen load, thereby preventing plants from accessing the benefits of AMF associations. Most pathogenic fungi are members of the Ascomycota and Basidiomycota, both of which we were unable to capture with our primers (Doehlemann et al. Reference Doehlmann, Ökmen, Zhu and Sharon2017).
Figure 2. Relative abundances of arbuscular mycorrhizal fungal (AMF) families by (A) site (CO, IL) and (B) between invaded and uninvaded soils. Operational taxonomic units (OTUs) that could not be identified to family level are listed as NA. Significantly different mean relative abundances from Kruskal-Wallis tests are represented by asterisks: *P < 0.05; **P < 0.01; ***P < 0.001.
AMF OTU composition was not significantly different between invaded and uninvaded plots in both sites (CO: R2 = 0.007, P = 0.859; IL: R2 = 0.021, P = 0.638; Supplementary Figure S3). This suggests that legacy effects were not the result of shifts within the AMF community. This finding was unexpected, given the effects of invasive plants on AMF communities documented in previous studies (e.g., Zubek et al., Reference Zubek, Majewska, Blaskowski, Stefanowicz, Nobis and Kaputsa2016). However, plant-driven microbial changes can take months to years to develop, so it is possible that the L. vulgaris populations we sampled had not colonized the area long enough to influence plant–AMF associations (Wolfe et al. Reference Wolfe and Klironomos2005). This possibility is supported by the similar levels of soil pH, N, and P in invaded versus uninvaded population. Despite the overall similarity in AMF community composition between invaded and uninvaded plots, analyses revealed that members of the Glomeraceae were more abundant in uninvaded than invaded soil across both sites (χ2 = 4.64, P = 0.031). Experimental studies have shown that plants inoculated with members of Glomeraceae had greater AMF root colonization and shoot biomass compared with plants inoculated with other AMF families (Säle et al. Reference Säle, Palenzuela, Azcón-Aguilar, Sánchez-Castro, Alves da Silva, Seitz, Sieverding, van der Heijden and Oehl2021). Thus, it is possible that L. vulgaris and native plants were more productive in uninvaded than invaded IL soils owing to the positive effects of Glomeraceae colonizing the roots. If this is the case, any negative legacy effects may not be caused by pathogenic microbes, but a lack of effective mutualists. Even so, this does not explain why L. vulgaris performed better in its own soil from CO. Future research should evaluate the roles and interactions between rhizosphere bacterial and fungal communities to better understand how plant–microbial associations with L. vulgaris might influence its distribution.
Indicator species analysis revealed several OTUs that differed either by state (Supplementary Table S2) or by invasion history (Supplementary Table S3). However, the limited taxonomic resolution in our sequence data made it difficult to distinguish differences in less-abundant AMF taxa between treatments. We had to remove many low-quality base pairs important in taxonomic identification, which limited our ability to define and classify unique OTUs. Possibly, the use of broader primers (e.g., ITS; Taylor et al. Reference Taylor, Walters, Lennon, Bochicchio, Krohn, Caporaso and Pennanen2016) could be used to glean differences in fungal and bacterial community that may better explain the different legacy effects. Instead, we used the relative abundance of AMF families to compare the effects of L. vulgaris invaded and uninvaded areas, an approach that utilized all available information of the AMF community and has been used in previous studies to identify broad patterns of plant–AMF interactions (Egerton-Warburton et al. Reference Egerton-Warburton, Johnson and Allen2007).
Multivariate analysis linking edaphic factors to AMF community composition showed that certain soil factors were significant correlates of the AMF communities (pH, N), consistent with geographic differences in soil properties between sites and the effects of N on AMF communities (Egerton-Warburton et al. Reference Egerton-Warburton, Johnson and Allen2007; Table 3). However, edaphic factors did not seem to drive differences in AMF community between invaded and uninvaded soil in both sites when considered together or separately. These results followed the trends of both AMF community and edaphic differences between our sites, signifying how regional environmental conditions may influence microbial community assembly (Kivlin et al. Reference Kivlin, Hawkes and Treseder2011).
Table 3. Permutational multivariate ANOVA (PERMANOVA) of effects of edaphic factors and their interactions with site, invasion, and site and invasion together on arbuscular mycorrhizal fungal (AMF) community composition a
a Significant P-values are denoted in bold.
Legacy Effects
In CO soil, L. vulgaris presence (χ2 = 7.72, P = 0.005; Figure 3), shoot biomass (χ2 = 3.98, P = 0.046; Figure 4), and root biomass (χ2 = 4.38, P = 0.036; Figure 4) were significantly greater in invaded soil than uninvaded soil, irrespective of N-fertilization treatment. The increased number of L. vulgaris plants growing in previously invaded CO soil relative to uninvaded soil supports our hypothesis that L. vulgaris shapes soil characteristics in CO, where it is aggressive, to become more successful. We attribute more L. vulgaris plants growing in invaded soil to increased sprouting of rhizomes, as we did not observe any plants produce seed during the study. Further evidence of positive legacy effects in CO soil is that both root and shoot biomass were greater among L. vulgaris plants grown in invaded than uninvaded soil. Increased biomass in previously invaded soil could indicate that L. vulgaris is better able to access resources such as water and nutrients due to microbial associations (Reinhart and Callaway Reference Reinhart and Callaway2006). Greater root mass creates more surface area for microbial colonization, resource acquisition, and resource storage and could lead to increased clonal reproduction (Pauchard et al. Reference Pauchard, Alaback and Edlund2003).
Figure 3. Percent of pots with plant survival and growth represented by plant presence/absence (1, 0) in pots containing (A) CO soil (n = 80) and (B) IL soil (n = 80). Error bars represent mean + 1 SE. Significantly different mean relative abundances from Tukey post hoc comparisons are represented by asterisks: *P < 0.05; **P < 0.01; ***P < 0.001. Missing bars (IL) signify treatments without any plant presence.
Figure 4. (A) Shoot and (B) root biomass for Linaria vulgaris and native plants grown in soil from CO (n = 47). Error bars represent mean + SE. Significantly different mean relative abundances from Tukey post hoc comparisons are represented by asterisks: *P < 0.05; **P < 0.01; ***P < 0.001.
Conversely, L. vulgaris only grew in previously uninvaded soil from IL, which could indicate negative legacy effects following invasion (Figure 3). Negative legacy effects could be caused by suppressive microbial interactions, such as those in its native range (Klironomos Reference Klironomos2002). However, we lacked statistical power to compare plant presence and biomass across treatments, so we cannot ascertain whether poor rhizome sprouting was the result of legacy effects or genetics, disease, poor quality of rhizomes, or other variables. Because we grew rhizomes in their soil of origin and did not study their population genetics, we cannot rule out genetically distinct populations as the source of L. vulgaris’s increased growth in CO soil compared with IL soil. A comparison of L. vulgaris populations from five western states showed that populations were genetically distinct from one another, with evidence of multiple introductions to the United States (Ward et al. Reference Ward, Fleischmann, Turner and Sing2009). Different plant genotypes within a species can form different plant–microbial associations that, in turn, may influence the extent of legacy effects we observed (Brown et al. Reference Brown, Grillo, Podowski and Heath2020). To maximize genetic diversity, we collected rhizomes from 10 paired plots within each site, but distinct populations between sites could produce contrasting legacy effects. Still, the absence of the strong positive legacy effect we observed in CO soil shows that legacy effects across the species’ introduced range could influence the distribution of its invasion.
We found some evidence that legacy effects impacted native plants differently across the two sites. However, there was no relationship between invasion legacy and native plant presence or shoot biomass in either CO or IL soil. While native plants grown in CO soil had greater root biomass in uninvaded than invaded soil (χ2= 4.31, P = 0.037), there was no effect of invasion on native plant biomass from IL soil (Figure 5). This implies that L. vulgaris invasion has a slightly negative affect on native plants grown in CO soil, but not IL soil. The negative effect of invasion on native plant root biomass and not abundance could be more indicative of impacts on native plant productivity and resource allocation to roots. Still, decreased productivity could result in diminished survival and reproduction over time. Managers in CO may experience challenges in native species’ reintroductions following invasion, whereas IL management might be centered around removal of L. vulgaris with the expectation that native plants will be able to recover.
Figure 5. (A) Shoot and (B) root biomass for Linaria vulgaris and native plants grown in soil from IL (n = 25). Error bars represent mean + SE.
Apart from legacy effects, we must also consider the preinvasion characteristics that may have made invaded areas more suitable for L. vulgaris. Rather than plants conditioning the soil via legacy effects, the sites where L. vulgaris historically grew may have localized microbial communities or soil chemistry that are more suitable for L. vulgaris establishment. We used a paired-plot design to reduce variation between invaded and uninvaded plots, but microbial communities and soil chemistry can vary over small scales and be highly localized (Nacke et al. Reference Nacke, Goldman, Schöning, Pfeiffer, Kaiser, Castillo-Villamizar, Schrumpf, Buscot, Daniel and Wubet2016). Similarly, we were unable to determine how long L. vulgaris populations had been established, which could influence the strength of legacy effects (Collins et al. Reference Collins, Spasojevic, Pombubpa and Diez2023). Ideally, we would have used the same soil across treatments and conditioned invaded soil by growing L. vulgaris for an extended period (Bever et al. Reference Bever, Westover and Antonovics1997).
An important caveat is that our study did not set out to test whether AMF communities drive L. vulgaris legacy effects. Instead, we posit a correlation between AMF community composition to the distribution of L. vulgaris invasion as one possible mechanism driving invasion success. A robust test of whether site-specific AMF and/or microbial communities influence invasion success by region would require reciprocal inoculation of soils from each site with the other site’s microbes followed by measurements of plant performance. Because our study compares AMF communities between the two sites, our results should be approached with caution when considering the relationship between invasion and AMF communities.
Plant N Levels and NUR
Interestingly, we did not detect any effect of N fertilization on the number of plants growing or biomass in either L. vulgaris or native plants across our two sites (Supplementary Figures S1 and S2). If plants were N limited, we would expect fertilization to increase the number of plants growing and/or biomass. We initially hypothesized that L. vulgaris would be a better competitor for N than native plants in CO soil where N is limited, allowing it to spread rapidly. Because both L. vulgaris and native plants were unaffected by N addition, this suggests that the legacy effects we observed were not related to N limitation.
Further support for our conclusion that N availability does not explain the range of L. vulgaris’s invasion came from analyses of NUR. The NUR in both L. vulgaris and native plants did not differ significantly by invasion status or fertilizer treatments in either site (Supplementary Table S1). However, an insufficient number of L. vulgaris plants emerged in IL soil to make a robust statistical comparisons of L. vulgaris NUR across treatments. If N uptake was a strong determinant of plant success, we would expect to see higher NUR in treatments where sprouting and growth were highest. Because invasion legacy did not appear to impact NUR, we suggest that legacy effects cannot be attributed to increased or decreased N acquisition or utilization. Similar levels of NUR in plants from N-fertilized and unfertilized treatments further illustrates that plants were not N limited, implying that environmental N availability was not a direct driver for L. vulgaris’s success.
Overall, our results indicate that soil characteristics, biological or chemical, unique to our two sites likely influenced the distribution and magnitude of L. vulgaris’s invasion. In addition, AMF communities were different between our two sites, which could influence the extent of aggressive invasion and native plant susceptibility to legacy effects. However, AMF communities were not different between previously invaded and uninvaded soil, which suggests that legacy effects are not the result of direct changes to AMF community composition. Further, NUR did not differ among invasion histories or fertilizer treatments in soil from either site, demonstrating that legacy effects are not the result of changes to N availability or uptake efficiency. Taken together, these results suggest that L. vulgaris invasion may alter rhizosphere microbial communities and/or soil chemistry in CO soil, but we did not have enough evidence to conclude whether the same is true in IL soil. Managers in CO should be mindful that previously invaded sites may be susceptible to reinvasion and take longer for native species to recover, whereas previously invaded sites in IL may not always show any considerable long-term impacts of invasion.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1017/inp.2025.4
Acknowledgments
We would like to thank Jeremie Fant, Hilary Noble, and Emma Leavens for assistance with DNA preparation, and Jennie Reithel for assistance locating field sites.
Funding statement
We thank the Plant Biology and Conservation Program at the Chicago Botanic Garden and Northwestern University, the Rocky Mountain Biological Lab, the Francis Y. Kwong Memorial, and the Elmhurst Garden Club for providing financial support for this project.
Competing interests
The authors declare no conflicts of interest.
Open access
