Introduction

The rhizosphere is a unique biological and chemical soil environment where some soil microbial communities attracted by the carbon resources as root exudates and rhizodeposition [1, 2]. In the rhizosphere, microbial consortia comprising fungi, bacteria, and archaea can form complex biological interactions that affect plant growth [2]. Among rhizosphere microbes, mycorrhizal fungi are directly associated with the roots of ca. 80% of terrestrial plants [3]. Fungi receive plant-derived photosynthetic products and provide their hosts with soil nutrients and water [4, 5]. Other rhizosphere microbes, known as plant growth-promoting rhizobacteria, positively affect plant growth by enhancing nutrient uptake and improving disease resistance [6]. Some of these bacteria are called mycorrhizal helper bacteria (MHB) because they contribute to the maintenance and function of mycorrhizal associations [7, 8].

To date, many MHB have been reported to function in promoting mycorrhizal fungal activities, improving soil nutrient acquisitions, and inhibiting pathogenic infections [7,8,9]. These reports focused mainly on two major arbuscular mycorrhizal and ectomycorrhizal (ECM) plants, which are found in diverse plant species that are geographically widespread on Earth [3, 8]. Positive effects on plants due to the cooperation between ECM fungi and MHB have been reported in inoculation experiments [10]. Moreover, rhizosphere bacterial communities, even as small as single ECM root tips, vary depending on the presence or absence of mycorrhizal fungal colonization [11, 12]. These findings suggest that cooperative bacteria in plants may be selectively screened by mycorrhizal fungal metabolic products [1]. However, the subsequent shifts of the bacterial community shaped by mycorrhizal fungi have not been investigated, despite their importance for the maintenance of mycorrhization, nutrient acquisition, and against pathogenic microbes. The rhizosphere bacterial community may shift as a result of mycorrhizal fungal degeneration caused by root turnover. An ecological importance of historical impacts in community assembly was conceptualized [13], but a recent review noted the lack of standardized study frameworks in ECM-bacterial interactions [14]. One potential reason is that monitoring the colonization progress of ECM roots is hardly applicable visually owing to the formation of fungal mantles. Monitoring of mycorrhizal fungi-bacterial relationships along with root growth is essential to reveal the tripartite mechanistic connections to maintain plant growth, mycorrhizal formation, and bacterial habitat niches.

Arbutoid mycorrhiza is a type of mycorrhizal structures where fungi form hyphal coils inside epidermal root cells [15]. Pyrola japonica (Ericaceae) is categorized as an arbutoid mycorrhiza, lacking fungal mantles. Transparent epidermal cells allows visual assessment of the extent of mycorrhizal development [15, 16]. Their mycorrhizas formed on the plant are presumably attributable to Russulaceae fungi that are known as well-known ECM fungi [16,17,18,19,20,21,22]. Moreover, the mycorrhizal developmental progress varies with the position of the root system from peripheral to proximal positions. For the nutrient acquisition strategy, P. japonica is known to be a mixotrophic plant that obtains carbon sources partially from their associated mycorrhizal fungi other than autotrophic photosynthesis [16, 17, 23]. Because P. japonica has various sets of cells with different levels of fungal colonization, the plant is an ideal model for studying the assemblage patterns of rhizosphere bacterial communities along with mycorrhizal development.

This study aimed to (1) characterize bacterial taxa associated with different conditions of mycorrhizal development and (2) determine whether the rhizosphere bacterial community structure changes with the development of P. japonica. To achieve these, rhizosphere bacterial communities were analyzed at different stages of mycorrhizal development of P. japonica roots that were visually distinguished using amplicon sequencing of the MiSeq system.

Materials and Methods

Study Site

The study site was located in a secondary forest in Mie Prefecture, central Japan (34°44′04″ N, 136°23′53″ E, 130 m a.s.l.), on a gentle south-facing slope. The upper layer was covered by mature broadleaf species Quercus serrata and Quercus glauca, with sporadic Rhododendron kaempferi and Cleyera japonica in the middle layer. The forest floor supported P. japonica along with other herbaceous plants that grow under mature trees. In 2023, the annual mean temperature and annual precipitation were 16.3 °C and 1613 mm, as recorded at the closest weather station (Mie district meteorological station, 34°44′00″ N, 136°31′10″ E), located 11.1 km from the site. The soil was classified as Cambisols, as defined by the United States Department of Agriculture, Natural Resources Conservation Service.

