Does soil history decline in influencing the structure of bacterial communities of Brassica napus host plants across different growth stages?

Abstract Soil history has been shown to condition future rhizosphere microbial communities. However, previous experiments have also illustrated that mature, adult plants can “re-write,” or mask, different soil histories through host plant–soil community feedbacks. This leaves a knowledge gap concerning how soil history influences bacterial community structure across different growth stages. Thus, here we tested the hypothesis that previously established soil histories will decrease in influencing the structure of Brassica napus bacterial communities over the growing season. We used an on-going agricultural field experiment to establish three different soil histories, plots of monocrop canola (B. napus), or rotations of wheat-canola, or pea-barley-canola. During the following season, we repeatedly sampled the surrounding bulk soil, rhizosphere, and roots of the B. napus hosts at different growth stages—the initial seeding conditions, seedling, rosette, bolting, and flower—from all three soil history plots. We compared composition and diversity of the B. napus soil bacterial communities, as estimated using 16S rRNA gene metabarcoding, to identify any changes associated with soil history and growth stages. We found that soil history remained significant across each growth stage in structuring the bacterial bulk soil and rhizosphere communities, but not the bacterial root communities. This suggests that the host plant’s capacity to “re-write” different soil histories may be quite limited as key components that constitute the soil history’s identity remain present, such that the previously established soil history continues to impact the bacterial rhizosphere communities, but not the root communities. For agriculture, this highlights how previously established soil histories persist and may have important long-term consequences on future plant–microbe communities, including bacteria.


