Pebbles (10 cm in size) were placed at the bottom of each pot, which measured 15 cm in diameter and 20 cm in height. Approximately 1 kg of soil was then added. The soil mixture comprised 50% autoclaved sand, 40% autoclaved ginseng field soil, and 10% raw ginseng field soil to introduce native microbiota. The raw ginseng field soil had been used for ginseng monoculture for 9 years. Each pot contained one-year-old ginseng seedlings and was treated with 50 mM of different nitrogen sources [NH₄Cl, glutamate (Glu), aspartate (Asp), asparagine (Asn), and valine (Val)] and untreated control (distilled water) on the 0th and 10th days. The ginseng seedlings were grown under controlled light conditions (25°C for 16 h) and dark conditions (20°C for 8 h) for 20 days per cycle and each cycle had three biological replicates (n = 25). This study was designed to explore the potential implications of microbiota shifts in the field, particularly concerning their association with nitrogen sources and the occurrence of root rot disease.

Ginseng plants were selected for sampling in the field located at coordinates 36°56′31.7″ N and 127°45′04.8″ E, under the administration of the National Institute of Horticultural & Herbal Science in Korea. The selected ginseng plants ranged in age from 2 to 6 years since their initial planting in the field. The ginseng roots were first gently freed from bulk soil to collect rhizosphere soil. Subsequently, each ginseng root was placed in a sterile beaker containing 200 mL of pre-chilled phosphate-buffered saline (PBS) buffer. The PBS buffer was prepared by dissolving 200 mg of KCl, 8 g of NaCl, 245 mg of KH₂PO₄, and 1.44 g of Na₂HPO₄ in L of distilled water, resulting in a pH of 7.4. The beaker, containing the ginseng root in PBS, was then subjected to a sonication bath (Bandelin Electronic GmbH & Co. KG, Berlin, Germany). The rhizosphere soil, which detached from the root due to ultrasound treatment, was subsequently separated and collected through centrifugation at 3,000 × g for 20 min using a 1736R centrifuge (LaboGene, Seoul, Republic of Korea).

Metagenomic DNA from the rhizosphere was extracted using the FastDNA™ SPIN Kit for Soil (MP Biomedicals, Solon, Ohio, United States) according to the manufacturer’s protocol. Initially, 978 μL of sodium phosphate buffer and 122 μL of MT buffer were sequentially added to a Lysing Matrix E tube. Subsequently, the tube was subjected to homogenization using the FastPrep-24™ Classic Instrument (MP Biomedicals, Santa Ana, California, United States) set to a speed of 6.0 for 40 s. Following homogenization, the tube was centrifuged at 18,000 × g for 10 min. The resulting clear supernatant obtained after centrifugation was carefully transferred to a clean 2 mL tube, to which 250 μL of protein precipitation solution was added and gently inverted 10 times. After the precipitation of proteins and their subsequent removal by centrifugation, the supernatant was transferred to a 15 mL conical tube. The binding matrix settled at the bottom of the solution was resuspended, and 1 mL of this resuspended solution was added to the conical tube. The tube was then inverted for 2 min, followed by a 3 min period of settling to allow the binding matrix to collect. The top 500 μL of clear supernatant was removed without disturbing the matrix, and the matrix was resuspended by pipetting and subsequently collected through SPIN™ and centrifugation (18,000 × g for 10 min). The collected matrix was washed with ethanol, and the extracted DNA was eluted using 50 μL of DNase-free water.

We conducted 16S rRNA V4 library sequencing using the Illumina adapter-linked primers 515F (5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG GTGCCAGCMGCCGCGGTAA-3′) and 805R (5’-GTCTCGTGGGCTCGGAGATGTGTATA AGAGACAGGACTACHVGGGTATCTAATCC-3′). This sequencing was performed with metagenomic DNA extracted from the rhizosphere of ginseng cultivated in a field setting. The KAPA HiFi HotStart ReadyMix (Kapa Biosystems, Wilmington, Massachusetts, United States) was utilized for polymerase chain reaction (PCR) amplification. The PCR conditions employed were as follows: an initial denaturation step at 95°C for 3 min, followed by 25 cycles consisting of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s. A final extension step was carried out at 72°C for 5 min. Subsequently, the sequencing process was outsourced to Macrogen (Seoul, Republic of Korea). To attach linker and barcode sequences, a second PCR thermal cycling step was conducted with the following conditions: an initial denaturation at 95°C for 3 min, followed by 8 cycles consisting of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s. A final extension step was carried out at 72°C for 5 min. The obtained sequences were read on a MiSeq 2 × 300 bp platform (Illumina, San Diego, California, United States) to facilitate downstream analysis.

