Seeds of sugar beet (KWS1176 from KWS Company of Germany) were selected as plant materials in this experiment. KWS1176 is a pelletized seed which is film-coated and wrapped in a thin layer containing thiram to kill pathogens on the seeds after disinfection. In our previous studies, it was found to exhibit strong resistance to adversity and is more widely used in Northeast China. The soil was sourced from the black soil area in Hulan District, Harbin City, Heilongjiang Province (latitude and longitude: 46°00′14′′ E, 126°38′49′′ N). Soil samples included soils from a maize-beet rotation (previous crop was maize) and soils from 3 years of continuous sugar beet crop (2-year recrop soils), both collected from adjacent plots with similar agricultural management practices. We sowed the beet pelleted seeds after washing off the coating with sterile water. The seeds were sown in a plastic pot containing 0.7 kg of soil. In this experiment, each pot was sown with six seeds and was irrigated with 50 mL of Hoagland nutrient solutions every 10 days. The samples were distributed as follows: non-continuous cropping group and continuous cropping group. To ensure the normal growth of plants, only one seedling was selected in each pot after 5 days of planting. All pots were incubated in a laboratory with the following photoperiod: day (lights on) 7 am–9 pm, 24°C; night (lights off), 19°C.

Harvesting was carried out on the 30th day. Beetroots, rhizosphere soils, and bulk soils were processed as previously described (Bulgarelli et al., 2012). Each treatment group was randomly selected from uniformly growing sugar beet, and five replicates of each treatment were recorded as one replicate by taking three pots of samples and mixing them. Beet plants were manually collected from the plastic pots and bulk soil aggregates were separated by shaking plant roots. The bulk soil samples were sampled from the pots and passed through the sterile 2-mm sieve. The bulk soil was divided into two parts for subsequent analysis. One portion was air-dried for soil properties and soil enzyme activities, and the other portion was frozen in liquid nitrogen and stored at −80°C for the sequencing of microbes. To standardize the sampling, we used a sterile scalpel to dissect the roots under the beet seedling stems to maximize repeatability. Beetroots were pooled into a sterile 50 mL tube containing 25 mL of sterile Silwet L-77 amended phosphate-buffered saline solution (PBS) (per liter: 7 mM Na₂HPO₄, 3 mM NaH₂PO₄, 200 μL Silwet L-77, and pH 7.0) and were vigorously mixed with a vortex for 15 s to detach from the root surfaces to produce the rhizosphere soil. Then samples were subjected to centrifugation at 3,200 g for 15 min to get the rhizosphere soil. For beetroot DNA extraction, roots were transferred to a new 50 mL sterile tube containing 25 mL of PBS, and vortexed. Repeat the steps until the phosphate buffer in the centrifuge tube was no longer turbid. Then, the root was transferred into another sterile tube and sonicated for 10 cycles, consisting of 30 s bursts and rests of 30 s. The purpose of this step was to remove the microorganisms attached and further clean the outer surface of the roots. Finally, the roots were transferred to a fresh volume of 25 mL of PBS and the sonicated roots were defined as beetroot. The beetroot samples and rhizosphere soil were then stored at −80°C until DNA extraction. Through the use of scanning electron microscopy and microbial culture techniques, Xiao et al. (2017) confirmed that the microorganisms on the root surface could be removed by the above steps.

Five plants were randomly selected for plant height, root fresh weight, and dry weight. Five sugar beet plants of uniform growth were randomly selected from each treatment group and washed with deionized water for measuring root vigor. The root vitality of the sugar beet plants was measured according to the alpha naphthylamine method (Ramasamy et al., 1997). Root morphology and root area (five plants) were determined by WINRHIZO software (Regent Instruments Inc., Quebec, QC, Canada). The air dried bulk soil samples were passed through a 2-mm sieve as described by Bao (2005). Air-dried soils were used to determine soil physicochemical properties and soil enzyme activity (pH, EC, AN, AK, AP, sucrase, acid phosphatase, CAT, and urease). Specifically, the hydrogen potential (pH) and electrical conductivity (EC) of the soil were measured with a potentiometer after shaking for 30 min at the soil to distilled water ratio of 1:2.5. Soil available nitrogen (AN) was measured by changing available nitrogen into ammonia. Soil available phosphorus (AP) was determined by leaching with sodium bicarbonate solution and molybdenum-antimony colorimetric method. Soil available potassium (AK, extracted in ammonium acetate) was determined by the flame photometer. Soil catalase activity (CAT) was measured according to KMnO₄ titration. Soil urease was based on urea and the enzyme activity was determined based on the enzymatic product ammonia interacting with phenol-sodium hypochlorite to produce blue indophenol. Soil sucrase content was determined by a colorimetric method using 3,5-dinitrosalicylic colorimetry and the amount of reducing sugars was used to express the enzyme activity. Soil acid phosphatase activity was determined by a colorimetric method using sodium benzene phosphate to express enzyme activity as phenol content. All bulk soil samples were assayed in five replicates, and inorganic controls were required for each treatment for the soil enzyme assay.

