# Spatiotemporal Heterogeneity and Intragenus Variability in Rhizobacterial Associations with Brassica rapa Growth

Scott A Klasek, Marcus T Brock, W John Calder, Hilary G Morrison, Cynthia Weinig, Lois Maïgnien

► To cite this version:

Scott A Klasek, Marcus T Brock, W John Calder, Hilary G Morrison, Cynthia Weinig, et al., Spatiotemporal Heterogeneity and Intragenus Variability in Rhizobacterial Associations with Brassica rapa Growth. mSystems, 2022, 7 (3), pp.e00060-22. 10.1128/msystems.00060-22. hal-05143735

HAL Id: hal-05143735

https://hal.univ-brest.fr/hal-05143735v1

Submitted on 4 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Spatiotemporal Heterogeneity and Intragenus Variability in Rhizobacterial Associations with *Brassica rapa* Growth

Scott A. Klasek,ª,ᵇ* Marcus T. Brock,ᵇ W. John Calder,ᵇ Hilary G. Morrison,ª Cynthia Weinig,ᵇ,ᶜ,>d Lois Maïgnienª,ᵉ

**ABSTRACT** Microbial communities in the rhizosphere are distinct from those in soils and are influenced by stochastic and deterministic processes during plant development. These communities contain bacteria capable of promoting growth in host plants through various strategies. While some interactions are characterized in mechanistic detail using model systems, others can be inferred from culture-independent methods, such as 16S amplicon sequencing, using machine learning methods that account for this compositional data type. To characterize assembly processes and identify community members associated with plant growth amid the spatiotemporal variability of the rhizosphere, we grew *Brassica rapa* in a greenhouse time series with amended and reduced microbial treatments. Inoculation with a native soil community increased plant leaf area throughout the time series by up to 28%. Despite identifying spatially and temporally variable amplicon sequence variants (ASVs) in both treatments, inoculated communities were more highly connected and assembled more deterministically overall. Using a generalized linear modeling approach controlling for spatial variability, we identified 43 unique ASVs that were positively or negatively associated with leaf area, biomass, or growth rates across treatments and time stages. ASVs of the genus *Flavobacterium* dominated rhizosphere communities and showed some of the strongest positive and negative correlations with plant growth. Members of this genus, and growth-associated ASVs more broadly, exhibited variable connectivity in networks independent of growth association (positive or negative). These findings suggest host-rhizobacterial interactions vary temporally at narrow taxonomic scales and present a framework for identifying rhizobacteria that may work independently or in concert to improve agricultural yields.

**IMPORTANCE** The rhizosphere, the zone of soil surrounding plant roots, is a hot spot for microbial activity, hosting bacteria capable of promoting plant growth in ways like increasing nutrient availability or fighting plant pathogens. This microbial system is highly diverse and most bacteria are unculturable, so to identify specific bacteria associated with plant growth, we used culture-independent community DNA sequencing combined with machine learning techniques. We identified 43 specific bacterial sequences associated with the growth of the plant *Brassica rapa* in different soil microbial treatments and at different stages of plant development. Most associations between bacterial abundances and plant growth were positive, although similar bacterial groups sometimes had different effects on growth. Why this happens will require more research, but overall, this study provides a way to identify native bacteria from plant roots that might be isolated and applied to boost agricultural yields.

**KEYWORDS** rhizosphere, microbial community, plant growth promotion, 16S RNA, *Brassica rapa*, feature selection, microbial communities, rhizosphere-inhabiting microbes

**Editor** Christopher R. Anderton, Pacific Northwest National Laboratory

**Copyright** © 2022 Klasek et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license.

Address correspondence to Scott A. Klasek, scott.a.klasek@gmail.com.

*Present address: Scott Klasek, Department of Plant Pathology, University of Minnesota, Minneapolis, Minnesota, USA.

The authors declare no conflict of interest.

**Received** 20 January 2022

**Accepted** 17 April 2022

**Published** 16 May 2022

<footnote>ªJosephine Bay Paul Center, Marine Biological Laboratory, Woods Hole, Massachusetts, USA</footnote>

<footnote>ᵇDepartment of Botany, University of Wyoming, Laramie, Wyoming, USA</footnote>

<footnote>ᶜProgram in Ecology, University of Wyoming, Laramie, Wyoming, USA</footnote>

<footnote>ᵈDepartment of Molecular Biology, University of Wyoming, Laramie, Wyoming, USA</footnote>

<footnote>ᵉLaboratory of Microbiology of Extreme Environments, UMR 6197, Institut Européen de la Mer, Université de Bretagne Occidentale, Plouzane, France</footnote>
Millimeters in thickness, the rhizosphere is the zone of soil that surrounds plant roots and supports bacterial communities distinct from those in bulk soils (1, 2). Plants induce this "rhizosphere effect" by exuding organic compounds from roots, which recruit rhizosphere bacteria and change bacterial community structure (3, 4). These rhizobacteria can enhance plant growth directly or indirectly by mobilizing soil nutrients (5, 6), suppressing disease (7) and herbivory (8, 9), mediating hormone signaling (10), and increasing tolerance to physical stresses such as drought (11) and salinity (12). Leveraging these plant growth-promoting rhizobacteria (PGPR) and their microbial communities to increase agricultural output is of great interest and importance (13–15).

Rhizosphere microbial communities are complex and spatiotemporally variable and can interact with hosts and with each other to maintain plant health in ways that individual isolates cannot (16, 17). Although rhizosphere composition varies to some extent across host species (18), host genetic and epigenetic regulatory mechanisms, including the circadian clock and histone methylation, can affect rhizosphere recruitment (19, 20). Variations in root exudate profiles across *Arabidopsis* growth stages have been interpreted to change rhizobacterial community composition (21), while subsequent investigations have detailed how certain exudate compounds can recruit or inhibit specific bacterial taxa within communities, often in response to nutrient stress (22–24). Population-level spatiotemporal variability has also been observed, such as in a reporter assay where the activity of a root-associated *Pseudomonas* aligned with expected rhizodeposition patterns during root growth (25). Not only guided by host influences, the assembly of plant-associated bacterial communities is to an extent shaped by stochastic forces across small spatial or temporal scales (26–28).

Considering this, high-resolution sampling is required to disentangle specific drivers of rhizosphere community assembly.

Several strategies have been used to identify plant growth-promoting rhizobacterial community sequence data (29, 30). Another used a random forest approach to identify a set of 75 operational taxonomic units (OTUs) whose relative abundances could explain 31% of the difference in grain weight of foxtail millet (31). Although the human microbiome literature has employed a variety of machine-learning strategies for identifying microbial community features (taxa or genes) associated with health or disease classifications (32), fewer studies have used these strategies to identify associations between rhizosphere community members and continuous plant growth data.

To identify and characterize plant growth-promoting bacteria within the context of rhizosphere community dynamics, we grew *Brassica rapa* in a 2-week time series greenhouse experiment using two soil treatments. These consisted of an autoclaved soil matrix that was inoculated with a plant growth-promoting soil (South Brush Creek [SBC], inoculated treatment) or left to allow microbial community recolonization by atmospheric deposition (disrupted treatment). Rhizosphere communities were sampled between 2 and 14 days of growth across six replicates (spatial blocks) and across diel cycles. Using an original approach based on generalized linear and Bayesian multilevel modeling, we identified 16S rRNA amplicon sequence variants (ASVs), ASV linear combinations, or clusters of cooccurring ASVs (network modules) whose abundances best predict leaf area, plant biomass, and growth rates across early or late stages of growth. Furthermore, we characterized cooccurrence patterns of growth-associated ASVs that may reveal their potential for plant growth promotion in field studies.

## RESULTS

**Soil inoculum promotes *Brassica rapa* growth.** *B. rapa* was grown in two soil settings to examine plant growth promotion (PGP) under different modes of community assembly, in an autoclaved soil favoring a more random assembly through migration and colonization of air/water-sourced microorganisms (disrupted treatment) or in an autoclaved soil followed by inoculation of a soil slurry favoring a more deterministic assembly (inoculated treatment). Compared to plants grown in disrupted soil treatments, inoculation of the *B. rapa* accession Wisconsin FastPlant (FPsc) with SBC soil resulted in
<figure><img src="image_3.png" /><figcaption>**FIG 2** Rhizosphere bacterial community diversity. (A and B) Alpha diversity (number of observed ASVs and Shannon indices) (A) and PCoA ordination of Bray-Curtis dissimilarities (B) between raw soil inoculum (SBC) and rhizosphere communities by treatment. Separate ordinations of disrupted (C) and inoculated (D) communities are shaded by the number of days of plant growth, with harvest times indicated by shape. (E) Bray-Curtis dissimilarities across the time series. In each treatment, dissimilarities from communities harvested at each day are shown relative to those from day 2. (Dissimilarities at day 2 were calculated from day 2 samples only.) Raw read counts were transformed using cumulative sum scaling.</figcaption></figure>

measurements averaged on days 10 and 11 is shown in Fig. S2 at https://doi.org/10.6084/m9.f figshare.c.5787851.v1. Block differences in biomass and relative growth rates showed less consistent patterns. Concentrations of bacterial 16S rRNA gene copies increased with time until around day 8 but did not significantly vary across treatments (Fig. 1C). Plants showed no indication of disease.

