Bacteria common to rhizosphere communities of Asiatic bittersweet across a post-glacial landscape

Nakatsu, Cindy H.; Pavlovic, Noel B.; Byappanahalli, Muruleedhara N.

doi:10.3389/fevo.2022.999099

ORIGINAL RESEARCH article

Front. Ecol. Evol., 16 November 2022

Sec. Conservation and Restoration Ecology

Volume 10 - 2022 | https://doi.org/10.3389/fevo.2022.999099

Bacteria common to rhizosphere communities of Asiatic bittersweet across a post-glacial landscape

$\r\nCindy H. Nakatsu*&#x;$ Cindy H. Nakatsu^1*†

Noel B. Pavlovic^2†

Muruleedhara N. Byappanahalli²

¹Department of Agronomy, Purdue University, West Lafayette, IN, United States
²U.S. Geological Survey, Great Lakes Science Center, Lake Michigan Ecological Research Station, Chesterton, IN, United States

Invasive plants such as Asiatic bittersweet (Celastrus orbiculatus Thunb.) are a significant problem for land managers as they impact plant species composition, disrupt nutrient dynamics and structure of native ecosystems, and are difficult to eradicate. As a result of the increasing abundance of Asiatic bittersweet across the eastern U.S., we have been investigating underlying factors potentially contributing to the success of this plant. Recently, ecologists have been investigating the role of plant-soil-microbe interactions contributing to plant invasion. This work has led to question: are there certain microbes (e.g., bacteria, fungi) contributing to the success of some invasive plants? We hypothesize that despite differences in geologic age of soils where Asiatic bittersweet has established in the Indiana Dunes National Park there are sufficient common factors that led to common bacterial taxa in their rhizosphere. The objectives were to determine differences and commonalities in the soil chemistry, plant community and bacterial communities of Asiatic bittersweet plants. To achieve these objectives, bittersweet plants were collected at thirteen locations in the national park from soils ranging in geologic age from 150 to over 14,500 years. Surrounding soil chemistry, plant cover and the 16S rRNA gene amplicon sequences of rhizosphere soil bacterial communities of these Asiatic bittersweet were compared. Asiatic bittersweet coverage of sampling sites ranged from 2 to 77% averaging 52 ± 2%. There were statistically significant differences (p < 0.05) in alpha diversity (Shannon, Faith’s PD and Pielou’s evenness) and beta diversity (Bray Curtis, Jaccard, unweighted Unifrac, weighted Unifrac) among the samples when grouped by soil age or habitat. Despite these differences in the bacterial communities from different soil ages and habitats, some bacterial taxa (e.g., Bacillus, Streptomyces, Sphingomonas and Rhizobiales) previously found in other studies to be beneficial to plant growth were found in every rhizosphere community sampled. These microbes provide insight into a possible contributing factor to the success of this invasive plant at the Indiana Dunes National Park, and a strategy for future work to reduce the impact of Asiatic bittersweet establishment and offer some new strategies to manage this nuisance species.

Introduction

Various microbiota-related strategies are used by invasive plants to outcompete other plants. Often through chemical secretions from the roots, plants can discourage the growth of adjacent native species through competition or kill them outright (Callaway and Ridenour, 2004). In addition, these secretions may discourage or encourage growth of certain microbes (bacteria, fungi) in the adjacent soil, giving the aggressive invader a growth advantage over other native plant species (Pascale et al., 2020). Another strategy is that bacteria or fungi can provide nutrients to the plant that it could not acquire otherwise. For example, the interactions between beneficial microbes and plants (e.g., nitrogen-fixing bacteria, mycorrhizal fungi) have been greatly exploited in crop production (Macik et al., 2020) and agroforestry for a long time. Recently, ecologists have begun to investigate the role of plant-soil-microbe interactions in facilitating plant invasion (Coats and Rumpho, 2014; Zhang et al., 2019; Shahrtash and Brown, 2021). But these studies are limited in number leading to a dearth of knowledge on the impacts of different invasive species on the soil microbial communities.

Throughout the eastern US, the highly invasive liana, Asiatic bittersweet (Celastrus orbiculatus Thunb.) (hereafter bittersweet) is dominating many landscapes (Pavlovic and Leicht Young, 2011; Leicht-Young and Pavlovic, 2012; Leicht-Young et al., 2013). Asiatic bittersweet is especially problematic due to its rapid growth, woody vine growth, ability to resprout from the roots when burned or damaged, hybridization with native American bittersweet (Celastrus scandens L.), and strangling trees, making them more susceptible to falling over. Rapid expansion of Asiatic bittersweet is a challenging problem for land managers due to their impacts on composition, nutrient dynamics, and structure of native ecosystems and because they are often difficult to eradicate. While confined to the boundary between the Dunes State Park and the town of Dune Acres in the 1980s within the Porter County in Indiana, monitoring of this species at the Indiana Dunes National Park (INDU) has shown a tenfold increase in frequency from 1998–2000 to 2010–2012 (Pavlovic et al., 2022). As a result of the increasing abundance of Asiatic bittersweet at INDU and across the eastern U.S., it is imperative that we gain a greater understanding of the underlying reasons for the success of this plant to develop strategies for its remediation before vulnerable ecosystems are lost. The long-term objective of this work is to identify methods to reduce and possibly minimize or eliminate this invasive plant from the park using existing (e.g., prescribed burning, mechanical removal) and emerging technologies such as microbial control of nuisance invasive plants (Kowalski et al., 2015; Trognitz et al., 2016; Shahrtash and Brown, 2021).

