The root and soil of healthy P. amurense were collected from Heilongjiang Province, China (45.7662°N, 126.6247°E). Five roots (at a depth of 10–30 cm) from one P. amurense were cut off and collected as one sample. In addition, bulk soil (at a depth of 10–20 cm) was collected at a distance of 100–200 cm away from the roots. After removing loose soil from the roots, only 2 mm of rhizosphere soil remained and was collected (Edwards et al., 2015; Hilton et al., 2021). These samples were labeled as R (root of wild P. amurense), BS (bulk soil), and RS (rhizosphere soil) with five replicates, respectively. All samples were stored at −80°C. Bulk and rhizosphere soil samples were then sieved by using a 2-mm sieve. The roots of P. amurense were carefully washed and sterilized by undergoing the following operations: 70% (v/v) ethanol for 3 min, 2.5% (v/v) sodium hypochlorite (NaClO) for 5 min, and sterile water four times (Wei et al., 2021). Half of each sample was used for the analysis of microbial characteristics, whereas the other half was employed for metabolite detection. Voucher specimens labeled as PAR20221-PAR20225, PABS20221-PABS20225, and PARS20221-PARS20225 were then deposited at the College of Life Sciences, Northeast Forestry University, China.

The collected root samples were dried at 60°C for 1 week, and, then, the dried materials were powered and sieved through a No. 60 mesh. An accurately weighed 0.2 g of sample of the powered materials was extracted ultrasonically with 70% aqueous methanol (2 mL) under 50°C for 1 h, with three replications. The extract was further centrifuged at 12,000 rpm for 10 min, and 1 mL of supernatant was filtered using a 0.22-µm nylon filter. As for the soil sample, 0.2 g of bulk and rhizosphere soil was extracted ultrasonically with 70% aqueous methanol (2 mL) under 50°C for 1 h, with three replications. The extract was centrifuged at 12,000 rpm for 10 min. Supernatant (1 mL) was evaporated (40°C) to dryness in vacuo, and the residue was redissolved in 0.2 mL of 70% methanol. Then, this solution was filtered using a 0.22 nylon filter.

After optimizing the experimental conditions, a method was developed using ultrahigh-performance liquid chromatography (UPLC) coupled with quadrupole time-of-flight mass spectrometry (MS) for the nontarget metabolomic profiling of P. amurense and its associated soil. Seven reference standards of phellodendrine, magnoflorine, magnocurarine, jatrorrhizine, tetrahydropalmatine, palmatine, and berberine were dissolved and then diluted with methyl alcohol to prepare a series of standard solutions of different concentrations. To ensure the robustness of the analytical method, the quality control samples were inserted into the sequence every five samples. The samples were injected into a Kinetex C18 column (2.6-μm particle size, 4.6 mm × 150 mm; Phenomenex, USA) maintained at 40°C. The gradient eluent consisted of mobile phase A (0.1% formic acid in water) and mobile phase B (acetonitrile). The gradient elution program was as follows: 0–1 min, 10% B; 1–10 min, 10% to 95% B; 10–12.3 min, 95% B; 12.3–13 min, 95% to 10% B; and 13–15 min 10% B. The flow rate used for separation was 0.5 mL/min, with an injection volume of 1 μL. The MS was equipped with an electrosprayionization (ESI) ion source and operated in the positive ion mode, in the full scan range from mass-to-charge ratio (m/z) of 100 to 2,000 with fragmentation of the five most intensive signals. The following MS conditions were used: ion spray voltage, 5,500 V; curtain gas, 25 psi; ion source gas 1, 40 psi; ion source gas 2, 60 psi; and temperature, 550°C.