The non-continuous cropping sugar beet bulk soil (Sn), its rhizosphere soil (Rn), its beetroot (Bn), and continuous cropping samples (Sc, Rc, and Bc) were extracted as research materials. DNA was extracted from bulk soil, rhizosphere soil, and root samples (0.5 g per sample) by the E.Z.N.A^®. Soil DNA Kit (Omega Bio-tek, Norcross, GA, United States). DNA quality and integrity were examined using 1% agarose gel electrophoresis. 16S rRNA genes were selected for the V3–V4 region and fungal ITS amplification was selected for the ITS1–ITS2 region (Huang et al., 2020). This 16S and ITS high throughput sequencing was performed by BIOZERON Ltd. (Shanghai, China).

The MiSeq platform PE250 (Illumina, Inc., CA, United States) was used for double-end sequencing and the PANDAseq software was used for splicing to obtain long reads with highly variable regions. QIIME was used to analyze the raw sequence data, and the 250 bp reads were truncated, receiving a low-quality score (≤20), and those shorter than 50 bp were discarded (Caporaso et al., 2010). Removal of chimeras using UCHIME software (Edgar et al., 2011). OTU clustering based on 97% sequence similarity using the UPARSE plug-in (Edgar, 2010). Species annotation of OTU representative sequences against the Greengenes database was performed using the RDP classifier (80% confidence interval) (Wang et al., 2007). OTUs with annotated results for chloroplasts and mitochondria or OTUs with only one sequence (singletons) were removed. Given the variation in sequencing depth between samples, the OTUs were normalized for subsequent analysis using the least square method.

The rarefaction was calculated using the software Mothur (Schloss et al., 2009) and conducted to reveal the alpha diversity, including the Chao1 and Richness indices (Schloss et al., 2011). To test the statistical significance of the structural similarity between communities with different sampling treatments, UniFrac-based hierarchical clustering (Lozupone et al., 2011) was performed using the community ecology package (Wang et al., 2012). Displacement multivariate analysis of variance (PERMANOVA) and non-metric multidimensional scaling analysis (NMDS) were used to examine and visualize differences in the structure of the bacterial and fungal communities in the three compartments, respectively. To identify biomarkers of important microbial taxa, LEfSe (linear discriminant analysis effect size) analysis was performed as previously reported (Segata et al., 2011). The relationship between changes and dissimilarities was analyzed by the Kruskal–Wallis sum-rank test, and the threshold was set at 0.05. For LDA analysis, the threshold for the size effect of each OTU was set at the threshold for abundant taxa (Ijaz et al., 2018). To further evaluate the relationship between microbial community composition and nine factors of physicochemical properties, the spearman rank correlation and mantel tests were visualized using the “ggcor” software package (Yuan M. M. et al., 2021). Functional annotation of 16S rRNA bacterial gene sequences from the SILVA database was performed using the Tax4Fun software package and the FAPROTAX package. Analysis of variance, Tukey’s test, and Duncan’s test were performed using IBM SPSS 19.0 software. Box line plots were drawn using PRISM (Version 7.0) software. The rest of our analyses were performed with the aid of R software (Version 4.0.0) and the BIOZERON cloud platform¹.

While the soil pH and AK were decreased, EC, AP, and AN were increased in the continuous cropping group compared to that of the non-continuous cropping group (Table 1). In addition, the continuous cropping group exhibited higher Catalase and urease activities, while sucrase and acid phosphatase activities were higher in the non-continuous cropping group (Table 1). We found that the plant height (Height), fresh weight (FW), root dry (DW), root vitality (RA), and root surface area (RSA) in the site with the continuous cropping group were significantly lower than the group with the non-continuous cropping group (Table 2 and Supplementary Figure 1).

We evaluated the alpha diversity of different compartments using the Chao1 and Richness indices to measure microbial community diversity and richness, respectively (Figure 1 and Supplementary Table 1). The continuous cropping group decreased the bacterial richness index (Chao1 and Richness) in three belowground compartments compared to the non-continuous crop group. However, the continuous cropping group had different effects on the fungal richness indices of the three underground compartments. The continuous cropping group significantly decreased the fungal diversity of the bulk soil and increased the fungal diversity of the beetroot compared to the non-consecutive cropping group. The results of the two-way ANOVA (Supplementary Table 2) showed that compartment had a significant effect on alpha-diversity (P < 0.01). This indicated that continuous cropping can affect the alpha diversity in plants and soils.

