The microbial inoculants used in this experiment were chosen for being able to maintain a high level of viability of microorganisms on granular fertilizer. The inoculants were Éra boost® (EB) or EVL Coating® (EVL). EB is composed of five Bacillus strains that were screened for their plant-growth promoting effects. As all strains in the EB consortia can sporulate, EB microorganisms can resist hostile environments such as granular fertilizer surface. EVL is composed of a consortium containing two Bacillus, one Lactobacillus, one Pseudomonas, and one Saccharomyces strains all combined in a bioreactor along with a proprietary substrate. The EVL consortium is added to a formulation designed for coating granular fertilizer (Garcia-Fraile et al., 2015).

A field experiment was conducted on a fine sandy soil located in St-Thomas de Joliette (46°00′01.8″N 73°21′02.2″W) during 2019. Three treatments were investigated: (1) un-inoculated control (mineral fertilizer only), (2) EB or, (3) EVL inoculants. Each product was coated onto granular starter fertilizer (GSF) applied at seeding. For all treatments, the rate of GSF was 1,050 kg ha⁻¹ with an NPK equivalent of 12–14–10 kg T⁻¹ of GSF (Synagri, St-Hyacinthe, QC, Canada). The microbial inoculant was applied at 1 L T⁻¹ of fertilizer for EB and 2 L T⁻¹ of fertilizer for EVL, following the manufacturer's recommendations. The fertilizer was banded as a 5 × 5 × 5 cm application, in such a way that fertilizer was placed 5 cm to each side of the seed tuber and both strips were 5 cm below the potato seed tuber. All treatments were replicated seven times and the experiment was structured following a randomized complete block design; four of the seven blocks were used to sample bulk soil and potato roots with closely adhering soil.

Each bulk soil sample consisted of four initial samples taken at the extremities of a 2 × 2 m² quadrant, pooled together. Bulk soil samples were collected directly on the potato row, between the seed tuber and the fertilizer rows. Within each quadrant, a 2.5 cm-diameter probe was used to sample soil to a depth of 30 cm. Samples were collected at four time points: (1) 4 days before seeding (May 27th), (2) 7 days after seeding (June 7th), (3) 21 days after seeding at shoot emergence stage (June 21st), and (4) 56 days after seeding at mid-flowering stage (July 26th). At the fourth time point, rhizospheric samples were also collected. Each rhizospheric sample consisted of fine roots and closely adhering soil sampled from three plants randomly selected from each plot. Samples of bulk and rhizospheric soils were immediately placed on ice until the arrival at the laboratory.

For each composite bulk soil sample, a representative subsample of about 100 g was prepared and stored at −20°C. Around 400 mg of the 100 g were collected in 2 mL collection tubes before DNA extraction. For each rhizospheric composite sample, a subsample of homogenized roots and closely adhering roots of 1 g was ground with liquid nitrogen. After grinding, 200 mg were collected and stored at −20°C in 2 mL collection tubes, for DNA extraction.

Analysis of soil chemical properties was conducted by Laboratoire Géosol (St-Hyacinthe, QC, Canada). P, K, Mg, Ca, Al, Zn, Mn, Cu, Fe, and B were extracted with Mehlich III solution (Mehlich, 1984) and quantified by Inductively Coupled Plasma Optical Emission Spectroscopy (ICP-OES, Perkin Elmer Optima 4,300 DV, Shelton, CT). Soil pH_water, pH_SMP and organic matter content were estimated with a 1:1 dilution in water (Hendershot et al., 1993), SMP (Shoemaker et al., 1961) and Walkley-Black solutions (Walkley and Black, 1934), respectively.

All DNA extractions were performed using a FastDNA Spin Kit for Soil (MPBio, Irvine, CA, USA). Bulk soil and rhizospheric subsamples were added to tubes containing 1.4 g of the bead matrix E and 1 mL of the lysis buffer supplied with the kit. Microbial inoculant samples consisted of the centrifuged cells contained in 5 mL of each commercial product. DNA extraction was performed according to the manufacturer's instructions. The resulting DNA pellet was suspended in 100 μL of sterile molecular-grade water.

