Front. Microbiol. Frontiers in Microbiology Front. Microbiol. 1664-302X Frontiers Media S.A. 10.3389/fmicb.2021.675352 Microbiology Original Research Soybean Roots and Soil From High- and Low-Yielding Field Sites Have Different Microbiome Composition Bandara Ananda Y. 1 Weerasooriya Dilooshi K. 1 Trexler Ryan V. 1 Bell Terrence H. 1 2 Esker Paul D. 1 * 1Department of Plant Pathology & Environmental Microbiology, The Pennsylvania State University, University Park, PA, United States 2Intercollege Graduate Degree Program in Ecology, The Pennsylvania State University, University Park, PA, United States

Edited by: Diogo Neves Proença, University of Coimbra, Portugal

Reviewed by: Angel Valverde, Consejo Superior de Investigaciones Científicas, Spanish National Research Council (CSIC), Spain; Juntao Wang, Western Sydney University, Australia

 *Correspondence: Paul D. Esker, pde6@psu.edu

This article was submitted to Microbe and Virus Interactions with Plants, a section of the journal Frontiers in Microbiology

 30 11 2021 2021 12 675352 03 03 2021 03 11 2021 Copyright © 2021 Bandara, Weerasooriya, Trexler, Bell and Esker. 2021 Bandara, Weerasooriya, Trexler, Bell and Esker

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The occurrence of high- (H) and low- (L) yielding field sites within a farm is a commonly observed phenomenon in soybean cultivation. Site topography, soil physical and chemical attributes, and soil/root-associated microbial composition can contribute to this phenomenon. In order to better understand the microbial dynamics associated with each site type (H/L), we collected bulk soil (BS), rhizosphere soil (RS), and soybean root (R) samples from historically high and low yield sites across eight Pennsylvania farms at V1 (first trifoliate) and R8 (maturity) soybean growth stages (SGS). We extracted DNA extracted from collected samples and performed high-throughput sequencing of PCR amplicons from both the fungal ITS and prokaryotic 16S rRNA gene regions. Sequences were then grouped into amplicon sequence variants (ASVs) and subjected to network analysis. Based on both ITS and 16S rRNA gene data, a greater network size and edges were observed for all sample types from H-sites compared to L-sites at both SGS. Network analysis suggested that the number of potential microbial interactions/associations were greater in samples from H-sites compared to L-sites. Diversity analyses indicated that site-type was not a main driver of alpha and beta diversity in soybean-associated microbial communities. L-sites contained a greater percentage of fungal phytopathogens (ex: Fusarium, Macrophomina, Septoria), while H-sites contained a greater percentage of mycoparasitic (ex: Trichoderma) and entomopathogenic (ex: Metarhizium) fungal genera. Furthermore, roots from H-sites possessed a greater percentage of Bradyrhizobium and genera known to contain plant growth promoting bacteria (ex: Flavobacterium, Duganella). Overall, our results revealed that there were differences in microbial composition in soil and roots from H- and L-sites across a variety of soybean farms. Based on our findings, we hypothesize that differences in microbial composition could have a causative relationship with observed within-farm variability in soybean yield.

 soybean spatial yield variation microbiome metagenomics fungal ITS bacterial 16S rRNA metabarcoding R2018-P09 Project PEN04660 and Accession 1016474 Project PEN04651 and Accession 1016233 Pennsylvania Soybean Board10.13039/100016706 U.S. Department of Agriculture10.13039/100000199 Introduction

Spatial heterogeneity in soil characteristics within agricultural fields is a frequently observed phenomenon that can influence yield and other crop traits (Khakural et al., 1996; Jaynes and Colvin, 1997; Paz et al., 1998; Kravchenko and Bullock, 2000; Vollmann et al., 2000; Batchelor et al., 2002; Shahandeh et al., 2005). Spatial field heterogeneity typically implies the occurrence of soil fertility gradients, owing to trends in different soil parameters, such as soil physical properties (Becher, 1995), element concentration (Berndtsson and Bahri, 1995), organic carbon content (Ball et al., 1993), and water and nitrogen content (Mulla et al., 1992; Wade et al., 1996). Spatial field heterogeneity can be closely linked with spatial variability in crop yields. In fact, understanding observed yield variability within corn and soybean fields is one of the most captivating problems in production research (Batchelor et al., 2002).

The site-scale spatial heterogeneity of various soil characteristics and their relationships with soybean yields have been previously reported in several locations including Iowa (Kaspar et al., 2004; Rogovska et al., 2007), Mississippi (Cox et al., 2003), and Michigan (Jiang and Thelen, 2004). These studies indicated that site-scale heterogeneity in soil factors, such as pH, P, K, slope of the site, site elevation, and soil texture explained the observed variability in soybean yield. Mann et al. (2019) studied bulk soil samples from farmer-designated “good” vs. “poor” sites across 34 farms in Maritime Canada and used Cornell’s Comprehensive Assessment of Soil Health (CASH), which integrates soil texture, available water capacity, surface and subsurface hardness, wet aggregate stability, autoclaved-citrate extractable soil protein, and soil respiration (Moebius-Clune et al., 2016), combined with phospholipid fatty acid analysis (PLFA), and conventional soil chemical analysis to assess the collected samples. Although CASH and PLFA assays did not statistically differentiate farmer-designated “good” and “poor” sites, differences between good and poor sites were evident when focusing on certain soil nutrients, pH, and organic matter. For example, the authors found that farmer-identified good sites contained significantly greater levels of B, Ca, Mg, and higher pH, compared to poor sites.

In order to understand the causes behind the occurrence of high- and low-yielding field sites (=within-farm-spatial-variation of soybean yields) across, we conducted a soil survey by using 14 farms within the Pennsylvania soybean on-farm network (Bandara et al., 2020b). Bulk soil samples were collected and a number of soil chemical, physical, topological, and biological factors were assessed. Samples from two site types did not exhibit significant differences based on the assessed factors. Our findings indicated that the causes behind within-farm-spatial-variation of soybean yields across selected farms are complex and investigations using additional variables such as soil and plant-associated microbial composition may be useful in understanding within-site variability in yield.

Although the potential association of fungal and prokaryotic composition with within-site spatial variability in crop yield has not been previously demonstrated, differences in fungal and prokaryotic composition have been reported as a function of space-related variables. For example, Srour et al. (2017) found significant differences in prokaryotic and fungal composition between soils collected from healthy areas and areas where sudden death syndrome symptoms were visible. Kasel et al. (2008) observed that sites differing in land use within the same location also differed in soil fungal composition. Furthermore, substantial differences in microbial composition have been observed in relation to different agronomic/cultural practices, such as tillage (Yin et al., 2010; Dorr de Quadros et al., 2012; Mathew et al., 2012), crop rotation (Yin et al., 2010; Dorr de Quadros et al., 2012; Mathew et al., 2012), pesticide/herbicide use (Crouzet et al., 2010; Lo, 2010; Banks et al., 2014), and organic/conventional agriculture (Moeskops et al., 2010; Sugiyama et al., 2010; Schmid et al., 2011).

The composition of microorganisms in the rhizosphere, which is the region of the soil in contact with plant roots, can affect plant performance. If host plants can capitalize on the microbial services presented by the microbiome, agricultural productivity could be augmented through more fully harnessing beneficial functions. Microbial services may include the production of phytohormones (Ping and Boland, 2004), improved tolerance to abiotic and biotic stresses (Redman et al., 2002), induction of the plant innate immune response mechanisms (Jain et al., 2011), delivery of nutrients (Janos, 2007), modification of plant functional traits (Friesen et al., 2011), or tissue chemistry (Larsen et al., 2006). All of these services are directly related to microbial composition. Therefore, spatial differences in microbial structure could potentially contribute to differential plant performance, manifested as yield variation either at a narrow geographic scale (e.g., site-to-site variation within the farm) or at a broader geographic scale (e.g., farm-to-farm variation). In this study, our objective was to investigate the association between site-to-site spatial variation of soybean yield within a farm (hereafter referred to as within-farm-spatial-variation of soybean yields) and fungal/prokaryotic composition in bulk soil, rhizosphere soil, and soybean roots using an amplicon-sequencing approach.

 Materials and Methods Experimental Locations, Sample Collection, and Sample Processing

Eight farms located in eight Pennsylvania counties (Cambria, Bedford, Bucks, Butler, Lancaster, Mercer, Northumberland, Tioga) within the “Pennsylvania Soybean On Farm Network” were randomly chosen for this study. Figure 1 summarizes the sampling strategy that was used for sample collection. At each farm, farmer-designated high- and low-yield sites were defined and the Global Positioning System (GPS) coordinates were recorded, as well as the sampling date (Supplementary Table 1). Within each site type, five soybean plants were randomly selected and surface debris around the plant was removed. The aboveground plant organs were aseptically removed. Using a shovel, plants were removed with the root ball intact. The size of each root ball was ∼15–20 cm3. Entire root balls were transported to the Esker Lab at The Pennsylvania State University on dry ice on the same day that they were sampled. Sampling was performed at two soybean growth stages (SGS) (V1 = one set of unfolded trifoliate leaf is visible, and R8 = 95% of the pods have reached their mature color). These two growth stages were selected to represent early vegetative and late reproductive stages of the soybean lifecycle.

Strategy adopted for sample collection. Bulk soil, rhizosphere soil, and roots were collected from one farmer designated high and one low yield site at V1 (one set of unfolded trifoliate leaf is visible) and R8 (95% of the pods have reached their mature color) soybean growth stages from eight soybean farms in Pennsylvania.