communities of Brassica napus host plants across different growth stages?
Authors: Andrew J.C. Blakney1*, Marc St-Arnaud1, Mohamed Hijri1, 2* List of Supplementary Tables Table S1.Environmental Data Table S2.Mock Community Composition Table S3.Primers Table S4.PERMANOVA List of Supplementary Figures Figure S1.Flower Supplementary Methods DNA extraction from Test Phase Brassicaceae root and rhizosphere samples For bulk soil and rhizosphere samples, ~500 mg was used for the NucleoSpin Soil gDNA Extraction Kit (Macherey-Nagel, Germany), while ~130 mg of seeds and roots were used with the DNeasy Plant DNA Extraction Kit (Qiagen, Germany) (Lay et al., 2018;Blakney et al., 2022).A no-template extraction negative control was used with the bulk soil, rhizosphere and root extractions, and included with the Test Phase samples (Fig. S2), to assess the influence of the extraction kits on our sequencing results, and the efficacy of our lab preparation.All extracted DNA samples were quantified using the Qubit dsDNA High Sensitivity Kit (Invitrogen, USA), and qualitatively evaluated by mixing ~2 µL of each sample with 1 µL of loading dye containing Gel Red (Biotium) and running it on a 0.7 % agarose gel for 30 minutes at 150 V.The no-template extraction negative controls were confirmed to not contain DNA after extraction, where the detection limit was > 0.1 ng (Qubit, Invitrogen, USA).16S rRNA gene amplicon generation and sequencing to estimate the bacterial community First, all DNA samples were diluted 1:10 into 96-well plates using the Freedom EVO100 robot (Tecan, Switzerland).To assess potential bias caused by lab manipulations, sequencing and downstream bioinformatic processing, a commercially available 16S rRNA mock community, of known composition (Table S2), was included on each plate (Fig. S2) following the manufacturer's instructions (BEI Resources, USA).The mock community contained DNA of 20 bacterial species (Table S2) in equimolar counts (10 6 copies/µL) of 16S rRNA genes.A no-template PCR negative control was also included on each plate, to assess the influence of the PCR reaction, and the efficacy of lab preparation on sequencing (Fig. S2).Four µL of each negative control was mixed with 1 µL of loading dye containing Gel Red (Biotium) and visualized on a 1% agarose gel after 60 minutes at 100 V. None of the negative controls for the DNA extractions nor the PCR reactions contained detectable amounts of DNA prior to submission.
The prepared plates of the Test Phase DNA samples were submitted to Génome Québec (Montréal, Québec) for 16S rRNA amplicon generation and sequencing (Bell et al., 2016;Lay et al., 2018;Blakney et al., 2022).PCR amplification used the S-D-Bact-0341-b-S-17 forward and S-D-Bact-0785-a-A-21 reverse primers, commonly referred to as 341F and 805R, respectively, to generate a 416 bp fragment from the V3-V4 region of the 16S rRNA gene (Klindworth et al., 2012).These amplicons were then prepared for paired-end 250 bp sequencing using Illumina's MiSeq platform (Génome Québec, Montréal) (Bell et al., 2016;Lay et al., 2018;Blakney et al., 2022).We estimated this should provide a mean of 60 000 reads per sample, which is in line with previous studies that describe bacterial communities (Bell et al., 2016;Lay et al., 2018;Blakney et al., 2022).Estimating ASVs from MiSeq 16S rRNA gene amplicons The 16S rRNA gene amplicons generated by Illumina MiSeq were used to estimate the diversity and composition of the bacterial communities present in the bulk soil, seed, rhizosphere and roots from across the Test Phase B. napus growth stages.The integrity and totality of the 16S MiSeq data downloaded from Génome Québec was confirmed using their MD5 checksum protocol (Roy et al., 2018).Subsequently, all data was managed, and analyzed in R (4.0.3 R Core Team, 2020), and plotted using ggplot2 (Wickham, 2016).
The dada2 package (Callahan et al., 2016a) was first used to filter and trim all 11 010 728 raw reads, forward and reverse, from the 16S rRNA gene amplicon data generated from the control samples, the mock communities, and the Test Phase B. napus samples, using the filterAndTrim function (Fig. S2), as described in Blakney et al., 2022.Filtered and trimmed reads were then processed through DADA2 for ASV inference (Fig. S3).Default settings were used throughout the DADA2 pipeline, except the DADA inference functions dadaF and dadaR which used the pool ='pseudo' argument, to increase the likelihood of identifying rare taxa.Consequently, the chimera removal function removeBimeraDenovo included the method ='pooled' argument (Callahan et al., 2016b).
ASVs identified from the 16S rRNA gene amplicon data were assigned taxonomy-to the species level where possible-using the Silva database release 138, which adopts the standardized Genome Taxonomy Database taxonomy structure (Yilmaz et al., 2013;Parks et al., 2018).The quality of the data was assessed using the included controls (Fig. S4); any ASVs identified as chloroplasts, or mitochondria were subsequently removed from the data, as were off-target archaeal ASVs.Rarefaction curves confirmed that we captured the majority of the bacterial communities in both the bulk soil, seed, rhizosphere, and roots (Fig. S5).Test Phase Brassicaceae 16S rRNA sequencing data was subsequently re-analysed independently following the described protocol to avoid any biases from the no-template negative controls and the mock communities.These are the Test Phase B. napus ASVs which are reported hereafter.

