The soil sampling was carried out in Jinhua, situated in the central region of Zhejiang Province, China. The sampling commenced on October 9, 2022, and continued for 10 days, with each sampling site’s process being completed within a single day. Jinhua is located approximately between 119°14′–120°47′ east longitude and 28°32′–29°41′ north latitude, spanning 129 km north to south and 151 km east to west. The orchard area covers 26,290 hectares, with 3,809 hectares dedicated to citrus orchards. This region falls within the subtropical zone, characterized by a warm and humid climate, with hot summers and mild winters. In 2022, the total annual sunshine hours were 1701.2, with a precipitation of 1300.3 mm and an average temperature of 19.4°C.

To ensure the similarity of sampling environments and eliminate other factors’ influence on soil microorganisms, we specifically selected family farms adjacent to rivers for environmental sampling. Four villages with lax management were chosen as sampling sites based on their geographic location. Lax management is defined as a management practice where, apart from essential watering, the intervention frequency (including pesticide and fertilizer usage) is equal to or less than twice a year; family farms with higher intervention frequencies were not considered. Among these, priority was given to farms that either abandoned management or had lower management frequencies.

Village A located on the river downstream of the center of the city and closely adjacent to the city. It is situated in a densely populated area and is strongly influenced by urban population mobility. On the other hand, villages B, C, and D are situated farther from the city and are next to the river but upstream of the city (Figure 1A). Fourteen family farms were designated as individual sampling site (Supplementary Table S1), where five sampling points were selected: one sampling point was in the center of the farm, and one was at each corner. In order to assess disparities in composition of the soil in the rhizosphere and the soil outside of the rhizosphere, two additional sampling points were strategically placed outside the rhizosphere (peripheral soil). Specifically, these points were positioned at the midpoint between two fruit trees, where no root growth was observed.

We conducted sampling using five-point sampling method during the ripening stages of Jinhua’s native citrus. Rhizospheric soil samples were gathered from the rhizosphere, which was considered to be the area around the root system. Five citrus trees aged 8–10 years were selected at each sampling point, and soil samples were collected around the roots of each tree to a depth of approximately 20 cm. We did not sample soil from the roots of citrus trees showing signs of disease or root damage. For the two sampling points located at the midpoint between two fruit trees, peripheral soil samples were collected to the same depth (20 cm), but without including the rhizospheric soil. During soil excavation, areas within a 1-m radius of the soil surface containing evident contaminants and decayed matter were excluded.

Before sampling, we removed the topsoil using a shovel. For sampling of rhizospheric soil, we gently removed soil that was loosely attached to the roots. Every 10 mL of soil from the five locations within each sampling point was placed in separate sterile plastic bags and thoroughly mixed to create a homogenized composite sample. Approximately 15 mL of the homogenized soil was then transferred into a labeled centrifuge tube and stored at −20°C during transportation. Upon arrival at the laboratory, the samples were promptly transferred to a −80°C ultra-low temperature freezer for long-term storage. The detailed locations of the sampling points and relevant environmental information can be found in Supplementary Table S1.

Finally, we systematically collected a total of 98 soil samples from four villages. Within this dataset, villages A (urban), C (rural), and D (rural) were each represented by 21 samples. Village B (rural) contributed a comprehensive set of 35 soil samples due to its broader and more diverse citrus cultivation landscape.

DNA extraction from soil samples was carried out using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. To assess the quality and concentration of the extracted DNA, a spectrophotometer (Nanodrop 2000; Thermo Fisher Scientific, Waltham, MA, United States) was employed. The extracted DNA was then stored at −80°C until PCR analysis.

The primers used to amplify the ITS1 region and the 16S V3–V4 region gene were: ITS-1F, 5′-CTT GGT CAT TTA GAG GAA GTA A-3′ and ITS-2R, 5′-GCT GCG TTC TTC ATC GAT GC-3′; 16S-338F, 5′-ACT CCT ACG GGA GGC AGC AG-3′ and 16S-806R, 5′-GGACTACHVGGGTWTCTAAT-3′.

