A total of 66 soil rhizosphere and bark samples from 11 representative regions (Zunyi, Panzhou, Wangcang, Heishui, Cili, Liuyang, Shuanglong, Shiyan, Linbao, Qimen, and Longshi) of E. ulmoides in China were used in the study (Supplementary Figure 1). Three healthy E. ulmoides, approximately 10 m from one another, were selected in each region, and bark samples from 1.5 m above the ground of each tree were collected (Dong et al., 2020). Sub-samples for chemical analysis were transported to the laboratory, dried, and crushed, followed by the measurement of pinoresinol diglucoside, aucubin, and geniposidic acid by tandem liquid chromatography quadrupole time of flight mass spectrometry (LC-QTOF MS/MS). Sub-samples for DNA extraction were frozen in dry ice during transportation and stored at −80°C until further processing. Subsequent chemical analysis and DNA extraction were performed as described previously Dong et al. (2020).

Similarly, three healthy E. ulmoides trees were chosen from each region, and rhizosphere soil sample was collected from each tree. For each tree, rhizosphere soil samples from three vertical directions (approximately 1 m away from the trunk) at a depth of 20 cm were collected, and were uniformly mixed into a sample (Xu et al., 2018). Soil samples were used to extract the total genomic DNA and measure physiochemical factors.

Total genomic DNA of E. ulmoides rhizosphere soil samples was extracted according to FastDNA Spin Kit for Soil (MP, United States). In this study, DNA concentration between 40 and 60 ng/μL was considered qualified. According to universal primer pair 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′), the 16S rRNA gene was amplified. The conditions and mixtures for PCR amplification referred to our previous method (Dong et al., 2021b). PCR products were sequenced by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China) based on an Illumina MiSeq PE300 platform (Illumina, San Diego, United States). Most of these analysis were performed on BMK Cloud.¹ The raw data was uploaded to the NCBI database (BioProject ID: PRJNA771480; BioSample: SAMN22311478).

The bark sampling method was according to our previous study (Dong et al., 2021b), and the current work is a follow-up study to explore the construction mechanism of the rhizosphere soil and bark microbiota and the function of the core microbiota. The bark sample was corresponding to the rhizosphere soil sample. The diversity of bark bacterial communities was characterized based on universal primer pair 338F and 806R. The sequencing data of E. ulmoides bark (BioProject ID: PRJNA771480; BioSample: SAMN22311479) was collated and compiled based on our previous study (Dong et al., 2021b).

The approach of genotype characteristics, active compounds, and soil physicochemical properties were described in our previous study (Dong et al., 2021b). The genetic characteristics of E. ulmoides were analyzed based on the chloroplast gene trnH-psbA (forward: 5′-CGCGCATGGTGGATTCACAATCC-3′; Reverse: 5′-GTTATGCATGAACGTAATGCTC-3′) and trnL-F (forward: 5′-ATTTGAACTGGTGACACGAG-3′; Reverse: 5′-CGAAATCGGTAGACGCTACG-3′), and these analyses were performed by the software on BioEdit v7.0.9.0, MAFFT v 7.037, MEGA6, SequenceMatrix1.7.8, and DnaSP version 6. The content of active compounds was determined based on LC-QTOF MS/MS. The determination of soil physicochemical properties mainly refers to “Physical and Chemical Analysis of Soil Properties.” In this study, the data of host genotypes, active compounds, and soil physicochemical properties were placed in Supplementary Tables 1, 2.

Core microbiota was defined with reference to the method of Dong et al. (2021a). The operational taxonomic units (OTUs) with all samples were screened out to construct a co-occurrence network. Highly significant (r > 0.6 or r < −0.6; p < 0.05) spearman correlations between bacterial co-occurrence network were visualized by Gephi, version 0.9.2. Similarly, network topology parameters (e.g., degree, betweenness centrality, closeness centrality, hub, clustering, and modularity) were calculated in the software Gephi, version 0.9.2. In the study, the threshold for core microbiota identification was degree ≥ 20, and closeness centrality ≥ 0.29. In addition, OTUs with relative abundance above 0.1% were placed within the dominant taxa (Xiong et al., 2021), while rare taxa refer to OTUs with relative abundance below 0.01% (Pascoal et al., 2021).

The habitat, energy source, and metabolic functions of the core microbiota were annotated based on METAGENassist² (Arndt et al., 2012). The phenotypic information of 11,000 microorganisms and the full-sequence of approximately 1,800 microorganisms were included in METAGENassist. Moreover, the phenotypic information contains approximately 20 categories for each microorganism, such as metabolism, habitat, energy source(s), oxygen requirements, preferred temperature range.