Each root ball was fractionated into bulk soil, rhizosphere soil, and roots. For bulk soil, loose soil was manually removed from the root ball by kneading and shaking the soil with sterile gloves on to a sterile work surface. Accumulated soil was homogenized manually and approximately 20 g was placed into a 50 ml sterile Falcon tube. Subsequently, the left-over soil around the root system was removed by patting roots with a sterile spatula to reach rhizosphere soil. Soil aggregates that extended up to 1 mm from the root surface were carefully collected into a 50 ml Falcon tube (approximately 5–10 g) and these were designated as rhizosphere soil. From the remaining root system, secondary roots were excised using a sterile scalpel and placed in a sterile 50 ml Falcon tube containing ice cold 25 ml phosphate buffered saline (doi: 10.1101/pdb.rec8247 Cold Spring Harb. Protoc. 2006). Only secondary roots were collected across all root balls for consistency. Tubes were vortexed at maximum speed for 20 s, which released most of the leftover soil particles from the roots and turned the buffer turbid. The turbid buffer was then decanted and filled with fresh, ice cold PBS (25 ml) and vortexed as indicated above. The buffer was decanted, and this process was repeated five times, although by the end of the fourth cycle, the PBS solution was typically clear, meaning that there were no visible soil particles in the wash buffer. After the last cycle, the remaining buffer was removed from root sample using Kimwipes. Roots were then placed in a new sterile 50 ml Falcon tube. Tubes containing bulk soil, rhizosphere soil, and roots were stored at −80°C.

 DNA Extraction, PCR Amplification, and High-Throughput Sequencing

All samples were lyophilized for 48 h prior to DNA extraction. Lyophilization is particularly important for efficient DNA extraction from root samples. However, for the sake of consistency, soil samples were also subjected to lyophilization. DNA was extracted from 0.25 to 0.30 g of bulk soil, rhizosphere soil, and roots using the Macherey-Nagel NucleoSpin® 96 Soil DNA Isolation Kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany).

Initial PCR reactions were performed targeting the ITS1 region of the ribosomal gene using universal primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and 58A2R (5′-CTGC GTTCTTCATCGAT-3′) (Gardes and Bruns, 1993; Martin and Rygiewicz, 2005). As described by Bell et al. (2016), both primers were modified to include overhangs needed for index attachment (ITS1F-Illu = TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCT TGGTCAT TTAGAGGAAGTAA; 58A2R-Illu = GTCTCGTGG GCTCGGAGATGTGTATAAGAGACAGCTGCGTTCTTCATC GAT). Initial 16S rRNA gene PCR reactions were performed according to the 16S Metagenomic Sequencing Library Preparation guide (part no. 15044223 rev. B) with some modifications (Howard et al., 2017). The universal prokaryotic primers 515F (5′-GTGYCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACNVGGGTWTCTAAT-3′) were used for 16S rRNA gene amplifications (Herlemann et al., 2011). Similar to the ITS primers, both 16S primers were also modified to include overhangs needed for index attachment (515F-Illu = TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGT GYCAGCMGCCGCGGTAA; 806R-Illu = GTCTCGTGGGCT CGGAGATGTGTATAAGAGACAGGGACTACNVGGGTWTC TAAT).

PCR reactions occurred in 30 μl volumes comprised of 0.3 μl (=1.5 U) of HotMaster Taq DNA Polymerase (QuantaBio, Beverly, MA), 3 μl of 10x Taq Buffer with 25 mM Mg2+ (QuantaBio, Beverly, MA), 1.5 μl of each primer at 10 μM, 0.6 μl of 10 mM dNTP mix (New England Biolabs Inc., Ipswich, MA), 20.6 μl of H2O, and 2.5 μl of template DNA. PCR cycling for fungal ITS amplicons was performed using the following protocol: 94°C for 3 min; 35 cycles of 94°C for 20 s, 45°C for 30 s and 72°C for 45 s; with a final elongation at 72°C for 5 min. For 16S rRNA gene amplifications, the reaction mix was the same as described above. PCR cycling for 16S rRNA gene amplicons was performed using the following protocol: 94°C for 2 min; 25 cycles of 94°C for 20 s, 55°C for 20 s, and 72°C for 30 s; and 72°C for 5 min for the final elongation.

PCR products were run on E-Gel™ 96 Gels with SYBR™ Safe DNA Gel Stain (Thermo Fisher Scientific, Waltham, MA) with E-Gel™ Low Range Quantitative DNA Ladder (Thermo Fisher Scientific, Waltham, MA) to visualize/confirm the PCR products/bands (fungal ITS1 ∼ 380 bp, prokaryotic 16S rRNA ∼ 420 bp).

The initial amplicon cleanup was carried out in clear 96-well plates using the Mag-Bind® TotalPure NGS Kit (Omega Bio-Tek Inc., Norcross, GA). Cleaned PCR amplicons were then subjected to a second round of PCR to attach Illumina Nextera-compatible barcode and adaptors to the amplicons. PCR reactions occurred in 25 μl volumes in 96-well plates and were comprised of 5 μl of PCR product (from the first round of PCR), 2.5 μl of forward and reverse primers (i5 and i7, each at 10 μM) with designated barcodes for overhang attachment, 0.3125 μl (=1.5625 U) of HotMaster Taq DNA Polymerase (QuantaBio, Beverly, MA), 3.125 μl of 10x Taq Buffer with 25 mM Mg2+ (QuantaBio, Beverly, MA), 0.625 μl of 10 mM dNTP mix (New England Biolabs Inc., Ipswich, MA), and 10.9375 μl of water. Conditions for index attachment PCR cycling were as follows: 98°C for 1 min; 8 cycles of 98°C for 15 s, 55°C for 30 s, and 72°C for 20 s; and 72°C for 5 min for final elongation.

All barcoded amplicons were normalized (∼25 ng) using the SequalPrep Normalization Kit (Thermo Fisher Scientific, Waltham, MA). From each normalized sample, 20 μl was combined into separate pools for fungal ITS fragments and 16S rRNA genes, concentrated using a SpeedVac, and resuspended in 90 μl of molecular grade water. From these suspensions, 45 μl was decanted, mixed with 15 μl loading dye, and run separately on 1.2% agarose gels (100 ml) precast with 3 μl of SYBR™ Safe DNA Gel Stain (Thermo Fisher Scientific, Waltham, MA). Target bands were cut out from the gel and cleaned using the PureLink™ Quick Gel Extraction Kit (Thermo Fisher Scientific, Waltham, MA), for a final pool volume of 30 μl. 16S rRNA gene and ITS pools were separately sequenced at the Cornell Genomics Facility (Ithaca, NY) on the Illumina MiSeq. A 500-cycle MiSeq Reagent Kit v.2 was used for both fungal ITS pool and prokaryotic 16S rRNA pool.

 Sequence Processing and Analysis

Initial sequence processing was conducted using DADA2 v.1.18 pipeline (Callahan et al., 2016). The taxonomic assignment was performed using the UNITE General FASTA release v.8.3 database (Kõljalg et al., 2005) for ITS gene sequences and Silva version 138.1 (Quast et al., 2012) database for prokaryotic 16S rRNA gene sequences.

All global analyses, including co-occurrence networks, linkage clustering, beta diversity, alpha diversity, and relative abundance were performed in R v. 3.5.3. All global analyses were conducted using rarefied read counts (based on the minimum available reads per sample among all samples). Rarefying helped in obtaining even numbers of reads by sample, thus normalizing inter-sample comparisons. In the case of ITS sequences, bulk soil, rhizosphere soil, and roots were rarefied to 11,113, 7,663, and 16,565 reads. With respect to 16S rRNA gene sequences, bulk soil, rhizosphere soil, and roots were rarefied to 3,754, 3,993, and 6,078 reads.

The co-occurrence networks containing both Fungal and prokaryotic taxa were constructed and analyzed using the SpiecEasi and Igraph packages in R (Csardi and Nepusz, 2006; Kurtz et al., 2015). Networks were constructed with ASVs that were present in 50% of samples or more. Network sparsity and stability were examined using SpiecEasi. The “spiec.easi” function was executed with the following specifications to construct all the networks: (i) method = “mb,” (ii) lambda.min.ratio = 1e–2, (iii) nlambda = 50, (iv) sel.criterion = “stars,” (v) pulsar.select = TRUE, (vi) rep.num = 99. As outlined in Agler et al. (2016), hub taxa were identified as those above the 90th percentile (1.3 standard deviations from the mean) of network ASVs for the measures of (i) betweenness centrality, (ii) hub scores (eigenvector centrality), and (iii) degree, for both Fungi and prokaryotes in that specific network. Following network construction in SpiecEasi and hub identification, networks were visualized with the ggnet2 function of the Ggally package of R (Schloerke et al., 2018). In order to test whether the layouts of the experimental networks were consistently and significantly different from stochastically created, scale-free networks with the same number of nodes as experimental networks, 100 stochastic networks were generated using Barbasi-Albert model of the “sample_pa” function in the igraph package of R. The degree distributions of random networks were also compared to those of experimental networks using the non-parametric two sample Kolmogorov-Smirnov using the “ks.test” function in the stats package of R. The “sample_pa” and “ks.test” functions have been previously used in a number of articles including (Longley et al., 2020).

The Bray-Curtis dissimilarities between samples were used to perform principal coordinates analyses (PCoA) to investigate β-diversity (OTU diversity between samples). Permutational multivariate analysis of variance (PERMANOVA) was performed using Bray-Curtis dissimilarity to identify the factors/explanatory variables (=Site type, soybean growth stage) that significantly contributed to total observed variation in PCoA plots. Significance was assessed from 999 permutations. PERMANOVA was performed separately for bulk soil, rhizosphere soil, and roots.

Alpha diversity was estimated for each sample using the Inverse Simpson index and Chao1 richness within the vegan R package. Subsequently, analysis of variance (ANOVA) was performed using generalized linear mixed model (GLMM) approach in SAS version 9.4 (SAS® Institute, 2017) to test the main/simple effects (α = 0.05) of site type, soybean growth stage, and sample type (fixed factors) on two alpha diversity measures. Location was considered as a random factor. As both measures were count data, modeling was performed with the negative binomial distribution and using integral approximations to the likelihood. The following specifications were used for modeling: (1) link function = Log; (2) variance component estimation method = Maximum Likelihood; (3) degrees of freedom method = Residual; (4) non-linear parameter optimization method = Newton-Raphson with Ridging; (5) overdispersion fixation method = Laplace. The inverse link function was used to create means and associated standard errors at the data scale.