The first round of PCR amplification was performed using the following cycling conditions: 3 min at 95°C; 25 cycles of 30 s at 95°C, 30 s at 55°C, and 30 s at 72°C; and a final extension step at 72°C for 5 min. Each 25 μL PCR mixture contained 12.5 μL of 2× KAPA HiFi HotStart ReadyMix, 1 μL of each primer (1 μM), and 12.5 ng of template DNA. Post-PCR purification was achieved using AMPure XP beads.

The second round of PCR amplification was performed using the following cycling conditions: 3 min at 95°C; 8 cycles of 30 s at 95°C, 30 s at 55°C and 30 s at 72°C; and a final extension at 72°C for 5 min. Each 50 μL PCR mixture contained 25 μL of 2× KAPA HiFi HotStart ReadyMix, 5 μL of each primer, 10 μL of water, and 5 μL of template DNA from products purified in the previous step. A second round of AMPure XP beads purification was performed.

Following purification, the amplicons were pooled in equimolar proportions and subjected to paired-end sequencing on an Illumina MiSeq platform (Illumina, San Diego, United States), according to the standard protocol provided by Majorbio Bio-Pharm Technology Co. Ltd. The number of sequencing reads obtained can be found in Supplementary Table S2. The sequencing reads have been deposited in the NCBI BioProject under ID PRJNA1007597.

All raw data underwent filtering using Trimmomatic (Bolger et al., 2014) (version 0.39) to eliminate adapters, primers, and low-quality sequences; the parameters were SLIDINGWINDOW:4:15, LEADING:3, TRAILING:3, and MINLEN:80. Subsequently, the processed data were imported into the Quantitative Insights Into Microbial Ecology version 2 (QIIME2) pipeline (Bolyen et al., 2019) (version 2021.11.0) for quality control. The 16S rRNA gene sequences were merged using vsearch (Rognes et al., 2016) (version v2.15.0) and then denoised via deblur using default parameters (Amir et al., 2017) to generate representative sequences. The representative sequences of ITS region gene were inferred using default parameters with the DADA2 plugin (Callahan et al., 2016) plugin.

Taxonomy assignment was performed on all representative sequences after training the species classifier using the q2-feature-classifier (Bokulich et al., 2018) plugin. The training process involved utilizing primers to extract target sequences, followed by filtering out corresponding taxonomic information using the RESCRIPt (Robeson et al., 2021) plugin. Finally, the feature-classifier plugin was used for training. The results of taxonomic analyses indicated that 98.0% of the bacteria were successfully identified at the phylum level using the 16S rRNA gene, whereas only 67.6% of fungi were successfully identified. Visualization of the taxonomy results was carried out on QIIME2 View¹. The amplicon package² of R was used to illustrate species compositions and to create Venn diagrams and circle plots.

Following alignment of the representative sequences exhibiting high abundance using MAFFT (Katoh et al., 2019) (version 7.490), the construction of the maximum likelihood (ML) tree was performed using IQTREE (Minh et al., 2020) (version 2.2.0.3). To determine the best DNA model, the Edge-linked Partition Model (Chernomor et al., 2016) was employed, and branch supports were assessed through ultrafast bootstrap (Hoang et al., 2018).

Functional prediction of the 16S rRNA gene representative sequences was performed with the picrust2 (Douglas et al., 2020) plugin within the QIIME2 pipeline. For BugBase (Ward et al., 2017) phenotypic predictions, submissions were made via the online platform³. The Linear Discriminant Analysis Effect Size (LEfSe) algorithm was processed using the format2lefse function within the amplicon package of R, followed by submission to ImageGP (Keighley et al., 2022) for further analysis.

For functional classification of the ITS region gene, FUNGuild (Nguyen et al., 2016) was employed to parse fungal community datasets based on trophic mode, trait and growth form. Differences between two groups were performed with the Mann–Whitney U test.

