Applications of the indole-alkaloid gramine modulate the assembly of individual members of the barley rhizosphere microbiota

Soil substrate

The experiments were carried out in “Quarryfield” soil, an unfertilized Scottish agricultural soil collected at the site of The James Hutton Institute, Invergowrie, Scotland (UK) (56°27′5″N 3°4′29″W). This soil was used previously to grow barley and left unplanted at least 4 years before being used for the experiments. Physical and chemical characterization: Silt: 39.35%; Clay: 11.08%; Sand: 49.58%; pH in water: 5.8; C.E.C.: 15.05 [meq/100g]; Org. Matter DUMAS: 5.95%.

Gramine adsorption

Five grams of Quarryfield soil were mixed with 10 mL of 10 mmol L⁻¹ CaCl₂ solution containing gramine at the following concentrations: 0, 1.5, 5, 10, 15, 20, 30, 60, 100, 150, 200, 250, 300, 400, 600, 800, 1,000 mg L⁻¹. Five replicates were prepared for each concentration. Soil suspensions were shaken for 24 h at room temperature. Then, suspensions were centrifuged for 5 min at 5,000 g, the supernatant was collected, filtered (0.45 µm, Phenomenex) and analysed by liquid chromatography (HPLC) using the method reported in (Maver et al., 2020). Pure reagent grade gramine was used to prepare two stock solutions in water of 0.5 mol L⁻¹ and 10 mmol L⁻¹. Adsorption solutions have been prepared by dilution from the stock solution with milliQ water. The adsorption rate of gramine in soil was obtained by difference between the initial and final concentration measured in the supernatant by HPLC.

Adsorption isotherms of gramine were fitted applying several nonlinear models: the two-parameter Langmuir and Freundlich isotherms and the three-parameter Sigmoidal Langmuir, Redlich-Peterson and Sips isotherms (Limousin et al., 2007; Foo & Hameed, 2010).

Plant material

Seeds of elite barley genotypes, two-row malting Barke and six-row malting Morex were selected for this experiment. Both are well represented in barley studies: Barke as a parental donor in the development of a nested associated mapping (NAM) population (Maurer et al., 2015), and Morex for being the first sequenced barley genotype (The International Barley Genome Sequencing Consortium, 2012).

Growth conditions

Barley seeds were cleaned using deionized water and gently shaken for 1 min. After that water was discarded and the whole process repeated 3 times. For seed germination, seeds were placed on petri dishes containing a semi-solid 0.5% agar solution. After a week, seedlings displaying a comparable development were individually transferred into 12-cm pots containing approximately 500 g of Quarryfield agricultural soil (Robertson-Albertyn et al., 2017; Alegria Terrazas et al., 2020) previously sieved to remove stones and large debris. Plants were grown in a randomized design in a glasshouse at 18/14 °C (day/night) temperature regime with 16 h daylight that was supplemented with artificial lighting to maintain a minimum light intensity of 200 µmol quanta m⁻² s⁻¹. The stock solution of gramine was prepared by adding pure reagent grade gramine (Sigma-Aldrich, >99%), in water, sonicated for 20 min and then stored at 4 °C. After 4 days of growing in soil, two different final gramine concentrations, 24 µmol L⁻¹ and 46 µmol L⁻¹, were added directly on the soil of selected pots. Mock controls (i.e., gramine 0 µmol L⁻¹) were supplemented with sterilised water. Additional watering was performed every 2 days with the application of 50 mL of deionized water to each pot. For each gramine treatment (i.e., 0, 24 and 46 µmol L⁻¹) we used five replicates (i.e., five individual pots) per barley genotype and unplanted pots containing the same soil substrates used as ‘bulk’ soil controls (i.e., 45 pots in total). Individual replicated pots were maintained in the glasshouse for 4 weeks post-transplantation, when the tested genotypes reached early stem elongation, corresponding to Zadoks stages 30–35 (Tottman, Makepeace & Broad, 1979).

Rhizosphere fractionation and sampling of soil-grown barley plants and bulk soil

The preparation of material for amplicon sequencing was performed following established protocols (Robertson-Albertyn et al., 2017; Alegria Terrazas et al., 2020). Briefly, four-week-old barley plants were carefully removed from the soil, and the shoot and root separated. The shoot was dried at 70 °C for 48 h and the dry weight collected. The roots were gently shaken to remove loosely bound soil particles, and the resulting root system and tightly adhered soil, operationally defined as rhizosphere, was further sectioned to retain the uppermost six cm of the seminal root system of each sample. This root material was transferred in a sterile 50 mL falcon tube containing 15 mL of phosphate buffered saline solution (PBS). Samples were then vortexed for 30 s, the soil sedimented for 2–3 mins, and the roots transferred in a new 50 mL falcon tube with 15 mL PBS, in which the samples were vortexed again for 30 s to separate the remaining rhizosphere soil from roots. The roots were discarded, the two falcon tubes were combined in one single falcon tube, now containing the rhizosphere soil fraction, and then centrifuged at 1,500 g for 20 min. After centrifugation, the supernatant was discarded, and the pellet immediately stored at −80 °C. In the unplanted soil controls (i.e., the bulk soil pots), a portion of soil corresponding to the area explored by roots was collected with a spatula and processed as described for planted soils. Until DNA extraction, all the samples were stored at −80 °C.