Metagenome data underwent rigorous processing procedures using a high-performance workstation equipped with an AMD Ryzen Threadripper 3,970X processor and 64 GiB of DDR4 RAM, running Ubuntu 18.04. The raw data, including previously processed datasets and newly acquired data from field-cultivated ginseng samples, were integrated to optimize error rate prediction through the divisive amplicon denoising algorithm (DADA) (Callahan et al., 2016) and minimize fragmentation. For the ginseng serial cultivation experiment, involving 40 pot samples, and the nitrogen source treatment, involving 30 samples, as well as the field study comprising 15 samples, the 805R Illumina adapter (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGGACTACHVHHHTWTCTAAT-3′) was employed. Primer sequences were meticulously removed using Cutadapt (version 2.10; Martin, 2011). The resultant primer-free raw sequences, obtained from the MiSeq platform with a 2 × 300 bp configuration, were subsequently imported into R (version 4.0.4). Sequences with an average accuracy of less than 99.9% at the 3′ end were trimmed. The DADA2 package (version 1.16.0) in R was employed to model the error rate for the trimmed sequences. Forward and reverse reading sequences were merged following sequence modification by DADA. Initially, clustering was performed based on perfectly identical sequences. Chimeric and non-bacterial amplicon sequence variant (ASV), classified using IDTAXA (Murali et al., 2018) referencing SILVA 138 SSU during the first clustering step, were removed to enhance alignment accuracy and to eliminate an additional 2 bp introduced by the 515F-806R amplicon. Following the removal of the 2 bp, a second clustering step based on completely identical sequences was executed. The ASVs generated during the second clustering were subjected to re-identification using the IDTAXA approach. The results of these procedures, including the number of reads, are presented in Supplementary Table S1.

Bacterial community analysis was conducted using a combination of software tools in the R programming environment. The analysis utilized phyloseq (version 1.32.0; McMurdie and Holmes, 2013), vegan (version 2.5–0; Oksanen et al., 2012), and DESeq2 (version 1.18.2; Love et al., 2014). Additional details and code for this analysis are available in the associated GitHub repository^{1, 2, 3}.

For predicting the pathways of each sample, we employed PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States; version 2.3.0-b) with KEGG and MetaCyc databases. This analysis was performed in Python (version 3.7.8) (Douglas et al., 2020).

To account for variations in α diversity between field and pot samples, we utilized the β-nearest taxon index (βNTI) (Stegen et al., 2012). The maximum likelihood tree required for βNTI calculations was generated using RAxML-NG (Randomized Axelerated Maximum Likelihood Next Generation, version 1.0.2) (Kozlov et al., 2019).

The R code used for these analyses is available in the following GitHub repository: https://github.com/gyeongjunCho/R-code-of-Ph.D.-Thesis (Cho, 2023).

The rarefaction curve demonstrates the effectiveness of ASV clustering, revealing a minimal likelihood of encountering new ASVs with increasing read depth (Supplementary Figure S1). Furthermore, we observed significant variations in diversity among different samples, including field samples, rhizosphere samples from serial cultivation in pots, and rhizosphere samples subjected to nitrogen source treatments. These differences were evident in the ASV counts, Shannon index, and Simpson index (Figure 1B). Notably, the α diversity indices followed a consistent order, with field samples exhibiting the highest diversity, followed by nitrogen treatment samples in pots and serial cultivation samples in pots. In the context of continuous ginseng replanting and harvesting in pots, the initial influencers displayed an initial suppression of root rot progression. However, their relative abundance in the rhizosphere gradually decreased, eventually leading to a peak in disease progression (Cho et al., 2023). When different nitrogen sources were applied to pots containing bulk soil from previous replanting pots, a significant increase in Pseudomonadaceae was observed. Interestingly, the disease was effectively suppressed in Asp (aspartate) treatment but significantly increased in Val (valine) treatment (Cho, 2023). This transition in the rhizosphere community dynamics during continuous ginseng cultivation in pots highlighted the transfer of disease suppression responsibility from the initial community influencers to Pseudomonadaceae. Consequently, distinct levels of disease occurrence were observed across these periods (Cho, 2023).