Representative sequences were processed using the EasyAmplicon (Liu Y. X. et al., 2021) pipeline. Normalization was conducted using the “otutab_rare.R” script⁴ with default parameters. The parameter “--depth 0” was utilized to automatically determine the minimum rarefaction depth. The diversity plugin of QIIME2 was used to calculate alpha and beta diversity. Alpha diversity was assessed using the Abundance-based Coverage Estimator (ACE) index, and differences were compared using the Kruskal–Wallis test. Principal Co-ordinates Analysis (PCoA) based on Bray–Curtis distance was performed. The resulting data were visualized using the amplicon package and pheatmap package of R.⁵ Differential analysis was conducted using “EdgeR” and results were visualized by using the ‘compare_volcano.R’⁶, ‘compare_manhattan.sh’ (), and ‘compare_heatmap.sh’⁷ scripts to generate volcano plots, Manhattan plots, and comparative heatmaps, respectively. A significance threshold of p < 0.05 was considered for all statistical analyses.

Among the successfully identified bacteria, the top seven abundant phyla were Proteobacteria, Acidobacteria, Actinobacteria, Chloroflexi, AD3, Gemmatimonadetes, and Firmicutes (Figure 1B). For fungi, the top seven abundant phyla were Ascomycota, Basidiomycota, Chytridiomycota, Mortierellomycota, Rozellomycota, Glomeromycota, and Blastocladiomycota. Notably, Ascomycota fungi were predominant across almost all of the soil samples, accounting for an average of 66.9% of the species (Figure 1C).

Rarefaction curve analyses of the sequencing results demonstrated stable patterns across all samples. From the rarefaction curves, it appeared that village A, which is urban, might have a significantly higher bacterial abundance as compared to the other three villages, which are considered rural (Figure 2A). These intriguing insights were confirmed by comparisons of the alpha diversity values. With regard to pairwise comparisons of bacterial alpha diversity, no statistically significant differences were observed between villages B, C, and D. In contrast, pairwise comparisons of the bacterial alpha diversity values from villages A identified significant differences with the alpha diversity values from villages B (p < 0.001), C (p < 0.001), and D (p < 0.001) (Figure 2C and Table 1). We conducted PCoA analysis using pairwise comparisons of Bray–Curtis distances between samples. This analysis confirmed the significant bacterial dissimilarity between urban soil and rural soil (Figure 2E). Furthermore, at the phylum taxonomic level, a noticeable increase in species from Acidobacteria-6 and decreases in Actinobacteria and Ktedonobacteria species were observed in urban soil (Figure 2F).

The most prevalent fungal genera in the samples were Talaromyces, Knufia, Fusarium, Chaetomium, Penicillium, Aspergillus, Paracamarosporium, Trichoderma, and Chaetomella. In terms of fungal alpha diversity, there were no statistically significant differences observed among the villages collectively (p = 0.787) (Figures 2B,D). The results of PCoA analysis revealed only subtle differences between urban soil and rural soil (Supplementary Figure S1A), indicating a generally similar fungal community structure between these two soil environments.

Notable patterns emerged among pathways represented in rural soil compared to urban soil, with nicotine degradation II (pyrrolidine pathway) (p < 0.001), toluene degradation I (aerobic) (via o-cresol) (p < 0.001), toluene degradation II (aerobic) (via 4-methylcatechol) (p < 0.001), nylon-6 oligomer degradation (p = 0.007), aromatic compound degradation via β-ketoadipate (p < 0.001), and 4-coumarate degradation (anaerobic) (p < 0.001) exhibiting significant increases. Conversely, urban soil exhibited distinctive shifts compared to rural soil, with formaldehyde assimilation II (RuMP Cycle) (p = 0.016), formaldehyde oxidation I (p = 0.017), toluene degradation V (aerobic) (via toluene-cis-diol) (p = 0.002), and anaerobic aromatic compound degradation (Thauera aromatica) (p < 0.001) experiencing a marked increase in relative abundance (Supplementary Table S11).

In Figures 3A,B, we quantified highly abundant bacteria and fungi (with average relative abundances exceeding 0.2%) using maximum likelihood trees to depict the bacterial and fungal community structures and phylogenetic relationships within the rhizospheric and peripheral soils. This approach provided insights into the community composition and evolutionary connections of bacteria and fungi in these distinct soil environments.