Field Sampling Procedures

At the study site, eight undamaged individuals of P. japonica, with leaves free from insect damage, were collected in August 2023. The selected plants were carefully excised as 25 × 25 × 15 cm soil blocks to avoid damage to the root systems. Individual plants were separated at intervals of 20 m or more. Bulk soil (250 mg) for a partial 16S rRNA analysis was also collected from the surrounding area for each sample. Root systems were excavated from soil blocks, and root cell status was visually categorized into three levels—limited, full, and digested—based on mycorrhizal development under a stereomicroscope (SZX7, Olympus, Japan; Fig. S1) [16]. In brief, “limited” indicates that hyphal coils of mycorrhizal fungi have not yet formed within epidermal cells; “full” indicates that the coils are fully developed within the cells, and “digested” indicates that the coils were not confirmed in the cells seemingly disintegrated (Figs. 1 and S1). For each plant, a total of 3 cm of root fragment (5–10 fragment pieces) from different locations in root systems showing mycorrhizal development was collected for subsequent a partial 16S rRNA analysis (Fig. S1). The collected roots and bulk soils were stored in a − 80 °C freezer until DNA extraction.

Fig. 1
figure 1

Images of three different stages of mycorrhizal development of Pyrola japonica: a a distal root composed of transparent cells without fungal coils, referred to as the limited condition; b a middle position root formed fungal coils (indicated by an arrowhead) within epidermal cells, referred to as a full condition; and c a proximal root aged epidermal cells with degenerated fungal coils, referred to as a digested condition

Full size image

Molecular DNA Analysis

Genomic DNA was extracted using the DNeasy Plant Mini Kit (Qiagen, Tokyo, Japan) for the roots and the FastDNA SPIN Kit (MP Biomedicals, CA, USA) for bulk soils following the manufacturer’s instructions. Due to low bacterial DNA concentration, after the first PCR amplification of the 16S rRNA using primer sets 9f (GAGTTTGATCCTGGCTCAG, [24]) and 1541r (AAGGAGGTGATCCAGCCG, [25]), the second PCR amplification of a partial 16S rRNA, a V4 region, was conducted with the primer set 515 F and 806R [26] using TaKaRa ExTaq (TaKaRa, Ohtsu, Japan) following the manufacturer’s recommendations. Amplified PCR products were sent to Seibutsu Giken Co., Ltd. (Kanagawa, Japan) for amplicon sequencing. Briefly, PCR products were purified with AMPure XP and used for PCR amplification with TaKaRa ExTaq HS (TaKaRa, Ohtsu, Japan) to generate libraries for each sample, and 2 × 300 bp paired-end V4 sequencing was performed using a MiSeq sequencer (Illumina, San Diego, CA, USA).

Data Analysis

Chimeric sequences were removed using dada2 plug-in in Qiime2 ver. 2023.5 [27], and amplicon sequence variants (ASVs) were assigned to bacterial taxa inferred using Greengene ver.13_8 [28]. Raw sequence data were deposited with links to BioProject accession number PRJDB16666 in the DDBJ database. Initially, we removed cyanobacterial sequences from the obtained ASV dataset that were likely derived from the mitochondria. To avoid contamination and sequencing errors, less than 0.02% of all ASVs were excluded as low-frequency occurrences. To assess our sampling methods in representing bacterial communities, sample-based accumulation curves were generated for all treatments using “specaccum()” function in the “vegan” package [29].

Phylogenetic trees were constructed using the neighbor-joining method in MEGA software version 7 [30] to determine the phylogenetic positions of the three most abundant orders of rhizosphere bacteria: Actinomycetales, Acidobacteriales, and Rhizobiales (see “Results”). For tree construction, ingroup sequences were selected if query sequences matched 100% with deposited sequences derived from roots in the GenBank database. The reliability of the tree branches was evaluated using bootstrap analysis with 1000 re-samplings [31]. A phylogenetic tree was also constructed for all 776 ASVs, including sequences from bulk soil, in the same manner as described above, to obtain a newick file for calculating UniFrac distances, as described below.