To compare the field and the pot rhizosphere with β diversity, relative abundance was calculated (Figure 1A) and principal coordinates analysis (PCoA) with Bray-Curtis distance was performed at the ASV level (Figure 2A). The result indicated that there was a significant difference between the pots and fields, as well as that pots could not describe fields’ results. However, Bray-Curtis distance from PCoA calculated with PICRUSt2 pathway prediction showed high similarity between field and NH₄Cl treated pot samples (Figure 2B). As a result of analysis with permutational multivariate analysis of variance (PERMANOVA) with Bray–Curtis distance, the field and almost all pot samples were significantly different, but there was no significantly difference between field and NH₄Cl treatments (Figure 2B and Supplementary Table S2). The distribution between the nitrogen source treated rhizosphere and the late replanting rhizosphere overlapped the relationship because the design of these previous studies at pot was well explained (Cho, 2023). These results indicate that the composition and diversity of rhizosphere microbiota communities are significantly affected by different nitrogen sources. The α diversity of the field was significantly higher than other treatments, and field the β diversity in PCoA at the ASVs level was different. Additionally, PCoA indicated that the nitrogen source treatment study was the directly continuous study of the cycling study because the nitrogen source treatment experiment was generally far from the first to third cycling samples in the ginseng cycling and overlapped with the 4th to 10th cycling samples (Figure 2A).

Since it was difficult to compare the field with the ASV level alone, the phylogenetic tree of the 16S rRNA V4 region of ASV was calculated as maximum likelihood, and beta NTI was calculated together with the phylogenetic tree and the relative abundance (Figures 3A,B). Importantly, the similarity of microbiota community structure between the field and pot rhizospheres was very high (βNTI < −2). When measuring how similar the phylogenetic diversity (β diversity) was, it was found to be the highest at 84.79% when compared to the field itself, and the second highest at 74.67% when compared to the field treated with NH₄Cl treatment (Figure 3C). The treatment first cycle, Glu treatment, and the second cycle had the highest similarity to field treatment (56.67, 56.00, and 55.00%, respectively). The remaining similarity did not exceed 50%. Assuming that the difference in β diversity is primarily influenced by the difference in α diversity, it is expected that the difference will decrease when considering β diversity with the phylogenetic tree. Therefore, βNTI was calculated.

Therefore, it is reasonable to group the field rhizospheres and the NH₄Cl-treated rhizospheres when comparing the relative abundance. Since the phylogenetic beta diversity was similar between the field rhizosphere and the NH₄Cl-treated rhizosphere, it was decided that relative abundance could be compared in the microbial community.

In our investigation, we quantified and visualized the shared microbial taxa using a Venn diagram (Figure 4A). We identified 3,723 ASVs in the field rhizosphere, 1,548 ASVs in the rhizosphere of the serial cultivation study, and 2017 ASVs in the rhizosphere of the nitrogen source treatment. Remarkably, 593 ASVs were found to be shared among all three categories, with their relative abundance medians measuring 65% in the field, 73% in the serial cultivation, and 78% in the nitrogen source treatment (Figure 4B). At the family level, we observed 236 ASVs in the field, 173 in the serial cultivation, and 192 in the nitrogen source treatment (Figure 4C). The relative abundance median of the 146 shared families was consistently high, with values of 99% in the field, 99% in the serial cultivation, and 99% in the nitrogen source treatment (Figure 4D). This indicated that family-level abundance comparisons were valid. To offer an alternative approach for comparison, we considered the predicted pathways using PICRUSt2, as the functional genes of the bacteria were expected to exhibit limited variation across different bacterial taxa. Our PCoA of the predicted pathways revealed significant distinctions in the bacterial community composition of the field rhizosphere compared to the pot rhizosphere, except in the case of the NH₄Cl treatment (Figure 2B; Table 1; Supplementary Table S2). The analysis of phylogenetic β diversity and the PCoA results of the predicted pathways provided further evidence that the bacterial community characteristics in the field closely resembled those of the NH₄Cl treatment.