Within the rhizospheric zone, noticeable increases were observed in six bacterial taxa, including Xanthomonadales, Bacilli, and Actinobacteria (Figures 4A,C). Similarly, 27 fungal taxa exhibited heightened representation, encompassing Agaricomycetes, Sordariomycetes, and Eurotiomycetes (Figure 4B). To further elucidate the effects of urban environmental factors and citrus root systems on soil microbiota, we classified samples simultaneously based on different sampling locations and soil types (Supplementary Figure S2). In villages B, C, and D, compared to peripheral soil, the alpha and beta diversity of rhizospheric soil exhibited greater fluctuations, the trend observed in both bacteria and fungi. Specifically, we observed that rhizospheric alpha diversity from the same village had larger interquartile ranges in boxplots compared to peripheral soil. Correspondingly, in the PCA plots of beta diversity, the confidence ellipsoids of rhizospheric samples from the same locations were more inclined to encompass those of peripheral soil. Conversely, in village A, the boxplots of alpha diversity for both bacteria and fungi indicated larger interquartile ranges in peripheral soil. Simultaneously, in the PCA plots, the confidence ellipsoids of peripheral soil and rhizosphere samples exhibited roughly equal intersection. These findings suggest that citrus root systems do indeed play a role in shaping soil microbial community structure, albeit weaker than the influence of urban environmental effects.

Furthermore, we contrasted bacterial functional disparities between citrus orchard rhizospheric and peripheral soils (including both urban and rural), pinpointing a significant reduction in the CMP-pseudaminate biosynthesis pathway in the rhizospheric soil (p = 0.001) (Figure 5A). Importantly, the divergence of this pathway was unaffected by urban factors, as evidenced by its substantial reduction both urban rhizospheric soil as compared with urban peripheral soil (p = 0.020, Figure 5B) and in rural rhizospheric soil as compared with rural peripheral soil (p = 0.014, Figure 5C). The outcomes of a LEfSe indicated the potential prominence of Xanthomonadaceae within the rhizospheric zone (p = 0.007) (Figure 6). While functional differences were identified in bacterial communities in rhizospheric and peripheral soils, no statistically significant differences between the two soil types were identified using phenotype prediction across all of the following nine phenotype traits: aerobic, anaerobic, mobile element-containing, biofilm-forming, Gram-negative, Gram-positive, potentially pathogenic, and stress-tolerant (Supplementary Figures S3, S4).

The results from the FUNGuild analysis, which categorizes fungal communities based on trophic mode, traits, and growth forms, reveal noticeable differences between the fungal communities inhabiting the rhizospheric and peripheral zones of the soil. Specifically, in peripheral soils, the trophic mode classification identified a comparatively heightened prevalence of pathotrophs (p = 0.042), coupled with an increased occurrence of plant pathogens as per the Guild classification (p = 0.011). Additionally, analysis of the growth form classification demonstrated significantly elevated proportions of smut (p = 0.026) and corticioid (p = 0.020) (Table 2).

Ascomycota exhibited a substantial dominance within the citrus orchard soils. This pronounced presence of a specific taxon likely exerted a multifaceted impact on both the ecological landscape of citrus orchard soils and the agricultural productivity of the area. For example, the elevated prevalence of beneficial Trichoderma spp. (0.3%) would suppress plant pathogens. The principal mechanisms not only include mycoparasitism, antibiosis, and competition for resources and space but also involve the induction of resistance pathways in plants, leading to increased plant growth and nutrient uptake. Additionally, these fungi produce a diverse array of antifungal enzymes, such as chitinases and beta-1,3 glucanases. These enzymes exhibit synergistic effects with each other and with other antifungal enzymes and materials (Harman, 2006; Bae et al., 2011). Finally, rhizosphere-competent T. harzianum strain, such as T22, have been shown to promote root growth in a variety of plants (Harman, 2000). Conversely, the heightened abundance of potential pathogenic fungi, such as Fusarium spp. (0.4%), might increase the risk of disorders like root rot, consequently impeding plant growth and leading to diminished yield (Buttar et al., 2024; Kamble et al., 2024).

In terms of fruit quality, the prominent occurrence of Penicillium (0.4%) and Aspergillus spp. (0.4%). might lead to surface mold formation on citrus fruits, thereby influencing their visual appearance and gastronomic value. Moreover, certain strains within the Penicillium spp. could generate deleterious toxins, thereby posing latent risks to human health. Additionally, the augmented abundance of fungi similar to Chaetomium (0.4%) and Chaetomella spp. (0.3%), in turn, could substantively contribute to organic matter decomposition.