Estimating absolute abundance of bacterial communities by quantitative PCR
To identify any changes in abundance of the bacterial ASVs within the Test Phase B. napus growth stages, we estimated the absolute abundance, or size, of the bacterial communities in each Test Phase DNA sample by qPCR (Azarbad et al., 2018;Blakney et al., 2022).Given the technicallimitations of high-throughput sequencing in assessing abundance, estimating the absolute abundance can provide data to better interpret the bacterial communities (Gloor et al., 2017;Props et al., 2017;Harrison et al, 2020;Jian et al., 2020).As such, we quantified the number of 16S rRNA gene sequences present in each DNA sample via qPCR (as Ct values) and then estimated the corresponding community size as the 16S rRNA gene copy number/µg from a standard curve (Zhang et al., 2017;Azarbad et al., 2018, Fig. S6).
First, a standard curve of 16S rRNA gene copy numbers was constructed.Near full-length 1.5 kb 16S rRNA gene fragments were PCR amplified using the primers PA-27F-YM and PH-R (Bruce et al., 1992;Table S3) from DNA extracted from previously used soil samples (Lay et al., 2018).The 16S PCR reactions consisted of 11.5 µL dH2O, 5.0 µL of 10X Buffer (Qiagen, Canada), 2.5 µL of 10 µM PA-27F-YM forward and PH-R reverse primers (Alpha DNA, Montréal, Canada; Klindworth et al., 2012), 1.0 µL of dNTPs (Qiagen, Canada), 0.5 µL of T. aq polymerase (Qiagen, Canada), and 2 µL of template DNA, for a total volume of 25 µL.PCR amplification was run in an Eppendorf Mastercycle ProS (Germany) thermocycler, and consisted of an initial denaturation of 2 minutes at 95°C, followed by 30 cycles of 30 seconds denaturation at 95°C, 30 seconds annealing at 55°C, and 1 minute elongation at 72°C, before a final elongation of 5 minutes at 72°C (Bell et al., 2016;Lay et al., 2018;Blakney et al., 2022).The amplified 1.5 kb 16S rRNA gene was visualized by a 1% gel electrophoresis, as described above, quantified using the QuBit dsDNA High Sensitivity Kit (Invitrogen, USA), and then serially diluted to 10 -9 .One µL of each dilution was then used as template in a 10 µL qPCR reaction.
The 16S rRNA gene qPCR reactions consisted of 5.0 µL of Maxima SYBR Green/ROX qPCR Mix (ThermoFisher Scientific, Canada), 3.4 µL dH2O, 0.3 µL of 10 µM Eub 338 forward and Eub 518 reverse primers (Alpha DNA, Montréal, Canada; Fierer et al., 2005).All qPCR reactions were set-up in triplicate in 96-well plates using the Freedom EVO100 robot (Tecan, Switzerland), with a no-template negative control included on each plate.Reactions were run in a ViiA 7 Real-Time PCR System (Life Technologies, Canada) following the same cycling conditions as described previously for the 16S rRNA PCR amplification.The Eub338/Eub518 qPCR reaction amplified a 200 bp region of the V3 region (Muyzer et al., 1993;Nogales et al., 1999;Bathe & Hausner, 2006;Davis et al., 2009).The number of 16S rRNA gene copies present in the serially diluted standard were calculated using the formula (Godornes et al., 2007): Number of 16S rRNA gene copies µL -1 = Avogadro's Constant x DNA (g µL -1 ) Number of base pairs x 600 Daltons The standard curve for each diluted sample was plotted, with an R 2 value of 0.9938 and an amplification efficiency of -3.2013 (Fig. S6), falling within acceptable values (Fierer et al., 2005).
16S rRNA gene copy numbers were then estimated for each Test Phase sample by using 1 µL of a 1:10 dilution of DNA as template in the same 16S rRNA qPCR reaction and cycling conditions as described above for the standard curve.Melt curves generated by 0.5°C increments at the end of the qPCR programme confirmed amplicon specificity, and the 16S rRNA gene copy number was then determined from the standard curve.A correction to determine the absolute abundance of each ASVs inferred from the Test Phase samples was achieved by multiplying the total 16S rRNA gene copy number per µg, as estimated from the qPCR reaction, by the relative abundance matrix of ASVs identified (Azarbad et al., 2018;Bakker, 2018).