Following analyses were performed using R software version 4.2.2 [32]. The community data of the filtered bacterial ASVs were treated as a binary matrix for a conservative evaluation. For calculating the unweighted UniFrac distance of ASV sequences, ASV627 was identified as an outgroup using the “pick_new_outgroup()” function (https://github.com/joey711/phyloseq/issues/1539). This function is designed to prevent random rooting and allows exact genetic distances to be measured by selecting appropriate outgroups. Then, UniFrac distances were calculated in the “phyloseq” package [33].

In sequences derived from different mycorrhizal developments and bulk soil, the numbers of both ASVs and UniFrac distances were tested for normality of variance using the Shapiro–Wilk test [34]. As normality was confirmed for all conditions (p > 0.05), the numbers among treatments were compared using a generalized linear model, followed by Tukey’s HSD tests [35]. An UpSet plot was drawn to display the number of unique or shared ASVs among samples using the “ComplexUpset” package [36]. To visualize the variation in community composition among all samples, bacterial ASV communities were scattered on a non-metric multidimensional scaling (NMDS) plot based on the UniFrac distance in the “phyloseq” package [29, 33]. The significance of clustering in the NMDS plot was compared using the “pairwise.adonis()” function in the “pairwiseAdonis” package, which is a modified version of the “adnis2()” function in the vegan package [29], by performing 999 permutations pairwise PERMANOVA tests [37]. The hierarchical clustering of bacterial community similarities based on UniFrac distance is displayed using “hclust()” function in the “cluster” package [38]. The average linkage method was used as the clustering strategy.

The following community analyses were used with not binary data, but abundance read data. A null model analysis quantifying ecological processes for community assembly was conducted by referring to the R script developed by [39] (https://github.com/stegen/Stegen_etal_ISME_2013). The analysis evaluates the relative impact of different five processes: homogeneous selection (homogeneous environmental conditions lead to more similar structures among communities), heterogeneous selection (heterogeneous environmental conditions lead to more dissimilar structures among communities), dispersal limitation (a low dispersal rate leads to increase structural variations among communities), homogenizing dispersal (a high dispersal rate leads to little variation among community structures), and drift (random processes of population via, e.g., birth, death, and reproduction, stochastically alter community structures). Indicator species analysis was performed to determine whether any bacterial taxa were presented among mycorrhizal developments [40]. The analysis combines information about uniqueness in a certain environmental condition and occurrence (fidelity) within replicates in the same treatment; thus, a good indicator species is expected to always occur only in particular communities. The indicator value (IndVal) was outputted using with the “labdsv” package [41].

For all analyses, the significance level was set at p < 0.05, unless otherwise stated.

Results

PCR amplification was successful for all root samples from all mycorrhizal developmental stages, with bulk soil samples used as a control. In the 4,766,675 raw reads (22,570 ASVs) obtained, the phylum Cyanobacteria and 0.02% of the lowest relative abundance of ASVs were excluded, and then, 741,211 reads (776 ASVs) were used for subsequent analysis. As a result, 440, 287, 357, and 310 ASVs were found in the limited, full, digested conditions, and bulk soils, respectively (Fig. S2, Table S1). The sample-based accumulation curves for the bacterial communities observed in this study did not reach a plateau for all treatments, suggesting additional samples can find more ASVs (Fig. S2). After pooling all root conditions, the top three orders were Actinomycetales (11.6%, 76/656 ASV), Acidobacteriales (10.7%, 70/656 ASV), and Rhizobiales (9.0%, 59/656 ASV; Table S1).

The number of ASVs for a full condition (76.0 ± 3.4 ASVs, mean ± standard error) was significantly lower than those in other conditions (limited condition, 96.8 ± 4.8 ASVs; digested condition, 90.9 ± 4.4 ASVs; bulk soil, 95.5 ± 2.6 ASVs; p < 0.05; Fig. 2a). The UniFrac distance for the full condition (0.42 ± 0.014) was significantly lower than that for the other conditions (limited condition, 0.55 ± 0.013; digested condition, 0.48 ± 0.019; bulk soil, 0.50 ± 0.019; p < 0.05; Fig. 2b). The distance for the limited condition was significantly higher than that for the digested condition (p < 0.05) and was not significant in bulk soils (p = 0.3; Fig. 2b). In the UpSet plot, the unique ASVs in the limited, full, digested, and bulk soils were 160, 64, 98, and 120, respectively (Fig. 3). A total of 88 ASVs were common to all treatments (Fig. 3). The community structures of the ASVs were significantly clustered among the treatments (PERMANOVA, p < 0.05; Fig. 4a). In pairwise clustering comparisons, the differences among mycorrhizal developments were not significant; however, the differences between each mycorrhizal development and the bulk soil were significant (pairwise PERMANOVA, p < 0.01; Fig. 4a, Table S2). Furthermore, the community structures at different mycorrhizal development stages derived from the same P. japonica individuals showed a tendency to be similar (Fig. 4b). Quantification of assemblage patterns of bacterial communities indicated different community generation processes during mycorrhizal developments and bulk soil (Fig. 5). The communities in the limited condition and bulk soil were shaped by both deterministic processes (7.1% for homogeneous selection and 32.1% for heterogeneous selection in the limited condition; 50% for heterogeneous selection in bulk soil) and stochastic ways (28.6% for drift and 32.1% for dispersal limitation in the limited condition; 25.0% each for both drift and dispersal limitation in bulk soil; Fig. 5). However, stochastic processes dominated the formation of the communities in the full and digested conditions (42.9% for drift and 57.1% for dispersal limitation in both the conditions), without the selection of deterministic processes (Fig. 5). Indicator taxa analysis was performed for all treatments and was identified 3 ASVs for the limited condition, 7 ASVs for the full condition, 11 ASVs for the digested condition, and 58 ASVs for bulk soil (Table S3). Neighbor-joining phylogenetic trees were constructed for the top three dominant bacterial orders in the total number of ASVs detected from all three mycorrhizal developments, including Actinomycetales, Acidobacteriales, and Rhizobiales (Figs. S3-S5). In the Actinomycetales tree, five clades contain ASVs of both mycorrhizal developments and reference sequences derived from roots with ≥ 99% bootstrap values. Specifically, clade 1 included ASVs from all mycorrhizal developments, clade 2 was full and digested conditions, clade 3 was a full condition, clade 4 was a digested condition, and clade 5 was a limited condition (Fig. S3). In clade 3, two ASVs detected only in the full condition were assigned to two reference sequences (OP720531 and MG199992) of Streptomyces sp. derived from the roots. In the Acidobacteriales tree, eight clades were found to contain ASVs of both mycorrhizal developments and root-derived reference sequences with ≥ 99% bootstrap values. Clades 1 and 2 were composed of all mycorrhizal developments; clade 3 was limited and digested; clade 4 was digested; clades 5, 7, and 8 were limited, and clade 6 was full and digested (Fig. S4). No clade comprising ASVs derived from the full condition was identified. In the Rhizobiales tree, six clades were clustered with both our ASV sequences and root or nodule-derived reference sequences with ≥ 99% bootstrap values. Clades 1 and 3 involved all mycorrhizal development; clades 2, 5, and 6 were composed of a full condition, and clade 4 was a digested condition (Fig. S5). Clades 2, 5, and 6 consisted of ASVs only from the full condition with reference sequences of uncultured Alphaproteobacteria (OY971958), Rhizobium spp. (MK638406 and MW965275), and Mesorhizobium spp. (e.g., OR754792).

Fig. 2
figure 2

Comparisons among different mycorrhizal development stages of Pyrola japonica and bulk soil, based on a the number of observed amplicon sequence variants (ASVs) as species richness metric and b UniFrac distances as dissimilarity index. Data were mean ± standard error (n = 8 for all treatments). Color bars indicate mycorrhizal development stages of limited (aqua blue), full (yellow), digested (brown) conditions, and bulk soil (dark gray). Letters above bars represent significant differences among mycorrhizal development in each plot, as determined by a general linear model and Tukey-HSD tests (p < 0.05)

Full size image
Fig. 3
figure 3

An UpSet plot displaying numbers of shared and unique amplicon sequence variants (ASVs) among mycorrhizal development stages of Pyrola japonica and surrounding bulk soils according to an intersection matrix. Connected dots in the bottom panel represent intersections of shared ASVs. The vertical bars in the upper panel showed intersection sizes among bacterial communities in the different mycorrhizal developments and the bulk soils. Numbers above the bars represent the counts of unique or shared ASVs in the mycorrhizal development stages or bulk soils marked by dots

Full size image
Fig. 4
figure 4

Similarities of bacterial community composition based on UniFrac distances detected from Pyrola japonica and surrounding bulk soil displayed by a a non-metric multidimensional scaling (NMDS) plot and b a hierarchical clustering plot. a Markers in the NMDS plot indicate limited (aqua blue circles), full (yellow triangles), and digested (brown diamonds) conditions of mycorrhizal development stages, and bulk soils (dark gray squares). Stress value = 0.15. Polygons showed clusters among the treatments. b The clustering plot is generated using the group averaging method. Each color corresponds to different treatments as indicated in (a). Sample names are labeled with plant ID and the treatments

Full size image
Fig. 5
figure 5

Ratio of ecological assemblage processes regulating different communities detected from mycorrhizal development of Pyrola japonica and surrounding bulk soil, examined by null model analysis

Full size image

Discussion

This study demonstrated progressive changes in the rhizosphere bacterial community structures of different mycorrhizal developments in P. japonica root systems. To the best of our knowledge, this is the first study to clarify the bacterial community structure associated with arbutoid mycorrhizas. Moreover, progressive changes in the community along with mycorrhizal development in situ have not been reported for any mycorrhizal associations. Here, we discuss the tripartite interactions among plants, mycorrhizal fungi, and rhizosphere bacteria. We also highlighted the potential of a novel nutrient acquisition pathway via bacteria, which has not been well considered in previous mixotrophic studies.

Diversity of Bacterial Community Changes Along with Mycorrhizal Development

We demonstrated that rhizosphere bacterial diversity varied among different levels of mycorrhizal development. The number of ASVs and UniFrac distances in the full condition were significantly lower than those in the other conditions and bulk soil. In full conditions, mycorrhizal fungal hyphae are filled within epidermal cells, indicating active nutrient movement considering the mixotrophic nature of P. japonica [16, 42]. In contrast, other mycorrhizal developments, i.e. limited and digested, are unlikely to function because mycorrhizal fungi have not yet colonized under limited conditions, or hyphal coils have disintegrated in the digested condition [16]. Because the predominance of Russulaceae ECM fungi was found in the mycorrhizal fungal communities of P. japonica [16,17,18,19, 22], they may be abundantly colonized under full conditions. If this is the case, extramatrical mycelia can extend from epidermal cells into the soil, forming hyphal networks for carbon acquisition [17]. Russulaceae is a common ECM fungal taxon that is ubiquitously distributed in mature forests [20, 21], and mycorrhizal associations in this family are closely related to MHB belonging to Streptomycetaceae, Bacillaceae, Rhizobiaceae, and Burkholderiaceae [8, 43]. Previous MHB studies on woody ECM plants have indicated that the presence or absence of fungal mantles can have a significant impact on rhizosphere bacterial communities [11, 12]. For P. japonica mycorrhizas, no fungal mantles but fungal coils within epidermal root cells were formed, and growth of extramatrical mycelium was observed [16, 17]. Thus, rhizosphere bacteria in P. japonica mycorrhizal associations are likely to be involved in nutrient exchange between plants and mycorrhizal fungi. Lower α- and β-diversities in the full condition may be due to convergence into certain bacterial taxa in the associations. In contrast, higher α- and β-diversities in the limited and digested condition may be an admix effect of both opportunistic associated and MHB-like selected bacteria. In the case of the limited condition, the root tip of P. japonica would have encountered soil bacteria by chance due to root growth. The intermediate level of β-diversity in the digested condition can be ascribed to the recruitment of new bacterial taxa as aging roots of full condition. Detached or decayed root tissues can provide the rhizosphere with root rhizodeposition, attracting neighboring saprophytic bacteria [44]. Therefore, the full condition may have bacterial community assemblies due to mycorrhizal colonization compared with the limited and digested conditions.

The Assembly Pattern of Bacterial Communities Along with Progressive Mycorrhizal Developments

Contrasting assembly patterns of bacterial community structures were found under limited and full conditions, whereas community compositions were consistent regardless of mycorrhizal development stages. The number of unique ASVs was significantly different, being highest in the limited condition and lowest in the full condition (Fig. 3). Physically, continuum habitat transitions from bulk soil to the rhizosphere allow for the sharing of common bacteria [2]. This was evidenced by 11.3% of all ASVs (88/776 ASVs) being shared among all treatments (Fig. 3). Thus, the unique ASVs in each condition were considered to be the remaining shared ASVs and minor members of the community with a low detection probability. The high number of unique ASVs under limited conditions suggests non-selective bacterial gathering owing to opportunistic associations. In contrast, under full conditions, a low number of unique ASVs may indicate the need for elaborate screening systems for bacterial requests during the development of mycorrhizal associations. The clustering pattern of the bacterial communities, as in the NMDS plot, suggested significant discrimination among the treatments (Fig. 4a). However, community compositions among mycorrhizal developments were highly similar, based on pairwise treatment comparisons (bulk soil vs. limited condition, p < 0.01; bulk soil vs. full condition, p < 0.01; bulk soil vs. digested condition, p < 0.01; Table S3). Thus, clustering was significantly different between the bulk soil and pooled data for all mycorrhizal developments. In this respect, community composition was stable regardless of mycorrhizal development. A nested analysis, incidentally, slightly suggested the community structure under the full condition is likely a subset of other conditions based on. However, owing to possible stochastic PCR and sequencing biases in the NGS short-reading approach, no statistical significance was found (NODF = 8; data not shown). The bacterial community composition remained highly similar irrespective of mycorrhizal development, but bacterial diversity decreased only under full conditions. This phenomenon can be explained as follows: First, the growing root tip encounters soil bacteria in a non-selective manner, leading to diversification under limited conditions. Second, under full conditions, a bacterial community that was not affected by mycorrhizal fungi was maintained. Finally, in the digested condition, certain members of the bacterial community from the previous full condition were almost completely retained, and additional bacteria that fed on rhizodepositions were recruited from the outside rhizosphere.

The explanation was also supported by the quantification of ecological processes in community assembly that showed different patterns between the limited condition and post-full conditions (Fig. 5). Notably, the bacteria communities in the limited condition and bulk soil were partially regulated by deterministic processes, although the communities in the full and digested conditions were shaped only by stochastic ways (Fig. 5). A similarity of the processes between the limited condition and bulk soil is additional evidence that the bacterial community in the limited condition is likely to be established contingently via surrounding soil bacteria (Figs. 2, 3, and 5). Only the community in the limited condition was partially regulated by homogeneous selection, where environment is uniform and community variation is less obvious (Fig. 5). This result indicates that certain bacteria groups may be selected in the rhizosphere during the formation of root tips (Fig. 5). Additionally, the homogeneous selection process was not predicted in subsequent mycorrhizal developmental stages, indicating that the phylogenetic selection of rhizosphere bacteria may only occur at an early stage of root growth (Fig. 5). Although both deterministic and stochastic ecological processes involve the regulation of bacterial communities with varying degree of contributions, the communities are more likely to be governed by the former process under extreme environmental situations (e.g., high pH) [45]. The rhizosphere environment in the limited condition and bulk soil may be less stable than those in the full and digested conditions owing to lacks of mitigation by the existence of plants or mycorrhizal fungi [1, 2]. Furthermore, dispersal limitation accounted for ca. 60% in the community assembly even after mycorrhizal fungi have degenerated (Fig. 5), inferring that bacterial communities were maintained largely unchanged in the rhizosphere. High similarities in bacterial communities among mycorrhizal developments within the same individuals—plants 1, 5, and 8—also supported the importance of historical contingency in rhizosphere community assembly (Fig. 4b). The historical impacts were conceptualized as “priority effect” in community ecology that occurs when the order and timing of species arrivals influence community assembly [13, 46]. In the case of P. japonica, roots remain functioning as conduits even after the mycorrhizal degeneration and are connected to newly formed mycorrhizal parts. Therefore, the historical continuity of rhizosphere bacterial communities may be a reasonable strategy for persistence in both nutrient acquisition and pathogenetic resistance. Surprisingly, the presence of mycorrhizal fungi did not determine bacterial communities as a biotic factor (Fig. 5); if mycorrhizal fungi attracted certain bacterial groups, homogeneous selection would have contributed to the community in the full and digested conditions. This result suggests that mycorrhizal fungi may function as an antagonist to pathogenic bacteria, rather than selecting bacteria to form a preferred bacterial community. Thus, some opportunistic bacteria can remain in the rhizosphere irrespective of mycorrhizal developments. This may result in diffused phylogenies and highly similar communities of bacteria in the full condition. These results indicate ecological importance that mycorrhizal fungi work not only as temporal drivers initiating the formation rhizosphere bacterial community but also as key founders exerting continuous influences to establish priority effects.

Potential Functions Characterized in Rhizosphere Bacterial Communities After Mycorrhizal Formation

In the phylogenetic trees of the top three taxa, four clades were composed of ASVs derived only from the full condition with reference sequences, such as Streptomyces spp. (ASV651 and ASV654, Fig. S3), uncultured Alphaproteobacteria bacterium (ASV378, Fig. S4), Rhizobium spp. (ASV1553, Fig. S5), and Mesorhizobium spp. (ASV295, Fig. S5). From our observations, epidermal cells in the full condition were harbored by healthy hyphal coils (Fig. 1b), indicating that the ASVs in the clades are likely to function in some roles with P. japonica and its associated mycorrhizal fungi. The reference sequences in three of the four clades were retained as MHB isolates, such as actinomycetes and rhizobia [8, 47]. Some mycorrhizal fungi in the genus Lactarius (Russulaceae) are symbionts of P. japonica and show co-occurrence relationships with Streptomycetaceae belonging to Actinomycetales [17, 22, 43]; therefore, it is not unexpected to list them as symbiotic candidates. For example, the candidate MHB Streptomyces sp. AcH 505, isolated from the ECM roots of Amanita muscaria, inhibits the growth of other fungi [48, 49]. Other ASVs were clustered with either Rhizobium spp. or Mesorhizobium spp., which belong to the Rhizobiaceae, known as nitrogen-fixing bacteria [50, 51]. However, they are not necessarily obligate symbionts and sometimes occur not only as free-living in the soil but also as MHB [52]. Considering the tripartite interactions among host plants, ECM fungi, and nitrogen-fixing bacteria [53], the two ASVs (ASV295 and ASV1553) detected in this study may have contributed to nitrogen acquisition by P. japonica.

Furthermore, the detection of indicator taxa revealed the potential contribution of bacterial community to host plant and mycorrhizal fungi (Table S3). Notably, four indicator ASVs selected from mycorrhizal developments belonged to Rhizobiales (ASV032 from the full condition; ASV012, ASV034, and ASV003 from the digested condition; Table S3). Among them, two ASVs were clustered with reference sequences derived from ECM mycorrhizosphere or nodules with bootstrap values ≥ 99% (Fig. S5). Moreover, these ASVs were detected in all mycorrhizal conditions (Fig. S5), suggesting their continuous association with P. japonica and mycorrhizal fungi. Combined with our understanding of the community analysis, rhizosphere bacterial communities associated with the plant may have been historically involved in nitrogen cycling.

Considering the mixotrophic nature of P. japonica, the nutrient acquisition process may function complementarily with rhizosphere microbes, with mycorrhizal fungi obtaining carbon and MHB to obtain nitrogen. Furthermore, MHB may indirectly contribute to carbon acquisition by assisting in mycorrhizal formation. Rhizosphere bacteria are potential new pathways for nutrient acquisition in mycoheterotrophic plants [54, 55]. The high nitrogen content found in mycoheterotrophic and mixotrophic plants has long been of interest [56]. Further understanding of mycoheterotrophic lifestyles would require clarification of the contribution of MHB. However, we acknowledge that the DNA approach with short-read sequencing provides a lower resolution for bacterial phylogenetic diversity. Functional evaluations using isolated bacterial strains in vitro are essential to examine the interactions between P. japonica and MHB.

As a caveat in the study, we focused on the mycorrhizal development as a factor regulating bacterial communities, while we did not fully consider plant-derived influences. The results of this study supported that bacterial shifts were caused by the mycorrhizal development but may be due to metabolic changes in root organs themselves along with plant growth. Since the only visual assessment of mycorrhizal developments was conducted in this study (Fig. S1), metabolic-functional aspects of how rhizosphere bacterial communities are regulated by plants and mycorrhizal fungi along with mycorrhizal development were remained to be unanswered.

Conclusion

The present study is the first to clarify the assemblage pattern of the rhizosphere bacterial communities of arbutoid mycorrhizal P. japonica, along with different mycorrhizal developments from the limited, full, and digested conditions of the associated fungi. As results, mycorrhizal developmental levels affected irreversibly of rhizosphere bacterial communities. Under the full condition, where mycorrhizal associations are likely to be active, both α- and β-diversity of the bacterial community were the lowest compared to other developments, and the community was historically maintained after the mycorrhizal fungal degeneration in the digested condition. Furthermore, some bacterial taxa known as mycorrhizal helper bacteria were characterized from the communities in the full and digested conditions. This study highlights the dynamic change of rhizosphere bacterial communities synchronized with the continuum of mycorrhizal developmental conditions within a single root system.