The quality and quantity of the DNA extracts were evaluated by spectrophotometry using a Biophotometer (Eppendorf, Mississauga, ON, Canada) with a G1.0 μCuvette (Eppendorf, Mississauga, ON, Canada) with readings at 260, 280, 230, and 320 nm. The V4 region of prokaryota (archaea and bacteria) rRNA 16S gene was amplified using 515FB and 806RB primers (Apprill et al., 2015; Parada et al., 2016). For the fungi, the eukaryotic (fungal) ITS1 gene was amplified (Bokulich and Mills, 2013). Both genes were amplified in a two-step dual approach PCR designed for Illumina instruments by Plateforme d'analyses génomiques (IBIS, Université Laval, Quebec City, QC, Canada). DNA sequencing was performed by IBIS on an Illumina MiSeq platform, following the methods of Jeanne et al. (2019).

Obtained sequences were demultiplexed based on the tag used. Sequence quality control and feature table construction was performed using QIIME2 (Caporaso et al., 2010) and DADA2 plugin (Callahan et al., 2016). Reference databases Greengenes 13.8 (Desantis et al., 2006) and UNITE version 8 (Kõljalg et al., 2013) were used for taxonomic identification of amplicon sequence variants (Callahan et al., 2017).

From the isolated DNA, 4 μL of DNA diluted at 1:20 in molecular-grade water was mixed with 6 μL of a master mix of SYBR green qPCR mix (Qiagen, Toronto, ON, Canada). Samples were prepared on 96-well plates with four water-filled wells acting as blanks. PCR conditions were 15 min at 95°C, followed by 40 cycles of 95°C for 1 min, 30 s at the annealing temperature and 72°C for 1 min. For prokaryotic DNA, the primer pair EUB-338/518 was used with an annealing temperature of 53°C (Fierer et al., 2005). For fungal DNA, the primer pair FF390/FR1 was used with an annealing temperature of 51°C (Emerson et al., 2015) all samples were replicated at least twice. Standard curves were generated from a known quantity of amplified DNA fragments, diluted over a 4-log range (efficiency of 89.1% and r² = 0.99 for prokaryotes; efficiency of 91.7% and r² = 1 for fungi). Fluorescence detection was performed on a CFX96 (Biorad, Hercules, CA, USA).

Amplification units (AU) were derived from the average CT values of samples using linear regressions designed by the Microbial Ecology Laboratory of Institut de Recherche et Développement en Agroenvironnement (IRDA) in which samples with known concentrations of DNA were used to fit the regressions. AU g⁻¹ of dry soil were obtained using moisture content and the weight of soil used during DNA extraction.

The entire analysis was performed using the R software with version 4.3 R, Core Team (2020). Before doing alpha and beta diversity analysis, normalization of the number of reads was performed using the rarefy-even-depth function from the phyloseq R package (Mcmurdie and Holmes, 2014).

Alpha diversity was measured by Shannon and Chao1 indices, evenness and number of observed ASVs. A three-way ANOVA was performed on alpha-diversity measurements of soil samples with microbial inoculants product, sampling time and block numbers as the three effects. A two-way ANOVA was performed on alpha-diversity measurements of rhizospheric samples with microbial inoculants product and block as effects. Bonferroni adjustment was used to determine significant differences during post-hoc testing.

To compare alpha-diversity measurements and relative abundances of certain microbial genera between rhizospheric and bulk soil communities at the same sampling date (July 26th), a paired T-test of rhizospheric against bulk soil samples was performed.

Beta-diversity was calculated using Bray-Curtis dissimilarity matrices with the ordinate function from the phyloseq package. Ordination of bulk soil and rhizospheric samples was determined using principal coordinate analysis (PCoA). Distance matrices were computed using the distance function from phyloseq. PERMANOVA was performed on these distance matrices using the Adonis function from the Vegan package with the explanatory variables set as microbial treatment, sampling date and treatment block number and 999 permutations. A post-hoc test was conducted using the pairwiseAdonis package with p-values adjusted using a Bonferroni adjustment.