Relative abundance of taxonomic groups (phylum and genus level) was determined using the transform_sample_counts () function in the R package phyloseq. Heatmaps and dendrograms based on average linkage clustering of ASV relative abundance were generated with the heatmap() function in R package Heatplus. For both fungal and prokaryote analyses, ASVs with >2% relative abundance were included in the analysis.

 Results Structure of Microbial Co-occurrence Networks

Microbial networks constructed for three sample types from high and low yield sites at two SGS differed in their network statistics (Table 1). Analyses performed with three sample types showed greater network size and total/positive/negative edges for high yield sites compared to low yield sites at both V1 and R8 growth stages, with the exception of rhizosphere soil at V1 stage (Table 1 and Figures 2A–D). All networks contained a greater number of prokaryotic than fungal nodes, except for root networks at V1 stage from both site types. Overall, networks had a diverse mix of bacterial and fungal phyla. In case of bacteria, bulk soil, rhizosphere soil, and root networks were dominated by Proteobacteria, Actinobacteriota, and Bacteroidota (Figures 2A–D, 3A–D, 4A–D). Fungal nodes were primarily from Ascomycota, Basidiomycota, and Mortierellomycota (Figures 2A–D, 3A–D, 4A–D) for bulk and rhizosphere soil networks. For root networks, Ascomycota and Glomeromycota were predominant in V1 networks while Ascomycota and Basidiomycota were predominant in R8 networks. When compared to 100 stocastic networks, each network except the rhizosphere soil-high sites-R8 and roots-low sites-V1 networks consistently had a significantly (p < 0.05) different network layouts than 100 stochastic networks (Table 1).

Attributes describing the bipartite microbial networks for different sample types collected from high and low yield sites across eight Pennsylvania farms at V1 (one set of unfolded trifoliate leaf is visible) and R8 (95% of the pods have reached their mature color) soybean growth stages.

Network # Nodes
 # Edges
 # Hubs
 #Modules Stability Sparsity P-value* P-value range±
Total Fungal Prokaryotic Total Positive Negative Total Fungal Prokaryotic
Bulk soil: low sites at V1 374 122 252 1,644 950 694 15 7 8 26 0.0494 0.0235 0.0019 9.86E-06–0.0140
Bulk soil: high sites at V1 419 123 296 2,066 1,119 947 24 8 16 22 0.0496 0.0235 0.0237 5.35E-05–0.0470
Bulk soil: low sites at R8 212 81 131 464 272 192 11 6 5 20 0.0479 0.0206 0.0241 0.0042–0.0330
Bulk soil: high sites at R8 313 96 217 1,159 648 511 14 5 9 23 0.0496 0.0237 0.0193 7.85E-04–0.0428
Rhizo. soil: low sites at V1 204 90 114 464 273 191 5 3 2 31 0.0480 0.0223 0.0048 0.0027–0.0226
Rhizo. soil: high sites at V1 193 81 112 385 214 171 8 6 2 28 0.0490 0.0207 0.0116 0.0077–0.0461
Rhizo. soil: low sites at R8 241 79 162 678 366 312 13 6 7 24 0.0495 0.0233 0.0006 0.0136–0.0365
Rhizo. soil: high sites at R8 274 82 192 890 488 402 12 2 10 28 0.0493 0.0237 0.0049 0.0434–0.2397
Roots: low sites at V1 96 53 43 52 34 18 3 0 3 46 0.0464 0.0113 0.1066 0.0174–0.1809
Roots: high sites at V1 110 57 53 56 35 21 2 1 1 54 0.0440 0.0093 0.0351 0.0036–0.0333
Roots: low sites at R8 163 35 128 264 155 109 5 1 4 18 0.0488 0.0199 0.1148 0.0096–0.0389
Roots: high sites at R8 177 37 140 293 171 122 9 3 6 26 0.0470 0.0187 0.0046 0.0054–0.0438

*P-value indicates the probability value obtained from two sample Kolmogorov-Smirnov tests to assess the significance of the node degree distribution between random networks and the of experimental network of interest. ±P-value range represent the range of computed p-values when comparing the experimental network against 100 random networks. All random networks were created with the same number of nodes as experimental networks.

Network of taxon associations for bulk soil samples. Co-occurrence of fungal and prokaryotic taxa for high yield sites at V1 soybean growth stage (A), low yield sites at V1 soybean growth stage (B), high yield sites at R8 soybean growth stage (C), and low yield sites at R8 soybean growth stage (D). Nodes represent exact sequence variants (ESVs). Colors represent specific phyla. Shapes represent specific domains. Red solid lines (edges/links) connecting nodes indicate statistically significant negative correlations, and solid blue lines, positive correlations between the connected taxa. Node sizes are proportional to the degree (=node connectivity). Nodes with degree ≥5 were included in the network for the purpose of visualization clarity.

Network of taxon associations for rhizosphere soil samples. Co-occurrence of fungal and prokaryotic taxa for high yield sites at V1 soybean growth stage (A), low yield sites at V1 soybean growth stage (B), high yield sites at R8 soybean growth stage (C), and low yield sites at R8 soybean growth stage (D). Nodes represent exact sequence variants (ESVs). Colors represent specific phyla. Shapes represent specific domains. Red solid lines (edges/links) connecting nodes indicate statistically significant negative correlations, and solid blue lines, positive correlations between the connected taxa. Node sizes are proportional to the degree (=node connectivity). Nodes with degree ≥5 were included in the network for the purpose of visualization clarity.

Network of taxon associations for root samples. Co-occurrence of fungal and prokaryotic taxa for high yield sites at V1 soybean growth stage (A), low yield sites at V1 soybean growth stage (B), high yield sites at R8 soybean growth stage (C), and low yield sites at R8 soybean growth stage (D). Nodes represent exact sequence variants (ESVs). Colors represent specific phyla. Shapes represent specific domains. Red solid lines (edges/links) connecting nodes indicate statistically significant negative correlations, and solid blue lines, positive correlations between the connected taxa. Node sizes are proportional to the degree (=node connectivity).

The co-occurrence analyses for bulk soils showed a reduction in network size and number of edges from V1 to R8 SGS for both site types (Table 1 and Figures 2A–D). However, for rhizosphere soil and roots from both site types, the network size as well as the number of edges increased from V1 to R8 (Figures 3A–D, 4A–D). Out of 118 hub ASVs detected in 12 networks (two site types, three sample types, and two SGS), 89 were restricted to a single network (Supplementary Table 1). Certain hub species were present in multiple networks. For example, the fungus Corynespora cassiicola was a hub in networks created for both site types with rhizosphere soil and roots from R8 soybean growth stage. Most prokaryotic hubs consisted of Proteobacteria and Actinobacteriota while fungal hubs were mainly comprised of Ascomycota and Basidiomycota (Supplementary Table 1). For bulk soil, high yield site network showed a greater number of hubs compared to that of low yield site network at both growth stages (Table 1). For, rhizosphere soil, although high yield site network showed a greater number of hubs compared to that of low yield site network at V1 growth stage, the opposite was observed at R8 growth stage. For, roots, high yield site network showed a lower number of hubs compared to that of low yield site network at V1 growth stage. However, the opposite was observed at R8 growth stage. Bradyrhizobium elkanii, which is the most common symbiotic nitrogen fixer associated with soybean in the continental United States, was present in networks created for all sample types (bulk/rhizosphere soil and roots). However, it was not detected as a hub in any of the networks (Supplementary Table 1).

 Beta-Diversity Analysis of Fungal and Prokaryotic Composition

Principal coordinates analyses (PCoA) of fungal data using Bray-Curtis dissimilarity values showed that the percent variation explained by the first two principal coordinates was greatest for root samples (33.8%) followed by bulk soil (31.4%) and rhizosphere soil samples (30.9%) (Figures 5A–C). Neither bulk soil nor rhizosphere soil samples were appeared to cluster based on site type (high/low yield) or soybean growth stage (Figures 5A,B). However, a clear grouping was observed among root samples based on soybean growth stage but not based on site type (Figure 5C). PERMANOVA showed non-significant effect of site type on the observed total variation for all sample types (bulk soil: R2 = 0.03, P = 0.6753; rhizosphere soil: R2 = 0.02, P = 0.9540; roots: R2 = 0.02, P = 0.8621). The effect of soybean growth stage on the observed total variation was non-significant for bulk soil samples (R2 = 0.03, P = 0.7173) and rhizosphere soil samples (R2 = 0.04, P = 0.1459). However, a significant growth stage effect was evident for root samples (R2 = 0.16, P < 0.0001).

Principal Coordinates Analysis (PCoA) of ASVs based on Bray-Curtis dissimilarity for ITS-bulk soil (A), ITS-rhizosphere soil (B), ITS-roots (C), 16S rRNA-bulk soil (D), 16S rRNA-rhizosphere soil (E), and 16S rRNA-roots (F). Marker shape and color indicate site type (high/low yield) and soybean growth stage (V1 = one set of unfolded trifoliate leaf is visible; R8 = 95% of the pods have reached their mature color), respectively.

For prokaryotic data, the percent variation explained by the first two principal coordinates was similar for bulk and rhizosphere soil samples (17.5 and 20.4%, respectively), while it was greater for root samples (56.1%) (Figures 5D–F). Bulk and rhizosphere soil samples did not cluster based on site type or soybean growth stage (Figures 5D,E). A clear grouping was observed among root samples based only on soybean growth stage but not based on site type (Figure 5F). PERMANOVA revealed non-significant effect of site type on the observed total variation for all sample types (bulk soil: R2 = 0.03, P = 0.7952; rhizosphere soil: R2 = 0.03, P = 0.8591; roots: R2 = 0.02, P = 0.8741). The effect of soybean growth stage on the observed total variation was non-significant for bulk soil samples (R2 = 0.04, P = 0.0769). However, a significant growth stage effect was evident for rhizosphere soil (R2 = 0.06, P = 0.0020) and root samples (R2 = 0.36, P < 0.001).

 Heatmap Analysis With Average Linkage Clustering