**Inoculated rhizospheres assemble deterministically.** Inoculated rhizosphere communities showed higher richness and evenness compared to disrupted ones, as measured by numbers of observed ASVs and Shannon's diversity indices (Fig. 2A). SBC inoculum communities were even more diverse. The top 100 ASVs from SBC soil inoculum communities included 15 bacterial classes in near-equal proportions, consisting of around 30% of sequence reads (Fig. S3A at https://doi.org/10.6084/m9.f figshare.c.5787851.v1). In contrast, 2 days after germination, the 100 most abundant ASVs in rhizosphere communities of both treatments consisted of only four classes (*Actinobacteria*, *Bacilli*, *Bacteroidia*,
<table><caption>TABLE 1 Co-occurrence network summary statistics for six sample groupsª</caption><thead><tr><th>Sample group</th><th>Count</th><th>No. of vertices</th><th>Mean degree</th><th>Mean strength</th><th>No. of edges</th><th>Mean wt</th><th>No. of modules</th></tr></thead><tbody><tr><td>All inoculated</td><td>144</td><td>718</td><td>19.98</td><td>0.798</td><td>7,173</td><td>0.040</td><td>56</td></tr><tr><td>All disrupted</td><td>144</td><td>548</td><td>9.26</td><td>0.620</td><td>2,537</td><td>0.067</td><td>54</td></tr><tr><td>Early inoculated</td><td>60</td><td>614</td><td>4.78</td><td>0.108</td><td>1,468</td><td>0.022</td><td>222</td></tr><tr><td>Late inoculated</td><td>60</td><td>787</td><td>13.88</td><td>0.652</td><td>5,463</td><td>0.047</td><td>70</td></tr><tr><td>Early disrupted</td><td>60</td><td>423</td><td>5.35</td><td>0.405</td><td>1,132</td><td>0.076</td><td>83</td></tr><tr><td>Late disrupted</td><td>60</td><td>469</td><td>6.49</td><td>0.315</td><td>1,522</td><td>0.048</td><td>59</td></tr></tbody></table>