A distance-based redundancy analysis (dbRDA) was performed using the “capscale” function from the Vegan package. Prokaryotic and fungal ASV features table was used with all soil chemical variables measured. PERMANOVA of soil chemical variables revealed statistically significant variables. These variables were represented as vectors on a constrained ordination for both prokaryotic and fungal communities.

Differentially abundant ASVs between the first block and the other three blocks were identified for bulk soil and rhizospheric soil. To identify these ASVs, the DESeq2 package was used (Love et al., 2014). To convert the biom files in DESeq files, the phyloseq_to_deseq2 function from phyloseq was used.

Potato tuber yield expressed as T ha⁻¹ were analyzed following a two-way ANOVA with microbial inoculant and block number as explanatory variables. Physicochemical properties of soil such as CEC, organic matter content, pH and P, Al, K, Mg, Ca, Zn, Mn, Cu, Fe, and B contents were analyzed following a three-way ANOVA with microbial inoculant, block number and sampling time as explanatory variables. The Tukey post-hoc test was used to differentiate significantly different pairwise comparisons for yield and physicochemical values.

Potato tuber yield (Table 1) was not significantly impacted by microbial inoculation of granular fertilizer (P = 0.396). The block number of each plot had a significant effect on tuber yield (P = 0.001). Tukey post-hoc test revealed that block 1 had significantly lower yields than block 2 (P_adj = 0.004), 3 (P_adj = 0.002), and 4 (P_adj = 0.004) while blocks 2, 3, and 4 were not statistically different from each other.

The HTS analysis of EVL bacterial and fungal mixed inoculant provided 17,110 and 16,311 reads, respectively, for prokaryotic and fungal communities. EVL inoculant was quite diverse with 30 ASVs encompassing 15 prokaryotic and 15 fungal ASVs. In EVL inoculant the major prokaryotic genera were Lactobacillus and Clostridium, and the major fungal genus was Saccharomyces (Figures 1A,B). The HTS analysis of EB bacterial inoculant provided 21,950 reads for the prokaryotic community. The eight ASVs identified in EB all belonged to the order Bacillales with seven Bacillus and one Paenibacillus (Figure 1A).

The average number of reads for soil samples was 30,068 ± 1,905 for prokaryotes and 22,717 ± 2,341 for fungi. Standardization of the number of reads resulted in 16,685 and 12,302 reads for prokaryote and fungal groups, respectively.

Prokaryotic and fungal genera present in EVL and EB inoculants were not significantly enriched in bulk soil microbial communities following microbial inoculation at seeding or in comparison to un-inoculated plots. Although the Bacillus genus was present in bulk soil samples, the relative proportions of this genus were not significantly different in EB treated plots compared to those in control plots. The predominant prokaryotic genus in the EVL consortium, Lactobacillus, was not detected in any soil samples. Similar results were observed for fungal genera present in the EVL consortium, in that the genera Saccharomyces and Wickerhamomyces were not observed in bulk soil samples (Figures 1C,D).

Microbial inoculation of fertilizer did not have a significant effect on any alpha-diversity indices, observed ASVs or evenness of soils at any time point (Supplementary Table 1). Sampling date had a significant effect on Chao1 (P_prok = 0.003, P_fung < 0.001) index and evenness (P_prok = 0.01, P_fung = 0.02) for prokaryote and fungal populations, respectively (Figure 2). Significant differences were observed in prokaryotic diversity between May 27th and June 21st for Chao1 index and between May 27th and July 26th for Chao1 index and evenness (Figure 2). Fungal Shannon and Chao1 indices and observed ASVs were significantly different between July 26th and all other sampling dates, while there were no significant differences in evenness for May 27th and July 26th.

The two bacterial phyla showing highest relative proportions were Proteobacteria and Actinobacteria across all sampling time points, while the fungal classes showing highest relative proportions were Sordariomycetes, Leotiomycetes, and Eurotiomycetes, with the latter having a noticeable enrichment at the final sample date (Figures 3C,D).