Based on ASVs >2% abundance criterion, 103, 115, and 115 fungal ASVs were included in the heatmap analysis for bulk soil, rhizosphere soil, and roots, respectively. Similarly, 34, 65, and 63 prokaryotic ASVs were included in the analysis for bulk soil, rhizosphere soil, and roots, respectively. Analysis with fungal as well as prokaryotic ASVs did not cluster bulk soil samples into two distinct groups based on either site type or soybean growth stage (V1 and R8) (Figures 6A–F). Similar results were observed for rhizosphere soil samples. Although not distinctly grouped into two groups, root samples showed some degree of clustering based on soybean growth stage with both fungal and prokaryotic ASV data (Figures 6C,F).

The heat map analysis of ASVs based on average linkage clustering for ITS-bulk soil (A), ITS-rhizosphere soil (B), ITS-roots (C), 16S rRNA-bulk soil (D), 16S rRNA-rhizosphere soil (E), and 16S rRNA-roots (F). The purple, red, green, and blue color boxes indicate the low yield site, high yield site, V1 soybean growth stage (one set of unfolded trifoliate leaf is visible), and R8 soybean growth stage (95% of the pods have reached their mature color), respectively. Note that, based on both ITS and 16S rRNA OTU data sets, there is no clear clustering of samples (bulk soil, rhizosphere soil, and soybean roots) based on either site type or soybean growth stage.

 Alpha-Diversity Analysis of Fungal and Prokaryotic Composition

Here we measured diversity within samples (i.e., how different are ASVs within samples) using Inverse Simpson index, which accounts for both the number of ASVs observed and evenness in the relative abundance of ASVs, where higher index values imply greater diversity and vice versa. We also used Chao1 index, which is a measure of the number of different ASVs present in a given sample. For fungal ASVs, ANOVA showed a significant main effect of soybean growth stage (P = 0.0016) and sample type × site type interaction effect (P = 0.0412) on mean Inverse Simpson index. The mean Inverse Simpson index for fungal ASVs was significantly greater (P = 0.0016) at V1 growth stage (9.94 ± 1.12) compared to that of at R8 stage (13.27 ± 1.45). Although the mean Inverse Simpson index for rhizosphere soil and roots from two site types did not significantly differ, bulk soil from high yield sites showed a significantly greater index than that from the low yield sites (Figure 7A). Per the ANOVA, sample type × growth stage interaction effect was significant (P < 0.0001) on the mean Chao1 richness. At V1 growth stage, Chao1 richness was greater in bulk soil compared to rhizosphere soil (P < 0.0001) and roots (P < 0.0001) while the Chao1 richness of rhizosphere soil was significantly greater than roots (P < 0.0001) (Figure 7C). At R8 growth stage, the Chao1 richness of bulk soil (P < 0.0001) and rhizosphere soil (P < 0.0001) was greater compared to that of roots (Figure 7C).

Bar graphs depicting the impact of studied variables on mean α-diversity measures. Designation of the significant mean difference in panel (A) was based on Fisher’s least significant difference (LSD). Same for panels (B,C) were based on the adjustment for multiple comparisons using Tukey-Kramer test at the 5% level of significance (=5% experiment wise error rate). Error bars represent standard errors. Means followed by a common letter within each letter type (i.e., uppercase, lowercase, and italic) are not significantly different. V1 = one set of unfolded trifoliate leaf is visible; R8 = 95% of the pods have reached their mature color.

For prokaryotic ASVs, the sample type × growth stage interaction effect was significant (P < 0.0001) on Inverse Simpson index. At both V1 and R8 growth stages, the Inverse Simpson index was greater in bulk soil compared to rhizosphere soil (P < 0.0001) and roots (P < 0.0001) while the same of rhizosphere soil was significantly greater than roots (P < 0.0001) (Figure 7B). The sample type × growth stage interaction effect was also significant (P < 0.0001) on the mean Chao1 richness. At V1 growth stage, Chao1 richness was greater in bulk soil compared to rhizosphere soil (P < 0.0001) and roots (P < 0.0001) while the same of rhizosphere soil was significantly greater than roots (P < 0.0001) (Figure 7C). At R8 growth stage, the Chao1 richness of bulk soil (P = 0.0002) was greater than that of roots (Figure 7C).

 Relative Abundance of Fungal and Prokaryotic Taxa

Phylum-level fungal abundance analysis showed that Ascomycota was the predominant phylum in all sample types at both growth stages and site types (Figures 8A–C). For bulk and rhizosphere soils, the relative abundance of Ascomycetes appeared to be relatively similar at both SGS for both site types; nonetheless, it tended to increase in roots from V1 to R8 at both site types. Further, the Ascomycota abundance in roots from high yield sites was greater than that of low yield sites at both growth stages. On the contrary, the abundance of Mortierellomycota and Basidiomycota in roots from high yield sites was lower than that of low yield sites at both growth stages. Although Zoopagomycota was prominent in both soil types, it was less prominent in roots. Glomeromycota was very prominent in roots, particularly at V1 growth stage.

Stacked bar plots showing the mean relative abundance of major fungal phyla for (A) bulk soil, (B) rhizosphere soil, and (C) soybean root samples collected from high and low yield sites at V1 (one set of unfolded trifoliate leaf is visible) and R8 (95% of the pods have reached their mature color) soybean growth stages across eight soybean farms in Pennsylvania.

Genus level fungal abundance analysis showed that the relative abundance of Apodus, Cladosporium, Fusarium, Podospora, and Septoria, was greater in bulk soil from low yield sites compared to high yield sites at both SGS (Figure 9A). Abundance of Clonostachys, Corynespora, Exophiala, Humicola, Lophotrichus, Metarhizium, Solicoccozyma, and Trichoderma was greater in bulk soil from high yield sites compared to low yield sites at both growth stages.

Stacked bar plots showing the mean relative abundance of major fungal genera for (A) bulk soil, (B) rhizosphere soil, and (C) soybean root samples collected from high and low yield sites at V1 (one set of unfolded trifoliate leaf is visible) and R8 (95% of the pods have reached their mature color) soybean growth stages across eight soybean farms in Pennsylvania.

For rhizosphere soil, low yield sites contained greater abundance of Cladosporium, Fusarium, Mrakia, Podospora, and Septoria than high yield sites at both growth stages (Figure 9B). The abundance of Conioscypha, Corynespora, Metarhizium, Preussia, and Solicoccozyma was greater in high yield sites than low yield sites at both growth stages.

Soybean roots from low yield sites had greater abundance of Ilyonectria, Macrophomina, Mortierella, and Septoglomus compared to those in high yield sites at both growth stages while the abundance of Fusarium and Rhizophagus was greater at R8 (Figure 9C). Acrocalymma, Corynespora, and Podospora were more abundant in roots from high yield sites than low yield sites at both growth stages, while Cercophora, Clonostachys, Diaporthe, and Hymenoscyphus were present in greater abundance at R8. For high yield sites, abundance of root associated Corynespora increased from V1 to R8, but it decreased in roots from low yield sites.

Phylum-level prokaryotic abundance analysis showed that Proteobacteria was the predominant phylum in all sample types at both growth stages and site types (Figures 10A–C). Actinobacteriota and Bacteroidota were also abundant in all sample types. Although Acidobacteriota, Chloroflexi, Crenarchaeota, Gemmatimonadota, Planctomycetota, and Verrucomicrobiota were abundant in both soil types (across growth stages and site types), they were less conspicuous in roots. Cyanobacteria was prominent in roots, particularly at V1 stage.

Stacked bar plots showing the mean relative abundance of major prokaryotic phyla for (A) bulk soil, (B) rhizosphere soil, and (C) soybean root samples collected from high and low yield sites at V1 (one set of unfolded trifoliate leaf is visible) and R8 (95% of the pods have reached their mature color) soybean growth stages across eight soybean farms in Pennsylvania.

Genus-level prokaryotic abundance analysis showed that the relative abundance of major prokaryotic genera in both bulk and rhizosphere soils was similar between high and low yield sites at both growth stages (Figures 11A,B). Further, for majority of the genera, there was no marked difference between two site types at both growth stages. For roots, abundance of Niastella, Piscinibacter, Pseudomonas, and Rhizobacter from low yield sites was greater compared to high yield sites at both growth stages (Figure 11C). Abundance of Bradyrhizobium, Flavobacterium, Novosphingobium, and Rhodoferax was greater in roots from high yield sites compared to low yield sites at both SGS.

Stacked bar plots showing the mean relative abundance of major prokaryotic genera for (A) bulk soil, (B) rhizosphere soil, and (C) soybean root samples collected from high and low yield sites at V1 (one set of unfolded trifoliate leaf is visible) and R8 (95% of the pods have reached their mature color) soybean growth stages across eight soybean farms in Pennsylvania.

 Discussion

Microbial composition in soil and roots can influence plant physiological performance and yield. However, little is known about the association between the spatial variation of microbial composition and the corresponding spatial variation of crop yields at the site scale. Since many soil- and root-inhabiting fungal and prokaryotic taxa are uncultivable using conventional methods, a metabarcoding approach offers insights into populations only known by sequencing (Prosser, 2015). Our results, among the first using this approach to address such questions, showed that the abundance of key taxa and network structures of fungi and bacteria are different between high and low yield sites.

Network analysis is a powerful tool for investigating and recognizing patterns in large, complex datasets, which can be more challenging to detect using common methods in ecology like standard alpha/beta diversity analyses (Proulx et al., 2005). Our results showed greater network sizes and total/positive/negative edges for high yield sites compared to low yield sites at both growth stages for all sample types, except for the V1-rhizosphere soil network. These results indicated that microbial networks from high yield sites are more complex compared to those of low yield sites. A previous study by Zheng et al. (2021) showed that the healthy tobacco plants had more complex microbial networks than the plants with bacterial wilt caused by Ralstonia solanacearum. Although we did not visually observe any disease conditions in soybean plant from low yield sites, it is possible that less complex microbial networks contributed to less healthy plants that will eventually translate into poorer yields. More in-depth investigations are required to gain insights into mechanistic understanding on how greater network complexity and soybean yields.