and *Gammaproteobacteria*) and contained roughly 80% of sequences (Fig. S3B at https://doi.org/10.6084/m9.fishshare.c.5787851.v1). Neither metric of alpha diversity at the community level varied notably across time. Treatment accounted for 22% of the variation in rhizosphere community structure (permutational multivariate ANOVA [PERMANOVA], *P* < 0.001), an observation reflected by principal coordinate analysis (PCoA) ordinations of Bray-Curtis dissimilarities (Fig. 2B). Days of growth across the 14-day time series explained 5.8% and 16.8% of community variation in disrupted and inoculated communities respectively (PERMANOVA, *P* < 0.001 for both), and separation between early and late stages of growth was more evident in inoculated communities (Fig. 2C and D). In rhizospheres intensively sampled at five times per day during two 48-h cycles (3 to 4 and 13 to 14 days), harvest times did not influence community structure across either treatment (Fig. 2C and D).

Several measures point toward deterministic community assembly in rhizospheres of inoculated treatments. Rhizosphere bacterial communities from disrupted treatments exhibited high but consistent dissimilarity throughout the time series, while inoculated communities showed low but increasing dissimilarity (relative to day 2 samples) until day 10 (Fig. 2E). Abundant ASVs were more highly prevalent (present in a higher percentage of communities) across inoculated communities than disrupted ones (Fig. S4 at https://doi.org/10.6084/m9.fishshare.c.5787851.v1). The core community of the inoculated treatment consisted of 46 ASVs present across all 144 samples and accounting for 51% of all sequence reads, while only 8 ASVs (36% of reads) were present in all 144 disrupted communities. Cooccurrence networks were constructed from rhizosphere communities belonging to the same treatment and time stage (here referred to as sample groups). Early groups correspond to rhizospheres at 3 to 4 days of growth and late groups to 13 to 14 days of growth, chosen to represent the increased sampling frequency within these 48-h intervals. Networks from all inoculated rhizospheres contained more connections, and vertices (ASVs) were more highly interconnected than in all the disrupted ones (Table 1). The mean number (degree) and strength of connections also increased conspicuously between early and late networks of inoculated communities only. Numbers of modules, or clusters of cooccurring ASVs in each sample group, ranged from 54 to 222 (Table 1).

Rhizosphere communities from both treatments exhibited spatial variation across the six greenhouse blocks, with block membership explaining 14% and 17% of variance among inoculated and disrupted communities, respectively (PERMANOVA, *P* < 0.001 for both). Differential abundance analysis (DESeq2 [33]) identified 33 ASVs whose abundances varied significantly across 4 of 6 blocks in inoculated communities (Fig. S5 at https://doi.org/10.6084/m9.fishshare.c.5787851.v1) and 49 ASVs across 5 of 6 blocks in disrupted communities (Fig. S6 at https://doi.org/10.6084/m9.fishshare.c.5787851.v1). These block-specific biomarkers made up highly variable proportions of their respective communities and included several dominant ASVs, particularly in the genera *Flavobacterium* and *Paenibacillus*. In both treatments, block 6 had highest numbers and abundances of block-

<footnote>ªCounts display the number of communities used for constructing networks, and vertices are the number of ASVs in each network. Degree and strength are attributes of vertices that represent the number of connections and sum of edge weights from each vertex. Edges are the number of co-occurrences in each graph, and weights are attributes of edges that represent cooccurrence strength. Modules of co-occurring ASVs are shown for each graph after removing negative edges.</footnote>
<table><caption>TABLE 2 Associations of rhizosphere communities and plant growth data tested in the GLM<sup>a</sup></caption><thead><tr><td>Time stage</td><td>Response</td></tr></thead><tbody><tr><td>Early (days 3 and 4)</td><td>Leaf PA, day 3</td></tr><tr><td>Early</td><td>PA residuals from days 3 and 4</td></tr><tr><td>Early</td><td>Biomass residuals from days 3 and 4</td></tr><tr><td>Late (days 13 and 14)</td><td>Avg PA, days 10 and 11</td></tr><tr><td>Late</td><td>PA, day 13</td></tr><tr><td>Late</td><td>PA residuals from days 13 and 14</td></tr><tr><td>Late</td><td>Biomass residuals from days 13 and 14</td></tr><tr><td>Late</td><td>RGR</td></tr><tr><td>All, from 3 days onwards</td><td>RGR</td></tr></tbody></table>

specific biomarkers (Fig. S5 and S6 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). ASVs whose abundances varied by growth stage (early or late, harvested before or after 7 days of growth) made up similar proportions of inoculated and disrupted rhizosphere communities (Fig. S7A at https://doi.org/10.6084/m9 figshare.c.5787851.v1). However, higher total numbers of growth-stage biomarker ASVs in inoculated communities likely reflect the higher alpha diversity (Fig. 2A). ASVs whose abundance peaked at one particular day of the time series made up smaller proportions of both treatments than growth-stage biomarker ASVs (Fig. S7B) but were more abundant in disrupted rhizospheres, recapitulating their higher degree of stochastic assembly. All ASVs that varied in differential abundance across treatment, block, treatment × time (growth stage or specific day), or treatment × block are presented as a table of logical values (Table S1 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). In both treatments, *Flavobacterium* represented the most diverse and the most abundant genus within rhizosphere communities, containing 630 ASVs and over 40% of reads. Overall, reads assigned to *Flavobacterium* increased in abundance with time in disrupted communities but decreased in inoculated ones (Fig. S8 at https://doi.org/10.6084/m9 figshare.c.5787851.v1).

**Rhizobacterial features associated with plant growth.** We used generalized linear modeling (GLM) to identify associations between normalized abundances of ASV or modules (here features) and numerical plant growth measurements while controlling for block heterogeneity. A least absolute selection and shrinkage operator (LASSO) was applied to refine the number of important, nonredundant predictive features. Rhizosphere communities were classified into sample groups based on inoculation treatment and time (early, 3 or 4 days growth; late, 13 or 14 days growth). GLM runs consisted of ASV or module abundances from these communities paired with contemporaneous plant growth data (Table 2). Temporal variation in RGR values between days 3 and 4 (Fig. S1 at https://doi.org/10.6084/m9 figshare.c.5787851.v1) prevented us from associating RGR with features in early-stage communities. After feature selection, Bayesian multilevel modeling (of individual features or linear combinations thereof) reduced the number of growth-associated ASVs or modules to only those that were most unambiguously associated with plant growth (i.e., nonzero slopes of feature abundances versus plant growth; Table S2 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). Further information is provided in Materials and Methods. This statistical framework allowed us to identify a total of 43 unique ASVs from rhizosphere communities that were good predictors of at least one measurement of plant growth across different time stages and treatments (here growth positive or growth negative) (Fig. 3). Twenty-six were associated with inoculated communities, 18 with disrupted communities, and only one (ASV117) with both. Abundances of individual ASVs explained from 5% to 40% of the variance in plant growth data (Table S2), while models containing linear combinations of multiple ASVs explained from 18% to 84% of variance in growth (conditional R<sup>2</sup>) (Table S3 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). Seven growth-associated ASVs were associated with more than one growth metric, four of

<footnote><sup>a</sup>These associations were run for both inoculated and disrupted soil treatments and with ASV and module abundances. RGR measurements were consistent across the time series (see Fig. S1 at https://doi.org/10.6084/m9 figshare.c.5787851.v1), allowing us to include communities from across the time series (excluding those from day 2, when only one PA measurement was available).</footnote>
<figure><img src="image_4.png" /><figcaption>**FIG 3** ASVs associated with plant growth in both treatments. ASVs were associated with projected area (PA), residuals of aboveground biomass controlling for days grown, residuals of PA controlling for days grown and sample collection date, and relative growth rates (RGR). Measurements from early and late columns denote rhizosphere communities from plants grown 3 to 4 or 13 to 14 days, respectively, while All denotes all plants grown >2 days. R² values show proportions of variance explained by Bayesian models of abundances of each ASV versus plant growth, with blue showing positive associations and red being negative. ASVs are labeled by most specific taxonomy information and shown in bold if also included in best multivariate Bayesian models. PA, residual, and RGR values represent averages from 3 to 7 plants (early stage) or 2 plants (late stage).</figcaption></figure>
A

Inoculated communities

<figure><img src="image_5.png" /></figure>

B

Disrupted communities

<figure><img src="image_6.png" /><figcaption>**FIG 4** Most abundant growth-associated ASVs across time. (A and B) Relative abundances of growth-associated ASVs from inoculated (A) or disrupted (B) communities harvested at times where they were positively or negatively associated with any metric of plant growth (early, days 3 and 4; late, days 13 and 14). ASVs with abundances of less than 0.2% of communities, on average, by treatment and time were omitted. Most specific taxonomies are shown. The y axis for positive growth-associated ASVs in panel B was truncated because ASV10 exceeded 15% abundance within three disrupted communities, with a maximum abundance of 36.6%.</figcaption></figure>

which were assigned to the genus *Flavobacterium* (Fig. 3). In general, ASVs associated with RGR explained less variance than those associated with other growth measurements. More traditional biomarker detection techniques, such as DESeq2 (33), applied to rhizosphere communities classified as growth-promoting, possibly growth-promoting, or non-growth-promoting categories did not return any significant results, justifying our GLM approach (for more detail, see Materials and Methods).

Eight ASVs belonging to the genus *Flavobacterium* were associated with growth across different treatments and time stages, several of which demonstrated particularly strong associations with growth (Fig. 3) and were highly abundant members of rhizosphere communities (Fig. 4). While most of these were growth positive, ASV37 was negatively associated with PA, biomass residuals, and PA residuals in early inoculated communities (Fig. 3). Flavobacterium ASV sequences did not cluster phylogenetically with growth associations in an obvious pattern, although three of the five unambiguously growth-positive ASVs grouped into a distinct subclade based on 16S rRNA V4V5 sequences (Fig. S9 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). Despite the high correlation (R² = 0.85) between PA values from inoculated plants at day 3 and residuals of PA from days 3 and 4, 11 additional ASVs were associated with PA residuals. Interestingly, ASV117, a *Flavobacterium* with a greater than 99% sequence identity to ASV37, was the only sequence that was associated both positively and negatively with growth and associated with growth in both inoculated and disrupted treatments. Five other genera (*Arcticibacter*, *Bacillus*, *Dyadobacter*, *Mucilaginibacter*, and *Paenibacillus*) also contained multiple ASVs that showed either positive or negative associations with growth (Fig. 3). These were identified across different treatments, stages, and measurements of plant
growth but also sometimes within the same pairings of sample groups and plant growth data (*Flavobacterium* in early-stage inoculated treatments, *Bacillus* with RGR in disrupted treatments).

Overall, relative abundances of ASVs associated with any measurement of growth were higher in inoculated than disrupted communities, and relative abundances of all growth-positive ASVs were higher than those of growth-negative ones in both treatments (Fig. 4). In addition to *Flavobacterium*, these dominant growth-associated ASVs included members of *Cellvibrio*, *Cytophaga*, *Adhaeribacter*, and *Pedobacter*. All dominant growth-associated ASVs in disrupted communities were associated with RGR throughout the entirety of the time series. Relative abundances of growth-positive ASVs showed high variability in disrupted communities, with *Flavobacterium* ASV10 ranging from less than 0.1% to 36% of reads and *Fluviicola* ASV58 from 0 to 4.6% (Fig. 4).

We then searched for growth-associated modules as described for ASVs. By condensing thousands of ASVs into dozens or hundreds of modules, we were able to include ASVs with low abundances or prevalences that would otherwise have been omitted from the GLM (see Materials and Methods for more detail). We identified 10 module associations with growth, where all modules consisted of only one or two ASVs (Fig. 5). Five of these consisted of low-prevalence or low-abundance ASVs omitted from the ASV-centric approach, while others included highly abundant singleton modules, such as *Flavobacterium* ASV6 and *Adhaeribacter* ASV74. Three of the four modules consisting of higher-abundance singleton ASVs revealed associations to PA or PA residuals similar to those shown in Fig. 3, and the two *Arcticibacter* ASVs in module 9 (ASVs 68 and 149) (Fig. 5B) both showed negative associations with RGR as individuals (Fig. 3).

**Co-occurrence networks.** Co-occurrence networks of growth-associated ASVs, showing ASVs displaying strong direct positive or negative connections to them, revealed differences in connectivity across treatments (Fig. 6). The inoculated network (Fig. 6A) consisted of 328 edges connecting 155 ASVs from six modules that made up an average of 39% of the reads within communities, while the disrupted network (Fig. 6B) contained 76 edges connecting 70 ASVs from six modules that made up 19% of reads. Several dominant growth-associated ASVs from both treatments (particularly *Cellvibrio gandavensis* ASV14 and *Adhaeribacter* ASV74 in inoculated and disrupted communities, respectively) showed no strong connections and were removed from network graphs. Certain modules within both networks contained both growth-positive and growth-negative ASVs, although growth-associated ASVs within these modules were generally weighted toward one direction. Some ASVs in both networks showed positive co-occurrences with both growth-positive and growth-negative ASVs. These comprised 50 ASVs across 26 genera in inoculated communities and two ASVs from the genus *Paenibacillus* in disrupted rhizospheres. Considering all (not just strong) co-occurrences, growth-associated ASVs from inoculated communities were more highly connected than those from disrupted communities and non-growth-associated ASVs from inoculated communities (Fig. 6C). Nevertheless, growth-positive and growth-negative ASVs showed no differences in connectivity in either network, and neither did ASVs that were kept versus omitted in multivariate Bayesian models.

DISCUSSION

The use of two soil treatments, consisting of a sterilized soil matrix that was either amended with a raw soil inoculum (inoculated) or exposed to colonization by atmospheric deposition (disrupted), allowed us to examine the plant growth-promoting potential of a native soil bacterial community relative to a recent assemblage sourced from greenhouse air and/or water. While 16S rRNA gene copy numbers are not indicative of live bacterial biomass, similar increases in copy numbers over time in both treatments allow us to discount biomass as a confounding factor. Trial experiments demonstrated no discernible difference in PA or aboveground biomass between plants grown in disrupted or freshly autoclaved soils, supporting the role of an intact microbial community
<figure><img src="image_7.png" /><figcaption>**FIG 5** Modules associated with plant growth in both treatments. Modules were associated with PA (A) and RGR (B). Columns denote time stages when rhizosphere communities were sampled (early, 3 to 4 days; all, > 2 days). R² values show proportions of variance explained by Bayesian models of each individual module's abundance versus plant growth, with blue showing positive associations and red being negative. Modules consisting of only one ASV show its most specific taxonomy, and modules indicated with an asterisk denote ASVs that were omitted from the ASV-centered workflow based on low abundance or prevalence. Module labels are shown in boldface if also included in best multivariate Bayesian models. PA, residual, and RGR values represent the averages from 2 to 7 plants.</figcaption></figure>

in *B. rapa* growth promotion (data not shown). Maximum percent yield increases of 28 to 39% (leaf PA and aboveground biomass, respectively) in inoculated plants were within range of several inocula from agricultural and forest soils that largely increased *Arabidopsis* shoot and root biomass relative to axenic controls (34). Temporal changes in rhizobacterial community composition across stages of plant growth have been described previously (35, 36) and specifically attributed to changes in root exudate chemistry (37) or microbial succession independent of host growth stage (38). Although some exudates in *Arabidopsis* and other members of *Brassicaceae* (notably glucosinolates) exhibit circadian fluctuation, many do not (39, 40). Although glucosinolate

<figure><img src="image_8.png" /><figcaption>**FIG 6** Subsetted networks showing growth-associated ASVs and their strong co-occurrences. Networks show inoculated (A) and disrupted (B) communities, subsetted to include only growth-associated ASVs and those they show strong co-occurrences with. Positive co-occurrences are shown as gray edges for weights exceeding the third quartile of all weights from unsubsetted networks, and negative co-occurrences are shown as red edges for weights below the first quartile. Nodes labeled in blue are ASVs positively associated with growth, and nodes labeled in red are negative. Unlabeled nodes are connected to one or more growth-associated ASVs. Node size corresponds to the square-rooted average relative abundance, and node color indicates module membership. Edge thickness corresponds to connection weight. Node sizes and edge widths across networks are not to scale. (C) Degree (number of connections) from each ASV in unsubsetted networks of inoculated or disrupted communities, separated by whether they were associated with growth. \*, P < 0.05; \***, P < 0.001.</figcaption></figure>
exudation patterns in FPsc are uncharacterized, differences in growing conditions (such as between greenhouse and growth chamber experiments) may partially explain the consistency we observed in rhizosphere community composition across diel cycles. The separation in community beta-diversity in both treatments between weeks 1 and 2 (Fig. 2C and D) occurred during thinning, which may be a result of density-dependent changes in root exudation between day 6 and later time points.

We interpret the contrasting patterns of beta-diversity and network connectivity of rhizosphere communities from disrupted and inoculated treatments as consequences of stochastic versus deterministic assembly. Although rhizosphere microbial community assembly is driven by both processes (41–43), the high dissimilarity in disrupted community beta-diversity over time (Fig. 2E) has also been noted in soils after wildfire disturbance (44). The corresponding reduction in plant growth associated with disturbance (sterilization) and recolonization highlights the importance of deterministic processes of rhizosphere community assembly for agricultural contexts.

In addition, our observation of several ASVs with differential abundance patterns across greenhouse blocks in both treatments suggests that random or dispersal forces also shaped rhizosphere community assembly, as previously noted in *Arabidopsis* phyllospheres and across meter-scale field plots in rhizospheres of several *Brassica* cultivars (26, 27). Alternatively, subtle and unmeasured variations in temperature or other factors across different areas of the greenhouse may play a role (particularly for block 6, which sat directly underneath a heater; see Fig. S5 and S6 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). It is important to note that whole-rhizosphere sampling as described here may obscure fine-scale spatial distributions of bacterial populations along the root axis, such as those recently linked with root flavone exudation during nitrogen deprivation (22).

Roles of growth-associated rhizobacteria may be direct (such as aiding nutrient uptake via N, P mineralization, or in competing with plants for these nutrients) or indirect, such as through hormonal suppression of root growth-inhibiting bacteria by *Variovorax* (10). Random forest associations have been used to identify bacterial OTUs, genera, and taxonomic orders positively and negatively associated with grain weight in foxtail millet (31), growth of tomato plants (45), and soybean yields (46). In these studies, the proportion of variance (R²) attributable to combinations of these identified features ranged from 0.31 to 0.38, which was on the lower end of the variance explained by our multivariate Bayesian models containing multiple features (mean R² = 0.4; Table S3 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). Although our decision to select features at the ASV level risks dispersing slightly different 16S rRNA gene sequences that may originate from the same genome into separate entities (47), this approach addresses strain-level heterogeneity that may provide insight into variability in rhizosphere functional potential (48–50). Connectivity patterns allow us to at least partially resolve this issue: despite sharing 99.2% sequence identity, two *Flavobacterium* ASVs (6 and 37) that negatively cooccur with each other and have different effects on plant growth (Fig. 6A) likely represent different strains, while two positively cooccurring *Arcticibacter* ASVs (68 and 149, with 99.3% sequence identity) that belong to the same module (Fig. 6B) likely arise from the same genome.

While our approach primarily focused on abundant members of rhizobacterial communities, rare members of the rhizosphere can contribute to maintenance of plant health, such as in producing compounds to suppress fungal pathogens (51). The higher statistical power inherent to more prevalent features (those detected in higher percentages of communities) likely allowed them to be more easily identified in GLM LASSO runs. Nevertheless, clustering ASVs into modules allowed us to detect several associations between plant growth and lower-abundance rhizosphere members.

Many growth associations we found were specific to early or late stages of plant growth only, as others have observed in PGPR isolates from banana rhizospheres (52). In particular, we identified seven early-stage biomarkers associated with growth within the time stages in which they were most abundant (ASVs 6, 37, 99, 113, 105, and 249).
Contemporaneous associations in late stages or disrupted treatments were less clear and identified lower-abundance ASVs overall, suggesting that early stages of plant growth are particularly influenced by members of diverse or stable rhizosphere communities. Despite the considerable spatial heterogeneity across rhizobacterial communities (Fig. S5 and S6), only four growth-associated ASVs varied significantly across greenhouse blocks, of which only ASV13 (genus *Pedobacter*) exceeded 0.2% average abundance. The overall lack of variability across spatial blocks in growth-associated features appears to reflect deterministic assembly processes (even in disrupted communities), which may be driven by selection in response to root exudation patterns.

Imperfect correlations between different growth metrics may account for differences in features with which they were associated. In particular, PA from inoculated plants at day 3 and PA residuals at days 3 and 4 showed an $R^2$ of 0.85 and were associated with ASVs 6 and 37, yet several additional ASVs were associated with PA residuals only (Fig. 3). This may highlight the sensitivity of the GLM in selecting features from similar or highly correlated data types, making the choice of growth data a critical consideration. Although most growth-associated features detected in univariate Bayesian models were also significant in multivariate models, the ones that were significant only in a multivariate context (Table S2 at https://doi.org/10.6084/m9 figshare.c.5787851.v1) may only appear so when controlling for the abundances of other ASVs or may require each other to synergistically affect growth.

Members of *Flavobacterium* have been associated with increased potato root and shoot growth in one of the first high-throughput sequencing studies of the rhizosphere (29). Plant-associated members of this genus are capable of gliding motility, possess enzymes to degrade glycans from plant cell walls, and can synthesize phytohormones to stimulate plant growth (53). Interestingly, a recent study associated abundances of *Flavobacterium* in wheat rhizospheres with fungal pathogen infection despite noting one outlier *Flavobacterium* OTU associated with healthy plants (54). While saprophytic members of *Cellvibrio* are capable of degrading diverse plant cell wall polysaccharides, at least one member can fix nitrogen (55, 56). Intriguingly, genera with multiple growth-negative ASVs identified in this study (*Adhaeribacter* and *Pedobacter*) have been associated with growth promotion in other rhizosphere systems (54, 57, 58). Within and across studies, the variable influences of rhizobacterial genera on plant health and growth appear to encompass a diversity of context-dependent host-bacterium interactions at finer taxonomic scales.

While treatment dramatically shaped co-occurrence patterns in rhizosphere communities, we detected no differences in connectivity between growth-positive or growth-negative ASVs, in contrast to Jin et al. (31). Different patterns of co-occurrences between growth- and non-growth-associated ASVs across treatments (Fig. 6C) suggests that higher connectivity is associated with growth promotion at the community level, although growth-positive and growth-negative ASVs from both treatments that did not show strong co-occurrences appear to be exceptions to this rule. We suggest that growth-positive ASVs with few strong co-occurrences represent promising candidates for isolation and inoculation to improve plant yields. In contrast, bioinoculation with features showing strong co-occurrences may require characterizing their mechanisms of interaction, designing synthetic communities, and/or directly inoculating seeds with growth-promoting soils.

In summary, we present differing effects of stochastic and deterministic rhizobacterial community assembly processes during a time series of *B. rapa* inoculated with a live or a knocked down soil treatment, noting spatial and temporal heterogeneity. By using a regularized regression approach, we identify correlations between relative abundances of ASVs or modules and several measures of plant growth in early or late stages of *B. rapa* development. Although the associations we report are correlative, different tests identified many of the same features. These associations provide hypotheses for laboratory and field studies on whether and how these rhizobacteria promote plant growth and associate with other members of rhizosphere communities. Because
many rhizobacteria remain uncultured or resistant to cultivation, culture-independent approaches such as genome-resolved metagenomics and metatranscriptomics can be used to identify adaptation strategies of recalcitrant populations in plant-associated systems (59). Further investigations should include comparative genomic investigations to characterize variability in growth association patterns among members of influential rhizosphere genera and link changes in rhizosphere communities with plant gene expression and root exudation patterns. A more complete understanding of spatial, temporal, and taxonomic heterogeneity in the rhizosphere will provide a valuable bioengineering foundation for a diversity of agricultural contexts.

## MATERIALS AND METHODS

**Greenhouse experiments.** To explore the timing of rhizosphere formation and its influence on plant growth, we initiated a time course greenhouse experiment at the University of Wyoming in December of 2018 (Laramie, WY; GPS location, 41.319787, -105.558253; elevation, 2,149 m). We planted 204 replicate pots of each of two treatments (disrupted control versus inoculated) with the *Brassica rapa* self-compatible Wisconsin FastPlant (FPsc; Rick Amasino and Scott Woody, University of Wisconsin). Pots in both treatments were filled with a soil matrix consisting of 50% (vol/vol) fritted clay (Profile porous ceramics greens grade; Profile, Buffalo Grove, IL, USA) to 50% sphagnum potting mix (Redi-Earth; Sungro Horticulture, Agawam, MA, USA) amended with powdered organic alfalfa meal (3.5 mL/liter) and adjusted to a pH of 7.0 with powdered lime. The soil matrix was sterilized with two 1-h autoclave cycles separated by 24 h, which has been shown to prevent amplification of extracted DNA (26). For the control (disrupted) treatment, autoclaved bags of soil matrix were opened to atmospheric inoculation in the greenhouse 72 h prior to planting seeds, allowing for colonization by the ambient microbial community. The day before planting, the inoculated treatment was initiated by mixing 2 mm sieved soil from a Wyoming site (South Brush Creek, or SBC; GPS location, 41.327404, -106.502784; elevation, 2,570 m) with sterile autoclaved soil matrix (5%, vol/vol, SBC to matrix). This site is a disturbed montane forest plot where *Boechera stricta* (*Brassicaeae*) is common. Inoculation experiments with *Boechera stricta* and *Brassica rapa* have previously demonstrated significant plant growth promotion using this microbial community (C. J. Hubbard, R. McMinn, M. T. Brock, and C. Weinig, unpublished data).

Pots (8.9 by 8.9 cm, 500 mL) were filled with either control or inoculated soil. The following day, FPsc seeds were sterilized by vortexing for 1 min in 70% ethanol (EtOH), 10 min in 10% bleach, and 5 rinses in sterile water, followed by planting in a 4-by-4 grid (for pots randomly assigned to early sampling days) or a 3-by-3 grid (for pots randomly assigned to late sampling days) and topped with autoclaved vermiculite. Pots with a maximum of 16 seedlings were necessary to attain sufficient rhizosphere mass at early growth stages. Pots were assigned unique locations in alternating cells of 3-by-6 trays, which were arranged into 6 spatial blocks. Greenhouse ambient temperatures were set to 20°C and 15.5°C (day/night), and supplemental LED lights were set to 14-h days (starting at 0500). Germination was first noticed 4 days after planting, which we used as time zero (70; 55% and 33% of seeds germinated 4 and 5 days after planting, respectively).

Six pots of each treatment (1 randomly selected replicate per treatment from each of six spatial blocks) were sampled in the afternoon (1400 h) of days 2, 3, 4, 6, 8, 10, 13, and 14 (Fig. 1). Nested within this 2-week sampling period, we also performed 48-h diel sampling events of 6 pots per treatment at early (days 3 and 4) and late stages (days 13 and 14) at the following additional time points over both days: 0900, 1800, 2200, and 0400. Seedlings within a pot were split between harvesting the rhizosphere (described below) or harvesting root tissue for transcriptome sequencing (RNA-seq) (not presented here). Pots destined for sampling later in the experiment were thinned down after 7 days to a final density of two plants per pot.

To harvest rhizospheres, we gently divided pots in half and pulled intact seedlings and roots from the bulk soil. Roots and rhizospheres were gently shaken to remove residual bulk soil, and then larger soil clumps of >2 mm were removed by hand. Gloves were changed and sterilized with 70% EtOH between handling pots. Roots plus rhizospheres were then placed in microcentrifuge tubes containing DNA/RNA Shield (Zymo Research Corporation, Orange, CA, USA) and placed on a vortexer for gentle disruption of the rhizosphere and mixing with DNA/RNA Shield. Samples were placed in a -80°C freezer until 16S amplicon preparation.

**Plant growth measurements.** To estimate average plant leaf area throughout the experiment, trays of pots were photographed (almost daily throughout the time series) using a dSLR camera at a fixed height on a copy stand. Resulting images were corrected for lens distortion in Photoshop CS6 (Adobe, San Jose, CA, USA) and then cropped to individual pots using custom Perl scripts and ImageMagick (60). Images of individual pots were batch processed with Easy Leaf Area (61) to extract the area of green pixels as our proxy for plant size after dividing green pixels per pot by the number of seedlings present. The relationship between this calculated projected area (PA) and the true leaf area should diminish over time, as plant architecture and leaf self-shading cannot be corrected from a single photo; however, correlations between average PA and average biomass remained highly significant on day 14 (r = 0.67; *P* < 0.0001). In addition to PA, we retained aboveground tissue when rhizospheres were harvested throughout the experiment. Tissue samples from each pot were dried for 48 h and weighed, and seedling number was tallied, allowing us to determine average aboveground biomass. Biomass and PA values were reported as means from all plants per pot. We calculated residuals for aboveground biomass
and PA to control for differences between days 3 and 4 or 13 and 14 among treatments. Relative growth rates (RGR) were calculated by log transforming PA measurements and calculating slopes from consecutive intervals on a per-pot basis.

**DNA extraction and sequencing.** Genomic DNA templates were prepared using the ZymoBIOMICS DNA/RNA Miniprep kit protocol according to the manufacturer's instructions (Zymo Research Corporation) using a parallel purification to extract RNA. We followed 16S rRNA gene amplification methods as previously described (62, 63). The starting master mix for the PCR contained 1× SuperFi buffer, 200 µM deoxynucleoside triphosphates, 2.5 U of SuperFi polymerase, and molecular biology-grade water. Before adding template DNA, we added individually indexed fusion primers (0.3 µM) to amplify the bacterial V4-V5 region of the 16S gene (*Escherichia coli* nucleotides [nt] 518 to 926, 408 nt). Primer sequences are described in reference 62. This volume was mixed by pipetting and 25 µL removed to serve as the negative control. One to 10 ng of template DNA was added to the remaining volume and mixed in by pipetting. The 100 µL was divided into triplicate reaction mixtures. Amplification conditions were an initial denaturation at 94°C for 3 min; 30 cycles of 94°C for 30 s, 57°C for 45 s, and 72°C for 1 min; and a final extension at 72.0°C for 2 min. Triplicate PCRs were pooled and the products cleaned using 0.75 volume of AMPure XP beads (Beckman Coulter) to 1 volume of sample.

Amplicon concentrations were quantified using the PicoGreen assay (ThermoFisher Scientific) and pooled in equimolar amounts based on PicoGreen results. Pools were size selected using BluePippin (Sage Scientific) to a range of 425 to 625 bp on a 1.5% gel cassette. The size-selected products were cleaned again using AMPure XP at a 1:1 ratio and the final pools visualized on an Agilent Bioanalyzer using the DNA1000 assay. The pool concentration was calculated using quantitative PCR (KK4835; KAPA library quantification kit) and the libraries sequenced on an Illumina MiSeq using the version 3 protocol to generate paired-end 300-nt reads.

**Sequence processing and diversity analyses.** DADA2 (64) was used to filter raw 16S fastq reads, infer ASVs, and remove chimeric sequences in a workflow similar to that of Callahan et al. (65). Taxonomy was assigned using version 138 of the SILVA nonredundant (SSU Ref NR99) 16S database (66). ASVs assigned to Eukaryota, mitochondria, and chloroplasts, as well as those that could not be classified at the domain level, were removed using the subset_taxa function in phyloseq (67). Before removal, mitochondria amounted to 2.5% of reads and chloroplasts to 0.7%. Next, contaminant ASVs were removed using the “combined” method of the R package decontam (68), which identified 77 sequence variants that were more abundant across three PCR blanks and in prepooled libraries with lower amplicon concentrations. After removal, the three PCR blanks (2,653 to 4,828 reads each) were pruned using phyloseq, leaving 20,961 ASVs in 293 communities, each sequenced at a depth of 25,438 reads. Alpha diversity calculations were performed using built-in phyloseq functions. Bray-Curtis dissimilarities calculated from cumulative sum scaling transformations of raw read counts (69) were used for PERMANOVA tests and ordinations implemented with the vegan R package (70). A maximum-likelihood phylogenetic tree was constructed in phangorn (71) using the Jukes-Cantor model and 100 bootstrap replications.

**Differential abundance analysis.** A first-order attempt to identify ASVs associated with plant growth promotion yielded no differentially abundant sequences across rhizosphere communities from inoculated-treatment plants based categorically on whether they were much larger, slightly larger, or no larger than their disrupted-treatment counterparts across the entire time series. DESeq2 (33) was used to identify differentially abundant ASVs across soil inoculation treatments and blocks using a threshold *P* value of 0.05 adjusted for multiple comparisons. To search for differentially abundant ASVs across inoculated plants (categorized as much larger, slightly larger, or no larger than their disrupted-treatment counterparts), we calculated aboveground biomass and PA at harvest values of inoculated plants as a percentage of their corresponding disrupted-treatment time point. This yielded three categorical groups: promoted (greater than 107%, the mean normalized score of the inoculated group, *n* = 23), possibly promoted (100 to 107%, *n* = 101), and not promoted (below 100%, *n* = 19). No differentially abundant ASVs were detected across these three groups (comparing communities from each individual group to those in the other two, using DESeq2 with an adjusted *P* value cutoff of 0.05).

**Network analysis.** We constructed cooccurrence networks of the most abundant and prevalent ASVs for six sample groups using SPIEC-EASI (72). These groups, split by treatment, consisted of all samples across the time series or those from days 3 and 4 (early) or 13 and 14 (late). ASVs were included for network analysis if they were detected in at least four communities and their mean abundance was greater than 0.5% within their corresponding sample group (e.g., early disrupted). This resulted in networks that contained 423 to 787 ASVs (Table 1). Associations between covarying ASVs were inferred with SPIEC-EASI using the “glasso” probabilistic inference method and a lambda.min.ratio of 0.01, and significant covariances were converted to correlations with the function cov2cor in a method similar to that of Lemonnier et al. (73). To define modules of highly correlated ASVs, negative edges of networks were removed and the Louvain algorithm (74) was implemented using the cluster_louvain function in igraph (75). igraph was also used to plot and subset networks (without negative edges removed) within customized R functions. Subsetted networks shown in Fig. 6 contain edges with weights below the first quartile value (strongly negative co-occurrences) or edges with weights above the third quartile (strongly positive co-occurrences). First and third quartile values were -0.073 and 0.152 for the inoculated network and 0 and 0.119 for the disrupted network, respectively.

**Feature selection.** We used regularized regression to identify associations between ASV abundances in rhizosphere communities and corresponding plant growth data. This consisted of a generalized linear model (GLM) with a penalized maximum likelihood paired with a least absolute shrinkage and selection operator (LASSO) used to perform variable selection and minimize spurious correlations and overfitting (76). Due to the sparsity in read count data and the low statistical power associated with low-
prevalence features, we restricted the GLM input to include only the 200 most abundant ASVs in each sample group that were detected in at least 30% or 50% of communities for disrupted or inoculated treatments, respectively. These prevalence cutoffs were determined heuristically after inspecting rank versus prevalence and cumulative abundance plots of most abundant ASVs across all sample groups (see Fig. S4 at https://doi.org/10.6084/m9 figshare.c.5787851.v1). While imposing different prevalence cutoffs to vary the number of input features, we noticed that a minority of runs returned different num- bers of features at different prevalence thresholds. For this reason, we additionally imposed stricter prevalence cutoffs of 50% or 90% for disrupted or inoculated sample groups, respectively, which reduced the number of input features from 200 to 120 to 130 (ASVs with prevalence below 70% were almost never identified, regardless of treatment). The top 200 ASVs in each sample group comprised 78% to 90% of the reads in these communities. Modules of cooccurring ASVs were run with the same growth data as the ASV data without subsetting by preva- lence or abundance, because each sample group consisted of fewer than 200 modules. (In the one exception, early inoculation, we omitted 22 modules that were each present in less than 27% of samples and totaled an av- erage of 0.2% of the reads in these communities).

The GLM was implemented with the cv.glmnet function of the R package glmnet (76) using a Gaussian distribution family within a customized R function. The customized R function was used to input raw read counts and corresponding plant growth data from samples (from the otu_table and sample_data of a phy- loseq object) as well as parameters to include features based on rank and prevalence and in specifying the number of lasso runs. Read counts were transformed using a centered log ratio (77), and reads of low- prevalence or low-abundance ASVs were summed together to preserve data compositionality.

We found slight variation among the ASVs identified with the GLM LASSO. We suspect that the highly correlated ASV abundances make each run sensitive to initial model conditions. To account for the variation among model runs, we ran each GLM LASSO model 100 times. Features with nonzero coefficients in at least 80 of 100 runs were considered "selected" and further evaluated in our Bayesian hierarchical models (described below). For both ASVs and modules, an optimal number of features corresponding to lambda values that minimized mean squared error during cross-validation was found to be 7 or 8 (for PA and RGR measurements) or 2 or 3 (for PA or biomass residuals).

Using the brms R package (78), Bayesian models of plant growth were constructed using abundan-
ces of individual features or combinations of features identified by the GLM. Our rationale for using an
additional Bayesian model after the GLM LASSO was due to the constraints and benefits of glmnet: the
glmnet function allowed us to efficiently evaluate different ASVs on plant growth but does not allow us
to include the block structure in our growth experiments. Bayesian models allowed us to incorporate
the block structure of our growth experiments using a hierarchical linear model with a random intercept
and slope for each block the plants grew in. A Gaussian family distribution was used to fit response (plant
growth) data, and features whose regression coefficient posterior distributions (95% credible intervals) did
not include zero were considered selected. Bayesian multivariate models were compared and refined using
a leave-one-out cross-validation approach with the loo R package (79). In the ASV data set, 89 ASV growth
associations from 73 unique ASVs were reduced to 59 after univariate Bayesian modeling, yielding 30 posi-
tive associations with growth and 29 negative. Likewise, the number of module-growth associations was
winnowed from 12 to 10 after these respective steps, resulting in six positive and four negative associations.
Twelve ASVs and one module were identified as significant by multivariate Bayesian modeling but not indi-
vidually; on the contrary, 23 ASVs and one module were detected with univariate Bayesian modeling only
(Fig. 3 and 5). All features and corresponding sample groups and growth metrics identified across different
steps are included in Table S2 (https://doi.org/10.6084/m9 figshare.c.5787851.v1). Univariate Bayesian models
preferentially upheld associations between ASVs and growth data identified by the GLM using both prevalence
thresholds (14 of 22, 63.6%) compared to those found at only one prevalence (25 of 68, 36.7%). Multivariate
Bayesian models explained anywhere from 9% to 80% of the variance in growth (marginal R² values), varying
with the number of features used in the model. Several models accounting for relative growth rates explained
less variance than other growth measurements. The percentages of variance explained by different models did
not vary depending on whether a feature was a module or an ASV or by sample group. Model summary statis-
tics are included in Table S3 (https://doi.org/10.6084/m9 figshare.c.5787851.v1).

To examine the specificity of the GLM at identifying spurious features, we randomly scrambled asso-
ciations between plant growth data and feature abundances and reran all tests. In doing so, we identi-
fied 41 associations compared to the 141 identified by the GLM in the unscrambled data. After following
up with univariate Bayesian modeling, this number was reduced to 33 (compared to 69 in the
unscrambled data). When scrambled, four tests that identified features in unscrambled data identified
more than one feature, corresponding to a false discovery rate of over 1% (Table S4 at https://doi.org/10
.6084/m9 figshare.c.5787851.v1). Thus, we conservatively decided to remove all associations from these
tests that were detected in the unscrambled data.

**Data availability.** Supplemental figures and tables are publicly available as a FigShare collection at https://doi.org/10.6084/m9 figshare.c.5787851.v1. Reproducible code for all analyses, including customized R functions and plant growth data files, are freely available as an archived GitHub repository at https://loimai.github.io/rhizosphere_rapa_16S. Raw fastq sequences of 16S amplicons were deposited to the NCBI sequence read archive (SRA) database under BioProject accession number PRJNA765038.

ACKNOWLEDGMENTS

We thank Nicole Robichaud and Victoria Frazier for conducting droplet digital PCR
measurements and Kirsten Peramba for DNA extraction and sequencing. Joshua G.
Harrison assisted with statistical interpretations.
This study was funded by NSF grant IOS-1444571. NSF had no role in study design, implementation, data collection or interpretation, or the decision to publish the manuscript.

M.T.B., C.W., and L.M. designed the study. M.T.B. grew plants, obtained plant growth data, and sampled rhizospheres. H.G.M. oversaw DNA sequencing and sequence data curation. W.J.C. contributed code and statistical analysis, and S.A.K. contributed code, conducted analyses, and wrote the manuscript, with assistance from M.T.B., C.W., H.G.M., and L.M.

REFERENCES

1. Marilley L, Aragno M. 1999. Phylogenetic diversity of bacterial communities differing in degree of proximity of Lolium perenne and Trifolium repens roots. Appl Soil Ecol 13:127-136. https://doi.org/10.1016/S0929-1393(99)00028-1

2. Bulgarelli D, Rott M, Schlaepki K, Ver Loren van Themaat E, Ahmadinejad N, Assenza F, Rauf P, Huettel B, Reinhardt R, Schmelzer E, Peplies J, Gloeckner FO, Amann R, Eickhorst T, Schulze-Lefert P. 2012. Revealing structure and assembly cues for Arabidopsis root-inhabiting bacterial microbiota. Nature 488:91-95. https://doi.org/10.1038/nature11336.

3. Sasse J, Martinioa E, Northen T. 2018. Feed your friends: do plant exudates shape the root microbiome? Trends Plant Sci 23:25-41. https://doi.org/10.1016/j.tplants.2017.09.003.

4. Harbort CJ, Hashimoto M, Inoue H, Niu Y, Guan R, Rombolà AD, Kopriva S, Voges M, Sattely ES, Garrido-Oter R, Schulze-Lefert P. 2020. Root-secreted coumarins and the microbiota interact to improve iron nutrition in Arabidopsis. Cell Host Microbe 28:825-837. https://doi.org/10.1016/j.chom.2020.09.006.

5. Rengel Z. 2015. Availability of Mn, Zn and Fe in the rhizosphere. J Soil Science Plant Nutrition 15:397-409. https://doi.org/10.4067/S0718-95162015005000036.

6. Gunina A, Kuzyakov Y. 2015. Sugars in soil and sweets for microorganisms: review of origin, content, composition and fate. Soil Biol Biochem 90:87-100. https://doi.org/10.1016/j.soilbio.2015.07.021.

7. Santhanam R, Luu VT, Weinhold A, Goldberg J, Oh Y, Baldwin IT. 2015. Native root-associated bacteria rescue a plant from a sudden-wilt disease that emerged during continuous cropping. Proc Natl Acad Sci U S A 112: E5013-E5020.

8. Berendsen RL, Pieterse CMJ, Bakker PAHM. 2012. The rhizosphere microbiome and plant health. Trends Plant Sci 17:478-486. https://doi.org/10.1016/j.tplants.2012.04.001.

9. Hubbard CJ, Li B, McMinn R, Brock MT, Maignien L, Ewers BE, Kliebenstein D, Weinig C. 2019. The effect of rhizosphere microbes outweighs host plant genetics in reducing insect herbivory. Mol Ecol 28:1801-1811. https://doi.org/10.1111/mec.14989.

10. Finkel OM, Salas-González I, Castrillo G, Conway JM, Law TF, Teix-Eira PJPL, Wilson ED, Fitzpatrick CR, Jones CD, Dangl JL. 2020. A single bacterial genus maintains root growth in a complex micro-biome. Nature 587: 103-108. https://doi.org/10.1038/s41586-020-2778-7.

11. Vejan P, Abdullah R, Khadiran T, Ismail S, Nasrulhaq Boyce A. 2016. Role of plant growth promoting rhizobacteria in agricultural sustainability—a review. Molecules 21:573. https://doi.org/10.3390/molecules21050573.

12. Santos SS, Rask KA, Vestergård M, Johansen JL, Priemé A, Frøslev TG, González AMM, He H, Ekelund F. 2021. Specialized microbiomes facilitate natural rhizosphere microbiome interactions counteracting high salinity stress in plants. Environ Exp Bot 186:104430. https://doi.org/10.1016/j.envexpbot.2021.104430.

13. Haskett TL, Tkacz A, Poole PS. 2021. Engineering rhizobacteria for sustainable agriculture. ISME J 15:949-964. https://doi.org/10.1038/s41396-020-00835-4.

14. Hutchins DA, Jansson JK, Remais JV, Rich VI, Singh BK, Trivedi P. 2019. Climate change microbiology—problems and perspectives. Nat Rev Microbiol 17:391-396. https://doi.org/10.1038/s41579-019-0178-5.

15. Hatfield JL, Antle J, Garrett KA, Izauralde RC, Mader T, Marshall E, Nearing M, Philip Robertson G, Ziska L. 2020. Indicators of climate change in agricultural systems. Clim Chang 163:1719-1732. https://doi.org/10.1007/s10584-018-2222-2.

16. Molina-Romero D, Baez A, Quintero-Hernández V, CastañedaLucio M, Fuentes-Ramírez LE, Bustillos-Cristales MR, Rodríguez-Andrade O, Morales-García YE, Munive A, Muñoz-Rojas J. 2017. Compatible bacterial

mixture, tolerant to desiccation, improves maize plant growth. PLoS One 12:e0187913. https://doi.org/10.1371/journal.pone.0187913.

17. Pascale A, Proietti S, Pantelides IS, Stringlis IA. 2019. Modulation of the Root Microbiome by Plant Molecules: the Basis for Targeted Disease Suppression and Plant Growth Promotion. Front Plant Sci 10:1741.

18. Fitzpatrick CR, Copeland J, Wang PW, Guttman DS, Kotanen PM, Johnson MTJ. 2018. Assembly and ecological function of the root microbiome across angiosperm plant species. Proc Natl Acad Sci U S A 115:E1157-E1165.

19. Hubbard CJ, Brock MT, van Diepen LT, Maignien L, Ewers BE, Weinig C. 2018. The plant circadian clock influences rhizosphere community structure and function. ISME J 12:400-410. https://doi.org/10.1038/ismej.2017.172.

20. Lv Z, Dai R, Xu H, Liu Y, Bai B, Meng Y, Li H, Cao X, Bai Y, Song X, Zhang J. 2021. The rice histone methylation regulates hub species of the root microbiota. J Genet Genom https://doi.org/10.1016/j.jgg.2021.06.005.

21. Chaparro JM, Badri DV, Vivanco JM. 2014. Rhizosphere microbiome assembly is affected by plant development. ISME J 8:790-803. https://doi.org/10.1038/ismej.2013.196.

22. Yu P, He X, Baer M, Beirinckx S, Tian T, Moya YAT, Zhang X, Deichmann M, Frey FP, Bresgen V, Li C, Razavi BS, Schaaf G, von Wirén N, Su Z, Bucher M, Tsuda K, Goormachtig S, Chen X, Hochholdinger F. 2021. Plant flavones enrich rhizosphere Oxalobacteraceae to improve maize performance under nitrogen deprivation. Nat Plants 7:481-499. https://doi.org/10.1038/s41477-021-00897-y.

23. Voges M, Bai Y, Schulze-Lefert P, Sattely ES. 2019. Plant-derived coumarins shape the composition of an Arabidopsis synthetic root microbiome. Proc Natl Acad Sci U S A 116:12558-12565. https://doi.org/10.1073/pnas.1820691116.

24. Shimasaki T, Masuda S, Garrido-Oter R, Kawasaki T, Aoki Y, Shibata A, Suda W, Shirasu K, Yazaki K, Nakano RT, Sugiyama A. 2021. Tobacco root endophytic arthrobacter harbors genomic features enabling the catabolism of host-specific plant specialized metabolites. mBio 12:397-409.

25. Herron PM, Gage DJ, Arango Pinedo C, Haider ZK, Cardon ZG. 2013. Better to light a candle than curse the darkness: illuminating spatial localization and temporal dynamics of rapid microbial growth in the rhizosphere. Front Plant Sci 4:323. https://doi.org/10.3389/fpls.2013.00323.

26. Maignien L, DeForce EA, Chafee ME, Eren AM, Simmons SL. 2014. Ecological succession and stochastic variation in the assembly of Arabidopsis thaliana phyllosphere communities. mBio 5:e00682-13.

27. Klasek SA, Brock MT, Morrison HG, Weinig C, Maignien L. 2021. Soil microsite outweighs cultivar genotype contribution to Brassica rhizobacterial community structure. Front Microbiol 12:645784. https://doi.org/10.3389/fmicb.2021.645784.

28. Walsh CM, Becker-Uncapher I, Carlson M, Fierer N. 2021. Variable influences of soil and seed-associated bacterial communities on the assembly of seedling microbiomes. ISME J 15:2748-2715. p. https://doi.org/10.1038/s41396-021-00967-1.

29. Manter DK, Delgado JA, Holm DG, Stong RA. 2010. Pyrosequencing reveals a highly diverse and cultivar-specific bacterial endophyte community in potato roots. Microb Ecol 60:157-166. https://doi.org/10.1007/s00248-010-9658-x.

30. Lay C-Y, Bell TH, Hamel C, Harker KN, Mohr R, Greer CW, Yergeau É, St-Arnaud M. 2018. Canola root-associated microbiomes in the Canadian prairies. Front Microbiol 9:1188. https://doi.org/10.3389/fmicb.2018.01188.

31. Jin T, Wang Y, Huang Y, Xu J, Zhang P, Wang N, Liu X, Chu H, Liu G, Jiang H, Li Y, Xu J, Kristiansen K, Xiao L, Zhang Y, Zhang G, Du G, Zhang H, Zou H, Zhang H, Jie Z, Liang S, Jia H, Wan J, Lin D, Li J, Fan G, Yang H, Wang J, Bai Y, Xu X. 2017. Taxonomic structure and functional association of foxtail millet root microbiome. GigaScience 6:1-12. https://doi.org/10.1093/gigascience/gix089.
32. Marcos-Zambrano LJ, Karaduzovic-Hadziabdic K, Loncar Tu-Rukalo T, Przymus P, Trajkovik V, Aasmets O, Berland M, Gruca A, Hasic J, Hron K, Klammsteiner T, Kolev M, Lahti L, Lopes MB, Moreno V, Naskinova I, Org E, Paciência I, Papoutsoglou G, Shigdel R, Stres B, Vilne B, Yousef M, Zdravevski E, Tsamardi-Nos I, Carrillo de Santa Pau E, Claesson MJ, Moreno-Indias I, Truu J. 2021. Applications of machine learning in human microbiome studies: a review on feature selection, biomarker identification, disease prediction and treatment. Front Microbiol 12:634511. https://doi.org/10.3389/fmicb.2021.634511.

33. Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15:550. https://doi.org/10.1186/s13059-014-0550-8.

34. Badri DV, Zolla G, Bakker MG, Manter DK, Vivanco JM. 2013. Potential impact of soil microbiomes on the leaf metabolome and on herbivore feeding behavior. New Phytol 198:264-273. https://doi.org/10.1111/nph.12124.

35. Shi S, Nuccio E, Herman DJ, Rijkers R, Estera K, Li J, da Rocha UN, He Z, Pett-Ridge J, Brodie EL, Zhou J, Firestone M. 2015. Successional trajectories of rhizosphere bacterial communities over consecutive seasons. mBio 6:e00746-15.

36. Ma Z, Yi Z, Bayar K, Fu Y, Liu H. May 2021. Community dynamics in rhizosphere microorganisms at different development stages of wheat growing in confined isolation environments. Appl Microbiol Biotechnol 105: 3843-3857. https://doi.org/10.1007/s00253-021-11283-1.

37. Zhalnina K, Louie KB, Hao Z, Mansoori N, da Rocha UN, Shi S, Cho H, Karaoz U, Loqué D, Bowen BP, Firestone MK, Northen TR, Brodie EL. 2018. Dynamic root exudate chemistry and microbial substrate preferences drive patterns in rhizosphere microbial community assembly. Nat Microbiol 3:470-480. https://doi.org/10.1038/s41564-018-0129-3.

38. Dibner RR, Weaver AM, Brock MT, Custer GF, Morrison HG, Maig-Nien L, Weinig C. 2021. Time outweighs the effect of host devel-Rhizobacteria associated with *Brassica rapa* growth optimal stage on microbial community composition. FEMS Microbiol Ecol 97:efiab102. https://doi.org/10.1093/femsec/fiab102.

39. Goodspeed D, Liu JD, Chehab EW, Sheng Z, Francisco M, Klieben-Stein DJ, Braam J. 2013. Postharvest circadian entrainment enhances crop pest resistance and phytochemical cycling. Curr Biol 23:1235-1241. https://doi.org/10.1016/j.cub.2013.05.034.

40. Badri DV, Loyola-Vargas VM, Broeckling CD, Vivanco JM. May 2010. Root secretion of phytochemicals in arabidopsis is predominantly not influenced by diurnal rhythms. Mol Plant 3:491-498. https://doi.org/10.1093/mp/ssi004.

41. Goss-Souza D, Mendes LW, Rodrigues JLM, Tsai SM. 2020. Ecological processes shaping bulk soil and rhizosphere microbiome assembly in a long-term amazon forest-to-agriculture conversion. Microb Ecol 79:110-122. https://doi.org/10.1007/s00248-019-01401-y.

42. Shen Y, Ryser ET, Li H, Zhang W. 2021. Bacterial community assembly and antibiotic resistance genes in the lettuce-soil system upon antibiotic exposure. Sci Total Environ 778:146255. https://doi.org/10.1016/j.scitotenv.2021.146255.

43. Trivedi P, Leach JE, Tringe SG, Sa T, Singh BK. 2020. Plant-microbiome interactions: from community assembly to plant health. Nat Rev Microbiol 18:607-621. https://doi.org/10.1038/s41579-020-0412-1.

44. Ferrenberg S, O'Neill SP, Knelman JE, Todd B, Duggan S, Bradley D, Robinson T, Schmidt SK, Townsend AR, Williams MW, Cleveland CC, Melbourne BA, Jiang L, Nemer gut DR. 2013. Changes in assembly processes in soil bacterial communities following a wildfire disturbance. ISME J 7:1102-1111. https://doi.org/10.1038/ismej.2013.11.

45. Gu Y, Dong K, Geisen S, Yang W, Yan Y, Gu D, Liu N, Borisjuk N, Luo Y, Friman VP. 2020. The effect of microbial inoculant origin on the rhizosphere bacterial community composition and plant growth promotion. Plant Soil 452:105-117. https://doi.org/10.1007/s11104-020-04545-w.

46. Chang H-X, Bowen CR, Haudenshield JS, Hartman GL. 2017. Metagenome-wide association study and machine learning prediction of bulk soil microbiome and crop productivity. Front Microbiol 8:519.

47. Větrovský T, Baldrian P. 2013. The variability of the 16S rRNA gene in bacterial genomes and its consequences for bacterial community analyses. PLoS One 8:e57923. https://doi.org/10.1371/journal.pone.0057923.

48. Burghardt LT, Epstein B, Guhlin J, Nelson MS, Taylor MR, Young ND, Sadowsky MJ, Tiffin P. 2018. Select and resequence reveals relative fitness of bacteria in symbiotic and free-living environments. Proc Natl Acad Sci U S A 115:2425-2430. https://doi.org/10.1073/pnas.1714246115.

49. Kumar A, Kumar A, Pratush A. 2014. Molecular diversity and functional variability of environmental isolates of Bacillus species. SpringerPlus 3: 312. https://doi.org/10.1186/2193-1801-3-312.

50. Maheshwari D, Agarwal M, Dheeman S. 2014. Trends and prospects of microbial diversity in rhizosphere, p 1-22. In Bacterial diversity in sustainable agriculture. Springer, Cham, Switzerland. https://doi.org/10.1007/978-3-319-05936-5_1.

51. Hol WHG, Garbeva P, Hordijk C, Hundscheid MPJ, Gunnewiek PJKAK, van Agtmaal M, Kuramae EE, de Boer W. 2015. Non-random species loss in bacterial communities reduces antifungal volatile production. Ecology 96:2042-2048. https://doi.org/10.1890/14-2359.1.

52. Gamez RM, Ramirez S, Montes M, Cardinale M. 2020. Complementary dynamics of banana root colonization by the plant growth-promoting rhizobacteria Bacillus amyloliquefaciens Bs006 and Pseudomonas palleroniana Ps006 at spatial and temporal scales. Microb Ecol 80:656-668. https://doi.org/10.1007/s00248-020-01571-0.

53. Kolton M, Erlacher A, Berg G, Cytryn E. 2016. The flavobacterium genus in the plant holobiont: ecological, physiological, and applicative insights, p 189-207. In Microbial models: from environmental to industrial sustainability. Springer, Singapore. https://doi.org/10.1007/978-981-10-2555-6_9.

54. Yin C, Casa Vargas JM, Schlatter DC, Hagerty CH, Hulbert SH, Paulitz TC. 2021. Rhizosphere community selection reveals bacteria associated with reduced root disease. Microbiome 9:86. https://doi.org/10.1186/s40168-020-00997-5.

55. DeBoy RT, Mongodin EF, Fouts DE, Tailford LE, Khouri H, Emerson JB, Mohamoud Y, Watkins K, Henrissat B, Gilbert HJ, Nelson KE. 2008. Insights into plant cell wall degradation from the genome sequence of the soil bacterium Cellvibrio japonicus. J Bacteriol 190:5455-5463. https://doi.org/10.1128/JB.01701-07.

56. Suarez C, Ratering S, Kramer I, Schnell S. 2014. Cellvibrio diazotrophicus sp. nov., a nitrogen-fixing bacteria isolated from the rhizosphere of salt meadow plants and emended description of the genus Cellvibrio. Int J Syst Evol Microbiol 64:481-486. https://doi.org/10.1099/ijs.0.054817-0.

57. Lin H, Liu C, Li B, Dong Y. 2021. Trifolium repens L. regulated phytoremediation of heavy metal contaminated soil by promoting soil enzyme activities and beneficial rhizosphere associated microorganisms. J Hazard Mater 402:123829. https://doi.org/10.1016/j.jhazmat.2020.123829.

58. Morais MC, Mucha Â, Ferreira H, Gonçalves B, Bacelar E, Marques G. 2019. Comparative study of plant growth-promoting bacteria on the physiology, growth and fruit quality of strawberry. J Sci Food Agric 99:5341-5349. https://doi.org/10.1002/jsfa.9773.

59. Zhang Y, Wang E, Xu J, Wang N. 2020. Mechanisms underlying the rhizosphere-to-rhizoplane enrichment of Cellvibrio unveiled by genome-centric metagenomics and metatranscriptomics. Microorganisms 8:583. https://doi.org/10.3390/microorganisms8040583.

60. Team TID. 2021. ImageMagick. https://imagemagick.org.

61. Easlon HM, Bloom AJ. 2014. Easy Leaf Area: automated digital image analysis for rapid and accurate measurement of leaf area. Appl Plant Sci 2: 1400033. https://doi.org/10.3732/apps.1400033.

62. Huse SM, Young VB, Morrison HG, Antonopoulos DA, Kwon J, Dalal S, Arrieta R, Hubert NA, Shen L, Vineis JH, Koval JC, Sogin ML, Chang EB, Raffals LE. 2014. Comparison of brush and biopsy sampling methods of the ileal pouch for assessment of mucosa-associated microbiota of human subjects. Microbiome 2:5. https://doi.org/10.1186/2049-2618-2-5.

63. Newton RJ, McLellan SL, Dila DK, Vineis JH, Morrison HG, Eren AM, Sogin ML. 2015. Sewage reflects the microbiomes of human populations. mBio 6:e02574-14.

64. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. 2016. DADA2: high resolution sample inference from Illumina amplicon data. Nat Methods 13:581-583. https://doi.org/10.1038/nmeth.3869.

65. Callahan BJ, Sankaran K, Fukuyama JA, McMurdie PJ, Holmes SP. 2016. Bioconductor workflow for microbiome data analysis: from raw reads to community analyses. F1000Res 5:1492. https://doi.org/10.12688/f1000research.8986.1.

66. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glöckner FO. 2013. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41: D590-D596.

67. McMurdie PJ, Holmes S. 2013. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8: e61217. https://doi.org/10.1371/journal.pone.0061217.

68. Davis NM, Proctor DM, Holmes SP, Relman DA, Callahan BJ. 2018. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome 6:226. https://doi.org/10.1186/s40168-018-0605-2.
69. Paulson JN, Stine OC, Bravo HC, Pop M. 2013. Differential abundance analysis for microbial marker-gene surveys. Nat Methods 10:1200-1202. https://doi.org/10.1038/nmeth.2658.

70. Oksanen J, Kindt R, Legendre P, O'Hara B, Simpson GL, Solymos P, Henry M, Stevens H, Wagner H. 2008. The vegan package. https://www.researchgate.net/profile/Gavin_Simpson/publication/228339454_The_vegan_Package/links/0912f50be86bc29a7f000000/The-vegan-Package.pdf. Accessed 19 May 2020.

71. Schleip K. 2011. Phangorn: phylogenetic analysis in R. Bioinformatics 27: 592-593. https://doi.org/10.1093/bioinformatics/btq706.

72. Kurtz ZD, Müller CL, Miraldi ER, Littman DR, Blaser MJ, Bonneau RA. May 2015. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput Biol 11:e1004226. https://doi.org/10.1371/journal.pcbi.1004226.

73. Lemonnier C, Perennou M, Eveillard D, Fernandez-Guerra A, Ley-Naert A, Marié L, Morrison HG, Memory L, Paillard C, Maignien L. 2020. Linking spatial and temporal dynamic of bacterioplankton communities with ecological strategies across a coastal frontal area. Front Mar Sci 7:e376. https://doi.org/10.3389/fmars.2020.00376.

74. Blondel VD, Guillaume JL, Lambiotte R, Lefebvre E. 2008. Fast unfolding of communities in large networks. J Stat Mech 2008:P10008. https://doi.org/10.1088/1742-5468/2008/10/P10008.

75. Csardi G, Nepusz T. 2005. The igraph software package for complex network research. InterJournal 2005:1695.

76. Friedman JH, Hastie T, Tibshirani R. 2010. Regularization paths for generalized linear models via coordinate descent. J Stat Softw 33:1-22.

77. Sisk-Hackworth L, Kelley ST. 2020. An application of compositional data analysis to multiomic time-series data. NAR Genom Bioinform 2:elqaa079. https://doi.org/10.1093/nargab/lqaa079.

78. Bürkner PC. 2017. brms: an R package for Bayesian multilevel models using Stan. J Stat Softw 80:1-28. https://doi.org/10.18637/jss.v080.i01.

79. Vehtari A, Gelman A, Gabry J. 2017. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput 27:1413-1432. https://doi.org/10.1007/s11222-016-9696-4.