The NMDS with the unweighted UniFrac algorithm demonstrated that the bacterial and fungal communities from the different compartment samples were clustered separately (Figures 2A,B). It showed that these samples under continuous cropping group and non-continuous cropping group formed distinct clusters in the plotted ordination space. The bulk soil and rhizosphere soil were clustered together, while they were significantly separated from the beetroot. A heat map of the Beta diversity index was constructed (Figures 2C,D). The results revealed that the bulk soil shared the highest level of correlation with rhizosphere soil. It suggested that microbial community composition varies considerably across compartments. Moreover, the soil and beetroot microbial communities differed considerably between the two treatment groups.

This study also showed the composition of bacterial and fungal species from phylum to genus level in beetroot, rhizosphere soil, and bulk soil samples subjected to the different types of soil treatment alternatives. The proportional abundance of dominant taxa changed under different cropping systems and compartment conditions.

The operational taxonomic units (OTUs) of the different compartments of bacterial communities belonged to 39 phyla, 113 classes, 259 orders, 413 families, and 846 genera. The top 10 bacterial phyla in terms of compartment abundance included Proteobacteria, Actinobacteria, Bacteroidota, Acidobacteriota, Cyanobacteria, Chloroflexi, Patescibacteria, Gemmatimonadetes, Myxococcota, and Verrucomicrobiota (Figure 3A). The total of these phyla accounted for more than 98.2% of the bacterial sequences. Proteobacteria and Actinobacteria accounted for 66% of the bacterial community and were the two largest phyla. The relative abundance of Proteobacteria increased in the continuous cropping, respectively, compared to levels in the non-continuous cropping, and the abundance of rhizosphere soil and beetroot was greater than that of bulk soil (Figure 3C). Furthermore, the abundance of Actinobacteriota decreased. Pseudomonas, Flavobacterium, and Pedobacter were significantly higher in abundance in the continuous cropping samples. In contrast to bulk soil, the rhizosphere soil and beetroot harbored a higher portion of Novosphingobium.

The fungal OTUs belonged to 7 phyla, 25 orders, 71 families, 163 families, and 363 genera in the high-throughput sequencing results. The four phyla were Ascomycota, Olpidiomycota, Basidiomycota, and Mucoromycota (Figure 3B). Ascomycota was dominant in the continuous cropping group, especially continuous cropping beetroot (84.0%), and Olpidiomycota was dominated by non-continuous cropping beetroot (85.8%) and rhizosphere soil (76.7%). The relative abundance of Tausonia, Fusarium, and Gilbellulopsis increased in the continuous cropping, respectively, compared to levels in the non-continuous cropping. Moreover, the Tausonia abundance of rhizosphere and bulk soil was greater than that of beetroot, and the Fusarium abundance of beetroot is greater than that of other treatment groups (Figure 3D). Furthermore, the continuous cropping group Olpidium was reduced and had the lowest abundance compared to the non-continuous cropping group.

In addition, ternary plots showed that many phyla were present in similar proportions in the three compartments (bulk soil, rhizosphere soil, and beetroot), but that some were comparatively more abundant at a specific position (Figure 4). We found that a proportion of bacterial phyla were significantly enriched in beetroots, while bulk and rhizosphere soils shared the majority of bacterial phyla (Figure 4A). The ternary plots showed that the dominant bacterial phylum was similarly distributed in the different compartments, while the fungal phylum differed more microbially at the level of the fungi. The specific fungal phylum in the beetroots of the continuous cropping group was Ascomycota, whereas the specific fungal phylum in the non-continuous cropping group was Olpidiomycota in both rhizosphere soil and beetroots (Figure 4B).

Differences in compartment microorganism community were assessed using linear discriminant analysis effect size (LEfSe) analysis at a linear discriminate analysis (LDA) threshold of 3 (Figure 5). Across the bacterial community, there were more species of Actinobacteria, Acidobacteria, Bacteroides, Proteobacteria, and there were species differences between the six test groups (Figure 5A).

