The experiments were performed on Nova Scotia Haskap farm located in Belmont (45°26′04.4″N 63°22′56.5″W). A total of 40 Aurora variety Lonicera caerulea plants, also known as blue honeysucklebare, or haskap were used in the experiments. The cuttings were planted 1-gallon pots using peat base soil rich in organic matter on June 10, 2018. Prior to planting, the root samples (around 3 g) of 21 plants were taken for root and rhizosphere microbiome analysis (S0 and R0). A total of 20 soil samples were taken from the pots for soil microbiome analysis (S0). The pots were watered to settle the soil and remove air pockets.

The strain RF67 was isolated from wild blueberry roots (Vaccinium angustifolium Ait) (Dussault, 2019) on MM-NH₄ media (Somerville and Kahn, 1983) and annotated as Rhizobiales based on the 16S rRNA sequencing. The sequence is available in the NCBI GenBank under the accession number OM753896. About 10 ml of inoculum containing 10⁸ CFU/ml of RF67 was introduced into 20 pots and 20 more pots were left without inoculation. The pots were placed in full sun, and nutrition and watering regimes were maintained as needed along with weeding into fall. After 3 months of plant growth, the soil and root samples were taken from each pot. In May 2019, inoculated and non-inoculated plants were randomly planted in the field and labeled for future identification. In September 2019, the root samples from these plants were taken for root microbiome analysis.

The topsoil litter was removed to expose the surface, if necessary. About 5–10 cm depth soil samples were collected using a sterile spatula, placed in sterile bags, and transported to the laboratory on ice. The soil samples were sieved (2 mm) and stored at −80°C for DNA isolation. Around 3–5 g of roots was collected from each plant. The roots were placed in the sterile bags and transported in the laboratory on ice. The root and rhizosphere samples were processed as described previously (Yurgel et al., 2017). About 0.250 g of soil, rhizosphere, and root tissue was set aside for DNA isolation.

DNA extraction was carried out using the PowerSoil DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA) according to the manufacturer’s protocol. DNA quality and concentration were measured using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham). About 5 μl of DNA sample was sent to the Dalhousie University CGEB-IMR¹ for V6-V8 16S rRNA (16S) library preparation and sequencing. Samples were multiplexed using a dual-indexing approach and sequenced using an Illumina MiSeq with paired-end 300 + 300 bp reads. All PCR procedures, primers, and Illumina sequencing details were described as mentioned in the study of Comeau et al. (2017). The DNA was sequenced for prokaryotic V6-V8 16S (ACGCGHNRAACCTTACC forward primer, ACGGGCRGTGWGTRCAA reverse primer) (Srinivasan et al., 2015). All sequences generated in this study are available in the NCBI sequence read archive under the accession numbers PRJNA804723, PRJNA804559, and PRJNA804564.

The overlapping paired-end forward and reverse reads were stitched together using PEAR (Zhang et al., 2014) and exported into QIIME2 (Bolyen et al., 2019). The sequences were trimmed of their primers using QIIME2’s Cutadept plug-in (Martin, 2011; Comeau et al., 2017). Low-quality sequences were filtered from the dataset using QIIME2’s q-score-joined function. Using QIIME2’s Deblur plug-in, the sequences were organized into amplicon sequence variants (ASVs)—high-resolution genomic groupings (Amir et al., 2017; Callahan et al., 2017; Comeau et al., 2017). To account for potential MiSeq bleed-through between runs (estimated by Illumina to be less than 0.1%), ASVs which accounted for less than 0.1% of the total sequences were removed. Taxonomic classifications were assigned to the ASV using QIIME2’s naïve-Bayes scikit-learn function, referencing SILVA databases (Quast et al., 2013; Bokulich et al., 2018). Additionally, ASVs assigned to mitochondria and chloroplasts were filtered out (Comeau et al., 2017).