Metagenomic DNA extraction from rhizosphere and bulk soil specimens

Total DNA was extracted from the rhizosphere and unplanted soil samples using FastDNA™ SPIN kit for soil (MP Biomedicals, Solon, USA) following the instructions by the manufacturer. To assess the concentration and the quality (260/280 nm and 260/230 nm ratios) of the extracted DNA, a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA) was used and samples were stored at −20 °C until further analysis. Aliquots at a final DNA concentration of 10 ng mL⁻¹ were prepared for each sample using sterilized deionized water, and stored at −20 °C.

16S rRNA gene amplicon library construction

The amplicon library was generated via a selective PCR amplification of the hypervariable V4 region of the 16S rRNA gene using the PCR primers 515F (5′-GTGCCAGCMGCCG CGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) as previously described (Robertson-Albertyn et al., 2017; Alegria Terrazas et al., 2020). Briefly, PCR primer sequences were fused with Illumina flow cell adapter sequences at their 5′ termini and the 806R primers contained 12-mer unique ‘barcode’ sequences to enable the multiplexed sequencing of several samples in a single pool (Caporaso et al., 2012).

For each individual bulk and rhizosphere preparations, 50 ng of DNA was subjected to PCR amplification using the Kapa HiFi HotStart PCR kit (Kapa Biosystems, Wilmington, USA). The individual PCR reactions were performed in 20 µL final volume and contained:

• 4 µL of 5X Kapa HiFi Buffer

• 10 µg Bovine Serum Albumin (BSA) (Roche, Mannheim, Germany)

• 0.6 µL of a 10 mM Kapa dNTPs solution

• 0.6 µL of 10 µM solutions of the 515F and 806R PCR primers

• 0.25 µL of Kapa HiFi polymerase

Reactions were performed in a G-Storm GS1 thermal cycler (Gene Technologies, Somerton, UK) using the following programme: 94 °C (3 min), followed by 35 cycles of 98 °C (30 s), 50 °C (30 s) 72 °C (1 min) and a final step of 72 °C (10 min). For each 515F-806R primer combination, a no template control (NTC) was subjected to the same process. To minimize potential biases originating during PCR amplifications, individual reactions were performed in triplicate and 2 independent sets of triplicate reactions per barcode were performed.

To check the amplification and/or any possible contamination, prior to purification, 6 µL aliquots of individual replicates and the corresponding NTCs were inspected on 1.5% agarose gel. Only samples that display the expected amplicon size and no detectable contamination in NTCs on gel were retained for library preparation.

Individual PCR amplicons replicates were then pooled in a single plate, moving each sample to a specific position according to their barcode. They were purified using Agencourt AMPure XP Kit (Beckman Coulter, Brea, USA) with 0.7 µL AmPure XP beads per 1 µL of sample. Following purification, 6 µL of each sample was quantified using PicoGreen (Thermo Fisher Scientific, Watham, USA). Once quantified, individual barcode samples were pooled to a new tube in an equimolar ratio to generate amplicons libraries.

Illumina 16S rRNA gene amplicon sequencing

Amplicon libraries were supplemented with 15% of a 4 pM phiX solution and run at 10 pM final concentration on an Illumina MiSeq system with paired end 2 × 150 bp reads (Caporaso et al., 2012) as recommended, to generate FASTQ sequence files for processing and analysis.

Amplicon sequencing reads processing

Sequencing reads were processed using a customized bioinformatics pipeline as described before (Terrazas et al., 2019). Briefly, sequencing reads were subjected to quality assessment using FastQC (Andrews, 2010). Amplicon Sequencing Variants (ASVs) were then generated using DADA2 version 1.10 (Callahan et al., 2016) and R 3.5.2 (Team R Development Core, 2018) following the basic methodology outlined in the ‘DADA2 Pipeline Tutorial (1.10)’. Read filtering was carried out using the DADA2 filterAndTrim method, trimming 10bp of sequence from the 5′ of each reads using a truncQ parameter of 2 and maxEE of 2. The remainder of the reads were of high quality so no 3′ trimming was deemed necessary. The dada method was run to determine the error model with a MAX_CONSIST parameter of 20, following which the error model converged after 9 and 12 rounds for the forward and reverse reads respectively. The DADA2 method was then run with the resulting error model to denoise the reads using sample pooling, followed by read merging, using the default minOverlap parameter of 12 bases, then chimera removal using the consensus method. Taxonomy assignment was carried out using the RDP Naive Bayesian Classifier through the ‘assign Taxonomy’ method, with the SILVA database (Quast et al., 2012) version 138, using a minimum bootstrap confidence of 50. The DADA2 outputs were converted to a Phyloseq object (version 1.26.1) (McMurdie & Holmes, 2013). ASVs assigned to chloroplast and mitochondria, putatively representing host contamination, as well as ASVs previously identified as putative lab contaminants (Pietrangelo et al., 2018) were removed in silico from the original Phyloseq object. Likewise, ASVs lacking taxonomic classification at Phylum level (i.e., classified as “NAs”) were removed from the dataset.