The weighted UniFrac distance matrix was calculated to evaluate beta-diversity (β-diversity), and then non-metric multidimensional scaling (NMDS) was used for ranking analysis (Xiong et al., 2021). The significance of different factors on community dissimilarity was tested with permutational multivariate analysis of variance (PERMANOVA) or nested PERMANOVA using the “adonis” function of the VEGAN package in R 3.6.3 (Oksanen et al., 2015). Null model analysis was performed “picante” package (Zhao et al., 2021), Structural equation model was analyzed the “lavaan” package in R 3.6.3 (Rosseel, 2012). Traceability analysis refers to the fast expectation-maximization for microbial source tracking (FEAST), which is a ready-to-use scalable framework that can simultaneously estimate the contribution of thousands of potential source environments in a timely manner, thereby helping unravel the origins of complex microbial communities³ (Shenhav et al., 2019).

Non-metric multidimensional scaling and PERMANOVA showed that variation of bacterial communities in the rhizosphere soil (R² = 78.84%, p = 0.001) and bark (R² = 67.14%, p = 0.001) of E. ulmoides was predominantly determined by environmental conditions (Figures 1A,B and Table 1). The compartment niche did not have a significant effect on E. ulmoides-associated bacterial communities (R² = 1.73%, p = 0.28) (Table 1), but the rhizosphere soil and bark samples clustered clearly and independently (Figure 1C and Supplementary Figure 2). Similarly, different genotypes had no significant effects on the rhizosphere soil and bark bacteria of E. ulmoides (Figures 1D–F and Table 1).

According to the unweighted β-NTI statistics, the β-NTI values of the rhizosphere soil and bark bacterial communities of E. ulmoides mainly ranged from −1 to 9 (Figures 2A,B), suggesting that the assembly of bacterial community was affected by both stochastic and deterministic factors (β-NTI value > 2, then the deterministic process was dominant; β-NTI value < 2, then the stochastic process was dominant). SEM further suggested that the bacterial α-diversity in the rhizosphere soil and bark of E. ulmoides was affected to varying degrees by host active components, soil physicochemical properties, and climatic factors (Figures 2C,D). For example, pH (with an estimated value of 0.56, p < 0.001), aucubin (0.31, p < 0.001), and relative humidity (0.30, p < 0.001) had significant positive effects on the Shannon diversity index, while geniposide acid had a significant negative effect on the Shannon diversity index (−0.47, p < 0.001). Altitude (0.61, p < 0.001) and relative humidity (0.46, p < 0.001) had significant positive effects on the Chao index, while rainfall had a significant negative effect on the Chao index (Figure 2C).

The bark and rhizosphere soil microbiota respond similarly to environmental factors. pH (0.57, p < 0.001) and aucubin (0.38, p < 0.001) had significant positive effects on the Shannon diversity index; geniposide acid (−0.49, p < 0.001) and total phosphorus (−0.34, p < 0.001) had significant negative effects on the Shannon diversity index; and pH (−0.72, p < 0.001), altitude (−0.43, p < 0.001), and relative humidity (−0.31, p < 0.001) had significant negative effects on the Chao index (Figure 2D). Overall, the composition and structure of the E. ulmoides microbiota were driven by multiple environmental factors.

Although compartment niches were not the main factor affecting the microbiota of E. ulmoides, there were still similarities and differences in the composition of the microbial communities of rhizosphere soil and bark. At the genus level, 690 taxa were shared between rhizosphere soil and bark microbial communities, among which the top five in relative abundance were Enterobacter, an unidentified genus of Acidobacteria, an unidentified genus of Muribaculaceae, Sphingomonas, and Lactobacillus. Bark contained five unique taxa, which were all unidentified genera of Verrucomicrobiae, Actinophytocola, Planctomycetacia, Myxococcales, and Armatimonadeles. In contrast, rhizosphere soil contained 43 unique taxa, among which the top five relative abundances comprised unidentified genera of Bacteroidia, Clostridiales, and Saccharibacillus, Spirochaeta, and Erysipelotrichaceae (Figure 3A).

Fast expectation-maximization microbial source tracking (FEAST) suggested that 54.92% of the members of the bark bacterial community were derived from rhizosphere soil, and 45.08% of the members were derived from other unknown environments (Figure 3B).