For root microbial communities, continuous biomarkers included Flavobacterium (its phylum to genus), Streptomyces (its order to genus), Rhizobium (its order to genus), Comamonadaceae (family), Pseudoxanthomonas (genus), and Stenotrophomonas (genus), while the non-continuous cropping biomarkers mainly included Glycomyces (its phylum and order), Lechevalieria (its order to genus), and Steroidobacter (its order to genus). For the rhizosphere microbial community, continuous cropping rhizosphere soil biomarkers included Massilia (its order, family and genus), Pseudomonas (its order to genus), Pedobacter (by order to genus), Xanthomonadaceae (its order and family) and Proteobacteria (family), while the non-continuous cropping rhizosphere soil biomarkers mainly included Novosphingobium (its order to genus) and Burkholderia (its order to genus). In the bulk soil microbial community, Gemmatimonadaceae (its phylum to family), Intrasporangiaceae (its phylum to family), Frankiales (order), Gaiellales (order), Micrococcaceae (family), Arthrobacter (genus), and Pseudarthrobacter (genus) were significantly enriched in the continuous cropping group, while RB41 (its phylum to genus), Vicinamibacteraceae (its order to family), Nocardioidaceae (its order to family), Solirubrobacterales (its order and phylum) Microtrichales (order), MB_A2_108 (order), Chloroflexia (its phylum and order), and KD4_96 (order) were significantly enriched in non-continuous cropping group.

Among the fungal communities, Ascomycota and Basidiomycota species were more numerous, and there were species differences between the five test groups (Figure 5B). For the beetroot microbial community, the phylum level of continuous cropping beetroot biomarkers included Ascomycota and Basidiomycota, while the genus level included Alternaria, Plectosphaerella, Dactylonectria, and Fusarium, while Olpidiomycota (from phyla to genus) was significantly enriched in non-continuous cropping beetroot. For the rhizosphere soil microbial community, continuous cropping rhizosphere soil biomarkers included Trichocladium and Tausonia (from phylum to genus). In the bulk soil microbial community, Phoma (its family and genus), Pseudombrophila (by order to genus), and Gibellulopsis (genus) were significantly enriched in continuous cropping bulk soil, while Chaetomiaceae (its order and family), Helotiales (its class and order), and Mortierella (by order to genus) were significantly enriched in non-continuous cropping bulk soil.

The environmental drivers of changes in the microbiome of three belowground compartments were explored by calculating correlations between the microbial community composition and soil properties (Figure 6). To conclude which environmental factors caused changes in the composition of the microbial community in the compartments, we correlated differences in functional community composition with soil properties employing distance correction. Soil pH, catalase, urease, and AK were significantly and negatively correlated with AN, AP, EC, acid phosphatase, and sucrase. As shown in Figure 6, the environmental factors were highly correlated with both bulk soil and beetroot of the bacterial (Figure 6A) and fungal (Figure 6B) microbiomes (P < 0.05). The fungal microbiome of rhizosphere soil was significantly related to environmental factors (P < 0.05). However, the bacterial rhizosphere soil microbiome only had a significant correlation with EC. This suggested the diversity and composition of the rhizosphere fungal community were more closely related to environmental factors than to bacteria. Besides, environmental factors were significantly related to both the functional (based on FAPROTAX annotation) and fungal taxonomic composition of the compartments’ microbiome. The bacterial taxonomic composition had no significant correlation with the soil environmental factors (P > 0.05).

In the FAPROTAX database of microbial ecological function predictions, the annotated OTU was assigned to 61 predicted functional groups. Nevertheless, in the Kruskal–Wallis test, only 50 groups showed significant differences between the six compartments (P < 0.05). Therefore, they were plotted as a functional heatmap (Figure 7). Among these function predictions, nitrogen (10 groups), carbon (6 groups), sulfur (3 groups), and manganese (1 group) were involved in the geochemical cycle. The N cycle-related functions showed different performances in the three compartments, and the N cycle-related functions were basically positively correlated with the bulk soil compartment, such as denitrification, nitrate denitrification, nitrate respiration, and nitrous oxide denitrification. Continuous cropping of bulk soil enhanced the soil’s functional advantages of phototrophy, photosynthetic cyanobacteria, oxygenic photoautotrophy, and aliphatic non-methane hydrocarbon degradation. These functions, including dark oxidation of sulfur compounds, nitrogen respiration, nitrate respiration, and human pathogens, were significantly (P < 0.05) enhanced in the rhizosphere soil and root zone of the continuous cropping area compared to the non-continuous cropping area.

The first reference pathway that contained genes in the first five of the KEGG pathway analysis was metabolic, environmental information processing, genetic information processing, human diseases, and cellular processes. The number of metabolic pathway genes with the highest content was about 25 times higher than that of the second residual genes (Supplementary Figure 2A). The metabolic pathways of the second reference pathway containing the first five genes in KEGG pathway analysis were global and overview maps; carbohydrate metabolism, amino acid metabolism; energy metabolism; and membrane transport (Supplementary Figure 2B). The metabolic pathways in which the third reference pathway contained the first five genes in KEGG pathway analysis are metabolic pathways, biosynthesis of secondary metabolites, microbial metabolism in diverse environments, biosynthesis of antibiotics, and carbon metabolism (Supplementary Figure 2C).