As a secondary quality filtering approach (Bokulich et al., 2013), we applied an abundance filtering for any given ASVs to be retained in the final dataset of 20 reads in at least 11% of the samples, representing the number of replicates for a given sample type/condition. In this way, an ASV with overall low abundance due to its association with only one given condition would have been retained and analysed in the dataset. Upon completion of this additional filtering step, we retained 6,615,714 reads (min = 25,615; max = 305,408, mean = 147,015.9) representing over 93.5% of the input quality filtered, non-contaminant sequences. Upon completion of this additional filtering step, individual ASVs were agglomerated at Genus level in Phyloseq. Finally, to control for sample-to-sample differences exceeding a factor of ∼10X in sequencing depth (Weiss et al., 2017), we downsized the Phyloseq object at 25,000 reads per sample.

Statistical analysis

Data analysis was performed in R software v 3.5.2. The following R packages were used: Phyloseq v.1.26.1 (McMurdie & Holmes, 2013) for Alpha and Beta-diversity indexes; DESeq2 v1.22.2 (Love, Huber & Anders, 2014) for the differential analysis of microbial enrichment; ggplot2 v.3.3.2 (Wickham, 2016) for data visualization; Vegan v.2.5-6 (Oksanen et al., 2019) for statistical analysis of beta-diversity; PMCMR v.4.3 (Pohlert, 2018) for non-parametric analysis of variance.

The analysis of the microbiota data was performed on the filtered “Phyloseq object” described above and linking into the analysis the mapping file (metadata information). For the Alpha-diversity analysis, Chao1, Observed ASVs and Shannon indices were calculated using the function to estimate richness included in the Phyloseq package.

For the Beta-diversity analysis, the rarefied ASV table was used as input to compute a Bray–Curtis, dissimilarity matrices. This dissimilarity matrix was visualized using Principal Coordinates Analysis (PCoA) and Canonical Analysis of Principal coordinates (CAP) (Anderson & Willis, 2003). Beta-diversity dissimilarity matrices were tested by Permutational Multivariate Analysis of Variance (Permanova) using Adonis function in Vegan package over 5,000 permutations, to calculate the statistical significance.

Microbial enrichment analysis

The analysis of the microbial enrichment was performed using the DESeq2 package (Love, Huber & Anders, 2014), in order to identify the number of genera significantly enriched in pair-wise comparisons with an adjusted p value (False Discovery Rate, FDR p < 0.05). The microbial enrichment analysis was carried out between bulk soil and planted soil to evaluate the impact of gramine on the ‘rhizosphere effect’. The number of genera enriched in the rhizosphere samples subjected to the three different concentrations of gramine, in both Morex and Barke, was plotted using the package UpSetR (Conway, Lex & Gehlenborg, 2017). In parallel, we performed a series of pair-wise comparisons between Morex and Barke using the genera enriched in these latter specimens from unplanted soil controls. Differentially enriched genera were visualized using ternary plots as previously described (Bulgarelli et al., 2012).

Gramine application impacts on the alpha- and beta-diversity properties of the microbiota thriving at the barley root-soil interface

Gramine availability and mobility in Quarryfield soil are two prerequisites to fulfill its allelopathic ability towards target plants and inter-organismal relationships. Thus, we assessed its adsorption applying several two- and three-parameter non-linearized isotherms models, with the best fit based on the residual sum of squares (RSSs; Fig. S1). This choice was dictated by the fact that non-linear forms permit greater accuracy in predicting parameters compared to linear forms (Foo & Hameed, 2010). The model with the best fit, the Sigmoidal Langmuir isotherm model, resulted in a sigmoidal curve, describing a cooperative adsorption phenomenon, which is common for non-polar compounds (Sparks, 2003; Limousin et al., 2007). A K_foc, defined as organic-carbon normalized Freundlich distribution coefficient, estimated at 1390 (OECD, 2001) classifies gramine as slightly mobile in this particular soil according to FAO Mobility Classification (FAO, 2000), yet remaining available in the soil solution. In addition, no significant differences were recorded by comparing plant dry weight at sampling time (Fig. S2, ANOVA followed by Tukey HSD test, Barke p value = 0.432 and Morex p value = 0.9). These data indicate that Quarryfield soil is a suitable substrate to investigate the impact of the exogenous application of gramine and the root-soil interface and that this latter is not associated to obvious pleiotropic effect on barley growth.