In the realm of farm management, a comprehensive understanding of the ecological roles of each species remains paramount. This comprehension is instrumental in devising judicious soil management strategies and disease control measures, designed to harness the proactive potential of beneficial species while mitigating the detrimental effects of harmful ones.

In the realm of diversity analysis, distinct differences were evident in the alpha diversity of bacterial communities between the citrus orchard soil of village A and those of villages B (p < 0.001), C (p < 0.001), and D (p < 0.001). The substantial disparity in species diversity became apparent through analyses of rarefaction curve plots, beta diversity plots, and taxonomic classifications at the phylum level (Figure 1). Considering the alignment of sampling design and real-world environmental observations, we inferred that the proximity to urban transportation zones played a pivotal role in driving differences in soil bacterial diversity between urban soil and rural soils. This inference was corroborated by the results of the functional disparity analysis of bacterial communities between urban soil and rural soil.

In the urban soil, the significant increase in the relative abundance of Formaldehyde assimilation II (RuMP Cycle) and formaldehyde oxidation I possibly stems from pollution sources, including vehicular emissions near urban areas, prompting bacterial communities to adapt metabolically to formaldehyde (Klein et al., 2022). Similarly, the notable increase in the relative abundance of toluene degradation V (aerobic) (via toluene-cis-diol) might be attributed to the prevalence volatile organic compounds, such as toluene, in urban environments, leading to metabolic adaptation among bacterial communities (Gelencsér et al., 1997). The significant rise in the relative abundance of anaerobic aromatic compound degradation (Thauera aromatica) is consistent with the presence of substantial aromatic compounds in the urban environment and the adaptation of bacterial communities to degrade these compounds (Keyte et al., 2016; Zhang et al., 2022).

The urban orchards were closely situated alongside a motorized road, resulting in heightened vehicular activity during routine periods. Research on the impact of highway-related activities on soil microorganisms has shown that these activities exacerbate vegetation degradation along the road, significantly alter soil physicochemical properties, and cause heavy metal pollution, thereby affecting the diversity and community structure of soil bacteria. Furthermore, this disturbance exhibits a gradual increase in intensity with proximity to the highway (Liu Z. et al., 2021). In comparison of bacterial species in village A (urban) and the other villages (rural), Bacillus muralis, Bacillus fumarioli, and Bacillus flexus all showed significant increases in urban. Bacillus spp. possess heavy metal tolerance due to the presence of chromium reductase and nitroreductase enzymes, which are involved in the chromium (VI) reduction pathway (Luo et al., 2022; Shahraki et al., 2022; Ramli et al., 2023). This finding suggests that the changes in soil microbial community structure observed are influenced by traffic factors, consistent with the environmental factors and research results we have observed. It also suggests that in the study area, the application of B. muralis, B. fumarioli, or B. flexus for soil microbial remediation could be considered (Jebeli et al., 2017; Kulkova et al., 2023; Schommer et al., 2023; Yu et al., 2023).

Litter can increase the network complexity of soil microbial communities (Feng et al., 2022). Changes in litter and roots can affect soil microbial activities and nutrient cycling (Chen et al., 2022; Zhu et al., 2022). Potentially important physical differences between rural and urban areas were noticed during soil collection. We observed dispersed waste materials in rural orchards. These materials were primarily composed of discarded cigarette remnants, abandoned bags, garments, and packaging materials (Supplementary Figure S5), whereas similar materials were not observed in the urban orchards. Correspondingly, there was a significant increase in the relative abundance of bacteria involved in nicotine degradation II (pyrrolidine pathway), indicating a potential bacterial metabolic adaptation to nicotine degradation. The significant elevation in the relative abundance of nylon-6 oligomer degradation might be linked to waste materials prevalent in rural environments.

