Abstract
Low-disturbance wild tea plantations provide useful field systems for examining how natural habitat heterogeneity is associated with rhizosphere soil conditions and microbial community organization. In this study, rhizosphere soils were collected from three naturally contrasting habitats in a wild tea plantation in Hainan, China: open sites, slope-shaded sites, and rainforest-shaded sites. Soil physicochemical properties were measured, and shotgun metagenomic sequencing was used to characterize microbial taxonomic composition and KEGG Ortholog-based functional potential. Clear habitat-associated edaphic differentiation was observed. Rainforest-shaded soils had relatively higher pH, organic matter, available phosphorus, and available potassium than the other habitats, whereas nitrogen-related indices remained comparatively stable. Microbial alpha diversity did not differ significantly among habitats, but Bray-Curtis-based beta diversity showed a significant overall habitat effect, with the principal compositional gradient involving RS relative to OP and SS. This pattern indicates that habitat-associated variation was expressed more strongly through community compositional reorganization than through changes in overall richness. Redundancy analysis further indicated that pH, organic matter, and available phosphorus were the main edaphic variables associated with microbial community differentiation. The most abundant KOs showed habitat-associated patterns in relative representation. Open-site soils showed higher representation of metagenome-inferred KOs annotated to substrate transport and environmental sensing, slope-shaded soils showed higher representation of KOs annotated to redox balance and stress-related processes, and rainforest-shaded soils showed higher representation of KOs annotated to genetic information processing and cellular maintenance. These patterns reflect differences in encoded functional potential rather than realized microbial activity. Overall, these results suggest that natural habitat heterogeneity in low-disturbance wild tea systems is associated with coordinated variation in rhizosphere soil fertility and microbial community composition, together with descriptive variation in selected components of metagenome-inferred functional potential. These findings provide field-based evidence for aboveground–belowground linkage in wild tea soils and highlight the importance of edaphic gradients in structuring rhizosphere microbiomes.
Introduction
In tea (Camellia sinensis) plantation ecosystems, rhizosphere microbial communities play important roles in mediating interactions among perennial tea roots, soil nutrients, and local environmental conditions, thereby contributing to organic matter turnover, nutrient cycling, soil aggregation, and plant health (; ). Previous studies in managed tea plantations have shown that fertilization, tillage, mulching, and organic management can substantially alter soil microbial communities and associated soil functions (; ). Recent evidence from continuous-cropping systems has further shown that long-term habitat and management conditions can reshape rhizosphere microbial communities together with soil physicochemical properties (). Collectively, these findings highlight the sensitivity of plant-associated microbiomes to environmental context. However, compared with management-driven changes, less attention has been paid to how naturally occurring habitat heterogeneity is associated with rhizosphere soil conditions and microbial community organization in low-disturbance tea systems.
Low-disturbance wild tea plantations provide useful field systems for examining plant-soil-microbiome relationships under relatively natural conditions. In such systems, variation in canopy cover, understory vegetation, litter input, and topographic position can jointly create contrasting habitat conditions within the same landscape. These habitat differences may influence soil microclimate, organic matter accumulation, nutrient retention, and rhizosphere resource availability, thereby providing distinct edaphic contexts for microbial community assembly (; ). Unlike controlled agronomic treatments, these naturally occurring habitats represent integrated field conditions in which multiple environmental factors act simultaneously. Therefore, they offer an opportunity to explore how aboveground habitat complexity is linked to belowground soil and microbial variation.
Soil pH, organic matter, and nutrient availability are widely recognized as important environmental factors associated with soil microbial community structure. In tea plantation soils, microbial communities are often sensitive to soil acidification, organic inputs, and nutrient status (). More broadly, habitat heterogeneity has been shown to contribute to spatial variation in belowground communities across forest and agricultural ecosystems (; ). Nevertheless, most studies on tea systems have focused on intensively managed plantations or specific agronomic interventions, whereas field-based evidence from minimally managed wild tea plantations remains limited. In particular, few studies have simultaneously examined rhizosphere soil fertility, microbial taxonomic composition, and metagenome-inferred functional potential across naturally contrasting habitats.
In this study, we investigated rhizosphere soils from three naturally contrasting habitats within a low-disturbance wild tea plantation in Hainan, China: open sites, slope-shaded sites, and rainforest-shaded sites. These habitats differed in canopy structure, vegetation cover, and topographic setting within the same plantation landscape, and were treated as field-based ecological categories rather than controlled shading treatments. Soil physicochemical properties were measured, and shotgun metagenomic sequencing was used to characterize microbial taxonomic composition and KEGG Ortholog-based functional potential. The objectives of this study were to: (i) determine whether natural habitat heterogeneity is associated with distinct rhizosphere soil physicochemical patterns; (ii) assess whether habitat-associated variation is expressed more strongly in microbial community composition than in alpha diversity; and (iii) examine whether taxonomic differentiation is accompanied by differences in metagenome-inferred functional potential. This study provides a field-based analysis of edaphic gradients and rhizosphere microbiome differentiation in a low-disturbance wild tea system, with implications for understanding aboveground–belowground linkages in soil ecological processes.
Materials and methods
Study site and habitat design
The study was conducted in a naturally regenerated wild tea plantation within Hainan Tropical Rainforest National Park, Wuzhishan City, Hainan Province, China (109°40′E, 18°51′N). The region has a tropical monsoon climate, with a mean annual temperature of approximately 25 °C and annual precipitation ranging from 1,300 to 1,800 mm. The soil at the study site is locally classified as lateritic red soil according to the regional soil classification and soil survey information available for Hainan Island (; ) and is typically acidic under tropical conditions. The plantation has remained unmanaged for more than 20 years, with no fertilization, tillage, or pesticide application during this period.
Three naturally contrasting habitat types were identified within the same wild tea plantation landscape: open sites (OP), slope-shaded sites (SS), and rainforest-shaded sites (RS). Open sites were characterized by relatively low canopy cover; slope-shaded sites were located on shaded slopes with moderate vegetation cover; and rainforest-shaded sites occurred beneath closed rainforest canopy. These habitats differed in canopy structure, understory vegetation, and topographic setting, and were therefore treated as naturally occurring ecological categories rather than controlled shading treatments.
Within each habitat type, three replicate plots of 20 m × 20 m were established. Plots belonging to the same habitat type were separated by at least 50 m to improve spatial independence. All plots were located within the same plantation landscape and shared a comparable long-term management history and soil type.
Rhizosphere soil sampling
Within each plot, seven healthy tea plants were randomly selected along an S-shaped sampling path. Fine roots were excavated from the 0–20 cm soil layer, and soil tightly adhering to the roots was collected by gentle brushing and defined as rhizosphere soil. Rhizosphere soils from the seven plants within the same plot were pooled and thoroughly homogenized to generate one plot-level composite sample. Each composite sample was treated as one biological replicate in subsequent analyses. All samples were immediately transported to the laboratory on ice. Subsamples for DNA extraction were stored at −80 °C, whereas subsamples for soil physicochemical analyses were air-dried and sieved before measurement.
Soil physicochemical analyses
Air-dried soils were sieved (2 mm) prior to analysis. Soil pH was measured in a 1:2.5 (w/v) soil-to-deionized water suspension using a calibrated PHS-3C pH meter (INESA Scientific Instrument Co., Ltd., Shanghai, China) following . Soil organic matter (OM) was determined by potassium dichromate oxidation using the Walkley-Black method (). Total nitrogen (TN) was determined by Kjeldahl method (), and alkali-hydrolyzable nitrogen (AN) by the alkali diffusion method (). Total phosphorus (TP) and total potassium (TK) were determined after H2SO4-HClO4 digestion (). Available phosphorus (AP) was extracted using the Olsen method and quantified by molybdenum-antimony colorimetry (; ). Available potassium (AK) was extracted with 1 M NH4OAC and measured using a flame photometer (). CEC (cation exchange capacity) was determined using the [Co(NH3)6]Cl3 extraction-spectrophotometry method according to Chinese national environmental standard HJ 889-2017. Humus carbon (HC) was determined following humus fraction extraction according to Chinese forestry standard LY/T 1238-1999.
Soil DNA extraction, metagenomic sequencing and bioinformatic analysis
Genomic DNA was extracted from rhizosphere soil samples using the E.Z.N.A.® Soil DNA Kit (Omega Bio-tek, Norcross, GA, USA) following the manufacturer’s instructions. DNA quality and concentration were assessed by agarose gel electrophoresis, and DNA concentration and purity were measured using a OneDrop OD-1000+ microvolume UV–Vis spectrophotometer (OneDrop, Shanghai, China). Qualified DNA samples were used for library construction.
DNA was fragmented using a Covaris system, and paired-end libraries were constructed through end repair, A-tailing, adaptor ligation, size selection, and PCR enrichment. The libraries were sequenced on an Illumina NovaSeq 6,000 platform using a paired-end 150 bp strategy.
Raw shotgun metagenomic reads were filtered to remove adaptor sequences and low-quality reads. After quality filtering, 16.14–20.52 million clean read pairs were retained per sample, with an average of 18.16 million read pairs per sample, corresponding to 4.84–6.15 Gb of clean sequence data. Across all nine samples, 163.44 million clean read pairs and 49.03 Gb of clean sequence data were obtained. Q30 values ranged from 95.35% to 96.69%. The resulting quality-filtered reads were de novo assembled using MEGAHIT v1.2.9. Gene prediction was performed using Prodigal v2.6.3.
A non-redundant gene catalog was constructed using CD-HIT v4.8.1 with the following parameters: identity = 0.95, aS = 0.90, G = 1, r = 1, and g = 1. The longest sequence in each cluster was retained as the representative sequence. Representative sequences were indexed using BWA v0.7.17, and clean reads were mapped back to the non-redundant gene catalog using BWA-MEM for gene abundance estimation. Alignments were filtered to retain reads with sequence identity ≥95%, and genes supported by two or fewer reads were removed. Gene abundance was summarized as read counts, and transcripts per million (TPM) values were calculated for length- and sequencing-depth-normalized comparisons when needed.
Taxonomic annotation and abundance profiling were performed using Kraken2 v2.1.2 with a confidence threshold of 0.2. The Kraken2 PlusPF database was used as the reference database, and taxonomic lineage information was standardized using TaxonKit. Taxon abundance was estimated by summing gene-level read counts assigned to each taxon and converting them into relative abundance at different taxonomic ranks.
Functional annotation was performed against the KEGG database using DIAMOND v2.0.6 in blastp mode with an e-value cutoff of 1 × 10−5. The best hit was retained for KEGG Ortholog (KO) assignment, and KO abundance was calculated by summing the read counts of genes assigned to each KO. KO abundance profiles were normalized for downstream comparative analyses.
Data processing and analysis
Statistical analyses were performed using SPSS 26.0 and R (v4.4.2.) Data are presented as mean ± standard error (n = 3). Statistical significance was determined at P ≤ 0.05 unless otherwise stated.
Differences in soil physicochemical properties and microbial alpha-diversity indices among habitats were tested using the Kruskal–Wallis test, followed by Dunn’s post hoc test with Benjamini–Hochberg adjustment for multiple comparisons. Alpha-diversity indices were calculated from taxonomic abundance profiles.
Microbial community composition was evaluated at both the phylum and genus levels. Bray-Curtis dissimilarity matrices were calculated based on relative abundance data, and principal coordinate analysis (PCoA) was used to visualize differences in community composition among habitats. Permutational multivariate analysis of variance (PERMANOVA) with 999 permutations was used to test for overall differences in community composition among habitat types. Pairwise post hoc PERMANOVA comparisons were subsequently performed for OP versus SS, OP versus RS, and SS versus RS at both taxonomic levels. For each taxonomic level, the three pairwise P values were adjusted using the Benjamini–Hochberg procedure.
Redundancy analysis (RDA) was used to examine relationships between soil physicochemical variables and microbial community composition. RDA was performed separately using the relative-abundance matrices of the top 10 phyla and the top 10 genera. Prior to RDA, microbial community abundance data were Hellinger-transformed. Environmental variables were log-transformed and subsequently standardized using Z-score scaling, with each variable centered to a mean of zero and scaled to a standard deviation of one. Multicollinearity among environmental variables was evaluated using variance inflation factors. The variable with the highest VIF was removed iteratively until all retained variables had VIF values ≤10. Following VIF screening, pH, available phosphorus (AP), soil organic matter (OM), and cation exchange capacity (CEC) were retained in the final RDA models. Ordination analyses and visualizations were performed using the vegan and ggplot2 packages in R.
For functional analysis, KO abundance profiles were compared among habitats. Overall differences in KO relative abundance among the three habitats were assessed using the Kruskal–Wallis test, with Benjamini–Hochberg correction for multiple testing and epsilon-squared (ε2) reported as the overall effect size. Pairwise comparisons were performed using Dunn’s test, with log2 fold changes used to describe the direction and magnitude of pairwise differences. For each KO, the P values from the three pairwise comparisons were adjusted using the Benjamini–Hochberg procedure. The corresponding top-ranked KO results are provided in Supplementary Table S1. Heatmaps were generated in R using the pheatmap package without row-wise scaling. The 35 KOs with the highest mean relative abundances across all nine samples were displayed in the heatmap. All KO-by-habitat cells were mapped to a common color scale based on the minimum and maximum mean relative-abundance values across the complete displayed matrix. Final figure arrangement was performed in Adobe Illustrator.
Results
Habitat-associated edaphic differentiation in rhizosphere soils
Rhizosphere soil physicochemical properties differed among OP, SS, and RS (Figure 1). All soils were acidic, but RS showed significantly higher pH than OP and SS (P ≤ 0.001), indicating relatively lower soil acidity.
FIGURE 1
Habitat-associated differences were also evident for OM and mineral nutrient availability. OM was significantly higher in RS than in OP (P ≤ 0.01). Potassium-related variables showed a similar pattern, with RS having significantly higher TK (P ≤ 0.001) and AK (P ≤ 0.05). Phosphorus availability also varied among habitats: TP differed moderately (P ≤ 0.05), whereas AP showed a stronger habitat-associated difference (P ≤ 0.01), with the highest values in RS and the lowest values in SS.
In contrast, nitrogen-related indices were comparatively stable across habitats. Neither TN nor AN differed significantly among OP, SS, and RS. Similarly, CEC and HC showed no significant habitat-associated differences. Overall, these results indicate that natural habitat heterogeneity within the wild tea plantation was mainly associated with variation in soil acidity, OM accumulation, and P and K availability, rather than with broad changes in N pools or exchange-related properties.
Microbial community composition and functional potential
Microbial taxonomic composition varied among OP, SS, and RS, although dominant phyla were generally shared across habitats (Figure 2). At the phylum level, Pseudomonadota was the most abundant phylum across all three habitats, but its relative abundance was lower in RS than in OP and SS. In contrast, Acidobacteriota showed higher relative abundance in RS. Other phyla, including Verrucomicrobiota and Planctomycetota, were consistently detected across habitats but occurred at comparatively lower relative abundances.
FIGURE 2
At the genus level, Bradyrhizobium was the dominant genus across all habitats, but its relative abundance was lower in RS than in OP and SS. Meanwhile, the aggregated “Other” category accounted for a larger proportion in RS, indicating that genera outside the top 10 collectively contributed a larger share of the community in RS. Overall, these results indicate that habitat-associated differences in the rhizosphere microbiome were expressed mainly through shifts in relative abundance structure rather than through complete replacement of dominant taxa.
The 35 most abundant KOs showed habitat-associated variation in relative abundance patterns (Figure 3). OP showed higher representation of several KOs related to substrate transport and resource acquisition. SS showed higher representation of several KOs associated with transport, environmental response, and redox-related processes. In contrast, RS showed higher representation of several KOs related to transcription, genetic information processing, and cellular maintenance. These patterns represent descriptive variation in selected components of metagenome-inferred functional potential and should not be interpreted as direct evidence of gene expression, enzyme activity, or realized in situ microbial functions. Separate KO-level statistical comparisons are summarized in Supplementary Table S1, which provides the top-ranked KO results together with their functional annotations, mean relative abundances across habitats, overall epsilon-squared effect sizes, pairwise log2 fold changes, and adjusted P values. No KO remained significant after global FDR correction.
FIGURE 3
Habitat-associated compositional differentiation exceeded alpha-diversity changes
Alpha-diversity indices were generally stable among OP, SS, and RS (Figures 4A–H). Shannon, Simpson, richness-related indices, and ACE showed no significant differences among habitats (P > 0.05). Although Chao1 showed a marginal trend (P = 0.051), the overall pattern indicated that habitat-associated variation did not result in marked changes in within-sample microbial diversity. Good’s coverage exceeded 96% in all samples, indicating high sampling completeness.
FIGURE 4
In contrast, beta-diversity analysis revealed significant overall habitat-associated differentiation in microbial community composition. PERMANOVA based on Bray-Curtis dissimilarities showed significant overall differences among habitats at both the phylum level (F = 7.206, R2 = 0.706, P = 0.032) and the genus level (F = 3.864, R2 = 0.563, P = 0.039) (Figures 4I,J). PCoA and pairwise effect-size patterns consistently indicated that the principal compositional gradient involved RS relative to OP and SS. At the phylum level, OP versus RS and SS versus RS had larger effect sizes than OP versus SS. At the genus level, the largest effect size was observed for SS versus RS, followed by OP versus RS, whereas OP and SS were comparatively similar. However, none of the pairwise comparisons was significant after Benjamini–Hochberg adjustment (Supplementary Table S2).
Together, these results indicate that habitat-associated differences in the rhizosphere microbiome were expressed more strongly through overall community compositional reorganization than through changes in within-sample diversity.
Edaphic correlates of microbial community differentiation
RDA was used to examine the relationships between soil physicochemical variables and microbial community composition. At the phylum level, the first two canonical axes explained 95.73% of the constrained variation, with RDA1 and RDA2 accounting for 78.77% and 16.96%, respectively (Figure 5A). RS samples were positioned mainly on the positive side of RDA1, broadly corresponding to the directions of pH, AP, and OM. OP samples occurred mainly on the negative side of RDA1 and were positioned closer to the direction of CEC, whereas SS samples occupied primarily the lower portion of the ordination space.
FIGURE 5
A similar pattern was observed at the genus level (Figure 5B). RDA1 and RDA2 explained 89.95% and 8.67% of the constrained variation, respectively, accounting for 98.62% in total. RS samples were positioned predominantly on the positive side of RDA1, broadly corresponding to the directions of pH and AP and, to a lesser extent, OM. OP and SS samples occurred mainly on the negative side of RDA1, although they differed in their positions along RDA2.
Together, these ordination patterns indicate that the principal habitat-associated gradient in microbial taxonomic composition, particularly the contrast involving RS relative to OP and SS, was associated with variation in the retained edaphic variables. Among these variables, pH, AP, and OM showed prominent directional correspondence with the RS side of the ordination, while CEC showed a different directional pattern. These relationships represent ecological associations and should not be interpreted as evidence of direct causal effects.
Discussion
Habitat-associated edaphic differentiation
The observed differences in rhizosphere soil properties suggest that contrasting natural habitats within the wild tea plantation were associated with distinct edaphic conditions. RS generally showed the highest mean pH, OM, AP, and AK among the three habitats, indicating that the closed-canopy habitat was associated with lower soil acidity and greater nutrient availability. Because OP, SS, and RS differed not only in canopy cover but also in vegetation structure and topographic setting, these patterns should not be interpreted as the effect of shading alone. The vegetation context also differed among the habitats. OP represented a relatively open vegetation structure with low canopy cover, SS represented an intermediate vegetation context on shaded slopes, and RS occurred beneath a closed rainforest canopy with comparatively denser surrounding vegetation. Such structural differences may influence the quantity and composition of litter and root-derived inputs, while also modifying soil temperature, moisture, and nutrient retention. These vegetation-associated pathways could have contributed to the distinct edaphic conditions and microbial composition observed in RS and to the comparatively transitional characteristics of SS. However, vegetation composition and litter inputs were not quantified, and these relationships should therefore be regarded as plausible habitat-level explanations rather than direct evidence of vegetation-driven effects. Instead, they more likely reflect the combined influence of habitat-associated environmental conditions under field settings.
The higher OM in RS may be associated with its closed-canopy vegetation structure, potentially greater litter and root-derived inputs, and more stable microclimatic conditions. Closed or semi-closed canopy conditions may buffer soil temperature and moisture fluctuations, while increased litter and root-derived carbon inputs may contribute to OM accumulation and nutrient retention (; ). Compared with the more exposed OP and transitional SS habitats, the relatively stable soil environment in RS may favor OM accumulation and nutrient retention (; ). This interpretation is consistent with previous studies showing that litter inputs and vegetation structure can influence soil carbon dynamics and microbial processes in terrestrial ecosystems (; ).
The higher AP and AK in RS may also be linked to the OM-rich soil environment. OM decomposition can release low-molecular-weight organic acids and other compounds that promote P mobilization through protonation, chelation, and competitive desorption of Fe/Al-bound phosphate (; ; ). Similarly, higher OM may improve soil buffering capacity and contribute to K retention, thereby increasing AK availability and reducing leaching risk (; ; ; ). However, because litter input, soil moisture, enzyme activity, and nutrient transformation rates were not directly measured in this study, these explanations should be regarded as plausible ecological interpretations rather than confirmed mechanisms.
Compared with RS, the intermediate nutrient status of SS may reflect its transitional habitat characteristics. In addition to moderate canopy cover, slope position may influence the redistribution of soil particles, OM, and labile nutrients through erosion and deposition processes (; ; ). Such topographic effects may partly contribute to the intermediate edaphic pattern observed in SS. The lower AP observed in SS, including a significant difference relative to RS, may be associated with slope-mediated losses and redistribution of labile phosphorus through surface runoff, soil erosion, and downslope transport. In addition, the relatively lower OM content in SS may have provided less capacity for phosphorus retention and mobilization than in RS. However, because runoff, erosion intensity, phosphorus fractions, and phosphorus transformation processes were not directly measured, these explanations remain tentative.
In contrast to P and K availability, TN and AN remained relatively stable among habitats. This may indicate that bulk N pools were less responsive to habitat-associated variation at the spatial scale examined here. Because TN and AN are pool-size indicators measured at a single sampling time, they may not fully capture dynamic nitrogen transformation processes, such as mineralization, nitrification, immobilization, and microbial nitrogen turnover. Habitat-associated differences in N cycling may therefore be expressed more strongly through process rates, microbial functional genes, or seasonal dynamics rather than through bulk pool size alone (; ). Therefore, the absence of significant differences in TN and AN does not necessarily imply that N cycling was unchanged among habitats.
Overall, these results indicate that natural habitat heterogeneity was associated with coordinated variation in soil acidity, OM accumulation, and P and K availability. Such edaphic differentiation provides an important environmental context for interpreting the microbial community patterns observed in this study.
Stable alpha diversity but habitat-associated compositional reorganization
Microbial alpha diversity did not differ significantly among OP, SS, and RS, indicating broadly comparable levels of within-sample richness and evenness across the three habitats. In contrast, overall community composition differed significantly among habitats at both the phylum and genus levels. The PCoA and pairwise effect-size patterns further indicated that the principal compositional gradient involved RS relative to OP and SS, whereas OP and SS were comparatively similar.
The contrast between stable alpha diversity and a significant overall habitat effect on community composition is ecologically meaningful. Similar levels of richness and evenness can be maintained across habitats even when the identities and relative abundances of community members differ. Thus, habitat-associated microbiome differentiation in this wild tea system was expressed primarily through compositional reorganization rather than through a net gain or loss of microbial diversity.
This pattern is consistent with habitat-associated environmental filtering, whereby differences in local edaphic conditions may favor different microbial taxa without necessarily altering overall community richness or evenness. At the same time, functional redundancy may allow taxonomically distinct communities to retain overlapping functional capacities. The habitat-associated abundance patterns observed among the most abundant KOs suggest possible variation in selected components of metagenome-inferred functional potential. However, because no KO remained significant after global FDR correction, these patterns should not be interpreted as evidence of statistically confirmed overall functional differentiation among habitats. Soil microbial communities may maintain similar levels of richness across habitats while differing in the identities or relative abundances of dominant and subordinate taxa. Similar patterns have been reported in soil microbial studies, where community composition often responds more sensitively to local edaphic and habitat conditions than diversity metrics alone (; ). In the present study, the genus-level ordination pattern further suggests that finer taxonomic resolution may better capture ecological differentiation associated with habitat heterogeneity.
The higher relative abundance of Acidobacteriota in RS should be interpreted cautiously. Although Acidobacteriota are often considered oligotrophic or acid-associated taxa (), this phylum contains ecologically diverse subgroups that may respond differently to pH, carbon availability, and nutrient status (; ). Therefore, the higher relative abundance of Acidobacteriota in RS may reflect subgroup-level responses to the local soil environment rather than a simple phylum-level ecological pattern. By contrast, the relatively lower abundance of Pseudomonadota in RS suggests that the responses of major phyla to habitat-associated edaphic conditions were not uniform.
At the genus level, the lower relative abundance of Bradyrhizobium and the larger “Other” category in RS indicating that genera outside the top 10 collectively contributed a larger share of the community in RS. Genera may respond more sensitively than broad phyla to local differences in substrate availability, rhizosphere conditions, and microhabitat structure (). Previous metagenomic studies have also suggested that finer taxonomic resolution can track environmental gradients more closely than broader taxonomic categories (). Therefore, the genus-level results support the interpretation that habitat-associated ecological variation within the plantation was reflected mainly in community compositional differentiation.
The lack of significant alpha-diversity differences does not imply that the microbiomes were functionally equivalent across habitats. Soil microbial communities often show functional redundancy (; ), but shifts in community composition may still be associated with differences in potential ecological roles, including substrate utilization, nutrient transformation, stress response, and cellular maintenance (; ). Accordingly, the observed beta-diversity and taxonomic differentiation provide an important compositional context for interpreting the habitat-associated KO patterns, while remaining consistent with an interpretation centered on metagenome-inferred functional potential rather than direct functional activity.
Edaphic context for functional-potential variation
RDA indicated that pH, AP, and OM were closely associated with microbial community differentiation. These variables were consistently aligned with the separation of RS from OP and SS, suggesting that soil acidity and nutrient availability provided an important edaphic context for microbiome differentiation. Soil pH is widely recognized as an important factor structuring soil microbial communities because it influences nutrient solubility, metal availability, and broader soil chemical conditions (; ; ). OM and AP are also closely related to microbial resource availability and substrate use (; ). Therefore, the association of RS with higher pH, OM, and AP may partly explain its distinct microbial community composition.
The 35 most abundant KOs displayed in Figure 3 showed habitat-associated patterns in relative abundance, providing a descriptive view of variation in major components of metagenome-inferred functional potential. In OP, the higher representation of several KOs related to transport and resource acquisition may be consistent with greater encoded potential for environmental sensing and resource acquisition under comparatively exposed or lower-fertility conditions (). In SS, several KOs associated with transport, environmental response, and redox-related processes showed relatively higher representation, which may be consistent with encoded responses to fluctuating or heterogeneous local conditions, possibly including those associated with slope position and microenvironmental variation (). In RS, the relatively higher representation of several KOs related to transcription, genetic information processing, and cellular maintenance may be associated with greater encoded potential for cellular maintenance and information processing under relatively higher nutrient availability (). These interpretations are based on descriptive abundance patterns and should therefore be regarded as possible ecological explanations rather than evidence of statistically confirmed functional enrichment.
The pattern observed in SS requires additional caution because the taxonomic analyses did not demonstrate significant pairwise separation between SS and OP. Although the overall PERMANOVA detected a significant habitat effect, none of the pairwise comparisons was significant after Benjamini–Hochberg adjustment, and the PCoA and RDA results indicated that OP and SS were comparatively similar at the phylum and genus levels. Nevertheless, a non-significant pairwise taxonomic difference does not establish that the two communities were identical. Broad phylum- and genus-level profiles may not resolve variation at the species or strain level, differences in gene repertoires among closely related taxa, or contributions from less-abundant community members. Therefore, the higher relative representation of selected KOs in SS should be interpreted as a descriptive KO-specific abundance pattern that may reflect finer-scale taxonomic or gene-content variation, rather than evidence that the overall functional profile of SS was significantly or categorically distinct from that of OP.
More broadly, shotgun metagenomic data describe functional potential encoded in community DNA rather than gene expression, enzyme activity, or actual in situ process rates. Accordingly, the observed KO patterns are interpreted as variation in the relative representation of selected encoded functions among habitats rather than direct measurements of microbial activity.
Taken together, the soil, taxonomic, and KO results suggest that habitat-associated edaphic differentiation was linked to variation in microbial community composition and selected components of encoded functional potential in this low-disturbance wild tea system. The principal taxonomic compositional gradient involved RS relative to the comparatively similar OP and SS communities, whereas the KO analysis identified habitat-associated relative-abundance patterns among selected functions. These results do not indicate that all three habitats, or specifically OP and SS, possessed completely distinct taxonomic or functional profiles. The findings support the view that aboveground habitat complexity, including canopy structure, vegetation cover, and topographic setting, can be linked to belowground microbial community patterns through associated changes in soil conditions. Future studies combining metatranscriptomics, enzyme activity assays, microbial biomass measurements, soil microclimate monitoring, and direct nutrient transformation measurements would be valuable for testing whether these metagenome-inferred patterns translate into realized ecological functions.
Conclusion
This study showed that natural habitat heterogeneity within a low-disturbance wild tea plantation was associated with distinct rhizosphere edaphic conditions and microbiome differentiation. RS was characterized by relatively higher pH, OM, AP, and AK, whereas N-related indices remained comparatively stable. Microbial alpha diversity did not differ significantly among OP, SS, and RS. Although the overall PERMANOVA indicated a significant habitat effect, the principal compositional gradient involved RS relative to the comparatively similar OP and SS communities. None of the pairwise contrasts was significant after Benjamini–Hochberg adjustment. Thus, the beta-diversity pattern should not be interpreted as a distinct three-way separation among all habitats.
Together, the soil properties, taxonomic composition, and RDA patterns indicate that pH, OM, and AP provided an important edaphic context for microbial community differentiation in this wild tea system. The KO analysis additionally revealed descriptive habitat-associated relative-abundance patterns in selected components of encoded functional potential, rather than statistically confirmed overall functional differentiation. However, because OP, SS, and RS represented naturally occurring field habitats rather than controlled treatments, these results should be interpreted as ecological associations rather than direct evidence of single-factor effects. Overall, this study provides field-based evidence that aboveground habitat complexity is linked to belowground soil and microbial patterns in low-disturbance wild tea soils, and highlights the role of edaphic gradients in structuring rhizosphere microbiomes.
Statements
Data availability statement
The datasets generated and analysed during the current study are available in the Sequence Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under BioProject accession PRJNA1452514. The remaining data and materials are included in this article.
Author contributions
YZ: Writing – original draft, Data curation, Formal analysis, Funding acquisition, Investigation, Software, Validation, Visualization. DL: Writing – review and editing, Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision. LZ: Writing – review and editing, Investigation, Validation. DZ: Writing – review and editing, Funding acquisition, Methodology. All authors contributed to the article and approved the submitted version.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This research was funded by the Hainan province science and technology special fund (ZDYF2024KJTPY031, ZDYF2026SXLH023, and ZDYF2024KJTPY056), and basic scientific research business expenses of HAAS (HAAS2025KJCX04).
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontierspartnerships.org/articles/10.3389/sjss.2026.17010/full#supplementary-material
References
1
AdhikariK.OwensP. R.AshworthA. J.SauerT. J.LibohovaZ.RichterJ. L.et al (2018). Topographic controls on soil nutrient variations in a silvopasture system. Agro. Geosci. Environ.1 (1), 1–15. 10.2134/age2018.04.0008
2
AliramayeeR.Khaledi DarvishanA.ArabkhedriM. (2019). Investigating the hydrological response and nutrient loss in rainfed lands in northeast of Iran using rainfall simulator. Agric. For.65 (2), 99–112. 10.17707/agricultforest.65.2.08
3
AlvarezA. F.GeorgellisD. (2023). Environmental adaptation and diversification of bacterial two-component systems. Curr. Opin. Microbiol.76, 102399. 10.1016/j.mib.2023.102399
4
BaoS. D. (2000). Soil Agrochemical Analysis. 3rd ed. Beijing: China Agriculture Press. (In Chinese).
5
BremnerJ. M. (1996). “Nitrogen—total,” in Methods of Soil Analysis, Part 3: Chemical Methods. Editor SparksD. L. (Madison, WI: Soil Science Society of America), 1085–1122. 10.2136/sssabookser5.3.c37
6
ChenH.MaK.LuC.FuQ.QiuY.ZhaoJ.et al (2022). Functional redundancy in soil microbial community based on metagenomics across the globe. Front. Microbiol.13, 878978. 10.3389/fmicb.2022.878978
7
ClarholmM.SkyllbergU. (2013). Translocation of metals by trees and fungi regulates pH, soil organic matter turnover and nitrogen availability in acidic forest soils. Soil Biol. biochem.63, 142–153. 10.1016/j.soilbio.2013.03.019
8
DangR.LiuJ.LichtfouseE.ZhouL.ZhouM.XiaoL. (2024). Soil microbial carbon use efficiency and the constraints. Ann. Microbiol.74, 37. 10.1186/s13213-024-01780-9
9
FageriaN. K. (2012). Role of soil organic matter in maintaining sustainability of cropping systems. Commun. Soil Sci. Plant Anal.43 (16), 2063–2113. 10.1080/00103624.2012.697234
10
FeketeI.BéniÁ.BíróB.GáborV.KatalinJ.MariannaM.et al (2022). Variability in litter inputs affecting soil fungi and bacteria through moisture and carbon content in forest soil. Soil Sci. annu.73 (4), 1–11. 10.37501/soilsa/157106
11
FiererN.BradfordM. A.JacksonR. B. (2007). Toward an ecological classification of soil bacteria. Ecology88 (6), 1354–1364. 10.1890/05-1839
12
GarcíaL. V.Maltez-MouroS.Pérez-RamosI. M.FreitasH.MarañónT. (2006). Counteracting gradients of light and soil nutrients in the understorey of mediterranean oak forests. Web Ecol.6 (1), 67–74. 10.5194/we-6-67-2006
13
GuoZ.LiuC. A.HuaK.WangD.WanS.HeC.et al (2022). Temporal variation of management effects on soil microbial communities. Geoderma418, 115828. 10.1016/j.geoderma.2022.115828
14
JiaY.WhalenJ. K. (2020). A new perspective on functional redundancy and phylogenetic niche conservatism in soil microbial communities. Pedosphere30 (1), 18–24. 10.1016/S1002-0160(19)60826-X
15
Jibola-ShittuM. Y.HengZ.KeyhaniN. O.DangY.ChenR.LiuS.et al (2024). Understanding and exploring the diversity of soil microorganisms in tea (Camellia sinensis) gardens: toward sustainable tea production. Front. Microbiol.15, 1379879. 10.3389/fmicb.2024.1379879
16
KabengiN. J.ZuraykR. A.BaalbakiR. Z.RyanJ. (2003). Phosphorus dynamics and characterization under long-term rotation trial. Commun. Soil Sci. Plant Anal.34 (3-4), 375–392. 10.1081/CSS-120017827
17
KerfahiD.GuoY.DongK.WangQ.AdamsJ. M. (2024). pH is the major predictor of soil microbial network complexity in Chinese forests along a latitudinal gradient. Catena234, 107595. 10.1016/j.catena.2023.107595
18
KielakA. M.BarretoC. C.KowalchukG. A.van VeenJ. A.KuramaeE. E. (2016). The ecology of acidobacteria: moving beyond genes and genomes. Front. Microbiol.7, 744. 10.3389/fmicb.2016.00744
19
LiX.HuoS.XiB. (2020). Updating the resolution for 16S rRNA OTUs clustering reveals the cryptic cyanobacterial genus and species. Ecol. Indic.117, 106695. 10.1016/j.ecolind.2020.106695
20
LiR.SunR.ZhangQ.LiuQ.LuoR.YinC.et al (2025). Shifts in core taxa drive microbial community stability and soil multifunctionality during alpine ecological restoration. J. Appl. Ecol.63, e70262. 10.1111/1365-2664.70262
21
LiangJ. X. (1988). Major soil types of Hainan Island. Chin. J. Trop. Crops9 (1), 53–72.
22
LinB. B.RichardsP. L. (2007). Soil random roughness and depression storage on coffee farms of varying shade levels. Agric. Water Manag.92 (3), 194–204. 10.1016/j.agwat.2007.05.014
23
LiuQ.ChenZ.WangS.LiangT.GaoZ.DongY. (2024). Changes in soil hydrological retention properties and controlling factors on shaded and sunny slopes in semi-arid alpine woodlands. Forests15 (7), 1136. 10.3390/f15071136
24
LiuW.YanR.FangL.ZhangH.ZengH.ShangguanW.et al (2025). Distinct responses of rare and abundant microbial taxa to long-term soil acidification. Soil Ecol. Lett.7, 250294. 10.1007/s42832-025-0294-2
25
MaaßS.MiglioriniM.RilligM. C.CarusoT. (2014). Disturbance, neutral theory, and patterns of beta diversity in soil communities. Ecol. Evol.4 (24), 4766–4774. 10.1002/ece3.1313
26
MiaoR.MaJ.LiuY.LiuY.YangZ.GuoM. (2019). Variability of aboveground litter inputs alters soil carbon and nitrogen in a coniferous - broadleaf mixed forest of central China. For10 (2), 188. 10.3390/f10020188
27
MurphyB. (2015). Key soil functional properties affected by soil organic matter - evidence from published literature. IOP Conf. Ser. Earth Environ. Sci.25 (1), 012008. 10.1088/1755-1315/25/1/012008
28
NainarA.TanakaN.BidinK.AnnammalaK. V.EwersR. M.ReynoldsG.et al (2018). Hydrological dynamics of tropical streams on a gradient of land-use disturbance and recovery: a multi-catchment experiment. J. Hydrol.566, 581–594. 10.1016/j.jhydrol.2018.09.022
29
NavarreteA. A.VenturiniA. M.MeyerK. M.KleinA. M.TiedjeJ. M.BohannanB. J. M.et al (2015). Differential response of acidobacteria subgroups to forest-to-pasture conversion and their biogeographic patterns in the western Brazilian amazon. Front. Microbiol.6, 1443. 10.3389/fmicb.2015.01443
30
OlsenS. R.ColeC. V.WatanabeF. S.DeanL. A. (1954). “Estimation of available phosphorus in soils by extraction with sodium bicarbonate,” in USDA Circular 939. Washington, DC: U.S. Department of Agriculture, 1–19.
31
PowerS. C.VerboomG. A.BondW. J.PackerK. F.CramerM. D. (2021). The role of shade in maintaining alternative stable states between open- and closed-canopy vegetation. J. Ecol.109 (11), 3835–3848. 10.1111/1365-2745.13760
32
SchreckingerJ.MutzM.GessnerM. O.GerullL.FrossardA. (2025). Fundamental shifts in soil and sediment microbial communities and functions during 10 year of early catchment succession. Soil Biol. biochem.203, 109713. 10.1016/j.soilbio.2025.109713
33
ShaoS.LiY.LiZ.MaX.ZhuY.LuoY.et al (2024). Impact of tea tree cultivation on soil microbiota, soil organic matter, and nitrogen cycling in mountainous plantations. Agronomy14 (3), 638. 10.3390/agronomy14030638
34
ShuY.XieS.FanH.DuanC.LiuY.ChenZ. (2025). Tea cultivation: facilitating soil organic carbon accumulation and altering soil bacterial community- leishan county, Guizhou province, southwest China. PeerJ13, e18683. 10.7717/peerj.18683
35
SinsabaughR.HillB.ShahJ. (2009). Ecoenzymatic stoichiometry of microbial organic nutrient acquisition in soil and sediment. Nature462 (7274), 795–798. 10.1038/nature08632
36
SunR.WuZ. X.LanG. Y.YangC.FraedrichK. (2021). Effects of rubber plantations on soil physicochemical properties on Hainan island, China. J. Environ. Qual.50 (6), 1351–1363. 10.1002/jeq2.20282
37
TaggartM. G.BaahD. S.AllenS.KhanZ.ArnscheidtJ.JordanP.et al (2025). Fitting soil extracellular enzyme activity into the complex network of abiotic and biotic soil properties often associated with soil health. Front. Microbiol.16, 1638267. 10.3389/fmicb.2025.1638267
38
TriplerC. E.KaushalS. S.LikensG. E.Todd WalterM. (2006). Patterns in potassium dynamics in forest ecosystems. Ecol. Lett.9 (4), 451–466. 10.1111/j.1461-0248.2006.00891.x
39
WalkleyA.BlackI. A. (1934). An examination of the degtjareff method for determining soil organic matter, and a proposed modification of the chromic acid titration method. Soil Sci.37, 29–38. 10.1097/00010694-193401000-00003
40
WangG. F.SunX. H.FangY.ZhouJ.ShenQ. T.XiangJ. L.et al (2019). Effects of shading on microbial characteristics and enzyme activities in matcha tea garden soil. J. Tea Sci.39 (3), 355–363.
41
WaniA. K.RahayuF.AlkahtaniA. M.AlreshidiM. A.YadavK. K.ParnidiN.et al (2024). Metagenomic profiling of rhizosphere microbiota: unraveling the plant-soil dynamics. Physiol. Mol. Plant Pathol.133, 102381. 10.1016/j.pmpp.2024.102381
42
WuW.ZhangY.TurnerB. L.HeY.ChenX.CheR.et al (2025). Organic amendments promote soil phosphorus related functional genes and microbial phosphorus cycling. Geoderma456, 117247. 10.1016/j.geoderma.2025.117247
43
YakovlevaL. V.DanilovD. A.NikolaevaE. A. (2020). Effect of mineral and organic fertilizers on potassium leaching in sandy loam soils. IOP Conf. Ser. Mater. Sci. Eng.828 (1), 012032. 10.1088/1757-899X/828/1/012032
44
YangX.ChenX.GuoE.YangX. (2018). Path analysis of phosphorus activation capacity as induced by low-molecular-weight organic acids in a Black soil of northeast China. J. Soils Sediments19 (2), 840–847. 10.1007/s11368-018-2034-z
45
YangY.DaiZ.YeC.ZhangX.LiuX.DaiW.et al (2026). Continuous cropping reshapes the rhizosphere microbiome and soil properties in tobacco production systems. Land Degrad. Dev.37 (10), 4694–4711. 10.1002/ldr.70391
46
YuT.WenY.ZhangQ.YangJ.ZangH.ZengZ.et al (2025). Organic management improves soil multifunctionality by enhancing soil quality and enriching key microbes in tea plantation. Appl. Soil Ecol.213, 106260. 10.1016/j.apsoil.2025.106260
47
ZhangX.LiuW.ZhangG.JiangL.HanX. (2015). Mechanisms of soil acidification reducing bacterial diversity. Soil Biol. biochem.81, 275–281. 10.1016/j.soilbio.2014.11.004
48
ZhangS.WangY.SunL.QiuC.DingY.GuH.et al (2020). Organic mulching positively regulates the soil microbial communities and ecosystem functions in tea plantation. BMC Microbiol.20, 103. 10.1186/s12866-020-01794-8
Summary
Keywords
Camellia sinensis, edaphic gradient, functional potential, habitat heterogeneity, metagenomics
Citation
Zhou Y, Li D, Zhang L and Zhang D (2026) Habitat-associated rhizosphere microbiome differentiation in wild tea soils. Span. J. Soil Sci. 16:17010. doi: 10.3389/sjss.2026.17010
Received
26 May 2026
Revised
27 July 2026
Accepted
19 August 2026
Published
31 August 2026
Volume
16 - 2026
Edited by
Raul Ortega Perez, University of Almeria, Spain
Updates
Copyright
© 2026 Zhou, Li, Zhang and Zhang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Dongliang Li, lidongliang@hnaas.org.cn
† These authors have contributed equally to this work and share first authorship
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.