Next, we generated a 16S rRNA gene amplicon library from 45 bulk soil and rhizosphere samples exposed to different concentration of gramine. Upon processing of the sequencing reads in silico (Material and Methods) with a protocol comparable to previous studies conducted in the same soil type (Robertson-Albertyn et al., 2017; Alegria Terrazas et al., 2020), we failed to identify a significant effect of the treatment on the alpha-diversity parameters of the tested communities (Fig. 1, Kruskal–Wallis test followed by Dunn’s post-hoc test, Observed genera p value = 0.5267, Chao1 p value = 0.6657 and Shannon p value = 0.6002, upon Bonferroni correction). Conversely, we observed a clear impact of the external application of gramine on the bacterial communities thriving at the root soil interface regardless of the applied concentration: both a Canonical Analysis of Principal Coordinates (CAP, Fig. 2) and a Principal Coordinates Analysis (PCoA, Fig. S3) built on a Bray–Curtis distance matrix revealed a partition of the microbiota associated to the applied treatment. Of note, the effect of gramine appeared more pronounced on bulk soil samples than rhizosphere specimens. Congruently, a permutational analysis of variance computed on both matrices indicated a significant effect of the individual microhabitat, bulk soil, Morex and Barke rhizosphere, respectively (R2 ∼38.6%, Adonis test p value = 0.00020; 5,000 permutations; Table 1) and the interaction term between gramine application and microhabitat (R2 ∼6.6; Adonis test p value = 0.01540; 5,000 permutations; Table 1). Conversely, the significance of the impact of gramine per se failed marginally to pass the threshold imposed (R2 ∼4.9%; Adonis test p value = 0.05559; 5,000 permutations; Table 1). Interestingly, we obtained strikingly similar results when this calculation was performed using individual ASVs (Figs. S4 and S5), suggesting that in the impact of gramine applications on the bacteria proliferating at the root-soil interface is coded by taxonomically conserved portions of their genomes.

The abundance of individual members of the barley microbiota is affected by gramine application

The observation that rhizosphere profiles tend to converge on the computed ordinations (see Figs. 2 and S3), suggested that the two barley genotypes evolved the capacity of reverting, at least in part, the selective pressure of gramine on soil bacteria. To quantify this phenomenon, we implemented pair-wise comparisons between individual genera retrieved from unplanted soil and rhizosphere communities at the three levels of gramine tested. Congruently with the initial observation, we determined that the majority of genera enriched in the rhizosphere of either genotype are comparable across gramine treatments (Fig. 3, Wald test, individual p values < 0.05, FDR corrected). Conversely, 18 genera, belonging to 16 distinct higher taxonomic ranks, whose cumulative relative abundance represented ∼1.66% and ∼1.12% of the Morex and Barke rhizosphere communities respectively, were identified as gramine-responsive in a genotype-independent manner (Fig. 4, Wald test, individual p values < 0.05, FDR corrected). Interestingly, we noticed that the host genotype drives the enrichment of genera whose abundance in gramine-treated bulk soil samples was almost obliterated. When looking at the taxonomic affiliation of these rhizosphere- and gramine-responsive genera we observed a broad phylogenetic affiliation, encompassing members of the phyla Acidobacteria, Actinobacteria, Bacteroidetes and Proteobacteria. Next, we inspected how gramine application impacted on the genotype-driven diversification of the rhizosphere bacterial microbiota, using the number of genera (a) enriched from soil and (b) differentially enriched between the tested genotypes as a readout for this analysis. This analysis revealed that plants exposed to no gramine or the lowest dosage displayed a differential enrichment between cultivars of 31 and 42 genera, respectively (Figs. 5A and 5B, Wald test, individual p values <0.05, FDR corrected). Conversely, at the highest gramine dosage the host genotype effect was limited to 12 differentially enriched genera (Fig. 5C, Wald test, individual p values <0.05, FDR corrected). When we inspected the taxonomic affiliation of these differentially enriched genera, we made two observations. First, members of the phylum Proteobacteria accounted for the majority of differentially recruited genera (14 out of 31 genera at gramine 0 µmol L⁻¹, 24 out of 42 at gramine 24 µmol L⁻¹ and 6 out of 12 at gramine 46 µmol L⁻¹, respectively). The second observation was that despite the number of differentially regulated genera at the intermediate dosage of gramine was comparable to mock-treated specimens, the magnitude of the host selection was modulated by the host genotype itself, as Barke enriched for almost three times the number of Proteobacteria upon gramine application.