Asiatic bittersweet is known to alter soil chemistry by increasing pH and soil concentrations of certain nutrients, such as potassium, calcium, and magnesium (Leicht-Young et al., 2009), likely influencing nutrient dynamics, successional pathways, and plant community structure. Nitrogen mineralization and litter decomposition rates were higher in soils impacted by the bittersweet, indicating elevated microbial activities. However, microbial community structure and its role in nutrient cycling was not directly assessed in that study. We hypothesize that despite differences in the soil ecology of Asiatic bittersweet establishment in the Indiana Dunes region, there are certain bacterial taxa common in their rhizosphere potentially enabling this invasive plant to outcompete adjacent vegetation. The objectives of this study were to (1) determine the range of different soil chemistries and plant vegetation where Asiatic bittersweet had established at INDU, (2) identify the differences in bacterial community composition and structure in the rhizosphere of these plants, and 3) determine if there were rhizosphere bacteria common to these invasive plants. Asiatic bittersweet plants were collected at thirteen locations in INDU from soils ranging in age from 150 to over 14,500 years. Surrounding soil chemistry, plant cover and the 16S rRNA gene amplicon sequences of rhizosphere soil bacterial communities of these Asiatic bittersweet impacted locations were compared to determine differences and commonalities among locations with Asiatic bittersweet establishment. This information will provide insight into potential underlying factors contributing to the success of this invasive plant and provide strategies for future work to control the spread of Asiatic bittersweet.

Materials and methods

Field sites

Asiatic bittersweet roots, bulk soil and adjacent vegetation were sampled on the lands of the Indiana Dunes National Park (41.63381°N, –87.05434°E), at the southern tip of Lake Michigan in Lake and Porter Counties, Indiana. The first herbarium collection of Asiatic bittersweet was collected in the Indiana Dunes in 1975 by Floyd Swink (Consortium of Midwestern Herbaria, 2022) and from 2000 to 2010 showed a tenfold increase in frequency across INDU (Pavlovic et al., 2022). Thirteen sites were selected across the Indiana Dunes incorporating sand dune and moraine habitats (Table 1). Sites encompassed a successional gradient from secondary dune-woodland transition to 14,500-year-old Glenwood dunes and two sites (Bailly and Chellburg Farm) on the older Lake Border Moraine [14.500 + years before present (ybp)] (Hansel et al., 1985; Argyilan et al., 2018). Habitats included secondary dunes transition, black oak savanna (recent and Tolleston dunes), woodland (Calumet dunes) and forest (Glenwood Dunes), and oak hickory forest. While the secondary dunes transition is sandwiched between grass dominated secondary dunes and oak savanna, these sites were on the edges of blowouts that formed approximately 150–200 years ago. However, blowout systems are undergoing a mosaic of disturbance occurring in the present time due to wind modifying sand movement across topographically complex landscape that produces patches of different successional stages. Consequently, we assigned woodland transitions to an age of 150 years rather than a younger age reflective of local blowout dynamics. Irrespective of the age assigned to the secondary dunes transition, they are younger than the other dunes/moraines and habitats. We sampled Asiatic bittersweet only and not American bittersweet (Celastrus scandens) because of concerns about its declining numbers and from the knowledge that sampling would kill the plants (Zaya et al., 2015).

TABLE 1

Table 1. Field sites, their habitats, geological phase, approximate system age, soil type.

Sample collection

Five independent Asiatic bittersweet plants were selected at each sampling site and each centered in 0.5 × 0.5 m quadrats for a total of 65 quadrats (13 sites × 5 plants). Samples were collected between August 24 and September 1 2017. The coverage and identity of herbaceous groundlayer vegetation surrounding the bittersweet plant was assessed within each quadrat. Representative bulk soil samples (multiple 10 cm × 1.0 cm soil cores) from each sampling location were collected from the quadrats for soil chemistry analysis. The soil core subsamples were pooled in the field and thoroughly mixed in the lab. From each of the five quadrats, one whole bittersweet plant was carefully excavated for rhizosphere analysis. Loose soil was shaken off the roots that were then placed into sterile bags for rhizosphere soil separation in the laboratory. All samples were placed on ice for transport to the lab. Bulk soils were stored at –20°C until processed. Root samples were processed after returning to the laboratory to collect closely adhered rhizosphere soils. The roots were placed in sterile centrifuge tubes containing phosphate buffer then vigorously shaken for 5 min then centrifuged. After centrifugation the buffer was poured off and the remaining rhizosphere soil was stored at –20°C until processed.

Soil chemistry