The mass spectrum data were analyzed by Analyst TF1.8 Software (https://sciex.com/products/software/analyst-software/). Using the metabolome software LibraryView (AB Sciex, Shanghai, China), the samples’ metabolites were analyzed qualitatively and quantitatively. The chemical structure of the compounds was characterized on the basis of the calculated accurate masses of the molecular ions, protonated molecules, fragment ions, and retention behavior. The observed constituent ion peaks [M]⁺/[M+H]⁺/[M+NH4]⁺ provide reliable information that can confirm the molecular weights and structures of the compounds. The relative contents (after taking log values) of compounds in different samples were displayed by R package pheatmap (v1.0.12).

Metagenomic DNA from root, rhizosphere soil, and bulk soil samples was extracted using the HiPure Soil DNA Kits (Magen, Guangzhou, China) according to the manufacturer’s recommendations. The DNA concentration and purity were accurately measured using the NanoDrop 2000 spectrophotometer prior to performing agarose (1% w/v) gel electrophoresis. The purified DNA was used to amplify either the fungal Internal Transcribed Spacer 2 (ITS2) region ITS3_KYO2 (5′-GATGAAGAACGYAGYRAA-3′) and ITS4 (5′-TCCTCCGCTTATTGATATGC-3′) (Guo et al., 2017) or the V3-V4 region of the bacterial 16S ribosomal RNA (rRNA) gene 341F (5′-CCT ACG GGNGGC WGC AG-3′) and 806R (5′-GGAC TAC HVGGG TAT CTA AT-3′) (Toju et al., 2012). The PCR-negative controls consisted of ddH₂O instead of DNA template. All samples were subjected to the following PCR conditions: initial denaturation at 95°C for 2 min, followed by 27 cycles of denaturation at 98°C for 10 s, annealing at 62°C for 30 s, extension at 68°C for 30 s, and a final extension at 68°C for 10 min. The paired-end approach was utilized for sequencing by using the Illumina Novaseq 6000 sequencing platform.

To obtain high-quality clean reads, raw reads were further filtered using FASTP (v0.18.0) (Chen S. F. et al., 2018). Paired-end clean reads were merged into raw tags using FLASH (v1.2.11) (Magoc and Salzberg, 2011) with a minimum overlap of 10 bp and mismatch error rates of 2%. The clean tags were clustered into operational taxonomic units (OTUs) of ≥ 97% similarity using UPARSE (v9.2.64) (Edgar, 2013), a high-accuracy and high-throughput OTU clustering pipeline. In each cluster, the tag sequence with the highest abundance was selected as the representative sequence.

The stacked bar plot of the community composition was visualized in the ggplot2 package (v2.2.1) (Wicham, 2011), and circular layout representations of species abundance at the genus level were created using circos (v0.69-3) (Krzywinski et al., 2009). To evaluate the diversity of species and community of the bacterial and fungal microorganisms in roots, rhizosphere soil, and bulk soil of P. amurense, alpha-diversity based on the diversity index and beta-diversity based on the weighted UniFrac distance matrix were conducted, respectively. Chao1, Abundance-based Coverage Estimator (ACE), Shannon, and Simpson index were calculated using QIIME (v1.9.1) (Caporaso et al., 2010), and PCoA (principal coordinate analysis) (Dixon, 2003) was generated using the Vegan package (v2.5.3) and ggplot2 package (v2.2.1) of R project. The biomarker features in bacterial and fungal microorganism communities of P. amurense were evaluated through Linear discriminant analysis Effect Size (LEfSe) (Segata et al., 2011).

WGCNA was performed to evaluate the correlation between the secondary metabolites in roots and their associated microbial communities. After obtaining the absolute contents of seven key BIAs in the bulk soil, rhizosphere soil, and roots of P. amurense, combined with the abundance of bacteria and fungi OTU in the bulk soil, rhizosphere soil, and roots of P. amurense, co-expression network modules were generated by WGCNA package in R. Eigengene value was calculated for each module to search the family-level microorganisms associated with BIA metabolism in P. amurense. The co-expression modules were obtained using automatic network construction function (blockwiseModules) with default parameters; apart from the soft threshold power of 18, networkType was unsigned, TOMtype was unsigned, mergeCutHeight was 0.1, and minModuleSize was 15. In addition, the Pearson correlation coefficient was also calculated. Microorganisms and metabolites with a correlation greater than 0.9 and P < 0.01 were selected, and the correlation network was drawn using Cytoscape (v3.9.1, https://cytoscape.org/).

To uncover the distribution characteristics of main active compounds of bulk, rhizosphere soil, and root of P. amurense, the targeted metabolomics was employed, which focuses on the analysis of specific categories of metabolites with more selectivity and sensitivity (Roberts et al., 2012). Obvious differences in the quantity and relative content of specialized metabolites were observed at the different compartments of P. amurense. A total of 31 metabolites were detected and identified (Table 1, Figure 1A). These metabolites include 28 alkaloids, two limonoids, and one quinic acid. Among the 28 alkaloids, 23 are BIAs, two are indole-type, two are quinoline-type, and one is another type of alkaloid. The 23 BIAs can be classified into three groups: three aporphine alkaloids, seven benzylisoquinoline alkaloids, and 13 protoberberine alkaloids (Supplementary Table S1, Supplementary Information). The number of BIAs is the greatest among the 28 alkaloids, accounting for 82% of all the BIA metabolites. In addition, 31 specialized compounds were characterized in the root of P. amurense, whereas 21 were shared in the rhizosphere soil. No specialized metabolites of P. amurense were detected in the bulk soil (Figure 1B). There was a significantly higher content of the 21 commonly specialized compounds in the roots compared with that in the rhizosphere soil.

To determine the major chemical compounds present in the rhizosphere soil and roots of P. amurense, we used UPLC-MS/MS to analyze the content of seven BIAs (phellodendrine, magnoflorine, magnocurarine, jatrorrhizine, tetrahydropalmatine, palmatine, and berberine), which have been shown to have well pharmacological activity. The highest contents of magnoflorine and berberine were found in the roots, whereas magnoflorine, magnocurarine, and phellodendrine were the main compounds in the rhizosphere soil (Figure 1C). The content of the seven BIAs in roots was significantly higher than that in rhizosphere soil, and these BIAs in roots and rhizosphere soil were, respectively, clustered in two groups, suggesting well sample repetition (Figure 1D).

To understand the differences in microbial composition among the bulk, rhizosphere soil, and root, statistical analyses were performed on different taxonomic levels. By using a cutoff of >97% sequence identity, high-quality reads were clustered into bacterial 15,950 OTUs, 16,471 OTUs, and 4,452 OTUs of bulk, rhizosphere soil, and root, respectively. Regarding fungi, high-quality reads were clustered into 3,359 OTUs, 3,526 OTUs, and 1,503 OTUs of bulk, rhizosphere soil, and root, respectively. Rarefaction curve analysis revealed that all samples were nearly parallel to the X-axis, suggesting that the obtained sequences adequately depict the overall bacterial and fungal diversity (Supplementary Figure S1, Supplementary Information). A taxonomic histogram showed that the composition of P. amurense root-associated microbial communities varied along with the compartments. At the bacterial phylum level, the microbial composition analysis revealed that Firmicutes, Acidobacteriota, and Proteobacteria exhibited the highest abundance in bulk soil, whereas Actinobacteriota, Proteobacteria, and Acidobacteriota were enriched in rhizosphere soil (Figure 2A). Notably, the root microbiota was characterized by a prevalence of Cyanobacteria, Proteobacteria, and Actinobacteriota. A marked increase in the proportion of the dominant Cyanobacteria species was observed in roots. At the genus level, the predominant bacterial genera in bulk soil were Bacillus (21.68%), Paenibacillus (11.52%), and RB41 (7.90%). In rhizosphere soil, Pseudomonas (6.23%), RB41 (3.63%), and Catenulispora (2.10%) exhibited the highest proportions. Similarly, Pseudomonas (4.76%), Steroidobacter (1.46%), and Bacillus (0.34%) were identified as the top three bacterial genera in the root (Figure 2B).

Regarding fungi, the phyla Ascomycota, Basidiomycota, and Mortierellomycota were identified as the most abundant, aligning with a globally observed soil pattern (Figure 2C). Analysis of the relative abundance of fungal genera indicated higher levels of Fusarium (BS, 1.71%; RS, 3.04%; and R, 27.71%), Clonostachys (BS, 0.50%; RS, 3.51%; and R, 5.15%), and Exophiala (BS, 0.51%; RS, 1.26%; and R, 3.89%) in the roots compared with that in the bulk soil and rhizosphere soil (Figure 2D). Conversely, the relative abundance of fungal genera Mortierella (BS, 8.48%; RS, 6.79%; and R, 0.52%), Tetracladium (BS, 11.04%; RS, 1.37%; and R, 0.68%), and Pseudogymnoascus (BS, 1.33%; RS, 7.90%; and R, 0.25%) demonstrated higher levels in the rhizosphere soil and bulk soil compared with that in the root.

The α-diversity index serves as a quantitative metric for assessing species richness within microbial communities, with higher diversity being considered advantageous for overall soil health. Analysis of α-diversity revealed significant variations among sampling compartments (bulk, rhizosphere soil, and root), demonstrating higher diversity in the bulk soil and rhizosphere soil compared with that in root. The Simpson, Chao1, ACE, and Shannon indices for both bacteria and fungi exhibited a sequential decrease in the order of bulk, rhizosphere soil, and root, with the α-diversity in root being the lowest (Figures 3A, B). Kruskal–Wallis analysis further confirmed α-diversity significant differences (p < 0.01) among bulk, rhizosphere soil, and root for both bacteria and fungi (Supplementary Figure S2, Supplementary Information). For β-diversity, PCoA using Bray–Curtis distance assessed the compositional relatedness of bacterial and fungal communities. Results indicated distinct separation of microbiome communities across different sampling compartments. In the case of bacteria, PCoA1 and PCoA2 accounted for 52.45% and 21.38% of the total variance, respectively, with the cumulative variance of PCoA contributing to 73.83% (Figure 3C). Similarly, for fungi, PCoA1 and PCoA2 explained 32.91% and 16.81% of the total variance, respectively (Figure 3D). The PCoA analysis underscored pronounced differences in the community composition of bacteria and fungi among bulk, rhizosphere soil, and root. Collectively, these findings emphasize the distinct microbial community compositions across the sampled compartments, underscoring the ecological heterogeneity within the studied soil system.

The LEfSe methods were employed to discern features exhibiting significant differential abundance among bulk, rhizosphere soil, and root, subsequently determining biomarker microbiota at the family level with a stringent criterion of linear discriminant analysis (LDA) score > 4 (Supplementary Figure S3). For bacteria (Figure 4A), the LDA scores for Bacillaceae and Paenibacillaceae were most elevated in bulk soil, whereas those for Micrococcaceae and Micropepsaceae peaked in rhizosphere soil. Notably, no bacterial biomarkers were annotated at the family level in roots with an LDA score exceeding 4. In the fungi (Figure 4B), the family-level LDA scores revealed distinct biomarkers. Helotiaceae and Mortierellaceae exhibited markedly higher scores in bulk soil, whereas Geastraceae and Pseudeurotiaceae dominated in rhizosphere soil. On the other hand, roots featured a higher abundance of Nectriaceae and Bionectriaceae. Consequently, we identified 17 biomarker microbial families, namely, Bacillaceae, Paenibacillaceae, Pyrinomonadaceae, Gemmatimonadaceae, Alicyclobacillaceae, Alicyclobacillaceae, Micrococcaceae, Micropepsaceae, Catenulisporaceae, Helotiaceae, Mortierellaceae, Leptosphaeriaceae, Geastraceae, Pseudeurotiaceae, Nectriaceae, Bionectriaceae, and Herpotrichiellaceae. These biomarkers contribute to the distinctive microbial signatures characterizing bulk, rhizosphere soil, and root in the studied soil ecosystem.

To unravel the correlation between microbiome communities and the principal bioactive chemical compounds in P. amurense, a total of 4,717 bacterial OTUs and 1,344 fungal OTUs from roots, rhizosphere soil, and bulk soil, alongside data on the seven main alkaloids, were utilized to construct a WGCNA network. This analysis identified 30 and 16 co-expression modules for bacteria and fungi, respectively, based on their similar expression patterns (Figures 5A, B). In bacteria, MEpurple is significantly related to Tetrahydropalmatine. Additonally, the module MEdarkgrey is significantly correlated with Magnoflorine; And the module MEturquoise is significantly related to Palmatine (Figure 5C). In fungi, MEmagenta is significantly related to Tetrahydropalmatine and the module MEred is significantly related to the Palmatine (Figure 5D). Subsequently, correlation analysis between key microorganisms within these modules and the seven BIAs of P. amurense was conducted, and the networks were visualized using Cytoscape. Regarding bacteria, all seven benzylisoquinoline alkaloids displayed significantly positive correlations (p < 0.01) with Streptosporangiaceae and Sphingobacteriaceae (Figure 5E). However, Chitinophagaceae showed a negative correlation with tetrahydropalmatine, whereas Paenibacillaceae exhibited negative correlations with tetrahydropalmatine, berberine, and jatrorrhizine. In fungi, berberine, palmatine, jatrorrhizine, phellodendrine, tetrahydropalmatine, magnocurarine, and magnoflorine exhibited significantly positive correlations (p < 0.01) with Nectriaceae and Cunninghamellaceae (Figure 5F). These findings provide insights into the intricate interplay between microbial communities and the biosynthesis of bioactive compounds in P. amurense.