Fruit tree root system can affect soil microbial community structure, but it was not the main factor. According to our analysis, we observed a significant difference (p = 0.001) in the CMP-pseudaminate biosynthesis pathway between the rhizospheric soil and the peripheral soil (Figure 5). Pseudaminic acid, a sialic acid-like sugar, is involved in the synthesis of bacterial cell surface glycans, which in turn affects pathogenicity and adaptability through bacterial adhesion, invasion, and immune evasion mechanisms (Knirel et al., 2003; Hsu et al., 2006). The reduction in the CMP-pseudaminate biosynthesis pathway may represent a plant-microbe interaction, where plant roots release certain compounds that influence bacterial metabolism and biosynthetic pathways in the rhizospheric soil, contributing to a protective effect. However, it is important to note that our study also revealed a relatively high false discovery rate (FDR = 0.415, Supplementary Table S10), which may indicate a certain degree of false positives in multiple comparisons. Further research is needed to validate these findings and explore the biological significance of the observed differences in more depth.

These variations in bacterial functional pathways support our hypotheses regarding the impact of urban environmental factors on citrus orchard soil. The exogenous factors such as urban traffic and rural litter did have a significant impact on the microbiome in citrus orchards. These factors played a significant role in the substantial differences in bacterial diversity, potentially further influencing citrus tree growth and fruit yield.

By elucidating the microbial dynamics in citrus orchard soil during the crucial fruit maturation period, our findings can inform management strategies aimed at enhancing citrus yield and fruit quality. The identification of beneficial soil microorganisms capable of promoting citrus growth and disease resistance can lead to the development of targeted microbial-based biofertilizers and biocontrol strategies (Scott et al., 2006; Das et al., 2022; Wang et al., 2023), reducing reliance on chemical inputs and fostering long-term sustainability (Ali et al., 2019; Li et al., 2019; Zhou et al., 2020).

We propose the implementation of physical barriers in rural areas to reduce the quantity of pollutants such as cigarette butts and nylon. Urban orchards in city areas should be located far away from densely populated regions with high traffic density (Hagler et al., 2012; Ranasinghe et al., 2019). However, considering the current pollution situation, microbial measures can also be employed for improvement. For example, the Trichoderma harzianum T22 strain can be cultured on potato dextrose agar medium, and the spores can be collected for seed coating (Harman and Shoresh, 2007; Alinc et al., 2021). Alternatively, soaking the seeds in a suspension containing Bacillus spp. for 2 h can help improve plant growth and crop productivity under various abiotic stresses, including heavy metals and drought (Anbuganesan et al., 2024). In addition, synthesizing a biodegradable carbon nanoparticle from bacterial biofilm is also an option (Li et al., 2024; Savadiya et al., 2024). A more direct approach is to consider inserting genes encoding antifungal protein internal chitinase or external chitinase into the citrus genome alone or in combination (Bolar et al., 2000, 2001).

By utilizing amplicon sequencing as a non-invasive and cost-effective method, we were able to gain insights into the microbial diversity and functional potential in citrus orchard soil at an unprecedented scale. Our study confirmed the beneficial application of metagenomics in addressing agricultural and environmental concerns. This included the identification of two pollutants that indeed affected microbial community functionality and proposed specific microbial methods for improving land heavy metal pollution and suppressing plant pathogens based on real-world microbial community structures. Consequently, future research could continue to utilize amplicon sequencing across different seasons and larger geographical scales, and shotgun sequencing methods could be employed to further validate the functionality of these microbes.

While our study has provided valuable insights into the microbial communities of citrus orchard soils, several limitations need to be considered. Firstly, our research solely relied on functional prediction methods to examine microbial community functional changes, which resulted in a lack of sufficient validation and interpretability. Future studies should incorporate more comprehensive functional analyses, such as employing shotgun sequencing methods to investigate specific differentially expressed genes. Moreover, our study did not include intervention studies using specific species, which could have strengthened the credibility of the proposed measures for soil health and agricultural productivity. Additionally, although we consider the fruit maturation stage as a crucial phase for investigation, it is essential to study soil in other growth periods to ensure the stability assessment of our findings. To maximize the comparability of sampling points, we applied specific restrictions, focusing on citrus orchards situated near rivers and lacking human management. Consequently, in urban settings, other suitable family farms meeting these criteria were not found. Further urban sampling and exploration across other soil environment categories are essential to enhance the universality of our conclusions. Finally, our study lacked a comprehensive assessment of environmental factors, which may have led to the omission of certain variables. Addressing these limitations in future research could further enhance our understanding of microbial communities in citrus orchard soils.