A portion of the pooled soil samples (for each of 65 locations) was sent to A&L Great Lakes in Fort Wayne, Indiana (N = 14, methods numbers follow (Soil Science Society of America, 2011) procedures based on (Page et al., 1982), for soil texture determination (hydrometer method) and a range of chemical parameters, including soil pH (Method 122.6), cation exchange capacity (CEC), base saturation of cations (% for Ca, Mg, Na, K Black, 1965), macronutrients (N, P (Method 24–5.1), K (Method 13–3.3.1), organic matter (Walkley-Black technique, Method 29–3.5.2), total and organic carbon, total N and micronutrients (S, Zn, Fe, Cu, B, Mn: various extractions and atomic absorption spectrometry).

Rhizosphere bacterial community analysis

Deoxyribonucleic acid extraction

A subsample (0.5 g) of rhizosphere soil was used for deoxyribonucleic acid (DNA) extraction using FastDNA™ kit (MP Biomedicals, LLC, Solon, OH, USA) following the manufacturer’s protocol. Agarose gels and spectrophotometry (Nanodrop 1,000 spectrophotometer; Thermo Scientific, Wilmington, DE) were used to ensure nucleic acid purity and quality. A NanoDrop 3,300 Fluorospectrometer (Thermo Fisher Scientific, Wilmington, DE) was used to quantify DNA that was stained with Hoechst dye.

Illumina 16S rRNA gene sequencing

Polymerase chain reaction (PCR) amplification of the 16S rRNA gene in soil samples was performed using V3-V4 region primers (343-forward TAC GGR AGG CAG CAG and 804-reverse CTA CCR GGG TAT CTA ATC C) (Liu et al., 2008; Nossa et al., 2010). PCR protocol and addition of primer tags were according to the manufacturer’s protocol (Illumina, San Diego, CA), with the exception that the first PCR reaction was performed for only 15 cycles and the second for 5 cycles to minimize PCR artifacts. PCR reactions were carried out using ∼10 ng of template DNA in Q5^® High Fidelity DNA Polymerase 2 × master mix (New England Biolabs, Ipswich, MA, USA). PCR amplicons were purified with AxyPrepMag PCR clean-up kit (Axygen Scientific, Union City, CA, USA) and subsequently quantified using a NanoDrop 3,300 Fluorospectrometer, after staining with the QuantiFluor dsDNA System (Promega). Equimolar amounts of amplicons from each sample were combined and sequenced (2 × 250 paired-end) using a MiSeq Illumina system at the Purdue Genomics core facilities, West Lafayette, IN, USA. Nucleotide sequences are available in NCBI-SRA database under accession number PRJNA826510.

16S rRNA gene amplicon sequence analysis

After primer tag removal using Illumina software, sequences were analyzed using the QIIME2 pipeline (version 2021.1) (Bolyen et al., 2019). DADA2 (Callahan et al., 2016) was used to check sequence quality, denoise, merge paired end reads and select amplicon sequence variants (ASV). ASVs with less than 0.1% frequency of the mean sample depth or found in only one sample were then filtered from the dataset. Taxonomic assignments were made using the Silva data set (version 138_99) (Yilmaz et al., 2014) then chloroplast, mitochondria and non-assigned sequences were filtered from the dataset. Datasets were rarified to use equal number of reads (12,600 reads) for alpha and beta community analyses. Bacterial alpha-diversity measurements calculated were Shannon diversity (richness and evenness), Chao1 index (richness), observed-features (ASVs), Faith’s phylogenetic diversity (PD) and Pielou’s evenness. Beta diversity metrics tested included both phylogenetic Unifrac distances (weighted and un-weighted) (Lozupone et al., 2011) and non-phylogenetic distances (Jaccard and Bray Curtis). Completeness of ASV representation in each sample was estimated by Good’s coverage, which was > 99.999% in all samples indicating that almost all the taxa were included in the analysis. Relative distribution of β-diversity measures among samples were visualized using Principal Coordinate Analysis (PCoA).

Statistical analyses

Prior to Principal components analysis (PCA) and cluster analysis, correlations among soil variables were calculated and variable were removed that had high correlations with other soil variables. Selected soil variables, P ppm, K ppm, pH, micronutrients ppm (Cu, Zn, Mn, Fe, S, B), total C and total N were standardized and centered in PCORD since variable units of measurement were different. PCA and cluster analysis was conducted in PAST (PAleontological STatistics) (Hammer et al., 2001). PCA was used to determine whether distinct clusters of chemistry were found across successional stages and ages and then were tested using discriminant analysis (SYSTAT, 2004) and a multiple response permutation procedure (MRPP) (McCune and Mefford, 2006).

Ground layer vegetation presence/absence and cover by growth form (forbs, graminoids, shrubs, trees, and lianas) was ordinated using non-metric multidimensional scaling (NMS) and NMS results were compared to soil chemistry variables.

For the microbial data most of the statistical comparisons were calculated for the samples grouped by geologic age and habitat. But note that habitat agglomerated geologic age groups within the oak woodland category (Table 2). Significant differences among alpha-diversity measurements were determined using Kruskal Wallis for all groups as well as pairwise comparisons followed by Benjamini-Hochberg FDR correction. Beta diversity differences among communities were determined using 999 permutations of perMANOVA (Anderson, 2001) of all groups and pairwise comparisons with Benjamini-Hochberg FDR correction. Then, permDISP (permutational analysis of multivariate dispersions) (Anderson et al., 2006) determined if dispersion could be contributing to significant differences. Differences were considered significant if p < 0.05 or q < 0.05 for pairwise comparisons.

TABLE 2

Table 2. Crosswalk between habitat and site age (ybp) for 13 sites sampled for rhizosphere bacteria, soil chemistry and vegetation.

Differential abundances in taxonomy across all habitats as well as geologic age were determined using Analysis of Composition of Microbiomes (ANCOM) (Mandal et al., 2015). ANCOM was used to identify differences in the log-transformed counts of taxa between habitat groups and geologic age groups. Bacterial taxa (genus and species level) common to all rhizospheres were determined using the “core feature” command in the feature tables plug-in in QIIME2. Taxa overlapping among the various habitats was illustrated using a Venn diagram program available from the University of Ghent Bioinformatics and Evolutionary Genomics group.¹

Canonical correspondence analysis (CCA) was used to find associations among the bacterial genera and differentiating soil and vegetation matrix factors using PAST and PCORD (McCune and Mefford, 2006). After removal of co-correlating factors, explanatory variables included in the preliminary model were dune age, slope, aspect, soil variables (pH, CEC,% organic matter, and ppm of B, Cu, Fe, K, Mn, P, S, and Zn), and vegetation (total vegetation, Shannon diversity (D), litter cover, bare ground cover, bittersweet and graminoids). CCA analysis in Vegan using the function ordistep were used to identify which soil and vegetation variables significantly contributed to the relationship with the bacterial genera (R Core Team, 2022).

Results

PCA of soil chemistry produced three components explaining, respectively, 42, 22, and 17% of the variation for a total of 81% (Figure 1A). The first axis separates soils with higher percentage sand from moraine soils. The second axis separates the sandy soils along the temporal chronosequence. The first principal component is strongly related to concentrations (mg/kg) of K, Mn, B, Zn, and total C and N. The second axis is an increasing pH gradient with negative loadings on P (mg/kg) and Fe (ppm). The third axis has highest loadings on Zn and Cu (ppm), an axis that is largely determined by one Calumet Dune sample that has high values of copper and zinc (10 ppm Cu, 20 mg/kg Zn). Sites with high soil iron included Glenwood 1 and 2, West Beach (> 50 mg/kg) whereas Central Ave., Kemil Road, and Portage Lakefront had low values (< 5 mg/kg). Soils were classified into four classes (discriminant Wilkes Lambda ₃₉ ₁₃₈ = 34.6, P < 0.001 with equal priors and 98% successful classification); i) Seventeen samples with high pH (5.1–7.9), sand content (94–96% and organic C/N ratio (11.3–82.1) representing young dune near Lake Michigan, ii) ten with high% silt,%clay, CEC, total C, total N, K ppm, Mn ppm and B ppm representing moraine soils (i.e., oak hickory), iii) 13 typical older sand dunes with low pH and high P, S mg/kg and Cu mg/kg, and lastly iv) 25 old dune samples with industrial impacts, low pH and high Zn and Fe (mg/kg). Soils were significantly different among geologic phase ages and habitats, but MRPP had a larger test statistic with dune age (MRPP dune age: T = –27.7, A = 0.67. P < 0.001, MRPP habitat: T = –25.1, A = 0.53, P < 0.001).