Here, we identified dominant taxa, core taxa and rare taxa (Figures 4A,B and Supplementary Figures 3, 4). A total of 201 OTUs (0.1% of total OTU number) were defined as the dominant microbiota, with Enterobacter, Actinobacteria, and Sphingomonas containing the most OTUs. The dominant microbiota accounted for 78.89% of the soil samples, and 70.61% of the bark samples (Supplementary Figure 4a). A total of 185 OTUs (0.01% of total OTU number) were defined as rare microbiota, of which the Proteobacteria Deltaproteobacteria, Burkholderiaceae, and Rhizobiales contained the most OTUs (Supplementary Figure 4b). Based on the network parameters, 18 OTUs were defined as core microbiota (Supplementary Figure 3), and these belonged to the genera Escherichia, Enterococcus, Methylobacterium, and Ramlibacter. The proportion of core microbiota was higher in bark samples and lower in soil samples compared to the dominant microbiota (Figures 4A,B).

To further characterize the impact of the core microbiota on other members of the host microbiota, potential interrelationships between the core microbiota and the dominant and rare microbiota were evaluated (Figures 4C,D). There was no significant difference in the complexity of the two microbial networks (average degree representing network complexity). Within co-occurrence networks, we found that correlations were usually positive (core-dominant taxa, 97.35%, n = 1057; core-rare, 99.1%, n = 783) (Figures 4C,D). Core microbiota was a highly interactive “hub” microbes connecting dominant and rare taxa, play an important role in determining microbiome structures.

Phenotypic and functional characteristics of the core microbiota of E. ulmoides were deciphered based on the METAGENassist platform and statistical analysis (Figure 5). The core microbiota with multiple habitat characteristics mostly exists in the soil rhizosphere, while the core microbiota related to the host predominantly exist in the bark of E. ulmoides. Based on energy source utilization, the core microbiota present in the rhizosphere soil were mainly methylotrophs, while the core microbiota existing in the bark were mostly autotrophs (Figure 5).

Metabolic functions of the core microbiota of E. ulmoides predominantly comprised ammonia oxidizer, dehalogenation, chitin degradation, sulfate reducer, and sulfide oxidizer. In addition, some metabolic functions closely related to the biological properties of E. ulmoides were worthy of attention, such as xylan degrader, nitrogen fixation, and stores polyhydroxybutyrate (Figure 5C). The correlation heatmap showed that Escherichia, Eubacterium, and an unidentified genus of the family Ruminococcaceae were significantly correlated with pinoresinol diglucoside, Lactococcus exhibited a significant correlation with aucubin, and Enterococcus and Lactococcus were significantly correlated with geniposidic acid (Figure 5).

Assembly of plant-associated microbiota is regulated by the host’s metabolism, genotype, soil physicochemical properties, climate, and niche (Wagner et al., 2016; Cregger et al., 2018; Sun et al., 2020; Xiong et al., 2021). In the current study, assembly of the rhizosphere soil and bark microbiota was predominantly shaped by environmental factors and not by genotype or compartment niche (Table 1). Environmental factors are generally accepted to be closely related to the composition of plant-associated microbiota (Breidenbach et al., 2015; Trivedi et al., 2020; Zhang et al., 2020; Bona et al., 2021). Relative humidity was significantly correlated with the diversity index of the soil and bark bacterial communities, which was consistent with previous studies (Dong et al., 2021c). As reported by Cruz-Paredes et al. (2021), environmental humidity within a certain range significantly impacts the respiration and growth of microorganisms. Moreover, pH could also significantly affect assembly of microbial communities in soil and bark (Figure 2). This may be because pH combines multiple soil chemical variables and is the result of many chemical reactions between organic (biological and abiotic) and inorganic soil components (Lopes et al., 2021). The bark microbiota was not directly in contact with soil pH but was still significantly affected by pH in the current study (Figure 2B). It was previously reported that the indirect impact of soil pH on the microbiota in other compartments of the plant is because microbes in the other compartments were derived from rhizosphere soil (Fitzpatrick et al., 2018; Lopes et al., 2021). Similarly, FEAST analysis suggested that 54.92% of the members of the bark bacterial community were derived from rhizosphere soil (Figure 3B).

Venn diagram analysis showed that 690 genera, such as Enterobacter, Sphingomonas, and Lactobacillus (Figure 3A), were shared between soil and bark communities. These genera were also found in previous studies on E. ulmoides (Dong et al., 2020, 2021a) and are beneficial microorganisms with potential ecological functions for medicinal plants (Quijada et al., 2018; Asaf et al., 2020; Webster et al., 2020). Highly abundant unique bacteria of the rhizosphere soil samples included Bacteroidia, Clostridiales, Saccharibacillus, Spirochaeta, and Erysipelotrichaceae, which have specific selective effects on the soil, predominantly exist in the soil, and play an important role in geo-biological-chemical circulation (Wilms et al., 2006; Namirimu et al., 2019). Most of the bacteria unique to the bark are those with relatively older evolutionary status, such as the Planctomycetes, which are important for the nitrogen cycle (Cavalier-Smith, 2002). Myxococcus, which is only detected in the bark of rhododendrons, has broad application potential in agriculture, biomedicine, and environmental protection due to its ability to produce a variety of natural products (Wang et al., 2021).