A principal coordinate analysis (PCoA) was used to visualize the community assemblages at the ASV level using Bray-Curtis dissimilarities (Supplementary Figure 1). The prokaryote community displayed a strong clustering of samples along the first PCoA axis, which explained 17% of the variation, with strong loading of sampling dates; the fourth sampling time date (July 26th) was segregated from the first three. The second axis, which explained 6.9% of the variation, showed a clear clustering of physical location indicated as block number. Fungal communities displayed a weak clustering effect of sampling time along the first axis, which explained 16.1% of the variation. The second axis does not represent any relevent biological variables and explained 9.4% of the variation (Supplementary Figure 1).

The PERMANOVA results indicated that microbial inoculants did not have a significant effect on bulk soil microbial community while sampling date and physical location (block number) significantly affected prokaryotic and fungal communities. The pairwise differences between sampling date revealed that May 27th prokaryotic composition was significantly different from that of June 21st (P = 0.01), and the July 26th prokaryotic composition was significantly different from the three other sampling dates (P = 0.006). Fungal composition was significantly different on July 26th from the three other sampling times (P = 0.006). The first block had a significantly different prokaryotic composition from the three other blocks (P = 0.006), as observed in the PCoAs (Supplementary Figure 1). Fungal community composition was significantly different between blocks 1 and 2 (P = 0.036).

Both alpha-diversity indices and beta-diversity showed a significant shift in prokaryotic and fungal community compositions between the third and fourth sampling dates. Many ASVs showed differential expression between June 21st and July 26th with the majority belonging to the classes Gammaproteobacteria, Actinobacteria, and Betaproteobacteria (Figures 4A,C). Additionally, many prokaryotic ASVs were unique to each block for bulk and rhizospheric prokaryotic communities (Figures 4B,D).

The microbial inoculation, sampling dates or sampling blocks for either prokaryotic or fungal communities did not significantly impact the number of prokaryotes and fungi detected in bulk soil. The average number of prokaryotes and fungi were, respectively, 3.57 × 10⁹ and 3.39 × 10⁷ AU g⁻¹ dry soil (Supplementary Table 1).

The physico-chemical properties of bulk soil sampled prior to seeding were used to determine if there are links with soil microbial communities (Supplementary Table 2). PERMANOVA following a dbRDA of prokaryotic and fungal communities' composition and soil chemical variables showed a significant effect of organic matter (P = 0.007), CEC (P = 0.001), Mg (P = 0.019), P (P = 0.001), and pH_SMP (P = 0.002) measurements on prokaryotic community and of CEC (P = 0.001) on fungal community (Figures 3A,B).

Rhizospheric samples had an average of 18,822 ± 1,382 and 14,718 ± 1,720 reads for prokaryotes and fungi, respectively. Similar to the reads obtained for the bulk soil samples, the number of reads normalized to 16,685 and 12,302 reads for prokaryote and fungal communities, respectively.

The highest relative proportion among potato rhizospheric prokaryotic phyla were detected for Proteobacteria, Cyanobacteria, and Actinobacteria. The fungal classes Sordariomycetes, Leotiomycetes, and Eurotiomycetes showed the highest relative proportions in rhizospheric soil (Figures 5C,D). The bacterial phylum Cyanobacteria is mostly composed of the Chloroplast genus which most likely come from lysed root cells containing plastids.

The microbial inoculation treatment showed no significant difference for Shannon and Chao1 indices, evenness and number of observed ASVs for either prokaryotic or fungal communities (Supplementary Table 1). The microbial inoculation treatment did not significantly impact the number of prokaryotes and fungi detected in rhizospheric samples. The average number of prokaryotes and fungi were, respectively, 1.27 × 10¹⁰ and 9.17 × 10⁷ AU g⁻¹ dry soil (Supplementary Table 1).

As observed with the bulk soil samples, the microbial inoculant genera were not enriched in rhizopheric microbial community. Lactobacillus, Saccharomyces, and Wickerhamomyces genera, found in the EVL inoculant, were not observed in rhizospheric samples. However, the Bacillus genus, which is the dominant bacterial genus in EB inoculant, showed a significant mean increase of 68% in rhizospheric samples compared to those detected in bulk soil (t₍₁₁₎ = 8.3, P < 0.0001) for inoculated and non-inoculated plots.