FIGURE 1

Figure 1. (A) Two dimensional principal components analysis (PCA) of Asian bittersweet bulk soil sample chemistry among the 13 sites. (B) Two dimensional non-metric multidimensional scaling ordination of vegetation presence/absence by species. Red biplots show the direction and magnitude of soil environmental variable corelates with vegetation composition. For both plots, each dot represents a sample and is colored by dune ages listed in the legend. Chemical abbreviations are B = boron, C = carbon, CU = copper, Fe - iron, K= potassium, Mn = manganese, N = nitrogen, P = phosphorous, and Zn = zinc. Samples were collected at Indiana Dunes National Park in 2017.

Vegetation composition

A total of 87 plant species were sampled including five exotic species: Asiatic bittersweet (Celastrus ortbiculatus), burning bush (Euonymus alatus), Asiatic bush honeysuckle (Lonicera japonica), and bouncing bet (Saponaria officinale), and 82 native species (nomenclature follows (Wilhelm and Rericha, 2017). All vegetation plots had Asiatic bittersweet ranging from 2 to 77% cover, averaging 52 ± 2% (± standard error (SE) and throughout) of the total cover of vegetation (7.5–100%). Mean cover of bittersweet was double in the black oak savanna (42 ± 7) and oak hickory forest (35 ± 7) compared to the secondary dune transition (15 ± 3) [F₍₅, ₆₀₎ = 28, P < 0.001]. All other comparisons were not significant. Average plant species richness per plot was 5 ± 1 ranging from 0 to 15 species per plot. Most common native species were Virginia creeper (Parthenocissus quinquefolia), jack-in-the -pulpit (Ariseama triphyllum), blue-stemmed goldenrod (Solidago caesia), sassafras (Sassafras albidum), and poison ivy (Toxicodendron radicans). Cover of native species was 21 ± 2 representing 45 ± 3% of the vegetation cover (0–93%). NMS of the presence absence vegetation data produced three ordination axes explaining the variation in vegetation composition (Figure 1B). Axis one was highly negatively correlated with percent sand (–0.55) and positively correlated with percent silt (0.55) and clay (0.53) representing a parent material productivity gradient with concentrations of manganese increasing along the axis. Zinc ppm (r = –0.51) was negatively correlated with axis 2. Silt loam oak-hickory sites were clustered in the upper right hand of the ordination graph but were not wholly distinct from the rest of the plots that had sandier soils (Figure 1B). Plants such as sugar maple (Acer saccharum), multiflora rose (Rosa multiflora), bedstraw (Galium aparine), wild lettuce (Lactuca canadensis) and gray dogwood (Cornus canadensis) distinguished the moraine soil sites from the sandy sites. NMS of growth form data produced a two-dimensional ordination (not shown) with the first axis negatively correlated with total vegetation cover (r = –0.69) and manganese ppm (–0.49) and positively correlated with bare ground (r = 0.42). Five growth form vegetation clusters were identified but they were based on growth form structure and were not related to habitat or dune age. Groups represented relative dominance of growth forms and one group stood out because of low vegetation cover compared to most of the plots and these were distributed among habitats and ages.

16S rRNA gene sequence analyses

PCR of one sample from the Central Avenue site (secondary dune) was unsuccessful and was excluded from all sequence analyses. Analysis of the remaining 64 samples produced 3,434,286 paired-end reads, with read numbers in individual samples ranging from 12,618 to 204,216 (mean reads 53,660 ± 35,941 S.D) after trimming, merging, and filtering for quality and non-bacterial reads. The read cutoff used for alpha and beta diversity analyses was 12,600. Goods coverage had a mean of 99.9% ranging from 99.9 to 100%.

Taxonomic annotations

Taxonomic annotations of the 16S rRNA gene sequences indicated there were 25 bacterial phyla represented in the collected rhizosphere soil samples. The phyla Proteobacteria (42.3 ± 1.2%, mean ± standard error), Actinobacteria (31.9 ± 1.2%), Acidobacteriota (15.3 ± 0.7%), Bacteroidota (3.5 ± 0.1%), Chloroflexi (2.1 ± 0.2%), Firmicutes (1.4 ± 0.2%), Myxococcota (1.3 ± 0.1%), Gemmatimonadota (0.9 ± 0.1%), and Methylomirabilota (0.5 ± 0.1%) cumulatively represent over 98% of the sequences (Figure 2). Other phyla and candidate phyla included, Bdellovibrionota, Cyanobacteria, Dependentiae, Desulfobacterota, Elusimicrobiota, Entotheonellaeota, Latescibacterota, MBNT15, Myxococcota, NB1-j, Nitrospirota, Patescibacteria, Planctomycetota, RCP2-54, SAR324_clade (Marine_group_B), Verrucomicrobiota, and WPS-2. These phyla represent 58 classes, 144 orders, 245 families, and 523 genera.

FIGURE 2

Figure 2. Comparison of relative abundances of phyla in Asiatic bittersweet rhizosphere soils collected from soils of increasing geologic ages (ybp). Only phyla with relative abundance > 1% in at least one soil age is included. Other phyla are represented by Bdellovibrionota, Cyanobacteria, Dependentiae, Desulfobacterota, Elusimicrobiota, Entotheonellaeota, Latescibacterota, MBNT15, Myxococcota, NB1-j, Nitrospirota, Patescibacteria, Planctomycetota, RCP2-54, SAR324_clade (Marine_group_B), Verrucomicrobiota, and WPS-2. Samples were collected at Indiana Dunes National Park in 2017.

Alpha diversity

Rarefaction curves reaching or approaching a plateau of all the alpha diversity metrics tested indicated 12,600 reads were adequate for making community comparisons (Supplementary Figure 1). Kruskal Wallis comparisons of alpha diversity metrics indicated there were significant differences for Shannon, Faith’s PD and Pielou’s evenness but not for the richness metrics Chao1 or observed features for samples grouped by geologic soil age or by habitat. Pairwise comparisons indicated there were some significant differences in the alpha diversity metrics for Faith’s PD and Pielou’s evenness between geologic age (Figures 3A,C) and habitat samples (Figures 3B,D). Faith’s PD and Pielou’s evenness values were significantly greater between the 14,000 + and 150 ybp samples. The Pielou’s evenness value was significantly greater for the oak-hickory forest habitat compared to secondary dune. Faith’s PD was significantly greater for both the oak-hickory forest and black-oak woodlands compared to secondary dune transition as well as to black oak forests.

FIGURE 3

Figure 3. Boxplots representing the medians and upper and lower quartile values of within sample (alpha) diversity metrics (A,B) PD whole tree (Faith’s phylogenetic diversity (PD), richness) and (C,D) Pielou’s evenness of samples from different soil geologic ages (A,B) and habitats (C,D). Kruskal Wallis indicates overall differences for all metrics tested are significant (p < 0.02) and pairwise-tests with 999 permutations with significant differences at q < 0.05 illustrated by different letters above bars. Outliers are represented by open circles. Samples were collected at Indiana Dunes National Park in 2017.

Beta diversity

PCoA of all the beta diversity metrics showed clear separation of samples by geologic age (representative example using unweighted Unifrac, Figure 4A). The unweighted Unifrac figure shows age separated mainly along PCoA1 but 14,500 ybp and 14,500 + ybp separated along PCoA2. PerMANOVA indicated the overall difference was significant for all the beta diversity metrics (p = 0.001) and pairwise comparisons also differed except for 2,000 compared to the 3,500–6,300 ybp samples for weighted Unifrac (Supplementary Table 1). However, permDISP was also significant for overall comparisons of all metrics (p < 0.05). The PCoA plot illustrates that there are differences in dispersion of samples among the different geologic age. But the 150 ybp samples clearly separate from the 14,500 and 14,500 + ybp samples along the first axis. There is overlap between the 2,000, 3,500–6,300, and 10,600–11,800 samples but clusters can be differentiated along axis 2.

FIGURE 4

Figure 4. Principal coordinate analysis of unweighted Unifrac beta diversity metrics of Asiatic bittersweet rhizosphere bacterial communities. Data points colored by (A) geologic age of soils (ybp- years before present), and (B) habitat type. perMANOVA comparisons of samples by age or habitat, overall differences significant p = 0.001 and all pairwise comparisons q < 0.02. Overall permDISP comparison also significant p < 0.03 for age and habitat. Pairwise analysis indicates that dispersion differences are significant (p < 0.05) between over 14,500 + y samples and all other ages except 14,000 y. For habitats, pairwise permdisp differences significant (p < 0.04) for oak hickory forest compared to secondary dune and to black oak forests, also black oak woodland compared to secondary dune and black oak forests. Samples were collected at Indiana Dunes National Park in 2017.

The differences noted in the PCoA plots of the geologic ages was also reflected when plots were differentiated by habitat (Figure 4B and Supplementary Table 2). As expected, the PCoA plot shows clear separation of secondary dune samples from the black oak and oak-hickory forest samples along PCoA1 and PCoA2. The differences are significant for habitat pairwise comparisons for all metrics except weighted Unifrac where black oak savanna did not differ from black oak woodland and secondary dune (Supplementary Figure 2 and Supplementary Table 2). Similarly, overall permDISP differed significantly but not for all pairwise comparisons.

Differentially abundant and common bacterial taxa

ANCOM identified similar bacterial taxa differing among the geologic age (Supplementary Table 3) and the habitat samples (Supplementary Table 4). There were six taxa that were highly significant based on volcano plots of the data (Supplementary Figure 3). In the oldest soils (14,500 + ybp) and oak-hickory forests there were significantly more uncultured Bacteriovoracaceae, Comamonas, Geobacter, and a Latescibacterota (Figure 5). In the 14,500 + ybp soils there were more Aurantisolimonas whereas there were more Starkeya found in the oak-hickory forest habitat (Supplementary Figure 4). In the youngest soils (150 ybp) and secondary dune transition there were more Promicromonospora.

FIGURE 5

Figure 5. Boxplot of medians and SD (standard deviations) of differentially abundant taxa identified using Analysis of Composition of Microbiomes (ANCOM) of microbial communities comparing the various soil geologic ages (ybp, years before present) along the y-axis. Taxa with W-values > 439 and clr > 34, (A) an uncultured Bacteriovoraceae (p__Bdellovibrionota), (B) Promicromonospora (p_Actinobacteriota), (C) Comamonas (p_Proteobacteria), (D) Geobacter (p Desulfobacterota), (E) an unnamed Latescibacterota and (F) Aurantisolimonas (p_Bacteoidota). W is number of time null hypothesis that were rejected and clr is the center log ratio.

More importantly for the objective of this study a Venn diagram comparison of the taxa present in each habitat indicated 290 of the 523 taxa (55.4%) identified at the genus level were found in at least one sample from every habitat (Supplementary Figure 5). The QIIME2 “core feature” program identified 23 taxa at the genus level found in every rhizosphere sample (Table 3). Analyses at the species are not included because none of the common genera were identified at the species level. A number of these core taxa are thought to be plant growth promoting bacteria, including seven taxa in the order Rhizobiales (Devosia, Hyphomicrobium, Labrys, Mesorhizobium, Bradyrhizobium, and two un-named Xanthobacteraceae genera), Bacillus, Streptomyces and Sphingomonas.

TABLE 3

Table 3. Bacteria at genus level of classification present in 64 samples collected at Indiana Dunes National Park in 2017.

Bacterial community association with soil chemistry and ground layer vegetation

Significant relationships existed between the bacterial genera distribution and soil chemistry and associated plant communities based on CCA (species- environmental correlation by axes 0 0.977–0.90–0.88, p < 0.001). The first and second canonical axes explained, respectively, 12 and 8% of the variation with a cumulative 20% and adding the third axis only increased the cumulative variation explained to 24%. The first axis represented the soil geologic age and nutrient gradient with secondary dunes transition bacterial communities on the left, black oak forest and oak hickory forest bacterial communities clustered on the right (Figure 6 and Supplementary Table 5). Oak hickory forest bacterial communities also clustered toward the right side of the plot but separated from the other samples along the second axis. The significant chemical factors relating to bacterial genera communities were pH (P = 0.005), total N (P = 0.005), dune age (P = 0.005), graminoid cover (P = 0.025), Zn (P = 0.03), total C (P = 0.04) and with Fe and slope included but did not individually contribute to the model (Figure 6, biplots and Supplementary Table 5). Zn is absent in Figure 6 since axis 3 is not shown (see factor loadings (Supplementary Table 5).

FIGURE 6

Figure 6. Canonical correspondence analysis (CCA) illustrating the relationship of bacterial genera (log transformed counts) to dune age, slope, soil chemical factors (pH, CEC, OM (% organic matter), concentrations of Fe, K, Mn, P, B ppm), % bare ground, litter cover, and graminoids. CCA1 and 2 explain 12.0 and 8.0% of the total constrained variation, respectively. Number of samples n = 64, overall p < 0.001 with 999 permutations. All variable tested in the model are listed in Supplementary Table 5.

Discussion

This study shows the broad range of soil conditions across a geologic age gradient in which the Asiatic bittersweet has successfully invaded in the Indiana Dunes region. These soils varying in geologic age differed in soil texture and productivity relative to substrate type, sand vs. clay. Within the sand dune habitats pH and total C and N varied in opposite directions, but otherwise the soils were fundamentally similar in physical structure and chemistry. Distinct microbiome communities are associated with the rhizospheres of Asiatic bittersweet. However, despite the differences there are also bacterial taxa common to all locations, indicating that these soils from different geologic ages each did not select for different bacterial communities. Indicating that these soils from different geologic ages did not differ substantially to select for entirely different microbial communities. Using alpha and beta diversity metrics to compare microbial communities in the rhizosphere of plants growing in soils of different ages and habitats, we show there are distinct differences in the composition and structure of rhizosphere microbiomes associated with Asiatic bittersweet along this soil chronosequence. Despite these microbiome differences, we found over half (290 of 523 total) of the bacterial genera were found in at least one sample from each soil age and 23 of them were found in every sample. These findings suggest that despite the recent entry of Asiatic bittersweet into these various habitats (since the 1970’s), there are certain bacterial communities, which have been shown by others to promote plant growth (e.g., Bacillus, Streptomyces), potentially associated with Asiatic bittersweet success in the park’s landscapes. But these taxa are very common soil bacteria (Delgado-Baquerizo et al., 2018) and unlikely newly introduced with the Asiatic bittersweet. Therefore, without further research we do not know why these taxa are common to these rhizospheres. Also, some groups, such as the Streptomyces, have both beneficial (Vurukonda et al., 2018) and pathogenic representatives (Loria et al., 2006).

The differences in soil chemical properties found in this study reflect trends reported over the last 50 years in the Indiana Dunes and Great Lakes sand dunes in general. Both Olson (1958) and Lichter (1998b) found that soil pH decreased from a nearshore peak between 7 and 8 to a low between 4 and 5 during the first 1,000 years of soil development. Other dune systems worldwide show a decrease in pH with time, but the magnitude of change and range varies (Nelson et al., 2020; Pollmann et al., 2020; Gaiero et al., 2021). Similarly, Lichter (1998a) at Wilderness State Park in lower peninsula of Michigan found that soil carbon and nitrogen asymptotically peaked between 500 and 1,000 years, respectively, at 40 and 2.3 mg/ha. The sometimes moderate to high concentrations of metals (Fe, Mn, Zn, Cu) are reflective of industrial steel production and human impacts on the Indiana Dunes soil (Cole et al., 1990; Esser et al., 1991; Lichter, 1998b). Despite this difference in soil chemical properties the Asian bittersweet has been able to rapidly establish across the Indiana dunes chronosequence (Pavlovic et al., 2022). Soil properties do not appear to influence invasion except for hydric soils where we have found that it has not been able to invade.

The ground layer vegetation composition alone did not separate into strong communities due to the dominance of invasive species and the small scale of the plots (Figure 1B. Soil ages and communities (not shown) are intermixed. Oak hickory understories on the moraine were distinctive compared to those on sandy soil that produce wide variation in composition with dune ages overlapping (upper right in Figure 1B). Bittersweet cover at the ground level was highest in the oak hickory forest and black oak savanna. The oak hickory samples were from secondary oak hickory forest at the edge of the older old growth oak hickory forest in the adjacent ravine, so canopy closure was less than in old growth forests. The black oak savanna was burned that resulted in higher densities of stems than otherwise would occur (Pavlovic et al., 2016 and pers. obs.). Graminoids were expected to be prevalent in the secondary dune transition, but their low abundance elsewhere may indicate how bittersweet shade impacts the native plant community in forested communities. In higher light levels (>70% PAR), bittersweet seedlings produced more leaves, thicker shoots, and longer and heavier roots compared to those grown in lower light levels (Greenberg et al., 2001).

The differences in microbiome composition along the soil chronosequence reported here have been also seen in other studies (Alfaro et al., 2017; Yeoh et al., 2017; Turner et al., 2019). Similar to other studies, pH appears to one of the main factors driving the differences in microbial community composition (Alfaro et al., 2017; Yeoh et al., 2017; Turner et al., 2019). The higher alpha diversity values in the oldest soils here was also reported in the Jurien Bay (Western Australia) dune chronosequence but this increase occurs up to soils in the 105 B.P. years and in much older soils the diversity decreases (Turner et al., 2019). Since these other studies were performed on dunal chronosequence soils in Australia and Peru these patterns of diversity differences are likely general to dunal chronosequences and not specific to specific geographic regions.

A model of the relationship among the soil factors, plant and microbial communities using CCA indicates that the main chemical factors associated with the different bacterial communities in the Indiana Dune soils of different geologic ages are pH, total N & C, and Zn content. Other studies have shown that soil chemical factors, especially pH, are major determinants of microbial communities (Lauber et al., 2009; Bahram et al., 2018). In soils of different ages plant succession and the members of the plant community (e.g., oak leaf litter) would contribute to the differences in soil pH (Olson, 1958) and as well as macronutrients such as carbon and nitrogen. Differences in bacterial communities potentially influenced by the plant communities are also evident in the CCA with grasses (Graminoids) more dominant in the younger soils making up the secondary dune transition with the shrubs (Liana) in the black-oak forest. These findings are similar to others along a soil chronosequence where soil pH declines with dune age (carbonate weathering, see Tombeur et al., 2020; Nelson et al., 2021), soil organic matter increases, and bacterial communities shift in composition as soil chemistry and vegetation changes (e.g., Cline and Zak, 2015; Alfaro et al., 2017; Yeoh et al., 2017). Leaf litter of Asiatic bittersweet has been shown to increase soil pH, and levels of potassium, calcium and magnesium (Leicht-Young et al., 2009). The increase in pH would counter the successional trend to decrease pH, but we do not have data to support that hypothesis from this study.

Despite the soil, vegetation and microbial differences expected along a soil chronosequence there were still many bacteria genera common among the soils. A set of core plant community associated rhizosphere bacteria was also found along a dunal chronosequence in the Great Sandy National Park on the Queensland coast of Australia (Yeoh et al., 2017). The authors thought that this core root microbiome had evolved over the 400-million-year history of the terrestrial plants in that park. While our dune system was less than 20 thousand years old, the native plant communities and bacterial communities might have coevolved for a longer time. But in our study the Asiatic bittersweet is a newly introduced species to these locations but were still able to out compete the native plant species. This is partially due to some of its growth strategies, such as, rapid growth with copious roots, and vine like habit enabling it to strangle trees (Fryer, 2011). But there is major lingering question. What factors do bacteria contribute to its rapid growth and successful invasion along this soil chronosequence?

Many of the bacteria genera common to the rhizosphere in all the samples are known to be involved in plant growth promotion and nitrogen acquisition. There are a number of taxa in the order Rhizobiales found in the rhizosphere samples. Although this group typically infects and forms nitrogen-fixing nodules in leguminous plants species, there is mounting evidence these bacteria are also beneficial to non-leguminous plants (Garrido-Oter et al., 2018; Fang et al., 2019). Laboratory growth experiments using a Rhizobium strain isolated from the non-leguminous invasive plant species Ageratina adenophora showed it was the best at improving plant growth (Fang et al., 2019). Unfortunately, we do not have direct evidence that any of the Rhizobiales associated with these bittersweet are fixing nitrogen or providing some other growth advantage to the plants. Rhizobia were also found as part of the core microbiome associated with the roots in the Australian dunal study (Yeoh et al., 2017). This suggests that although some bacterial taxa have possibly evolved along with the plants over geologic time, this invasive species potentially has mechanisms to out compete those plants to successfully invade diverse ecosystems. Since there is an apparent a core of microbes with plant growth promoting traits, we hypothesize that these Asiatic bittersweet plants do not need to enhance the growth of new microbes for invasion success but only some means of optimizing attraction of existing growth promoting bacteria to their rhizosphere. Since all the soil biota were not assessed in this study, it is possible that some other biota (e.g., fungi) are also contributing to the success of these plants. For instance, Asiatic bittersweet has been known to harbor endomychorrhizae in its roots possibly aiding its invasion ability (Leicht-Young et al., 2016). But it is likely that soil biota contributes to the success of these invasive plants to outcompete the indigenous plants.

The identification of some common bacteria from this study might provide a microbe-based means to control this plant species in diverse habitats. One possible means of reducing unwanted bacteria is to use taxa specific phages to reduce their numbers. Phage therapy is being used in medicine to treat human pathogens (Pirnay et al., 2022; Venturini et al., 2022) and more recently to control plant pathogens (Ye et al., 2019). Unfortunately, the bacteria associated here with Asiatic bittersweet can be beneficial to other plants in this ecosystem therefore it may not be reasonable to eliminate them. The Asiatic bittersweet is fairly widespread from Ontario and Quebec in Canada to mainly the states in the Atlantic seaboard and into Great Lakes region but its distribution is expanding (Fryer, 2011). There might be antagonistic microorganisms (bacteria, fungi) that can have parasitic relationships with plants rather than beneficial organisms that can be used to control these plants. Other possibilities would be application of nematodes to prey on the roots of the target plants (Wilschut and Geisen, 2021). Therefore, additional studies are needed to better understand the interactions between bittersweet and soil microbial communities, including fungi, to better understand the role of these communities in plant invasion and to identify potential microbiota for controlling this nuisance species. One weakness of this work is that funds were not available to determine the distribution of fungi in these samples. This study now paves the way for future studies to test microbial-plant interactions using controlled greenhouse experiments (Bickford et al., 2022). Elucidating the mechanisms used by invasive plants for a competitive advantage over native species also has broader implications for microbially influenced plant-plant interactions. For example, similar mechanisms may be used by weeds to provide them competitive advantage over crop plants. At this time, we do not know if there are common signals produced that give some plants a competitive advantage over others.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA826510. Soil and plant data are available in USGS data release https://doi.org/10.5066/P99W7Y45.

Author contributions

CN, NP, and MB conceived and designed the research. CN performed the bioinformatics and statistical analyses for the bacterial community data. NP collected the field data, plant identification, and statistical analyses of the plant and soil data. CN wrote the article with significant contribution from NP. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by the US Geological Survey. The USDA National Institute of Food and Agricultural Hatch grant IND01825 provided partial funding to CN. Any use of trade, product, or firm names was for descriptive purposes only and did not imply endorsement by the U.S. Government.

Acknowledgments

We thank Erin Argylan and Todd Thompson for assistance with dune formation dates, Andrew Davies and Justin Chapelle for collecting the field samples, and Stephanie Grebinoski for laboratory assistance.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.999099/full#supplementary-material

Footnotes

^ https://bioinformatics.psb.ugent.be/webtools/Venn/

References

Alfaro, F. D., Manzano, M., Marquet, P. A., and Gaxiola, A. (2017). Microbial communities in soil chronosequences with distinct parent material: the effect of soil pH and litter quality. J. Ecol. 105, 1709–1722. doi: 10.1111/1365-2745.12766