Core microbiota, an important subset of the host-associated microbiota, can not only promote the growth of host plants but also recruit functional species and block pathogens/pests (Banerjee et al., 2018; Hassani et al., 2018; Toju et al., 2018). In recent years, the core microbiota of medicinal plants has been accepted as a processing plant or reservoir for secondary metabolites (Koberl et al., 2013; Chen et al., 2018; Xu et al., 2021). Using microbial network analysis, 18 core bacteria of E. ulmoides belonging to the genera Escherichia, Enterococcus, Methylobacterium, Eubacterium, and Lactococcus were identified (Figure 4A). Among them, species of Escherichia, Eubacterium, and Ruminococcaceae were significantly related to pinoresinol diglucoside; a species of Lactococcus had a significant correlation with aucubin; and species of Enterococcus and Lactococcus were significantly correlated with geniposide acid (Figure 4D). Previous studies have shown that Escherichia sp. (Leonard et al., 2008), Enterococcus sp. (Müller et al., 2001), Eubacterium sp. (Newman and Cragg, 2020), and Lactococcus sp. (Golomb and Marco, 2015), Ruminococcaceae sp. (Yin et al., 2020) could promote plant growth or accumulate secondary metabolites. The phenotypic information indicated that the core microbiota mainly contained methylotrophic bacteria in the rhizosphere soil and mostly autotrophic bacteria in the bark (Figure 5B). This indicated that the core microbiota from the rhizosphere soil may play the role of consumers, while the core microbiota from the bark may predominantly play producer roles (Thakur and Geisen, 2019; El-Sayed et al., 2021). The metabolome of the 18 core bacteria was further analyzed, and in addition to some common metabolic functions of the plant microbiota (Xu et al., 2018; Trivedi et al., 2019), the functions of xylan degradation, nitrogen fixation, and polyhydroxybutyrate storage were worthy of additional attention. Xylan-degrading and nitrogen-fixing bacteria were reported to potentially exhibit direct or indirect effects on plant growth and metabolism (Gopalakrishnan et al., 2017; Jung et al., 2018). Polyhydroxybutyrate was a type of polyhydroxyalkanoate (PHA), a polyester polymer (Sirohi et al., 2021) that has similar structure and function to E. ulmoides rubber, such as strong plasticity (He et al., 2014; Sirohi et al., 2021). Therefore, the core taxa with polyhydroxybutyrate metabolism functions may be involved in the synthesis of E. ulmoides rubber. In general, these findings reveal that the core microbiota may participate in soil organic matter decomposition and metabolite synthesis in E. ulmoides.

Network analysis can be used to identify the core or hub microbiota at the center of the microbial network (Banerjee et al., 2018; Dong et al., 2021a). Core microbiota is composed by important taxonomic units that dominate host microbial community construction, influencing community structure through strong interactions with the host or other microbes rather than simply being present in high abundance (Banerjee et al., 2018; Hamonts et al., 2018; Trivedi et al., 2020; Dong et al., 2021c). A comprehensive investigation of the Arabidopsis phyllosphere microbiota found that the core microorganisms Albugo sp. and Dioszegia sp. regulate the entire microbial community by inhibiting the growth and diversity of other microbes (Agler et al., 2016). The current study found that the core microbiota was highly connected with the dominant microbiota, forming closely related microbial clusters (Figure 4C). Dominant microbiota and core microbiota have vital ecological roles in microbiota assembly and ecosystem functions (Banerjee et al., 2018). Core microbiota can selectively adjust and assist these dominant microbiota (Banerjee et al., 2018), such as by forming a highly specific microbiota with the dominant microbiota for biological nitrogen fixation (Fierer, 2017). In contrast with the dominant microbiota, rare microbiota was often overlooked. Recent studies confirmed that rare microbiota was instrumental in the construction, function, and stability of microbiota (Jousset et al., 2017; Dong et al., 2021c; Pascoal et al., 2021). Similarly, significant interactions between core and rare microbiota of E. ulmoides were detected in the current study (Figure 4D). This is congruent with previous findings that rare microbiota often interacts with other species although they occupy narrow ecological niches (Fierer et al., 2007). Together these data suggest the core microbiota not only has important functions but may also be a keystone species that determines the network structure of the E. ulmoides microbial community.