QIIME2’s diversity function was used to calculate Shannon indices (alpha diversity) and also UniFrac matrices (beta diversity) (Lozupone et al., 2011; Kim et al., 2017). These UniFrac matrices were then subjected to an ADONIS test through which their values were fitted to linear regression to determine what proportion of variance in community structure could be attributed to the treatment. Non-metric multidimensional scaling (NMDS) of bacterial communities was performed on Bray–Curtis matrices using the Vegan R package (Oksanen et al., 2016). Differential abundances bacterial taxa were determined using ALDEx2 (Fernandes et al., 2014) with Benjamini–Hochberg corrected p-value of the Kruskal–Wallis test (p < 0.05). The graphics were produced using ggplot2 (Wickham, 2016). The co-occurrence analysis was performed using the Compositionality Corrected by REnormalization and PErmutation (CCREPE) R package (Schwager et al., 2020) with 1,000 bootstrap iterations and default settings. To obtain comparable datasets from each treatment, 14 replicate samples from R1, R2, R3, and R4 datasets were randomly selected. The co-occurrence and co-exclusion patterns in the samples were scored. The results were filtered to remove non-statistically significant relationships. We generated the network based on the strong correlations with p-values < 0.01. The networks were visualized with Cytoscape (Shannon et al., 2003) and were represented as graphs with microbial functions as vertices or nodes and the edges as interaction types.