The prokaryotic and fungal communities' diversity for bulk soil samples collected on the same date as rhizospheric soil samples (July 26th) were analyzed as paired T-tests. Prokaryotic Shannon (P < 0.001) and Chao1 (P < 0.001) indices, evenness (P < 0.001) and number of observed prokaryotic ASVs (P < 0.001) were significantly smaller in rhizospheric soil compared to bulk soil. Fungal Chao1 index (P = 0.003) and number of observed fungal ASVs (P = 0.003) were significantly smaller for rhizospheric soil than bulk soil. Mean numbers of prokaryotes and fungi estimated as AU g⁻¹ dry soil were significantly greater in rhizospheric soil.

PERMANOVA analysis of rhizospheric and bulk soil samples collected on July 26th showed a significant effect (P < 0.001) of sample type (roots or soil) for prokaryotic and fungal communities. Microbial inoculation and block location did not have significant effects on the prokaryotic community (P_treatment = 0.43 and P_block = 0.15) and on fungal community (P_treatment = 0.76 and P_block = 0.36). Analysis of rhizospheric soil samples alone showed a significant effect of block location (P = 0.019) on prokaryotic community, but no significant effect (P = 0.415) on fungal community.

We hypothesized that microbial inoculation of granular fertilizer would directly influence soil microbial community by increasing relative abundances of inoculated microbial genera. We had assumed that since HTS technique and the resultant ASV can discriminate microbial genera with a higher level of precision thantechniques such as culture-dependent techniques or phospholipid-derived fatty acids, the use of HTS would detect changes in the relative abundance of the inoculated microbial genera. However, our HTS results showed that neither microbial communities of bulk soil nor rhizospheric samples were significantly affected by microbial treatments. Previous experiments using an HTS approach have successfully detected an impact from inoculants (Deng et al., 2019; Dong et al., 2019). An experiment on Fragaria × ananassa (strawberry) showed that VESTA®, a soil amendment composed of a broad spectrum of microbes, fermentation by-products and organic acids, applied at a rate of 75 L ha⁻¹ every month for 5 months, had an effect on bulk soil prokaryotic community (Deng et al., 2019). Inoculation of P. ginseng with a one time-application of 90–270 L ha⁻¹ of microbial inoculant consisting of Burkholderia and Rhizobium diluted in pig and cow manure as a biofertilizer significantly impacted the bulk soil prokaryotic community (Dong et al., 2019). An inoculation of 10⁹ spores L⁻¹ of B. subtilis in 12 different Chinese soils cultivated in a greenhouse with Cucumis sativus (cucumber) did have an effect on the rhizospheric prokaryotic community sampled 80 days after inoculation (Tian and Gao, 2014). The fertilizers coated with the microbial inoculants used in this experiment, were applied as narrow strips across the field at seeding, making their effect extremely localized within the soil. Additionally, the banded application rates of both microbial inoculants were 1–2 L ha⁻¹ (about 4–8 × 10¹¹ cells ha⁻¹) which is drastically lower than the application rates used in previous experiments. That may explain the lack of direct effects on bulk and rhizospheric soil microbial communities.

Relative abundances of the microbial inoculants of each commercial formulation in bulk soil and rhizospheric samples do not follow the same pattern. EB inoculant consisted almost exclusively of Bacillus genus. The ASVs present in EB were significantly more abundant in rhizospheric soil than bulk soil, even tough EB inoculation did not increase their abundance. This result indicates that microorganisms present in EB are recruited by plants, a sign of their potential for positive effects on plant growth (Islam et al., 2016; Backer et al., 2018). On the other hand, ASVs identified in EVL were never detected in bulk and rhizospheric samples. Their absence in our samples may indicate that they cannot compete with native soil microbial community (Khare and Arora, 2015).

In this experiment, bulk soil cores were collected from zero to 30 cm depth while the fertilizer was banded 5 cm on each side of the seed tuber and both strips were 5 cm below the potato seed tuber. A dilution effect on the microbial community caused by the bulk soil sampling depth compared to the depth of the inoculant may have impaired our capacity to detect the microbial inoculant genera. Most research has looked at the rhizospheric community instead of the bulk soil community. In this experiment, rhizospheric community samples were only taken the same day as the fourth bulk soil sampling. The root system was not well-developed at the time of the first and second samplings, and at the third sampling, rhizospheric soil was not sampled. In hindsight, analysis of rhizospheric microbial community at this stage of root development would be an ideal moment to investigate microbial inoculant effect on microbial community.

