SLRs enabled a taxonomically resolved assessment of the soybean rhizosphere microbiome
The present study employed SLR technology in tandem with avidity sequencing to explore the composition and structure of the soybean rhizosphere microbiome. For 18S-ITS amplicon sequencing targeting eukaryotes, this method generated a collective 852 M short reads assembled into 1.4 M SLRs averaging a length of 1,108.31 ± 663.99 bp (full and partial length) (Supplementary Fig. 1, Additional File 2; Additional File 3). Similarly, the full-length SLRs were 1,084.69 ± 671.98 bp. From these, 1,014 denoised ASVs were classified at the genus level using the Ensemble reference database, and 44.77% could be further classified at the species level. The mapped eukaryotic SLRs were 2,328.54 ± 213.7 bp. Conversely, the 16S (prokaryotic) dataset comprised 192 M short reads, resulting in 1.4 M SLRs with a mean length of 1,362.44 ± 291.46 bp (Supplementary Fig. 1, Additional File 2; Additional File 4). The full-length contigs were 1,476.06 ± 70.51 bp in length. Mapping these to the Silva138 database permitted the identification of 895 ASVs, all classified at the species level, with 207 (23.13%) achieving strain-level classification. The mapped prokaryotic SLR length was consistent with that of all full-length SLRs, being 1,477.64 ± 25.15 bp.
The soybean rhizosphere microbiome displayed a diverse taxonomic profile spanning seven kingdoms. The eukaryotic fraction consisted primarily of partly aquatic, saprotrophic fungi that exhibit perithecial fruiting bodies and filamentous mycelial growth forms, while the prokaryotic fraction largely comprised gram-negative, mesophilic, aerobic bacteria. Beyond bacteria and fungi, protist populations from five kingdoms were observed: Alveolata, characterized by membrane-bound sacs (alveoli) beneath the plasma membrane [103]; Apusozoa, flagellated unicellular eukaryotes [104]; Heterolobosa (i.e., Heterolobosea), protists with both amoeboid and flagellated stages [105]; Rhizaria, identified by their thread-like pseudopodia [106]; and Stramenopila, which encompasses diatoms, brown algae, and oomycetes, differentiated by their heterokont flagella [107].
Growth stage and spatial heterogeneity best explained microbiome structure
In assessing microbiome complexity and structure, the within-sample characteristics of ASV richness, diversity, and evenness were determined and compared across treatments, cultivars, and growth stages (Fig. 2A, B, C). ASV richness was measured with the non-parametric estimator Chao1 given its capacity to project undetected taxa based on the abundance of those rarely observed in a dataset. This extrapolated measure accounts for potential undersampling in high-diversity environments (e.g., soil), providing a more thorough assessment of community richness [49, 108]. Herein, the mean eukaryotic Chao1 index value was 1,033 ± 36, and was increased marginally in growth stage R6 vs R2 (p-value = 0.07, ME = 220.62) and in cultivar CZ4979X vs CZ4810X (p-value = 0.09, ME = 205.38) (Fig. 2C). A cultivar-growth stage interaction was also observed and explained the highest proportion of variance (marginalized R2 = 0.88), with CZ4979X:R6 being significantly less than the reference (p-value = 0.0023, ME = -533.75) (Fig. 2C). The prokaryotic Chao1 index was 85.44 ± 8.76 and demonstrated a significant increase in R6 vs R2 (p-value = 0.0099, ME = 50.28) and significant decreases in CZ4979X vs CZ4810X (p-value = 0.04, ME = -32.24) and in treatment control vs biostimulant (p-value = 0.04, ME = -31.72) (Fig. 2C). Moreover, growth stage explained the highest proportion of variance for prokaryotic Chao1 (marginalized R2 = 0.45). The baseline (V1) sampling showed no significant trends in the eukaryotic or prokaryotic Chao1 indices.
Shannon and Simpson indices were estimated to comprehensively assess α diversity. The Shannon index integrates ASV richness and evenness, with increasing values indicative of greater diversity and uniformity among ASV abundances [51]. Conversely, the Simpson index quantifies ASV dominance, where elevated values denote diminished dominance and heightened diversity [50]. The two indices displayed cooperative, statistically insignificant trends across fixed effect levels in the present study. Eukaryotic α diversity (6.3 ± 0.07 Shannon and 0.99 ± 0.002 Simpson) peaked at the R2 growth stage (evidenced by ME at V6 vs R2 and R2 vs R6), was reduced in CZ4979X vs CZ4810X, and was increased in control vs biostimulant, with growth stage explaining the highest proportion of variance (marginalized R2 = 0.73 and 0.33 for Shannon and Simpson indices, respectively) (Fig. 2C). The prokaryotic α diversity indices (3.82 ± 0.09 Shannon and 0.94 ± 0.007 Simpson) demonstrated an opposing trend, increasing with time and being reduced in control vs biostimulant (Fig. 2C). The direction of change between CZ4979X and CZ4810X was consistent with eukaryotic α diversity, and variance was most attributed to cultivar for both indices (marginalized R2 = 0.41 and 0.58 for Shannon and Simpson, respectively). Notably, statistically significant differences were detected in the baseline between control plots and those to which biostimulants would be applied, with the control plots displaying reduced eukaryotic and prokaryotic α diversity (eukaryotic Shannon p-value = 0.1, ME = -0.19; eukaryotic Simpson p-value = 0.05, ME = -0.001; prokaryotic Shannon p-value = 0.005, ME = -0.41; prokaryotic Simpson p-value = 0.002, ME = -0.04).
While α diversity incorporates both richness and evenness, relying solely on composite diversity indices might obscure their individual contributions to ecosystem function [109]. Pielou's evenness was therefore employed to quantify the count distribution across ASVs. The index was consistent between fixed effect levels and domains, with a mean of 0.91 ± 0.007 for eukaryotes, 0.90 ± 0.009 for prokaryotes, and no statistically significant trends observed (Fig. 2A, C). This observation was supported by a strong association between Chao1 richness and Shannon diversity (ρ = 0.95, p-value = 4.58e− 63) (Fig. 2B). Eukaryotic evenness was influenced predominantly by growth stage (marginalized R2 = 0.98), while no patterns were present for the prokaryotic dataset. Consistent with Shannon and Simpson diversity, the baseline control plots showed reduced prokaryotic evenness in comparison to biostimulant plots (p-value = 0.0004, ME = -0.089).
Compositional dissimilarity was assessed using CSS-normalized counts from mapped ASVs. The phylum-level composition was first visualized across individual samples and experimental conditions, revealing that the eukaryotic rhizosphere microbiome was largely composed of Ascomycota and the prokaryotic microbiome dominated by Proteobacteria (Fig. 3A, B, D, E). Dissimilarity matrices were constructed subsequently leveraging Bray-Curtis, Euclidean, and Jaccard distances, and compositional trends were inferred between fixed effect levels with PERMANOVA and β dispersion estimation. Findings were consistent across the matrices (Supplementary Table 1, 2, Additional File 2); therefore, representative ordinations were generated with Bray-Curtis dissimilarity (Fig. 3C, F). Notably, growth stage had the most significant impact on eukaryotic microbiome composition (Bray-Curtis R2 = 0.08, p-value = 0.0001) (Fig. 3G; Supplementary Table 1, Additional File 2), yet ANOVA suggested heterogeneous β dispersion across growth stage levels (Bray-Curtis p-value = 0.003) (Supplementary Table 1, Additional File 2). Both PCoA and NMDS ordinations implicated a strong spatial effect on eukaryotic microbiome composition as well (Fig. 3F). Prokaryotic microbiome composition showed neither strong statistical nor qualitative trends, albeit a treatment-growth stage interaction explained the highest proportion of variance (Bray-Curtis R2 = 0.06, p-value = 0.08) (Fig. 3G, H; Supplementary Table 2, Additional File 2). Lastly, while no compositional trends emerged from the baseline measurements, β dispersion did vary between treatment levels for Euclidean dissimilarity (p-value = 0.05) (Supplementary Fig. 2; Supplementary Table 3, Additional File 2).
Complementary to PERMANOVA, the Bray-Curtis indices were decomposed with the Similarity Percentage method [71] to discern ASVs most influential for pairwise similarity/dissimilarity between fixed effect levels. The five ASVs to which dissimilarity was most attributed were
identified for each comparison, and their lowest taxonomic classification retrieved. Interestingly, Cyathus stercoreus and an Acremonium species were identified for all five fixed effect comparisons with the eukaryotic Bray-Curtis matrix, followed by a Nectria species and Phallus rugulosus in four comparisons, and Neocosmospora falciformis in three. Plectosphaerella cucumerina was exclusive to comparisons between vegetative and reproductive growth stages, while Mortierella elongata and a Polymyxa species were exclusive to R6 vs R2 (Fig. 3I). Likewise, ASVs corresponding to Nitrosomonas europaea, Clostridium sporogenes, Bacteroides thetaiotaomicron, and Rhodospirillum rubrum F11 were most influential for all fixed effect comparisons of prokaryotic Bray-Curtis dissimilarity, with Bifidobacterium adolescentis identified in four comparisons and Bradyrhizobium japonicum CCBAU 15618 exclusive to R6 vs R2 (Fig. 3I). Given the vast overlap of these ASVs between pairwise comparisons, and that all demonstrated statistical insignificance (p-value > 0.05 based on 9,999 permutations), ASV identification more likely reflected high abundance/variation across the amplicon datasets than contribution to dissimilarity, which is a common (yet often misinterpreted) element of Similarity Percentage analysis [71].
Microbiota with agriculturally-relevant life strategies exhibited distinct membership trends across fixed effect levels
Core, unique, and differentially abundant taxa were identified at genus and species levels for eukaryotes and prokaryotes, respectively, across all growth stages (baseline and those succeeding biostimulant application). Unsurprisingly, 20/22 core taxa (defined as those with a ≥ 0.5 prevalence across all samples) were eukaryotic genera that mostly exhibit a partly aquatic, saprotrophic lifestyle with filamentous mycelial growth and perithecial fruiting bodies (Fig. 4A). Furthermore, 7/20 (35%) were annotated as plant pathogens (Fig. 4A). The prokaryotic core members included the nitrite-oxidizing bacterium Nitrospira japonica [110] and the scarcely reported bacterium Pseudolabrys taiwanensis [111] (Fig. 4A).
Eukaryotic genera/prokaryotic species unique to a fixed effect level encompassed a collective 217 taxa between treatments (30 eukaryotes, 187 prokaryotes), 296 between cultivars (88 eukaryotes, 208 prokaryotes), and 260 between growth stages (74 prokaryotes, 186 eukaryotes) (Fig. 4B). Those unique to biostimulant-treated samples were all prokaryotes and included four species in the symbiotic genus Bradyrhizobium [112], the additional rhizobia Mesorhizobium ciceri, Mesorhizobium plurifarium, Rhizobium grahamii, and Rhizobium massiliae, Pseudomonas fluorescens [113], three species of Bacillus [114], and 10 Streptomyces species [115] (Fig. 4B). Notable taxa exclusive to the control treatment were the arbuscular mycorrhizal/root-associated genus Paraglomus [116], six plant pathogens (including the soybean disease-causing oomycete genus Phytophthora [117]), and Bradyrhizobium stylosanthis (Fig. 4B).
Perhaps the most distinct trend was the exclusivity of plant pathogens to a particular cultivar. Eight pathogens were unique to CZ4979X (SDS-tolerant cultivar), including Phytophthora and the soybean-parasitizing fungal genus Septoria [118] (Fig. 4B). Other CZ4979X-specific taxa were Bradyrhizobium algeriense, Bradyrhizobium betae, Mesorhizobium plurifarium, Nitrospira multiformis [119], and Pseudomonas fluorescens (Fig. 4B). Samples from CZ4810X, which presumably have heightened susceptibility to the soybean disease SDS in comparison to CZ4979X, had 15 plant pathogens not present in CZ4979X samples, including the soybean disease-causing genera Diaporthe [120], Macrophomina [121], and Rhizoctonia [122] (Fig. 4B). Additional taxa unique to CZ4810X included the genus Paraglomus, Bradyrhizobium lupini, Bradyrhizobium stylosanthis, Mesorhizobium ciceri, Rhizobium cellulosilyticum, and nine species of Streptomyces.
Of the 260 taxa exclusive to a single growth stage, 48 were unique to V1 (18 eukaryotes, 30 prokaryotes), 68 to V6 (17 eukaryotes, 51 prokaryotes), 61 to R2 (22 eukaryotes, 39 prokaryotes), and 83 to R6 (17 eukaryotes, 66 prokaryotes) (Fig. 4B). Those corresponding to V1 included six plant pathogens (including Septoria), Nitrospira multiformis, and Bradyrhizobium algeriense (Fig. 4B). The V6 growth stage was characterized by three distinct plant pathogens, Paraglomus, Bradyrhizobium lupini, two Bacillus species, and nine Streptomyces species (Fig. 4B). Similarly, R2 possessed five unique plant pathogens (including Rhizoctonia), Bradyrhizobium betae, Mesorhizobium plurifarium, and Rhizobium cellulosilyticum (Fig. 4B). The final sampled growth stage contained 5 unique plant pathogens (including Macrophomina, Phytophthora, and Cercospora [123]), two Bacillus species, Bradyrhizobium stylosanthis, and Pseudomonas fluorescens (Fig. 4B). A comprehensive list of condition-specific taxa can be found in Additional File 5.
Differentially abundant taxa between fixed effect levels were determined by fitting CSS-normalized counts with a ZINB regression model. The greatest number of those statistically enriched/depleted (q-value < 0.25) was observed between levels of treatment (22 eukaryotes, 19 prokaryotes, total n = 41), with 18 taxa being enriched and 23 depleted in biostimulant vs control (Fig. 4C). Notably, this included the differential abundance of saprotrophic fungi, the depletion of five fungal pathogens, and the enrichment of Bradyrhizobium elkanii, Bradyrhizobium japonicum, Bradyrhizobium lablabi, and Mesorhizobium amorphae (Fig. 4C). The most depleted taxa in biostimulant-treated samples were the potential human/foodborne pathogenic bacteria Clostridium sporogenes and Escherichia coli [124, 125] (Fig. 4C). Twenty-seven taxa (13 eukaryotes, 14 prokaryotes) were differentially abundant between cultivars, 14 of which were enriched and 13 depleted in CZ4979X vs CZ4810X (Fig. 4C). Most of the identified eukaryotes displayed marginal depletion in CZ4979X (including 3/4 differentially abundant plant pathogens) (Fig. 4C). In contrast, the majority of prokaryotes were enriched, including the inorganic phosphate-solubilizing bacterium Bacillus acidiceler [126], Bradyrhizobium elkanii, Bradyrhizobium japonicum, and Escherichia coli (Fig. 4C).
As expected, an increase in differentially abundant taxa was associated with the temporal distinctiveness of compared growth stages. The greatest number was observed in the R6 vs V1 comparison (distance = 3 growth stages) (23 eukaryotes, 15 prokaryotes, total n = 38), followed by R6 vs V6 (distance = 2 growth stages) (16 eukaryotes, 13 prokaryotes, total n = 29), R2 vs V1 (distance = 2 growth stages) (19 eukaryotes, 10 prokaryotes, total n = 29), V6 vs V1 (distance = 1 growth stage) (14 eukaryotes, 13 prokaryotes, total n = 27), R2 vs V6 (distance = 1 growth stage) (15 eukaryotes, 10 prokaryotes, total n = 25), and R6 vs R2 (distance = 1 growth stage) (12 eukaryotes, 10 prokaryotes, total n = 22) (Fig. 4C). The eukaryotic dataset was defined by the enrichment of saprotrophic genera Acremonium and Chaetomium with the progression of time, which was particularly distinct between vegetative and reproductive growth, and a depletion in the saprotrophic genus Conocybe (Fig. 4C). The reproductive growth stages also displayed an overall enrichment in Microbacterium rhizosphaerae, Bradyrhizobium elkanii, and Bradyrhizobium japonicum, although the latter was depleted at all stages compared to V1 (Fig. 4C). To this end, the V1 growth stage showed unique enrichment of eukaryotic genera Ciliophora, Cyberlindnera, Podospora, and Setophoma, of prokaryotic species Bradyrhizobium japonicum and Pseudarthrobacter chlorophenolicus, and a depletion in the genus Neocosmospora and species Arthrobacter humicola, Bacillus megaterium, and Microlunatus panaciterrae (Fig. 4C). Differentially abundant taxa and associated metadata are provided in Additional File 6.
A treatment-cultivar interaction defined genus-level microbial co-occurrence network structure
Putative genus-genus associations were inferred with Spearman and Pearson correlation methods. The associations were determined from CSS-normalized absolute abundances to mitigate the concomitant limitations of compositionality bias and biases stemming from differential sampling efficiency of taxa [127, 128]. Moreover, the p-values of pairwise associations were corrected for multiple testing given the prevalence of Type I errors during microbial co-occurrence network construction [91]. Herein, the global Spearman co-occurrence network comprised 826 edges (associations) and 188 nodes (genera), all of which presented a low mean relative abundance (< 0.5%) across the dataset (Fig. 5A, C). The global Pearson network possessed 1,007 edges and 294 nodes, with evident variability in mean relative abundance observed (Fig. 5B, C). In addition, 462/1,371 (33.7%) of the edges were shared between the networks (Fig. 5C).
Condition-specific co-occurrence networks, defined by treatment-cultivar-growth stage combinations, were constructed with Spearman and Pearson associations as outlined for global networks. The Spearman networks presented an increase in network density with time (vegetative node n = 11.13 ± 2.54, edge n = 16.38 ± 4.88; reproductive node n = 22.88 ± 1.51, edge n = 37.50 ± 4.27), and a potential treatment effect at the V6 growth stage (control node n = 5.0 ± 1.0, edge n = 6; biostimulant node n = 21.0 ± 2.0, edge n = 35.50 ± 1.50) (Fig. 5D; Table 1). Given the density of the biostimulant-CZ4979X-V1 network (node n = 16, edge n = 27), which is the baseline sampling preceding biostimulant application, the variation observed at V6 may better reflect a
Treatment-cultivar interaction between biostimulant and CZ4810X (Fig. 5D; Table 1). In support of this notion, the greatest trend in Pearson network density was observed in biostimulant-CZ4810X networks at growth stages succeeding biostimulant application (control node n = 99.0 ± 14.42, edge n = 612.33 ± 86.81; biostimulant node n = 136.0 ± 4.58, edge n = 1,533.67 ± 67.85) (Fig. 5E; Table 2). Pearson networks were overall denser than Spearman networks (Spearman node n = 17.0 ± 2.08, edge n = 26.94 ± 4.15; Pearson node n = 107.63 ± 4.78, edge n = 791.50 ± 99.20), which was consistent with the global co-occurrence networks, and 257/10,665 (2.41%) of edges were shared between the association methods (Fig. 5E, G).
Condition-specific network topology was further defined by centralization degree (the extent to which a single node “controls” a network), cluster count (the number of separate, interconnected groups), connectance (the proportion of all possible links that are actual connections), giant component size (the size of the largest connected subgraph), hub count (the number of nodes with a Kleinberg hub centrality score > 0.2), KS test statistic (how closely degree distribution adheres to a scale-free topology), mean node degree (the average number of connections per node), and
Table 1
Topological properties for condition-specific Spearman networks
Network
|
Centralization degree
|
Cluster count
|
Connectance
|
Edge count
|
Giant component size
|
Hub count
|
KS stat
|
Mean degree
|
Modularity
|
Node count
|
V1-Bst-79X
|
0.11
|
4.00
|
0.23
|
27.00
|
6.00
|
6.00
|
0.44
|
3.38
|
0.62
|
16.00
|
V6-Bst-79X
|
0.06
|
4.00
|
0.22
|
37.00
|
6.00
|
6.00
|
0.23
|
3.89
|
0.71
|
19.00
|
R2-Bst-79X
|
0.05
|
7.00
|
0.13
|
55.00
|
6.00
|
12.00
|
0.28
|
3.67
|
0.80
|
30.00
|
R6-Bst-79X
|
0.04
|
7.00
|
0.10
|
24.00
|
4.00
|
4.00
|
0.36
|
2.18
|
0.84
|
22.00
|
V1-Ctrl-79X
|
0.00
|
2.00
|
0.40
|
6.00
|
3.00
|
6.00
|
NA
|
2.00
|
0.50
|
6.00
|
V6-Ctrl-79X
|
0.00
|
2.00
|
0.40
|
6.00
|
3.00
|
6.00
|
NA
|
2.00
|
0.50
|
6.00
|
R2-Ctrl-79X
|
0.06
|
6.00
|
0.13
|
31.00
|
5.00
|
5.00
|
0.36
|
2.82
|
0.79
|
22.00
|
R6-Ctrl-79X
|
0.13
|
6.00
|
0.16
|
36.00
|
7.00
|
7.00
|
0.36
|
3.27
|
0.63
|
22.00
|
V1-Bst-10X
|
0.00
|
2.00
|
0.40
|
6.00
|
3.00
|
6.00
|
NA
|
2.00
|
0.50
|
6.00
|
V6-Bst-10X
|
0.05
|
6.00
|
0.13
|
34.00
|
5.00
|
5.00
|
0.25
|
2.96
|
0.80
|
23.00
|
R2-Bst-10X
|
0.01
|
4.00
|
0.20
|
21.00
|
4.00
|
4.00
|
0.36
|
2.80
|
0.73
|
15.00
|
R6-Bst-10X
|
0.08
|
6.00
|
0.14
|
39.00
|
6.00
|
6.00
|
0.33
|
3.25
|
0.77
|
24.00
|
V1-Ctrl-10X
|
0.00
|
3.00
|
0.25
|
9.00
|
3.00
|
9.00
|
NA
|
2.00
|
0.67
|
9.00
|
V6-Ctrl-10X
|
0.00
|
1.00
|
1.00
|
6.00
|
4.00
|
4.00
|
NA
|
3.00
|
0.00
|
4.00
|
R2-Ctrl-10X
|
0.08
|
6.00
|
0.16
|
51.00
|
7.00
|
7.00
|
0.23
|
3.92
|
0.72
|
26.00
|
R6-Ctrl-10X
|
0.10
|
5.00
|
0.19
|
43.00
|
7.00
|
7.00
|
0.28
|
3.91
|
0.68
|
22.00
|
Column 1: Growth Stage-Treatment-Cultivar; Ctrl = Control; Bst = Biostimulant; 79X = CZ4979X; 10X = CZ4810X
|
modularity (the strength of division into distinct modules/communities) (Table 1, Table 2). Of these, centralization degree, cluster count, giant component size, mean degree, and modularity demonstrated a statistically significant, positive association (Spearman’s ρ statistic > 0, p-value < 0.05) with both edge count and node count across Spearman networks, while connectance presented a negative association (Fig. 5F; Table 1). The KS test statistic was negatively associated with edge count (Fig. 5F; Table 1). Conversely, only hub count and mean degree presented significant associations with Pearson network node count (both of which were positive), and centralization degree, connectance, hub count, mean degree, and modularity were associated
Table 2
Topological properties for condition-specific Pearson networks
Network
|
Centralization degree
|
Cluster count
|
Connectance
|
Edge count
|
Giant component size
|
Hub count
|
KS stat
|
Mean degree
|
Modularity
|
Node count
|
V1-Bst-79X
|
0.16
|
7.00
|
0.12
|
627.00
|
31.00
|
29.00
|
0.21
|
12.29
|
0.64
|
102.00
|
V6-Bst-79X
|
0.20
|
6.00
|
0.11
|
489.00
|
36.00
|
21.00
|
0.24
|
10.19
|
0.72
|
96.00
|
R2-Bst-79X
|
0.11
|
8.00
|
0.06
|
447.00
|
27.00
|
19.00
|
0.13
|
7.58
|
0.81
|
118.00
|
R6-Bst-79X
|
0.16
|
6.00
|
0.12
|
719.00
|
32.00
|
31.00
|
0.17
|
13.31
|
0.67
|
108.00
|
V1-Ctrl-79X
|
0.13
|
5.00
|
0.17
|
700.00
|
30.00
|
27.00
|
0.17
|
15.56
|
0.61
|
90.00
|
V6-Ctrl-79X
|
0.14
|
5.00
|
0.13
|
485.00
|
24.00
|
21.00
|
0.18
|
11.02
|
0.73
|
88.00
|
R2-Ctrl-79X
|
0.10
|
5.00
|
0.10
|
399.00
|
38.00
|
19.00
|
0.19
|
8.77
|
0.77
|
91.00
|
R6-Ctrl-79X
|
0.29
|
5.00
|
0.14
|
951.00
|
61.00
|
49.00
|
0.20
|
16.12
|
0.35
|
118.00
|
V1-Bst-10X
|
0.14
|
5.00
|
0.15
|
725.00
|
29.00
|
29.00
|
0.20
|
14.50
|
0.64
|
100.00
|
V6-Bst-10X
|
0.14
|
5.00
|
0.17
|
1,398.00
|
43.00
|
40.00
|
0.14
|
21.51
|
0.61
|
130.00
|
R2-Bst-10X
|
0.26
|
4.00
|
0.18
|
1,599.00
|
88.00
|
61.00
|
0.25
|
24.05
|
0.25
|
133.00
|
R6-Bst-10X
|
0.18
|
6.00
|
0.15
|
1,604.00
|
49.00
|
48.00
|
0.18
|
22.12
|
0.57
|
145.00
|
V1-Ctrl-10X
|
0.11
|
5.00
|
0.12
|
684.00
|
29.00
|
25.00
|
0.18
|
12.91
|
0.74
|
106.00
|
V6-Ctrl-10X
|
0.15
|
7.00
|
0.18
|
448.00
|
56.00
|
23.00
|
0.18
|
12.62
|
0.62
|
71.00
|
R2-Ctrl-10X
|
0.16
|
5.00
|
0.11
|
743.00
|
56.00
|
31.00
|
0.18
|
12.49
|
0.73
|
119.00
|
R6-Ctrl-10X
|
0.16
|
4.00
|
0.11
|
646.00
|
53.00
|
30.00
|
0.17
|
12.07
|
0.68
|
107.00
|
Column 1: Growth Stage-Treatment-Cultivar; Ctrl = Control; Bst = Biostimulant; 79X = CZ4979X; 10X = CZ4810X
|
with edge count (all positive associations except modularity) (Fig. 5F; Table 2). The discrepancy in network topology was reflected in the ranking of co-occurrence networks, as eight of the 10 metrics showed statistically significant differences in conditional ranking between the association methods (Friedman rank sum test p < 0.05) (Supplementary Fig. 3B, Additional File 2). Despite this, all nodes represented in Spearman networks were present in Pearson networks (Supplementary Fig. 3A, Additional File 2), and the hierarchical clustering of Pearson networks by node centrality score supported the overarching trend of treatment-cultivar interaction defining co-occurrence network structure (Fig. 5H).
Edaphic property dynamism and integration within phenotype-taxon networks
To contextualize microbiome dynamics, 24 edaphic parameters were measured for each soil sample, and differences between fixed effect levels and their interactions were determined using GLMMs. Of these, 18 parameters displayed significant variation (p < 0.05) between levels of one or more explanatory variables, with soil K and SOM affected most (n = 5). These were followed by ARS, B, Ca, and P (n = 3), and ALP, Ca/Mg, CEC, Fe, K/Mg, Na, and NO3- (n = 2). The parameters GBA3, Mg, NAG, S, and soil pH each displayed variation between levels of a single explanatory variable (Fig. 6A). Inversely, cultivar explained observed variation for 14 edaphic parameters, 12 of which were increased in CZ4979X vs CZ4810X (Fig. 6A). Ten parameters were increased in V6 vs R2, and six differed significantly in the R6 vs R2 comparison (two increased and four decreased) (Fig. 6A). The explanatory variables treatment, treatment-cultivar interaction, and cultivar-growth stage interaction each explained observed variation in a lesser number of edaphic parameters (Fig. 6A). Furthermore, rank-based measures of association were inferred between edaphic parameters with Spearman’s ρ statistic, rendering 109 significant associations (positive n = 73, negative n = 36) (Fig. 6B). Of these, the Ca/Mg ratio presented the most significant associations of the parameters (positive n = 6, negative n = 9; total n = 15) (Fig. 6B).
Significant microbial associations spanning one or both global networks were coupled with edaphic data to construct phenotype-taxon networks. Briefly, lasso regression was used to identify taxa putatively associated with each edaphic parameter given its propensity to assign coefficient penalties in instances when sample size is small relative to feature (node) count [89, 129]. A reduced GLM was deployed thereafter to provide directionality to phenotype-node associations and was overlaid with the global microbial association dataset to derive final networks. The phenotype-taxon networks comprised 285.63 ± 1.98 nodes (eukaryotic n = 165.71 ± 1.70, prokaryotic n = 119.92 ± 0.31) and 1,417.46 ± 24.33 edges, and a collective 306 nodes (183 eukaryotes, 123 prokaryotes) were represented in at least one network (Fig. 6C, D). Given the inherent zero inflation in taxon abundance (Fig. 6E), node prioritization was preceded by filtering to those with a ≥ 0.2 prevalence. A composite score was then computed with modularity measures and Kleinberg's hub centrality (see “Materials and Methods”), and the top 20 nodes by mean composite score were extracted. These nodes included 17 eukaryotic genera that are predominantly aquatic/partly aquatic with mycelial growth and a saprotrophic lifestyle, as well as the prokaryotic genera Bradyrhizobium, Lysobacter, and Mycobacterium (Fig. 6F, G). Six of the 20 (30%) were represented in the core microbiome. In addition, seven (35%) of the selected eukaryotic genera were annotated as plant pathogens.
Rank-based measures of association were determined thereafter between prioritized nodes, as well as between nodes and edaphic phenotypes. Thirty node-node associations were significant (positive n = 22, negative n = 8), with the greatest number of those encompassing the pathogenic/saprotrophic genus Fusarium (positive n = 6, negative n = 1, total n = 7) and the saprotrophic genus Neocosmospora (positive n = 5, negative n = 2, total n = 7) (Fig. 6F). The node-phenotype analysis rendered 90 significant associations (positive n = 34, negative n = 56), with the most represented nodes being the pathogenic Apiospora (positive n = 6, negative n = 5, total n = 11) and the saprotrophic Phallus (positive n = 4, negative n = 7, total n = 11), and the most represented phenotypes being B (positive n = 1, negative n = 6, total n = 7) and Buffer pH (positive n = 3, negative n = 4, total n = 7) (Fig. 6G).
Agronomic property dynamism and integration within phenotype-taxon networks
At the R8 growth stage, the agronomic parameters 100-seed weight, aboveground biomass, belowground biomass, pods/plant, and theoretical yield were determined for each of the 16 plots. Variations between levels of treatment, cultivar, and their interaction were then determined using GLMMs. Biostimulant-treated plots had significantly increased 100-seed weight (p-value = 0.05, ME = 0.83), and a marginal decrease was observed in CZ4979X vs CZ4810X (p-value = 0.08, ME = -0.73) (Fig. 7A, B). Expectedly, consistent directional changes were present for theoretical yield, with biostimulant demonstrating a marginal increase (p-value = 0.09, ME = 464.66) over the control, and CZ4979X being decreased in comparison to CZ4979X (albeit insignificantly) (Fig. 7A, B). Above- and belowground biomass showed opposing trends between cultivars, being statistically increased (p-value = 0.002, ME = 3.66) and insignificantly decreased, respectively, in CZ4979X vs CZ4810X (Fig. 7A, B). Both parameters were increased insignificantly in biostimulant vs control (Fig. 7A, B). Lastly, pods/plant was increased significantly in CZ4979X vs CZ4810X (p-value = 0.02, ME = 14.23) and increased insignificantly in biostimulant vs control (Fig. 7A, B). The cultivar explained the highest proportion of variance for aboveground biomass and pods/plant (marginalized R2 = 0.85 and 0.71, respectively), while variance in 100-seed weight, belowground biomass, and theoretical yield were best explained by treatment (marginalized R2 = 0.56, 1.0, and 0.95, respectively).
Phenotype-taxon networks were constructed and visualized as described for edaphic parameters, encompassing 128 ± 6.66 nodes (eukaryotic n = 75.20 ± 5.64, prokaryotic n = 52.80 ± 1.56) and 448.20 ± 35.11 edges (Fig. 7C, D, E). Furthermore, 148 nodes (89 eukaryotes, 59 prokaryotes) were present in one or more of the networks (Fig. 7C, D). Consistent with the previous analysis, node prioritization identified 17 eukaryotic genera that are predominantly aquatic/partly aquatic with mycelial growth and a saprotrophic lifestyle, seven of which are also annotated as plant pathogens (Fig. 7F, G). The remaining genera included the prokaryotes Bradyrhizobium, Burkholderia, and Lysobacter (Fig. 7F, G). Six of the 20 were represented in the core microbiome. In addition, 10 (50%) nodes were prioritized for both agronomic and edaphic networks. Rank-based measures of association were next inferred between the top 20 nodes and between nodes and agronomic parameters. There were 35 significant associations between nodes (positive n = 31, negative n = 4), and Neocosmospora was most represented (positive n = 8, negative n = 1, total n = 9) (Fig. 7F). Additionally, six significant associations were discerned between nodes and parameters, with four nodes positively associated with aboveground biomass (Naegleria, Setophoma, Neonectria, and Fusariella), Leptosphaeria positively associated with pods/plant, and Phaeosphaeriopsis associated negatively with theoretical yield (Fig. 7G).