Generating phylogenetic trees
In order to employ phylogeny-based analysis methods, including phylogenetic diversity (PD), and UniFrac distances, we assembled phylogenies for the Test Phase B. napus bacterial communities.Following the method described by Callahan et al., 2016, 16S rRNA gene sequences for each ASV inferred from the Test Phase data were aligned using a profile-to-profile algorithm (Wang & Dunbrack, 2004) with a dendrogram guide tree using the decipher package (Wright, 2016).With the phangorn package (Schliep, 2011), the maximum likelihood of each site was calculated using the dist.mlfunction using a JC69 equal base frequency model, before assembling phylogenies using the neighbour-joining method.An optimized GTR nucleotide substitution model was fitted to the phylogeny using the optim.pmlfunction.Phylogenies were subsequently added to the phyloseq object.a-diversity of the Brassica napus bacterial communities Faith's PD was calculated as an a-diversity index from the B. napus samples using the pd function from the picante package (sum of all branch lengths separating taxa in a community).Log transformed PD indices were confirmed to respect normality.Normality of the residuals was established with a Shapiro-Wilk test, shapiro.test,while the heteroscedascity of residuals was confirmed with using a Bartlett test, bartlett.test.Within each compartment, we assessed differences of the mean PD between growth stages for each soil history, and their interactions, using a Multi-Factor ANOVA and Tukey's Post-Hoc test for significant groups that respected the assumptions of normality.For significant ANOVAs, a post-hoc Tukey's Honest Significant Difference test, TukeyHSD, was used to determine which groups were statistically different.
Differential taxa clusters identified significantly enriched (i.e., abundant), using the nonparametric Kruskal test, followed by the post hoc pairwise Wilcox test, with an FDR correction.
Significantly enriched taxa, labelled in bold, were tested between each pair of growth stages and soil history.A significant indicator value is obtained if an ASV has a large mean abundance within a group, compared to another group (specificity), and has a presence in most samples of that group (fidelity).The fidelity component complements the differential abundance approach between taxa clusters, which only considers abundance.b-diversity of the Brassica napus bacterial communities PERMANOVA divides any variation in the ordinated data distance matrix among all the pairs of specified experimental factors.The UniFrac distance index incorporates the phylogenetic relationship of each ASV and their absolute abundance, as estimated by qPCR of the 16S rRNA gene.Determining community distances based only on the number of shared taxa does not account for evolutionary distances between taxa, which are often extremely diverse among microbes.Conversely, using a UniFrac index illustrates how bacterial community composition varies by phylogeny, which provides insight into how different community assembly mechanisms, including dispersal, drift, selection, and speciation, may be at work.
To further characterize what ecological mechanisms may be responsible for changes to the b-diversity, we partitioned it into turnover (i.e.species replacement) and nestedness (i.e.loss/gains that result in poor species richness being a subset of richer sites) components.To take advantage of our abundance data, where turnover and nestedness are analogous to balanced variation & abundance gradients, we first arranged the abundance data between all samples using the betapart.core.abundfunction.Then, to compare each sample to itself through time, we used the beta.pair.abundfunction, with Bray-Curtis dissimilarity distances.As nestedness remained below 0.1 throughout, it was not further analyzed (data not shown).As turnover among compartments was not normally distributed, the non-parametric Kruskal test was used to determine significance, followed by the post-hoc pairwise Wilcox test with an FDR correction.Within each compartment, we identified any significant differences of the mean turnover between growth stages for each soil history, and their interactions with a Multi-Factor ANOVA, as described above for a-diversity, as normality was respected.
Table S2.Bacterial strains included in the mock community (BEI Resources, USA) of known composition, was included on each plate (Fig. S2).The mock community contains DNA of 20 bacterial taxa in equimolar counts (10 6 copies/µL) of 16S rRNA genes.  .Conceptual design of the experiment.(A) Soil history was established during the growing season of 2018 (Conditioning Phase, t = 0), while the effect of soil history on Brassica napus bacterial rhizosphere and root communities, or associated bulk soil, was observed the following season, in 2019 (Test Phase, t = 1, 2, 3 … 5).Samples were harvested throughout the growing season at five growth stages; seed, seedling, rosette, bolting, and flower.(B) Field plan for the experiment.The experimental design was a split-plot replicated in four complete blocks.In the 'Conditioning Phase', three soil history treatments were randomly assigned: treatment 1 (green) another generation of canola (B.napus L., cv.L252LL) as a monocrop, treatment 4 (brown) spring wheat (Triticum aestivum cv.AAC Brandon), as the wheat phase of a two-year crop rotation with B. napus, and treatment 8 (black) barley (Hordeum vulgare cv.Canmore), as a three-year crop rotation with B. napus and pea (Pisum sativum L. cv AAC Lacombe).Plot numbers appear in subscript.