Aspects of sample collection, rather than microbial inoculation, had significant effects on microbial community composition, namely sampling date and sample location (block number). Changes in bacterial and fungal communities over the course of a growing season have been observed repeatedly on agricultural lands (Houlden et al., 2008; Lauber et al., 2013; Degrune et al., 2017). Even with minimal human inputs, microbial community composition changes are correlated with changes in soil temperature and/or humidity (Lauber et al., 2013). In addition to changes in environmental conditions, application of pesticides and in-season tillage activities influence the presence of microbial species in certain layers (strata) of the soil, as well as their relative abundance (Lo, 2010; Degrune et al., 2017; Sun et al., 2018).

In this experiment, a significant shift could be observed for both prokaryote and fungal communities between potato shoot emergence and mid-flowering, much more than between tuber seeding and shoot emergence. Examining crop management activities, between shoot emergence and mid-flowering, fungicides were applied at a weekly interval and soil was tilled once around potato rows, to promote tuber production and apply additional fertilizer. In comparison, between seeding and shoot emergence, when the first three sampling time points occurred, only seeding of seed tubers and application of banded fertilizer were carried out; there was no pesticide application or tillage. Therefore, the changes observed in microbial community composition seem to be more important during the mid-season than the early season in field potato production, around the same time when pesticides applications and tillage are frequently conducted.

Since this experiment followed an RCBD design, we had physical treatment replicates pooled into seven blocks across a field of <1 ha. Our hypothesis that microbial inoculation would have a stronger effect than the block effect on microbial community composition was rejected. Prokaryotic community composition and soil organic matter, P, Ca, Mn, Fe, and B contents were significantly different between the first block and other blocks. A db-RDA revealed that P content significantly explained prokaryotic community composition. The C:N:P ratio was found to be a strong predictor of bacterial diversity in previous research (Delgado-Baquerizo et al., 2017), further indicating that P content in bulk soil can be an important driver of soil bacterial community composition.

Interestingly, potato tuber yield was also significantly different between the first block and the three other blocks. Therefore, the changes in microbial community observed are associated with changes in potato yield. Previous work has shown that potato yield was correlated with microbial community indices (Jeanne et al., 2019) and that these indices could be used to estimate potato productivity. Here, we report ASVs that have relative abundances significantly different between two different yield productivity zones. ASVs with higher abundances in high yield regions (Figure 4C) could be used as potential microbial inoculants that could boost potato tuber yields.

Rhizospheric microbial community composition differed substantially from bulk soil microbial community composition (Figure 5). Enrichment of the rhizospheric community over bulk soil for Proteobacteria and in Eurotiomycetes was evident. Such differences in microbial community associated with potato roots have been reported in numerous other studies (Smalla et al., 2001; Berendsen et al., 2012; Xue et al., 2018).

According to Pfeiffer et al. (2016) rhizospheric microbial community consist of a core community along with a dynamically changing portion which varies among plant locations and growth stages. The dynamically changing portion of the rhizospheric community is not necessarily linked with changes in bulk soil community for potato (Inceoglu et al., 2011), soybean (Sugiyama et al., 2014), wheat (Donn et al., 2015), and wild oats (Shi et al., 2015). Instead, the rhizospheric community is shaped by plant development and location, in concordance with the results observed here and in other publications (Nahar et al., 2020). In this experiment, Actinobacteria were significantly enriched in rhizospheric soils of plots with a high tuber yield. This class of prokaryotes could be utilized as plant-growth promoting microbial inoculants, especially organisms belonging to genera Streptomyces and Mycobacterium (Sathya et al., 2017). Understanding the functions related to plant growth contributed by the core and the changing portions of rhizospheric microbiome could help in designing microbial inoculants that would provide beneficial effects under a wide variety of soil conditions.