Bradyrhizobium elkanii is reported as the major soybean nodulating rhizobia (Jordan, 1982). Although B. elkanii was present in all microbial co-occurrence networks constructed in the present study, it was not detected as a hub in any of the networks. This finding indicated that the potential interactions of B. elkanii with other prokaryotic/fungal members of the community is minimal, despite its importance in soybean cropping system. Poudel et al. (2016) discussed that the most important taxa for plant health may neither have any links to other taxa, nor have correlations with plant or pathogen performance components. Therefore, it was not surprising to see B. elkanii as a non-hub taxon in the current study.

Although the network association structure and related attributes tended to change depending on soybean growth stage, such changes did not appear to depend on site type. For instance, the network size and total/positive/negative edges increased from V1 to R8 for rhizosphere soil and roots for both high and low yield sites. Therefore, while there were temporally defined ecological rearrangements of microbial composition, as evidenced by the changes in topological attributes with the co-occurrence network, such changes appeared to be common between the two site types. As such, at least considering the seasonal scale, the temporal dynamics of microbial composition may not contribute to site-dependent spatial yield variation.

For both fungal and prokaryotic data, regardless of growth stage and site type, bulk soil showed the greatest Inverse Simpson index and Chao1 richness (measure of α-diversity), followed by rhizosphere soil and roots. As such, there was a diversity gradient from bulk soil to roots via the rhizosphere interface. Results were consistent with studies that have demonstrated the highest alpha diversity of both fungi and bacteria in the soil (Lebreton et al., 2019; Suárez-Moo et al., 2019).

As revealed by relative abundance analysis, fungal and prokaryotic composition varied based on site type (high/low), soybean growth stage (V1/R8), and sample type (bulk/rhizosphere soil, roots) at both phylum and genus level. Previous reports clearly demonstrated growth stage (Mougel et al., 2006; Houlden et al., 2008; Xu et al., 2009; Chaparro et al., 2014) and sample type (Sugiyama et al., 2014; Yurgel et al., 2017) dependent relative abundance variation of prokaryotic and fungal taxa.

Our findings showed a greater relative abundance of the fungal genus Fusarium in roots from low yield sites compared to high yield sites. Other than a single ASV that corresponds to an unclassified species of genus Fusarium, we detected F. acutatum, F. solani, F. cuneirostrum, and F. oxysporum in these root samples. Interestingly, all of these Fusarium species are known soybean pathogens (Daamen et al., 1991; Aoki et al., 2005; Zhang et al., 2010; Arias et al., 2013; Yan and Nelson, 2020; Degani and Kalman, 2021). We also found a greater relative abundance of the fungal species, Macrophomina phaseolina in roots from low yield sites compared to high yield sites. M. phaseolina is one of the most devastating pathogens in soybean causing charcoal rot disease (Mengistu et al., 2013; Bai et al., 2015; Bandara et al., 2020a). As such, soybean roots from low yield sites appeared to be under higher pathogen pressure and corresponding disease risks, although external disease symptoms such as rots, spots, and lesions were not visually observed during the sampling process.

Interestingly, the relative abundance of root-colonizing genus Corynespora (represented by a single species, C. cassiicola) was greater in high yield sites compared to low yield sites. This organism causes target spot on leaves, stems, roots and flowers of more than 280 plant species, including many economically important crops, such as soybean (Silva et al., 1995). However, Dixon et al. (2009) and Shimomoto et al. (2011) demonstrated greater pathogenicity variability among C. cassiicola isolates. It is possible that the pathogenicity of C. cassiicola from high yield sites is significantly lower, despite higher relative abundance. Alternatively, as an endophyte, C. cassiicola is reported to secrete an array of polyketide and fatty acid derivatives where their biological activities have not been tested (Chagas et al., 2015). These compounds can potentially be antifungal or antibacterial in nature. Given this possibility, greater relative abundance of C. cassiicola in roots from high yield sites could also be beneficial and contribute to greater yields. More investigations are essential to precisely determine the association between root colonizing C. cassiicola and soybean yields.

The bacterial genus, Pseudomonas was present in greater percentages, particularly in rhizosphere and root samples collected from low yield sites at both growth stages. Pseudomonads from rhizosphere soil was predominantly represented by P. umsongensis while P. chlororaphis and P. frederiksbergensis were the two species detected in roots. Most of the Pseudomonas species are pathogenic to plants (Young, 2010; Mansfield et al., 2012; Lamichhane et al., 2015). However certain Pseudomonas species such as P. fluorescens and P. protegens are useful as biocontrol agents (Ramette et al., 2011; Gull and Hafeez, 2012; Santoyo et al., 2012; Sivasakthi et al., 2014). Interestingly, the two species that we observed in greater quantities in roots from low yield sites are reported as antifungal and plant growth promoting (Bardas et al., 2009; Chatterjee et al., 2017). Although it seems counterintuitive to observe a greater number of beneficial Pseudomonads in roots from low yield sites compared to that of high yield sites, more in-depth studies could justifiably explain this phenomenon.

A noteworthy observation was the greater relative abundance of the genus Flavobacterium in roots from high yield sites compared to low yield sites at both SGS. Some Flavobacterium species are known as plant growth promoters (Sessitsch et al., 2004) while other species are reported to protect plants against microbial infections through root colonization (Sang and Kim, 2012) or bioremediation (Kang et al., 2013). Given these beneficial effects, greater Flavobacterium colonization in roots from high yield sites can in turn contribute to better physiological performance and greater yields.

The bulk soils from high yield sites at both SGS contained greater percentages of Trichoderma (predominantly T. spirale and an unclassified Trichoderma species) fungal genera. The beneficial mycoparasitic and other roles of various Trichoderma species are widely reported in relation to soybean production (Menendez and Godeas, 1998; Shovan et al., 2008; Bagwan, 2010; John et al., 2010; Khalili et al., 2016; Khaledi and Taheri, 2016; Zhang et al., 2016). Certain Trichoderma species can also control phytopathogenic nematodes like soybean root knot nematodes (Oyekanmi et al., 2007). In particular, the utility of T. spirale as a biocontrol agent against C. cassiicola is previously reported (Baiyee et al., 2019). Note that, as discussed above, the abundance of C. cassiicola was greater in high yield sites. Even if we assume that C. cassiicola isolates from high yield sites are pathogenic, the greater abundance of T. spirale in high yield sites can potentially suppress C. cassiicola and mitigate its negative impacts on soybean yields. Moreover, rhizosphere soils from high yield sites at both SGS contained greater percentages of genus Metarhizium (M. anisopliae and M. marquandii) compared to low yield sites. The beneficial entomopathogenic behavior of Metarhizium including M. anisopliae and M. marquandii is widely documented in relation to soybean production (Clifton et al., 2018; Baron et al., 2020; Lopes et al., 2020). Therefore, the greater occurrence of beneficial fungal genera such as Trichoderma and Metarhizium can potentially eases the soilborne pathogen and insect pressure. This could eventually translate into greater yields in high yield sites.

While the relative abundance of symbiotic nitrogen fixers like Bradyrhizobium elkanii was greater in roots from high yield sites at both V1 and R8 growth stages, free-living nitrogen fixers such as Klebsiella and Kosakonia were more abundant in roots from low yield sites compared to high yield sites. Symbiotic nitrogen fixation occurs under a narrow range of environmental conditions while free-living nitrogen fixation can occur under a wide range of environmental conditions (Smercina et al., 2019). In contrast to symbiotic fixation, free-living nitrogen fixation can take place even under non-optimal conditions such as greater oxygen concentration and lower carbon (energy) availability (Smercina et al., 2019). The greater abundance of free-living nitrogen fixers in low yield sites is therefore indicative of a less efficient nitrogen fixation via symbiotic fixers. Nevertheless, these symbiotic nitrogen fixers directly depend on the host (soybean plant in this case) to fulfill their carbon requirements. Given these observations, we hypothesize that soybean plants from low yield sites lose more carbon due to symbiotic nitrogen fixers while not getting enough fixed nitrogen, which in turn contributes to the poor yields. Furthermore, natural environments are ubiquitously inhabited by unproductive rhizobia strains that extort benefits without compensating costs and thus proliferate more efficiently than nitrogen-fixing cooperators (Gibson et al., 1975; Bottomley and Jenkins, 1983; Moawad et al., 1998; Burdon et al., 1999; Fujita et al., 2014). In the current study, the genus Allorhizobium (comprised of A. cellulosilyticum, A. daejeonense, A. mesosinicum, A. phaseoli, and two unclassified Allorhizobium species) was present in greater abundance in roots from low yield sites compared to that of high yield sites. The reports that reveal their use as efficient Nitrogen fixers are sparse. Therefore, we further hypothesize that different Allorhizobium species detected in the current study could be parasitic. The greater presence of potentially parasitic Allorhizobia in roots from low yield sites could use surplus energy for their own growth or for creating storage substances and in turn results in lower plant productivity and yields. A comprehensive comparison of Allorhizobia strains isolated from high and low yield sites in terms of their ability to affect soybean yields could provide additional insights into the relationship between site scale spatial variation of soybean yields and symbiotic nitrogen fixers.

 Conclusion