We conducted a comparative analysis between the field and NH₄Cl treatment group and another group comprising various rhizosphere samples using the DESeq2 R package. In the field and NH₄Cl treatment groups, we observed a significant reduction in the abundance of pathways that had previously shown negative correlations with root rot disease progression in a study involving different nitrogen sources (Wang et al., 2019; Ding et al., 2021). These pathways included ‘arginine, ornithine, and proline interconversion,’ ‘glucose degradation (oxidative),’ ‘pyridoxal 5’-phosphate biosynthesis I,’ ‘l-arginine degradation II (AST pathway),’ and ‘polymyxin resistance’ (Figure 5A). In contrast, the ‘terpenoid biosynthesis’ pathway, including the ‘super pathway of geranylgeranyl diphosphate biosynthesis I (via mevalonate),’ ‘mevalonate pathway I,’ and ‘isoprene biosynthesis II (engineered),’ which had shown negative correlations with F. solani population density in the previous nitrogen source treatment study (Griffin, 1970; Morgan and Timmer, 1984), remained unchanged. All of the identified pathways were subjected to ontology analysis at the MetaCyc ontology class level (Figure 5C), revealing notable differences in the carbon and nitrogen source circulation within the rhizosphere bacterial community between the field group and the other samples.

We further explored the DESeq2 results at family level, focusing on families that were instrumental in influencing the initial rhizosphere community to inhibit root rot progression. In the context of the serial cultivation study, these influential families included Rhizobiaceae, Beijerinckiaceae, Devosiaceae, Xanthomonadaceae, Sphingobacteriaceae, Sphingomonadaceae, and Comamonadaceae. Additionally, Pseudomonadaceae, which had been associated with pathways negatively correlated with root rot progression and F. solani population density in the nitrogen source treatment (Figure 5B). Within this analysis, we found that Comamonadaceae and Sphingomonadaceae were more abundant in the field and NH₄Cl treatment group, whereas Pseudomonadaceae and Sphingobacteriaceae were less abundant in this group. These findings were consistent with the results from the field where ginseng was grown for 2–6 years. The pathways negatively correlated with root rot progression, including ‘arginine, ornithine, and proline interconversion,’ ‘glucose degradation (oxidative),’ ‘pyridoxal 5’-phosphate biosynthesis I,’ ‘l-arginine degradation II (AST pathway),’ and ‘polymyxin resistance,’ were less abundant in the Field and NH4Cl treatment group. Notably, the density of F. solani, the root rot causative pathogen, which had shown negative correlations with pathways in previous studies, did not differ between the field and NH₄Cl treatment group and the other samples (Morgan and Timmer, 1984). Regarding the initial community influencers identified in the continuous study, we observed varying levels of abundance in the field and NH₄Cl treatment groups. However, Pseudomonadaceae, which had played a significant role in the eight identified pathways, was significantly less abundant in the field and NH₄Cl treatment group. This observation suggests that the rhizosphere field may be more susceptible to disease, independent of F. solani population density, and emphasizes the potential for disease suppression by Pseudomonadaceae due to its lower abundance in the field compared to other samples.

A previous field study extends and compares the findings from previously reported artificial continuous cultivation of ginseng and nitrogen source studies (Cho, 2023), both of which are associated with the rhizosphere community. This investigation focuses on the later stages of the continuous cultivation, providing an evaluation of the ginseng rhizosphere in a field setting. Our results demonstrate that the initial community in the field deteriorated within 2 years after planting and closely resembled the NH₄Cl treatment group at the time of sampling. This shift is attributed to the predominant use of fertilizers primarily composed of ammonium rather than amino acids during continuous cultivation. Furthermore, we observed a reduction in disease-correlated pathways and Pseudomonadaceae, which had played a crucial role in these pathways in the previous study, in the field of rhizosphere. This study underscores the susceptibility of the field rhizosphere to disease, corroborating previous research findings. Considering our findings, we propose the potential for root rot disease suppression in ginseng fields using Asp and Pseudomonadaceae.