Fig. S2
Figure S2.Organization of our lab workflow for the Test Phase Brassicaceae samples from harvest to generating amplicon sequence variants (ASVs).The Test Phase Brassicaceae samples were harvested in from 12 plots organized as split plots that were randomly assigned three soil history treatments (Fig S1B).Each plot was harvest at five growth stages (seed, seedling, rosette, bolting and flower); note the seed was all the same for each plot (n = 1).Bulk soil samples were taken from each plot at each growth stage (n = 60).In the field, each plant had its root systems and associated immediately flash-frozen in liquid nitrogen and kept on ice.In the lab, samples were sieved to remove debris (rocks, undecomposed straw, twigs etc…), and root systems were divided into rhizosphere and root samples (n = 48 each).Roots and seeds were then ground in liquid nitrogen, and DNA was extracted from all the Test Phase samples.No-template extraction controls were included to assess what contaminates, or biases, the extraction kits (Machery-Nagel Nucleospin Soil gDNA kit, at left in brown, and Qiagen Plant DNEasy kit in green) might impart.We also included no-template PCR negative controls, and confirmed by gel electrophoresis that none of the no-template extraction controls, nor the no-template PCR controls, contained DNA.To help identify sequencing biases, or bath effects, a replicate of the bacterial mock community (BEI Resources, USA) was included on each plate submitted for sequencing.All DNA samples were submitted to Génome Québec for 16S rRNA PCR amplification, library preparation, and paired-end 250 bp Illumina MiSeq sequencing.All reads were subsequently trimmed and processed through the DADA2 pipeline for ASV inference.Figure S7.Phylogenetic diversity was significantly different between compartments (p < 0.001) in bacterial communities throughout the Brassica napus growing season.Root and seed communities were significantly different from each compartment across growth stages.Bulk soil and rhizosphere communities were strikingly similar across growth stages, except for the flower and bolting stage, respectively, which were both more diverse than the rhizosphere communities at the seedling and rosette stages.Diversity across growth stages was tested with a Multi-Factor ANOVA, which confirmed that the compartments and the Test Phase B. napus growth stages were significant and did interact.Statistically significant groups were identified using Tukey's post-hoc test.Principal co-ordinate analysis illustrated that seed and root bacterial communities were phylogenetically distinct from bulk soil and rhizosphere communities.Bacterial root communities were more similar at the seedling and rosette stages, and appeared more diverse over time.Axis 1 and 2 captured 69.4% of the variability among the bacterial communities.(B) Bacterial communities from the bulk soil and rhizosphere were more phylogenetically similar at the seed and seedling stages, and increasingly diverse at following growth stages.Axis 1 and 2 captured 35% of the variability among the bacterial bulk soil and rhizosphere communities.Principal co-ordinate analysis were plotted using UniFrac distances weighted by absolute abundance, where phylogenetically similar communities were plotted closer together.

Fig. S9
Figure S9.Specific bacterial taxa were significantly different between adjacent compartments at specific Brassica napus growth stages.(A) Bacterial rhizosphere communities at the seedling stage were enriched in taxa compared to the bulk soil and root communities.(B) At the rosette stage, bulk soil and rhizosphere communities were relatively similar, while both the rhizosphere and root communities were enriched in different taxa.(C) The bacterial communities in the bulk soil and rhizosphere were also similar at the bolting stage, though the the rhizosphere and root communities both exhibited specific enrichments in different taxa.(D) Bacterial rhizosphere communities were enriched in taxa compared to both the bulk soil and root communities.The abundance of each taxonomic group was compared between compartments, using the using the non-parametric Kruskal test and the post-hoc pairwise Wilcox test, with the FDR correction.Taxa that were significantly (p.adj < 0.05) more abundant were highlighted.

Fig. S8 Figure S8 .
Fig. S8 Taxa have been provided to illustrate the level of comparison.

Table S3 .
Primers used in this study.

Table S4 .
PERMANOVA for all the sampled Test Phase communities identified compartment (bulk soil, rhizosphere, or root), growth stage, and soil history established in the Conditioning Phase, as significant experimental factors.PERMANOVA was calculated using a Bray-Curtis and Weighted Unifrac distance matrix, with 9999 permutations.