Spatial variation of crop yields exists even within small geographic regions (i.e., within a farm) is a commonly observed phenomenon. Several factors can contribute to such variation including soil microbiota. Recent advances in high-throughput sequencing have allowed us to explore microbial composition in much greater detail, as we can identify even poorly known phyla with limited phenotypic data. Using an amplicon sequencing approach, we found that high yield sites contain significantly fewer fungal genera that contain known pathogenic species and more fungal genera that contain species that are known to be mycoparasitic or nematophagous. Moreover, the number of root-colonizing bacterial genera that contain known plant growth promoting bacteria was also greater in high yield sites. The greater occurrence of free-living nitrogen fixing bacterial genera in roots from low yield sites indicates that prevailing edaphic conditions in low yield sites were not ideal for symbiotic nitrogen fixation. Additionally, results from microbial network analysis showed that the size and the total number of edges in networks from high yield sites is greater with a larger number of specific biological interactions among them. Therefore, our findings showed that the occurrence of high and low yield spots in soybean fields was associated with a consortium of microbial taxa suggesting the importance of certain microbial taxa to the establishment and stability of healthy soil. Taken together, our findings provide new insights into the relationship between fungal and prokaryotic composition to observed site-scale spatial heterogeneity in soybean yield. These findings will be helpful in devising future site-specific management practices of soil and root-associated microorganisms toward better soil health and increased agricultural production.

 Data Availability Statement

The datasets generated for this study can be found in the NCBI SRA Archive: PRJNA707296.

 Author Contributions

PE and TB secured the funding and edited the manuscript. PE conceptualized and supervised the project. AB and PE designed the experiments. PE and DW coordinated the sample collection. AB and DW processed the samples. AB, DW, and RT prepared the amplicon libraries for sequencing. AB analyzed the data, created the figures, and wrote the draft manuscript. All authors approved the final manuscript.

 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.

 Funding

This project was supported by the Pennsylvania Soybean Board (Project R2018-P09). This project was also supported by the USDA National Institute of Food and Federal Appropriations under Projects PEN04660 and PEN04651 and Accession numbers 1016474 and 1016233. Funders had no role in the design, collection, and analysis of data.