C-terminal enhanced green fluorescent protein (eGFP) fusion of CspA2-GFP cloned into pK19 mob sacB was used to label RF67 as described previously (Ogden et al., 2019). Arabidopsis thaliana seeds were surface-sterilized with 70% (v/v) ethanol for 5 min, air-dried, and sown on 1% (w/v) agar plus half-strength Murashige and Skoog’s basal salt medium (#M524, PhytoTech Laboratories), pH 5.7 plates. Following vernalized for 2 days at 4°C, seedling was grown vertically at 22°C with 16-h light–8-h dark at 100 to 140 μE/m². RF67 overnight culture was added on the top of 10-day-old seedlings and incubated for 24 h. Then, roots were rinsed in water, mounted in water, and imaged using Leica SP8 laser confocal scanning microscope equipped with 40 × 1.3NA oil immersion objective, 488 nm excitation, and 494-535 emission range.

Cells from 2- to 3-day-old plates were resuspended in Min-salt solution to OD600 = 0.5, and the cell suspensions were diluted with Min-salt solution at 1–10, 1–100, 1–1,000, 1–10,000, and 1–100,000 times in a 96-well microplate. Aliquots of these suspensions were then transferred using a sterile bolt replicator on the plates containing solid nutrient agar (Abdulkadir and Waliyu, 2012) with some modifications. More specifically, since RF67 was very sensitive to NaCl, it was removed from the recipe, and the pH for the media was adjusted to 5, 6, 7, or 8. The colony size was scored after 2–5 days. The endophytes isolated from apple roots and annotated as rhizobium, as well as Sinorhizobium meliloti strain 1021 (Galibert et al., 2001), were used as controls.

Since RF67 was isolated from wild blueberry root, which usually grows in high acidity soils (Korcak, 1989), we tested the ability of the strain to grow on the plates with a range of pH between 4 and 8. We did not detect any inhibitory effect of pH up to 7 on the growth of strain (Supplementary Figure 1), suggesting its potential tolerance to less acidic environments. To evaluate the ability of RF67 to establish symbiosis with plants, the strain was labeled with GFP protein expressed from cspA2 promoter, which is highly active in rhizobia (Ogden et al., 2019). Inoculation of A. thaliana roots resulted in the colonization of apoplast of root apical meristem (Figures 1A–C) and in the root differentiation zone where bacterial cells were detected at the root center proximally to xylem cells (Figures 1D–F). Reconstructing transverse (Z-plane) sections (Figure 1G) through the root differentiation zone also confirmed the presence of bacteria three cell layers below the root surface (Figure 1H). Individual bacterial cells in the apoplast were observed under higher digital zoom settings (Figures 1H–J).

The dataset retained a total of 12,916 features encompassed by 1,755,900 reads spread across 209 samples, with a mean frequency of 8,401 reads per sample and a median frequency of 5,643 reads per samples. For normalization purposes, the samples were rarefied to a depth of 2,660 reads per sample with total 462,860 reads. This process removed 35 samples with insufficient depth producing 174 samples which included: roots from 21 bare root cuttings (R0); roots from 17 inoculated and 19 non-inoculated plants after 3 months of growth in pots (R1 and R2, respectively); roots from 14 non-inoculated and 15 inoculated plants after 1 year of field growth (R3 and R4, respectively); rhizosphere from 14 bare root cuttings (RS0); rhizosphere form 16 inoculated and 12 non-inoculated plants after 3 months of growth in pots (RS1 and RS2, respectively); 13 samples of soil used for planting (in the time of planting, S0); and 18 inoculated and 15 non-inculcated soil samples from the pots after 3 months of plant growth (S1 and S2, respectively). The final dataset of 12,677 features containing Alphaproteobacteria, Bacteroidia, Gammaproteobacteria, Actinobacterian, and Polyangia as the most abundant bacterial classes represented by 27, 18, 16, 7, and 5% of total microbiome reads, respectively (Supplementary Figure 2).

We did not detect significant changes in alpha diversity (Shannon) in the root microbiome at the end of 3-month growth. However, after a year of growth in the field, the Shannon diversity of the microbiome (R4 group) was significantly lower compared to that of R0 and R2 microbiomes (Figure 2). The structure of the root microbiome was visually different when it was grouped based on the time of plant growth (Figure 3), with 47% of community dissimilarity explained by the period of plant growth (R0 vs. R2 vs. R4, R² = 0.47, p < 0.01; Table 1). We identified 59 bacterial classes differentially represented between the root microbiomes from the plants with different growth stages (Supplementary Table 1). These taxa represented by at least 1% of total 16S reads are shown in Figure 4A. Most notable classes included Rhizobiales and Polyangiales, which were underrepresented in R2 group; and Proteobacteria Burkholderiales, Caulobacterales, Xanthomonadales, and Pseudomonadales, Bacteroidetes Cytophagales and Flavobacteriales, and Sphingobacteriales, which relative abundances decreased over time (Figure 4A). The relative abundances of Actinobacteria Micromonosporales, Streptomycetales, and Propionibacteriales, as well as Bacillales, were increased over time in the root microbiome (Figure 4A).

Rhizosphere microbiome also underwent a significant transformation over the period of first 3-month plant growth, including an increase in Shannon diversity (Figure 2) and the significant dissimilarity in rhizosphere sample grouping (RS0 vs. RS2, R² = 0.27, p < 0.01; Table 1 and Figure 3). Additionally, bacterial classes Polyangiales and Gammaproteobacteria CCD24 were depleted in the rhizosphere of the plants after 3-month growth compared to the bare root plants (Figure 4B).

After 3-month plant growth, we detected an increase in the bacterial alpha-diversity in the pot’s soil (Figure 2), which was reflected in increased Shannon diversity in S0 compared to the S2 group. Visualization of dissimilarity between soil communities over time revealed a visible trend in beta diversity (Figure 3, S0 vs. S2), and the analysis of strength and statistical significance of sample groupings (ADONIS test) indicated that the time was a significant factor shaping bacterial community (S0 vs. S2, R² = 0.50, p < 0.01; Table 1). We also detected variations in the relative abundances of several bacterial taxa between initial (S0) and 3-month-old (S2) soils (Figure 4C). More specifically, Chitinophagales, Sphingomonadales, Xanthomonadales, Burkholderiales, Blastocatellales, and Gammaproteobacteria were overrepresented, whereas Polyangiales, Haliangiales, and Pedosphaerales were underrepresented in S0 compared to S2 soils (Figure 4C).

We used the RF67 16S rRNA sequence to detect a corresponding ASV, which was identified as 2280c05c4198790e14682350fb738135 and annotated as Alphaproteobacteria Allorhizobium–Neorhizobium–Pararhizobium–Rhizobium group. The ASV had a low relative abundance in the non-rarefied ASV table with an average of 0.22 reads per sample. One-way ANOVA showed a significant increase in the number of 2280c05c4198790e14682350fb738135 reads in the R1 group (Tukey, p < 0.05; Supplementary Table 2). Additionally, this ASV was detected in S1, RS1, and R3 groups. However, the ASV was also sporadically detected in non-inoculated R0 and R2 groups (Supplementary Table 3). More specifically, 3, 4, 2, and 1 reads corresponding to 2280c05c4198790e14682350fb738135 were detected in samples BF13(R0), BF5Rall(R0), AF22Rall(R2), and AF26Rall(R2), respectively.

The introduction of RF67 into pots’ soil resulted in the significant decrease of alpha diversity in root microbiome after 3-month plant growth (Figure 2). Shannon diversity in the inoculated roots (R1) was 7.3 vs. 8.4 in the non-inculcated roots (R2). NMDS plots showed a strong visual separation between R1 and R2 groups (Figure 3). This visual community separation was supported by the analysis of strength and statistical significance of sample groupings, indicating that the inoculation was a significant factor shaping bacterial community after 3-month growth (R1 vs. R2, R² = 0.23, p < 0.01). Additionally, families, such as Enterobacteriaceae and Rhizobiaceae, were overrepresented in inoculated roots (R1) compared to non-inoculated roots (R2) (Figure 5A).

Bulk soil microbiome responded strongly to the introduction of inoculum. The strength and statistical significance of sample groupings indicated that the inoculation was a significant factor shaping soil microbiome after 3-month growth (S1 vs. S2, R² = 0.14, p > 0.001) (Table 1), and some visual separations between S1 and S2 groups were detected in NMDS plots. Interestingly, a number of bacterial families had decreased relative abundances in inoculated soils compared to non-inoculated ones, including Alphaproteobacteria Xanthobacteraceae, Sphingomonadaceae, Reyranellaceae, Beijerinckiaceae, and Rhizobiaceae, Gammaproteobacteria Xanthomonadaceae and Burkholderiales SC-I-84, Verrucomicrobiae Chthoniobacteraceae, Thermoleophilia Gaiellales, Nitrospiria Nitrospiraceae, and Chloroflexi Gitt-GS-136 (Figure 5B).

In contrast to the root and bulk soil microbiomes, the inoculation with RF67 did not have detectable effect on rhizosphere microbiome. More specifically, we did not detect any significant differences in bacterial alpha diversity (Figure 2) with an estimated Shannon diversity in the inoculated rhizosphere (RS1) at 9.2 and in the non-inculcated rhizosphere (R2) at 9.3. The strength and statistical significance of sample groupings indicated that the inoculation was not a significant factor shaping rhizosphere microbiome (RS1 vs. RS2, R² = 0.05, p > 0.05) (Table 1), which was correlated with no visual separation between RS1 and RS2 groups in NMDS plots (Figure 3). Furthermore, no bacterial taxa were differentially represented between RS1 and RS2 groups.

We analyzed the co-occurrence pattern in root microbiomes. Based on the number of taxa associated with clusters representing strong positive interaction, root microbiome from non-inoculated roots exhibited less cooperation compared to inoculated roots (Figure 6). Then, 3 months post-inoculation, network from inoculated roots (R1) comprised of 150 nodes (taxa) with 420 edges (interactions) and maximum node degree (the number of edges connected to the node) 24, compared to 148, 242, and 14, respectively (Supplementary Tables 4, 5). Additionally, in the non-inoculated microbiome (R2), Firmicutes were the major taxa with strong cooperation among them (Supplementary Table 5), whereas Proteobacteria, Bacteroidota, Gemmatimonadota, and Actinobacteriota form a strong cooperation in inoculated roots (Supplementary Table 4).

After 1 year of field growth, we did not detect any differences in alpha-diversity between inoculated and non-inoculated roots (R3 vs. R4, Figure 2). There was no significant community variation explained by the inoculation (R3 vs. R4, R² = 0.06, p > 0.05) (Table 1), and NMDS analysis indicated strong visual similarity between the groups. Additionally, no bacterial taxa were differentially represented between inoculated and non-inoculated roots after 1 year of field growth. However, similar to network from 3 months post-inoculated roots, after a year post-inoculation, network from inoculated roots (R4) exhibited much stronger cooperation, compare to non-inoculated roots (R3). More specifically, after a year post-inoculation, R4 network comprised of 155 nodes with 460 edges and maximum node degree 22, compared to 164, 267, and 13 forms R3 network, respectively (Supplementary Tables 6, 7). In the inoculated roots, Actinobacteriota Gaiellales and Acidimicrobiia IMCC26256, Bacteroidota Niastella, Firmicutes Bacillales, Planifilum and Bacillus thermolactis, Gemmatimonadota Gemmatimonas, Proteobacteria Pseudolabrys, and Halanaerobiaeota Halocella were the most connected taxa with 17 or more degrees (Supplementary Table 6).