We thank our farmer cooperators and the following extension educators (the Pennsylvania State University) who contributed to site identification and sample collection from the Pennsylvania Soybean On-Farm Network: Adriana Murillo-Williams, Andrew Frankenfield, Anna Busch, Claire Coombs, Del Voight, Elizabeth Bosak, Jeff Graybrill, Justin Brackenrich, Nicole Santangelo, and Zach Larson. We also thank Ravin Poudel for providing a piece of R script for the original network analysis using igraph package in the draft manuscript. Note that the published version does not contain any results or inferences/interpretations from the original network analysis.

 Supplementary Material

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

 References Agler M. T. Ruhe J. Kroll S. Morhenn C. Kim S. T. Weigel D. (2016). Microbial hub taxa link host and abiotic factors to plant microbiome variation. PLoS Biol. 14:e1002352. 10.1371/journal.pbio.1002352 26788878 Aoki T. O’Donnell K. Scandiani M. M. (2005). Sudden death syndrome of soybean in South America is caused by four species of Fusarium: Fusarium brasiliense sp. nov., F. cuneirostrum sp. nov., F. tucumaniae, and F. virguliforme. Mycoscience 46 162–183. 10.1007/s10267-005-0235-y Arias M. D. Munkvold G. P. Ellis M. L. Leandro L. F. S. (2013). Distribution and frequency of Fusarium species associated with soybean roots in Iowa. Plant Dis. 97 1557–1562. 10.1094/PDIS-11-12-1059-RE 30716864 Bagwan N. B. (2010). Evaluation of Trichoderma compatibility with fungicides, pesticides, organic cakes and botanicals for integrated management of soil borne diseases of soybean [Glycine max (L.) Merril]. Int. J. Plant Prot. 3 206–209. Bai L. Cui J. Jie W. Cai B. (2015). Analysis of the community compositions of rhizosphere fungi in soybeans continuous cropping fields. Microbiol. Res. 180 49–56. 10.1016/j.micres.2015.07.007 26505311 Baiyee B. Pornsuriya C. Ito S. I. Sunpapao A. (2019). Trichoderma spirale T76-1 displays biocontrol activity against leaf spot on lettuce (Lactuca sativa L.) caused by Corynespora cassiicola or Curvularia aeria. Biol. Control 129 195–200. 10.1016/j.biocontrol.2018.10.018 Ball S. T. Mulla D. J. Konzak C. F. (1993). Spatial heterogeneity affects variety trial interpretation. Crop Sci. 33 931–935. 10.2135/cropsci1993.0011183x003300050011x 34798789 Bandara A. Y. Weerasooriya D. K. Murillo-Williams A. White C. M. Collins A. A. Bell T. H. (2020b). Relationship between soybean yield from high and low yielding field sites and selected soil characteristics. Agrosyst. Geosci. Environ. 3:e20126. 10.1002/agg2.20126 Bandara A. Y. Weerasooriya D. K. Bradley C. A. Allen T. W. Esker P. D. (2020a). Dissecting the economic impact of soybean diseases in the United States over two decades. PLoS One 15:e0231141. 10.1371/journal.pone.0231141 32240251 Banks M. L. Kennedy A. C. Kremer R. J. Eivazi F. (2014). Soil microbial community response to surfactants and herbicides in two soils. Appl. Soil Ecol. 74 12–20. 10.1016/j.apsoil.2013.08.018 Bardas G. A. Lagopodi A. L. Kadoglidou K. Tzavella-Klonari K. (2009). Biological control of three Colletotrichum lindemuthianum races using Pseudomonas chlororaphis PCL1391 and Pseudomonas fluorescens WCS365. Biol. Control. 49 139–145. 10.1016/j.biocontrol.2009.01.012 Baron N. C. de Souza Pollo A. Rigobelo E. C. (2020). Purpureocillium lilacinum and Metarhizium marquandii as plant growth-promoting fungi. PeerJ 8:e9005. 10.7717/peerj.9005 32518715 Batchelor W. D. Basso B. Paz J. (2002). Examples of strategies to analyze spatial and temporal yield variability using crop models. Eur. J. Agron. 18 141–158. 10.1016/s1161-0301(02)00101-6 Becher H. H. (1995). On the importance of soil homogeneity when evaluating field trials. J. Agron. Crop Sci. 174 33–40. 10.1111/j.1439-037x.1995.tb00192.x Bell T. H. Stefani F. O. P. Abram K. Champagne J. Yergeau E. Hijri M. (2016). A diverse soil microbiome degrades more crude oil than specialized assemblages obtained in culture. Appl. Environ. Microbiol. 82 5530–5541. 10.1128/AEM.01327-16 27371586 Berndtsson R. Bahri A. (1995). Field variability of element concentrations in wheat and soil. Soil Sci. 159 311–320. 10.1097/00010694-199505000-00004 Bottomley P. Jenkins M. (1983). Some characteristic of Rhizobium meliloti isolates from alfalfa fields in Oregon. Soil Sci. Soc. Am. J. 47 1153–1157. 10.2136/sssaj1983.03615995004700060019x Burdon J. Gibson A. Searle S. Woods M. Brockwell J. (1999). Variation in the effectiveness of symbiotic associations between native rhizobia and temperate Australian Acacia: within-species interactions. J. Appl. Ecol. 36 398–408. 10.1046/j.1365-2664.1999.00409.x Callahan B. J. McMurdie P. J. Rosen M. J. Han A. W. Johnson A. J. A. Holmes S. P. (2016). DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13 581–583. 10.1038/nmeth.3869 27214047 Chagas F. O. Caraballo-Rodriguez A. M. Pupo M. T. (2015). “Endophytic fungi as a source of novel metabolites,” in Biosynthesis and Molecular Genetics of Fungal Secondary Metabolites, eds Zeilinger S. Martín J. F. García-Estrada C. (New York, NY: Springer), 123–176. 10.1007/978-1-4939-2531-5_8 Chaparro J. M. Badri D. V. Vivanco J. M. (2014). Rhizosphere microbiome assemblage is affected by plant development. ISME J. 8 790–803. 10.1038/ismej.2013.196 24196324 Chatterjee P. Samaddar S. Anandham R. Kang Y. Kim K. Selvakumar G. (2017). Beneficial soil bacterium Pseudomonas frederiksbergensis OS261 augments salt tolerance and promotes red pepper plant growth. Front. Plant Sci. 8:705. 10.3389/fpls.2017.00705 28523010 Clifton E. H. Jaronski S. T. Coates B. S. Hodgson E. W. Gassmann A. J. (2018). Effects of endophytic entomopathogenic fungi on soybean aphid and identification of Metarhizium isolates from agricultural fields. PLoS One 13:e0194815. 10.1371/journal.pone.0194815 29566067 Cox M. S. Gerard P. D. Wardlaw M. C. Abshire M. J. (2003). Variability of selected soil properties and their relationships with soybean yield. Soil Sci. Soc. Am. J. 67 1296–1302. 10.2136/sssaj2003.1296 Crouzet O. Batisson I. Besse-Hoggan P. Bonnemoy F. Bardot C. Poly F. (2010). Response of soil microbial communities to the herbicide mesotrione: a dose-effect microcosm approach. Soil Biol. Biochem. 42 193–202. 10.1007/s11356-015-4797-8 26122568 Csardi G. Nepusz T. (2006). The igraph software package for complex network research. Inter. J. Complex Syst. 1695 1–9. 10.1186/1471-2105-12-455 22115179 Daamen R. A. Langerak C. J. Stol W. (1991). Surveys of cereal diseases and pests in the Netherlands. 3. Monographella nivalis and Fusarium spp. in winter wheat fields and seed lots. Neth. J. Plant Pathol. 97 105–114. 10.1007/bf01974274 Degani O. Kalman B. (2021). Assessment of commercial fungicides against Onion (Allium cepa) basal rot disease caused by Fusarium oxysporum f. sp. cepae and Fusarium acutatum. J. Fungi 7:235. 10.3390/jof7030235 33801156 Dixon L. J. Schlub R. L. Pernezny K. Datnoff L. E. (2009). Host specialization and phylogenetic diversity of Corynespora cassiicola. Phytopathology 99 1015–1027. 10.1094/PHYTO-99-9-1015 19671003 Dorr de Quadros P. Zhalnina K. Davis-Richardson A. Fagen J. R. Drew J. Bayer C. (2012). The effect of tillage system and crop rotation on soil microbial diversity and composition in a subtropical acrisol. Diversity 4 375–395. 10.3390/d4040375 Friesen M. L. Porter S. S. Stark S. C. von Wettberg E. J. Sachs J.L. Martinez-Romero E. (2011). Microbially mediated plant functional traits. Annu. Rev. Ecol. Evol. Syst. 42, 23–46. 10.1146/annurev-ecolsys-102710-145039 Fujita H. Aoki S. Kawaguchi M. (2014). Evolutionary dynamics of nitrogen fixation in the legume–rhizobia symbiosis. PLoS One 9:e93670. 10.1371/journal.pone.0093670 24691447 Gardes M. Bruns T. D. (1993). ITS primers with enhanced specificity for Basidiomycetes-Application to the identification of mycorrhizae and rusts. Mol. Ecol. 2 113–118. 10.1111/j.1365-294x.1993.tb00005.x 8180733 Gibson A. Curnow B. Bergersen F. Brockwell J. Robinson A. (1975). Studies of field populations of Rhizobium: effectiveness of strains of Rhizobium trifolii associated with Trifolium subterraneum L. pastures in South-Eastern Australia. Soil Biol. Biochem. 7 95–102. 10.1016/0038-0717(75)90005-x Gull M. Hafeez F. Y. (2012). Characterization of siderophore producing bacterial strain Pseudomonas fluorescens Mst 8.2 as plant growth promoting and biocontrol agent in wheat. Afr. J. Microbiol. Res. 6 6308–6318. Herlemann D. P. R. Labrenz M. Jurgens K. Bertilsson S. Waniek J. J. Andersson A. F. (2011). Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea. ISME J. 5 1571–1579. 10.1038/ismej.2011.41 21472016 Houlden A. Timms-Wilson T. M. Day M. J. Bailey M. J. (2008). Influence of plant developmental stage on microbial community structure and activity in the rhizosphere of three field crops. FEMS Microbiol. Ecol. 65 193–201. 10.1111/j.1574-6941.2008.00535.x 18616582 Howard M. Bell T. H. Kao-Kniffin J. (2017). Soil microbiome transfer method affects microbiome composition, including dominant microorganisms, in a novel environment. FEMS Microbiol. Lett. 364:fnx092. 10.1093/femsle/fnx092 28472491 Jain A. Singh S. Sarma B. K. Singh H. B. (2011). Microbial consortium mediated reprogramming of defense network in pea to enhance tolerance against Sclerotinia sclerotiorum. J. Appl. Microbiol. 112 537–550. 10.1111/j.1365-2672.2011.05220.x 22181349 Janos D. P. (2007). Plant responsiveness to mycorrhizas differs from dependence upon mycorrhizas. Mycorrhiza 17 75–91. 10.1007/s00572-006-0094-1 17216499 Jaynes D. B. Colvin T. S. (1997). Spatiotemporal variability of corn and soybean yield. Agron. J. 89 30–37. 10.2134/agronj1997.00021962008900010005x Jiang P. Thelen K. D. (2004). Effect of soil and topographic properties on crop yield in a north-central corn-soybean cropping system. Agron. J. 96 252–258. 10.2134/agronj2004.0252 John R. P. Tyagi R. D. Prévost D. Brar S. K. Pouleur S. Surampalli R. Y. (2010). Mycoparasitic Trichoderma viride as a biocontrol agent against Fusarium oxysporum f. sp. adzuki and Pythium arrhenomanes and as a growth promoter of soybean. Crop Prot. 29 1452–1459. 10.1016/j.cropro.2010.08.004 Jordan D. C. (1982). Transfer of Rhizobium japonicum Buchanan 1980 to Bradyrhizobium gen. nov., a Genus of slow-growing, root nodule bacteria from Leguminous plants. Int. J. Syst. Bacteriol. 32 136–139. 10.1099/00207713-32-1-136 Kang J. Y. Chun J. Jahng K. Y. (2013). Flavobacterium aciduliphilum sp. nov., isolated from freshwater, and emended description of the genus Flavobacterium. Int. J. Syst. Evol. Microbiol. 63 1633–1638. 10.1099/ijs.0.044495-0 22904226 Kasel S. Bennett L. T. Tibbits J. (2008). Land use influences soil fungal community composition across central Victoria, south-eastern Australia. Soil Biol. Biochem. 40 1724–1732. 10.1016/j.soilbio.2008.02.011 Kaspar T. C. Pulido D. J. Fenton T. E. Colvin T. S. Karlen D. L. Jaynes D. B. (2004). Relationship of corn and soybean yield to soil and terrain properties. Agron. J. 96 700–709. Khakural B. R. Robert P. C. Mulla D. J. (1996). “Relating corn/soybean yield to variability in soil and landscape characteristics,” in Precision Agriculture, eds Robert P. C. Rust R. H. Larson W. E. (Madison, WI: ASA, CSSA, SSSA), 117–128. 10.2134/1996.precisionagproc3.c12 Khaledi N. Taheri P. (2016). Biocontrol mechanisms of Trichoderma harzianum against soybean charcoal rot caused by Macrophomina phaseolina. J. Plant Prot. Res. 56 21–31. 10.1515/jppr-2016-0004 Khalili E. Javed M. A. Huyop F. Rayatpanah S. Jamshidi S. Wahab R. A. (2016). Evaluation of Trichoderma isolates as potential biological control agent against soybean charcoal rot disease caused by Macrophomina phaseolina. Biotechnol. Biotechnol. Equip. 30 479–488. 10.1080/13102818.2016.1147334 Kõljalg U. Larsson K. H. Abarenkov K. Nilsson R. H. Alexander I. J. Eberhardt U. (2005). UNITE: a database providing web-based methods for the molecular identification of ectomycorrhizal fungi. New Phytol. 166 1063–1068. 10.1111/j.1469-8137.2005.01376.x 15869663 Kravchenko A. N. Bullock D. G. (2000). Correlation of corn and soybean grain yield with topography and soil properties. Agron. J. 92 75–83. 10.1007/s100870050010 Kurtz Z. D. Müller C. L. Miraldi E. R. Littman D. R. Blaser M. J. Bonneau R. A. (2015). Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput. Biol. 11:e1004226. 10.1371/journal.pcbi.1004226 25950956 Lamichhane J. R. Messéan A. Morris C. E. (2015). Insights into epidemiology and control of diseases of annual plants caused by the Pseudomonas syringae species complex. J. Gen. Plant Pathol. 81 331–350. 10.1007/s10327-015-0605-z Larsen E. H. Lobinski R. Burger-Meÿer K. Hansen M. Ruzik R. Mazurowska L. (2006). Uptake and speciation of selenium in garlic cultivated in soil amended with symbiotic fungi (mycorrhiza) and selenate. Anal. Bioanal. Chem. 385 1098–1108. 10.1007/s00216-006-0535-x 16770577 Lebreton L. Guillerm-Erckelboudt A. Y. Gazengel K. Linglin J. Ourry M. Glory P. (2019). Temporal dynamics of bacterial and fungal communities during the infection of Brassica rapa roots by the protist Plasmodiophora brassicae. PLoS One 14:e0204195. 10.1371/journal.pone.0204195 30802246 Lo C. C. (2010). Effect of pesticides on soil microbial community. J. Environ. Sci. Health Part B 45 348–359. 10.1080/03601231003799804 Longley R. Noel Z. A. Benucci G. M. N. Chilvers M. I. Trail F. Bonito G. (2020). Crop management impacts the soybean (Glycine max) microbiome. Front. Microbiol. 11:1116. 10.3389/fmicb.2020.01116 32582080 Lopes R. B. Sosa-Gómez D. R. Oliveira C. M. Sanches M. M. de Souza D. A. Benito N. P. (2020). Efficacy of an oil-based formulation combining Metarhizium rileyi and nucleopolyhedroviruses against lepidopteran pests of soybean. J. Appl. Entomol. 144 678–689. 10.1111/jen.12787 Mann C. Lynch D. Fillmore S. Mills A. (2019). Relationships between field management, soil health, and microbial community composition. Appl. Soil Ecol. 144 12–21. 10.1016/j.apsoil.2019.06.012 Mansfield J. Genin S. Magori S. Citovsky V. Sriariyanum M. Ronald P. (2012). Top 10 plant pathogenic bacteria in molecular plant pathology. Mol. Plant Pathol. 13 614–629. 10.1111/j.1364-3703.2012.00804.x 22672649 Martin K. J. Rygiewicz P. T. (2005). Fungal-specific PCR primers developed for analysis of the ITS region of environmental DNA extracts. BMC Microbiol. 5:28. 10.1186/1471-2180-5-28 15904497 Mathew R. P. Feng Y. Githinji L. Ankumah R. Balkcom K. S. (2012). Impact of no-tillage and conventional tillage systems on soil microbial communities. Appl. Environ Soil Sci. 2012 548620. Menendez A. B. Godeas A. (1998). Biological control of Sclerotinia sclerotiorum attacking soybean plants. Degradation of the cell walls of this pathogen by Trichoderma harzianum (BAFC 742). Mycopathologia 142 153–160. 10.1023/A:1006910707804 Mengistu A. Arelli P. Bond J. Nelson R. Rupe J. Shanno G. (2013). Identification of soybean accessions resistant to Macrophomina phaseolina by field screening and laboratory validation. Plant Health Prog. 14:25. 10.1094/php-2013-0318-01-rs Moawad H. El-Din S. Abdel-Aziz R. (1998). Improvement of biological nitrogen fixation in Egyptian winter legumes through better management of Rhizobium. Plant Soil 204 95–106. 10.1007/978-94-017-2321-3_9 Moebius-Clune B. N. Moebius-Clune D. J. Gugino B. K. Idowu O. J. Schindelbeck R. R. Ristow A. J. (2016). Comprehensive Assessment of Soil Health – The Cornell Framework. Manual, 3rd Edn. Geneva, NY: Cornell University. Moeskops B. Buchan D. Sleutel S. Herawaty L. Husen E. Saraswati R. (2010). Soil microbial communities and activities under intensive organic and conventional vegetable farming in West Java, Indonesia. Appl. Soil Ecol. 45 112–120. 10.1016/j.apsoil.2010.03.005 Mougel C. Offre P. Ranjard L. Corberand T. Gamalero E. Robin C. (2006). Dynamic of the genetic structure of bacterial and fungal communities at different developmental stages of Medicago truncatula Gaertn. cv. Jemalong line J5. New Phytol. 170 165–175. 10.1111/j.1469-8137.2006.01650.x 16539613 Mulla D. J. Bhatti A. U. Hammond M. W. Benson J. A. (1992). A comparison of winter wheat yield and quality under uniform versus spatially variable fertilizer management. Agr. Ecosyst. Environ. 38 301–311. 10.1016/0167-8809(92)90152-2 Oyekanmi E. O. Coyne D. L. Fagade O. E. Osonubi O. (2007). Improving root-knot nematode management on two soybean genotypes through the application of Bradyrhizobium japonicum, Trichoderma pseudokoningii and Glomus mosseae in full factorial combinations. Crop Prot. 26 1006–1012. 10.1016/j.cropro.2006.09.009 Paz J. O. Batchelor W. D. Colvin T. S. Logsdon S. D. Kaspar T. C. Karlen D. L. (1998). Analysis of water stress effects causing spatial yield variability in soybeans. Trans. ASAE 41 1527–1534. 10.13031/2013.17284 Ping L. Y. Boland W. (2004). Signals from the underground: bacterial volatiles promote growth in Arabidopsis. Trends Plant Sci. 9 263–266. 10.1016/j.tplants.2004.04.008 15165555 Poudel R. Jumpponen A. Schlatter D. C. Paulitz T. C. Gardener B. B. M. S. Kinkel L. L. (2016). Microbiome networks: a systems framework for identifying candidate microbial assemblages for disease management. Phytopathology 106 1083–1096. 10.1094/PHYTO-02-16-0058-FI 27482625 Prosser J. I. (2015). Dispersing misconceptions and identifying opportunities for the use of ‘omics’ in soil microbial ecology. Nat. Rev. Microbiol. 13 439–446. 10.1038/nrmicro3468 26052662 Proulx S. R. Promislow D. E. Phillips P. C. (2005). Network thinking in ecology and evolution. Trends Ecol. Evol. 20 345–353. 10.1016/j.tree.2005.04.004 16701391 Quast C. Pruesse E. Yilmaz P. Gerken J. Schweer T. Yarza P. (2012). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41 D590–D596. 10.1093/nar/gks1219 23193283 Ramette A. Frapolli M. Fischer-Le S. M. Gruffaz C. Meyer J. M. Défago G. (2011). Pseudomonas protegens sp. nov., widespread plant-protecting bacteria producing the biocontrol compounds 2, 4-diacetylphloroglucinol and pyoluteorin. Syst. Appl. Microbiol. 34 180–188. 10.1016/j.syapm.2010.10.005 21392918 Redman R. S. Sheehan K. B. Stout R. G. Rodriguez R. J. Henson J. M. (2002). Thermotolerance generated by plant/fungal symbiosis. Science 298 1581–1581. 10.1126/science.1072191 Rogovska N. P. Blackmer A. M. Mallarino A. P. (2007). Relationships between soybean yield, soil pH, and soil carbonate concentration. Soil Sci. Soc. Am. J. 71 1251–1256. 10.2136/sssaj2006.0235 Sang M. K. Kim K. D. (2012). The volatile-producing Flavobacterium johnsoniae strain GSE09 shows biocontrol activity against Phytophthora capsici in pepper. J. Appl. Microbiol. 113 383–398. 10.1111/j.1365-2672.2012.05330.x 22563881 Santoyo G. Orozco-Mosqueda M. D. C. Govindappa M. (2012). Mechanisms of biocontrol and plant growth-promoting activity in soil bacterial species of Bacillus and Pseudomonas: a review. Biocontrol Sci. Techn. 22 855–872. 10.1080/09583157.2012.694413 Schloerke B. Crowley J. Cook D. Briatte F. Marbach M. Thoen E. (2018). Ggally: Extension to ggplot2. R Package Version. Schmid F. Moser G. Müller H. Berg G. (2011). Functional and structural microbial diversity in organic and conventional viticulture: organic farming benefits natural biocontrol agents. Appl. Environ. Microbiol. 77 2188–2191. 10.1128/AEM.02187-10 21278278 Sessitsch A. Reiter B. Berg G. (2004). Endophytic bacterial communities of field-grown potato plants and their plant growth-promoting and antagonistic abilities. Can. J. Microbiol. 50 239–249. 10.1139/w03-118 15213748 Shahandeh H. Wright A. L. Hons F. M. Lascano R. J. (2005). Spatial and temporal variation of soil nitrogen parameters related to soil texture and corn yield. Agron. J. 97 772–782. 10.2134/agronj2004.0287 Shimomoto Y. Sato T. Hojo H. Morita Y. Takeuchi S. Mizumoto H. (2011). Pathogenic and genetic variation among isolates of Corynespora cassiicola in Japan. Plant Pathol. 60 253–260. 10.1111/j.1365-3059.2010.02374.x Shovan L. R. Bhuiyan M. K. A. Begum J. A. Pervez Z. (2008). In vitro control of Colletotrichum dematium causing anthracnose of soybean by fungicides, plant extracts and Trichoderma harzianum. Int. J. Sustain Crop Prod. 3 10–17. Silva W. P. K. Multani D. S. Deverall B. J. Lyon B. R. (1995). RFLP and RAPD analyses in the identification and differentiation of isolates of the leaf spot fungus Corynespora cassiicola. Aust. J. Bot. 43 609–618. 10.1071/bt9950609 Sivasakthi S. Usharani G. Saranraj P. (2014). Biocontrol potentiality of plant growth promoting bacteria (PGPR)-Pseudomonas fluorescens and Bacillus subtilis: a review. Afr. J. Agric. Res. 9 1265–1277. Smercina D. N. Evans S. E. Friesen M. L. Tiemann L. K. (2019). To fix or not to fix: controls on free-living nitrogen fixation in the rhizosphere. Appl. Environ. Microbiol. 85:e02546-18. Srour A. Y. Gibson D. J. Leandro L. F. Malvick D. K. Bond J. P. Fakhoury A. M. (2017). Unraveling microbial and edaphic factors affecting the development of sudden death syndrome in soybean. Phytobiomes 1 91–101. 10.1094/pbiomes-02-17-0009-r Suárez-Moo P. D. J. Vovides A. P. Griffith M. P. Barona-Gómez F. Cibrián-Jaramillo A. (2019). Unlocking a high bacterial diversity in the coralloid root microbiome from the cycad genus Dioon. PLoS One 14:e0211271. 10.1371/journal.pone.0211271 30726265 Sugiyama A. Ueda Y. Zushi T. Takase H. Yazaki K. (2014). Changes in the bacterial community of soybean rhizospheres during growth in the field. PLoS One 9:e100709. 10.1371/journal.pone.0100709 24955843 Sugiyama A. Vivanco J. M. Jayanty S. S. Manter D. K. (2010). Pyrosequencing assessment of soil microbial communities in organic and conventional potato farms. Plant Dis. 94 1329–1335. 10.1094/PDIS-02-10-0090 30743621 Vollmann J. Winkler J. Fritz C. N. Grausgruber H. Ruckenbauer P. (2000). Spatial field variations in soybean (Glycine max [L.] Merr.) performance trials affect agronomic characters and seed composition. Eur. J. Agron. 12 13–22. 10.1016/s1161-0301(99)00042-8 Wade S. D. Foster I. D. L. Baban S. M. J. (1996). The spatial variability of soil nitrates in arable and pasture landscapes implications for the development of geographical information system models of nitrate leaching. Soil Use Manag. 12 95–101. 10.1111/j.1475-2743.1996.tb00965.x Xu Y. X. Wang G. H. Jin J. Liu J. J. Zhang Q. Y. Liu X. (2009). Bacterial communities in soybean rhizosphere in response to soil type, soybean genotype, and their growth stage. Soil Biol. Biochem. 41 919–925. 10.1016/j.soilbio.2008.10.027 Yan H. Nelson B. Jr. (2020). Effect of temperature on Fusarium solani and F. tricinctum growth and disease development in soybean. Can. J. Plant Pathol. 42 527–537. 10.1080/07060661.2020.1745893 Yin C. Jones K. L. Peterson D. E. Garrett K. A. Hulbert S. H. Paulitz T. C. (2010). Members of soil bacterial communities sensitive to tillage and crop rotation. Soil Biol. Biochem. 42 2111–2118. 10.1016/j.soilbio.2010.08.006 Young J. M. (2010). Taxonomy of Pseudomonas syringae. J. Plant Pathol. 92 S5–S14. Yurgel S. N. Douglas G. M. Comeau A. M. Mammoliti M. Dusault A. Percival D. (2017). Variation in bacterial and eukaryotic communities associated with natural and managed wild blueberry habitats. Phytobiomes 1 102–113. 10.1094/pbiomes-03-17-0012-r Zhang F. Ge H. Zhang F. Guo N. Wang Y. Chen L. (2016). Biocontrol potential of Trichoderma harzianum isolate T-aloe against Sclerotinia sclerotiorum in soybean. Plant Physiol. Biochem. 100 64–74. 10.1016/j.plaphy.2015.12.017 26774866 Zhang J. X. Xue A. G. Zhang H. J. Nagasawa A. E. Tambong J. T. (2010). Response of soybean cultivars to root rot caused by Fusarium species. Can. J. Plant Sci. 90 767–776. 10.4141/cjps09133 Zheng Y. Han X. Zhao D. Wei K. Yuan Y. Li Y. (2021). Exploring biocontrol agents from microbial keystone taxa associated to suppressive soil: a new attempt for a biocontrol strategy. Front. Plant Sci. 12:655673. 10.3389/fpls.2021.655673 33959142