Fertilizer regimes reshape microbial interaction networks without altering sugarcane rhizosphere diversity

Sugarcane (Saccharum spp.), an economically important crop in the food and bioenergy industries, has historically played a central role in the Cuban economy, shaping its agricultural landscape and international relations. Although production has declined in recent decades, sugarcane remains a strategic crop, with its byproducts contributing to national energy and industrial outputs. However, how long-term fertilization shapes soil microbial communities across varying edaphoclimatic conditions remains largely unknown. Here, we investigated a traditional sugarcane plantation to evaluate how distinct fertilizer formulations varying in nitrogen, phosphorus, and potassium affect the soil and rhizosphere microbiota. Using high-throughput sequencing of 16S rRNA (bacterial) and internal transcribed spacer (fungal) genes, we identified 421 bacterial and 471 fungal genera from 5741 amplicon sequence variants across different fertilization regimes. While microbial composition and diversity did not differ significantly between treatments, a global co-occurrence network revealed a highly modular and stable trans-domain architecture. Furthermore, taxonomic overlap analysis highlighted a nutrient-specific recruitment of accessory taxa, suggesting that fertilization regimes may influence rhizosphere functioning by modulating specialized ecological niches rather than through broad taxonomic turnover. The findings provide a detailed characterization of the rhizosphere microbiome of Saccharum spp. in brown sialitic soils (inceptisol), offering ecological insights into its bacterial–fungal associations. This highlights the importance of understanding how long-term fertilization regimes influence rhizosphere microbial dynamics, which is key to designing more sustainable soil management and fertilization practices in sugarcane production systems.

1 Introduction

Cultivated across tropical and subtropical regions for sugar and bioethanol production, sugarcane (Saccharum officinarum L.) is one of the most important industrial crops worldwide (Cherubin et al., 2018; Otto et al., 2022). Sugarcane is predominantly cultivated as a monoculture under intensive management systems that rely on long-term fertilization practices (Cherubin et al., 2018; Dinesh Babu et al., 2022; Khan et al., 2023) due to its high biomass productivity and nutrient demand. This results in cultivation that requires substantial inputs of nitrogen, phosphorus, and potassium (NPK) fertilizers (Cherubin et al., 2019). Sugarcane is primarily produced in Brazil, accounting for 38 % of global production, followed by India (23 %) and China (5 %) (Ritchie et al., 2023). Farmers worldwide have indiscriminately applied nitrogen fertilizers to ensure high crop yields, thereby emitting nitrous oxide into the environment (Tian et al., 2020). To increase nutrient availability, accumulation, and use efficiency in the current sugarcane crop system, new fertilization management strategies are needed. In several production systems, fertilizer application is one of the highest economic and environmental costs of sugarcane farming. This has strong implications for greenhouse gas emissions, soil acidification, and nutrient leaching (Wang et al., 2024). Therefore, assessing the interactions between fertilizer regimes and soil microbial communities in sugarcane agroecosystems has become a key priority (Liu et al., 2024). In Cuba, sugarcane cultivation has historically relied on fixed-dose fertilization schemes, primarily involving NPK (FAO, 2003). Even when these schemes account for soil heterogeneity and plant nutrient demands, sugarcane productivity in Cuba has declined due to inadequate attention to the crop, shortages of fuel and herbicides, reduced water availability, and other limitations (FAO, 2003, 2025). Fertilization regimes often do not consider microbial interactions, which can influence ecosystem functioning (Guo et al., 2020). Nutrient over- or under-application may occur, leading to reduced yields or environmental problems. This can include soil contamination through nutrient leaching and a loss of ecosystem services due to reduced microbial diversity (Berger et al., 2024; Du et al., 2025).

Closely related to fertilization systems, soil microbial communities are central regulators of terrestrial ecosystem functioning. They play essential roles in nutrient cycling, organic matter turnover, soil structure, and plant health (Du et al., 2025). Different ecosystems, including tropical agroecosystems, can host highly diverse microbial communities adapted to specific edaphoclimatic environments (Guo and Chen, 2022; Jiménez-Venegas et al., 2024; Ramos-Tapia et al., 2023; Tripathi et al., 2016; Zamora-Leiva et al., 2025). These microbial consortia are particularly relevant because they directly influence nutrient availability, plant productivity, and resilience to environmental stress (Antoszewski et al., 2022; Kumar et al., 2023; Núñez Salazar et al., 2020; Ramos-Tapia et al., 2022; Souza et al., 2025). Agricultural practices, including tillage, crop rotation, irrigation, and fertilization regimes, can have strong selective pressures on soil microbiomes. This can lead to shifts in microbial community composition and functionality (Liu et al., 2021b; Orrù et al., 2021; Sui et al., 2023). Climatic conditions, such as high rainfall and intense weathering, which are common in tropical ecosystems, often produce nutrient-poor, acidic, and highly leached soils with low cation-exchange capacity (Oishy et al., 2025). Therefore, understanding how soil microbial communities respond to different fertilization strategies is crucial to optimizing sustainable crop management and maintaining soil health.

Although continuous monocropping, combined with sustained fertilizer inputs, can strongly influence soil health and microbial community structure, the effects are not easily predictable. Long-term fertilization can produce contrasting outcomes, ranging from enhanced nutrient cycling to reduced microbial diversity and ecosystem resilience (Akter et al., 2025; Chen et al., 2023; Kracmarova et al., 2022; Pan et al., 2014). Nitrogen fertilization has been associated with decreases in microbial alpha diversity, community shifts toward copiotrophic taxa, and reductions in nitrogen-fixing microorganisms (Sui et al., 2023). Phosphorus and potassium amendments affect fungal communities, particularly arbuscular mycorrhizal fungi, which are highly responsive to changes in nutrient availability (Semenov et al., 2022; Yang et al., 2025). This indicates that fertilizer regimes can shape soil microbial assemblages in ways that alter nutrient cycling processes, soil–plant interactions, and ecosystem sustainability (de Castro et al., 2022; Mitter et al., 2021). However, soil microbial community responses to fertilization vary across ecosystems and are highly context-dependent (Kracmarova et al., 2022; Moneda et al., 2022). This response is influenced by several factors, including soil type (Huang et al., 2025), climate (Du et al., 2025), cropping history (Pan et al., 2014), and management intensity (Suman et al., 2022). In certain cases, microbial communities show considerable resilience, maintaining similar diversity and composition across a wide range of fertilization treatments (Philippot et al., 2021). This variability complicates extrapolating results across agroecosystems and underscores the need for site-specific studies (Souza et al., 2025).

Microbial communities are commonly studied through high-throughput sequencing techniques (Lahlali et al., 2021). These technologies have substantially expanded our ability to explore the structure and function of microbial communities in environmental samples, describing the taxonomic composition and relative abundance of microbial taxa across diverse ecosystems (Lahlali et al., 2021). However, while attention has been traditionally focused on how community composition and diversity shift under different environmental or management conditions, an equally important aspect lies in the relationships among microbial taxa (Zamkovaya et al., 2021). These interactions, whether cooperative, competitive, or neutral, shape the functional stability and resilience of microbial ecosystems (Freilich et al., 2011). Network-based approaches are powerful tools to explore these relationships, allowing the identification of co-occurrence patterns, ecological clusters, and potential keystone taxa that play central roles in ecosystem functioning (Shaw et al., 2020). Integrating network analysis with community profiling can provide a deeper understanding of how microbial communities respond to environmental factors, such as fertilization regimes (Gao et al., 2022).

In the context of Cuban agriculture, sugarcane is largely cultivated on brown sialitic soils, which are mineral-rich soils with relatively underexplored microbiological dynamics. Understanding how fertilization interacts with local edaphoclimatic conditions is essential for improving sustainable management strategies and maintaining soil function in tropical agroecosystems. However, traditional fertilization recommendations in Cuba, such as those from SERFE (Servicio de Recomendaciones de Fertilizantes), primarily focus on crop nutritional requirements and do not explicitly account for soil microbial processes. Notably, available yield records from the long-term fertilization trial, established in 1998 at the Sugar Cane Research Institute in Santiago de Cuba, indicate that while fertilized plots consistently outperform non-fertilized controls, sugarcane yield tends to plateau across a wide range of NPK fertilization regimes, with limited differences among moderate and high input treatments (Table S1 in the Supplement (Tamayo-Isaac et al., 2024)). This apparent decoupling between fertilizer dose and yield response suggests that factors other than nutrient availability alone may regulate system functioning under long-term management. We therefore hypothesized that sustained fertilization under tropical conditions leads to physicochemical homogenization of these soils, generating changes in microbial community composition and promoting nutrient-specific reorganization of microbial interactions within the rhizosphere. To test this hypothesis, we leveraged the long-term fertilization trial established in 1998 at the Sugar Cane Research Institute in Santiago de Cuba, which provides a rare opportunity to assess the legacy effects of sustained fertilization on soil microbial network structure. Using high-throughput sequencing of bacterial (16S rRNA gene) and fungal (internal transcribed spacer; ITS region) markers, we evaluated how gradients of NPK inputs shape both microbial composition and co-occurrence networks in rhizosphere and bulk soil samples.

2 Materials and methods

2.1 Experimental design and sample collection

2.1.1 Profile of the experimental site

The long-term experiment, encompassing different mineral fertilization rates of NPK, was established in 1998 on sialitic brown soil at the Estación Territorial de Investigaciones de la Caña de Azúcar (ETICA) Oriente Sur, and has maintained the same fertilization regime across the years. Crucially, fertilization regimes assigned to each plot have been strictly and continuously maintained without interruption since the trial began. The study reproduced the traditional sugarcane cultivation system used in Cuba, featuring manual planting and an annual harvesting scheme. Crop management followed consistent 5-year cycles, consisting of one initial plant cane harvest followed by four ratoon harvests, before the plots were entirely replanted. The station is located in the municipality of Palma Soriano in the province of Santiago de Cuba (Sieiro and Pablos, 2021). Sialitic brown soils have moderate organic matter content, high clay fraction, and a dominance of exchangeable Ca²⁺ and Mg²⁺. These properties contribute to a relatively stable pH and high cation exchange capacity. The area is characterized by a humid tropical climate, with an average annual rainfall of 1375 mm, and two well-defined periods: rainy (May–October) and dry (November–April). For the current study, sampling was conducted in March 2022, during the second ratoon crop on the ongoing cycle (which was planted in spring 2020 with cultivar C86-12). Because there were no gaps in management between 1998 and 2022, the soil properties and microbial communities evaluated herein reflect the robust, cumulative 24-year legacy of these continuous, treatment-specific fertilization practices.

2.1.2 Experimental design

The completely randomized block experiment included 11 treatments, with different co-concentrations of NPK, all expressed in kg ha⁻¹ (Table 1): (i) N0 (N0P25K50), no nitrogen application; (ii) N1 (N50P25K50), as a basal application of nitrogen, phosphorus and potassium; (iii) N2 (N100P25K50), medium nitrogen application; (iv) N3 (N150P25K50), high nitrogen application; (v) N4 (N200P25K50), very high nitrogen application; (vi) N5 (N50P0K0), only nitrogen application; (vii) P1 (N50P75K50), medium phosphorus application; (viii) P2 (N50P100K50), high phosphorus application; (ix) K1 (N50P25K100), medium potassium application; (x) K3 (N50P25K200), high potassium application; and (xi) K4 (N50P25K300), very high potassium application.

Table 1 Concentration of fertilization for the treatments.

Derived from the experimental design, three comparison groups were generated to evaluate different parameters (Table 2): (i) effect of the increase in nitrogen, (ii) effect of the increase in phosphorus, and (iii) effect of the increase in potassium.

Table 2 Fertilization regimes and comparison groups.

Two controls were considered: bulk soil without planting (blank) and rhizosphere without fertilization (control). Each treatment was repeated three times, and each replicate was 48 m². The chemical fertilizers applied were urea (46 % N), calcium superphosphate (12 % P₂O₅), and potassium chloride (60 % K₂O). All fertilizers were applied at the bottom of the furrow before planting and 20 d post-harvest, manually and without fractioning doses (INICA, 1990). The evaluated plot had received three fertilizer applications by the sampling date, with planting in spring 2020 and sampling in March 2022, corresponding to the initial application plus those made after each cutting/ratoon.

2.1.3 Soil sampling

Sugarcane root systems were extracted at a depth of 10 cm and gently shaken to remove the loosely bound bulk soil. The soil that remained firmly adhered to the roots was considered the rhizosphere. These root-adhering soil samples were then carefully collected into 50 mL Falcon™ tubes (Thermo Fisher Scientific, Waltham, MA, USA) until the volume was filled with roots. The tubes were immediately sealed to avoid environmental contamination, labeled, and stored at 4 °C during the laboratory transport (Ramos-Tapia et al., 2022). Plants were vortexed to dislodge soil particles adhering to the roots, and this soil was subsequently used for DNA extraction. Simultaneously, the soil removed from around the roots (within 3 cm) was collected as a composite soil sample for physicochemical analysis. The composite soil (control and blank) samples comprised five subsamples of approximately 200 g, which were combined and mixed to better capture the natural heterogeneity of the sampled area. Subsequently, these composite soil samples from each experimental plot were used for chemical analysis and soil respiration, developed by ETICA (Tamayo-Isaac et al., 2023, 2024). Sampling was conducted in March, one of the driest months of the year (Centella et al., 2007; Durán, 2012). For physicochemical parameters, pH in water and KCl was measured with potentiometric method (soil:solution ratio 1:2.5); organic matter with Walkley–Black colorimetric method, digestion with H₂SO₄, K₂Cr₂O₇ 1 N (ox.) and FeSO₄ x7 H₂O 0.5 N (red.), using Vam Bemmen's factor (1.724) to convert carbon into organic matter; assimilable phosphorus as P₂O₅ mg per 100 g with Oniani method, extraction with 0.1 N H₂SO₄, soil:solution ratio 1:2.5, determination by colorimetry; assimilable potassium as K₂O mg per 100 g with Oniani method, extraction with 0.1 N H₂SO₄, soil:solution ratio 1:2.5, by flame photometry; and exchangeable bases (K⁺, Ca²⁺, Mg²⁺, and Na⁺) as cmol(+) kg⁻¹ with extraction with 1N NH₄OAc at pH 7, determining calcium and magnesium by EDTA titration, sodium and potassium by flame photometry (INICA, 1990; Tamayo-Isaac et al., 2024). Soil physicochemical properties were analyzed using a one-way ANOVA (p<0.05) followed by Tukey's HSD post-hoc test.

2.1.4 DNA extraction, PCR, and sequencing

Microbial DNA was successfully extracted from 39 samples, comprising 33 fertilized rhizosphere samples (11 treatments with three replicates each), three non-fertilized rhizosphere controls, and three bulk-soil (blank) controls. For each sample, genomic DNA was extracted from 0.5 g of fresh soil using the DNeasy Powersoil DNA isolation kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. DNA concentration and purity were checked using a Qubit dsDNA HS Assay kit and Qubit 4.0 fluorometer (Thermo Fisher Scientific). The V3-V4 hypervariable regions of the 16S rRNA gene were amplified with the primers 5^′-CCTACGGGGGNGGCWGCAG-3^′ and 5^′-GACTACHVGGGTATCTAATCC-3^′ (Klindworth et al., 2013). Meanwhile, for the fungal ITS1 and ITS2 regions, the primers 5^′-CTTGGTCATTTAGAGGAAGTAA-3^′ and 5^′-GCTGCGTTCTTCATCGATGC-3^′ were used (Orgiazzi et al., 2012). Amplification and bulk sequencing were performed with Illumina MiSeq technology by Austral-Omics (Universidad Austral de Chile). The sequencing process was performed with the Miseq Reagent v3 600-cycle kit (2 × 250 pair-end). The sequences were imported as FASTQ files with forward and reverse for each sample.

2.1.5 Bioinformatic analysis

The FASTQ files were processed using QIIME 2 v2025.7 (Caporaso et al., 2010). Read quality was evaluated using FastQC v0.11.9 (Andrews, 2010). Based on this evaluation, adapters were removed, and the forward and reverse sequences were trimmed to 260 and 205 bp for bacteria and to 250 and 220 bp for fungi, respectively, using the DADA2 denoise-paired plugin (Callahan et al., 2016). Quality was maintained using PHRED 20 (Callahan et al., 2016). Taxonomic assignment was performed using the classify-sklearn plugin in the SILVA database version 138.2 (Quast et al., 2013) and UNITE database version 10.0 (Abarenkov et al., 2024), for bacterial and fungal classification, respectively, both preprocessed with RESCRIPt (Robeson et al., 2021). New unified Amplicon Sequence Variant (ASVs) with 99 % identity were aligned using MAFFT (Katoh and Standley, 2013). The maximum likelihood was inferred using FastTree v2.2.0 (Price et al., 2009). The abundance, taxonomy, phylogeny, and metadata of the ASVs were integrated into a Microeco object for subsequent analyses using the Microeco v1.15.0 R package (Liu et al., 2021a), in R v4.3.2 (R Core Team, 2023). ASVs assigned to archaea, mitochondria, chloroplasts, metagenomes, and ambiguous taxa were deleted, as well as fungal ASVs in bacterial analysis, and bacterial ASVs for fungal analysis. Samples with < 1000 reads were excluded; however, all samples met the criteria. ASVs that did not occur more than twice in at least 10 % of the samples were also excluded (Callahan et al., 2017). To minimize potential bias, unclassified ASVs were strictly retained in all analyses but excluded from taxonomic visualizations to ensure that the discussion focused on identifiable biological entities. The Microeco object was used to calculate alpha and beta diversities and for constructing relative abundance plots (Liu et al., 2021a). Alpha diversity was calculated using observed features – species richness on samples, Shannon – entropic information on abundance of observed ASVs, and Simpson (1-dominance) indices. Alpha diversity metrics were assessed using the Kruskal-Wallis test, and beta diversity differences were evaluated using Permutational Multivariate Analysis of Variance (PERMANOVA).

To determine beta diversity, all samples were compared across sampling points using the Bray–Curtis, Canberra, Jaccard, Robust Aitchison, Unifrac, and weighted Unifrac matrices. Significant differences were evaluated using Kruskal–Wallis tests and PERMANOVA. ASVs were transformed into proportions, and relative abundance was used to measure the composition of microbial communities at the phylum, class, order, family, and genus levels. All samples were rarefied according to the minimum sample size of 11 248 and 18 829 for bacteria and fungi, respectively. Differential taxon abundance was analyzed by implementing Linear Discriminant Analysis Effect Size (LEfSe, linear discriminant analysis [LDA] value > 2, α = 0.05), supported with an Analysis of Compositions of microbiomes with Bias Correction 2 (ANCOM-BC2) analysis.

For the co-occurrence network and overlapping community analyses, the bacterial and fungal datasets were merged at the ASV level. To account for differences in sequencing depth, the complete dataset was rarefied to a uniform depth of 11 248 reads per sample, using the trans_rarefy function from the Microeco R package. To evaluate microbial interactions, a single, unified global trans-domain network incorporating all 39 samples was constructed using the cal_network function in the Microeco R package (which implements igraph v2.1.4). Pairwise correlations between ASVs were computed using Spearman's rank with Benjamini-Hochberg false discovery rate (FDR) correction (Looney and Hagan, 2011), retaining only robust and significant connections (Spearman rho > 0.6, adjusted p < 0.01). Network modularity was defined as the extent to which a network is organized into cohesive subgroups of nodes (modules), with higher within-module connectivity than between modules. Modules were detected using the cluster_fast_greedy algorithm, and modularity values were computed with the modularity function in the igraph R package. Network hubs were identified based on within-module connectivity (Zi) and among-module connectivity (Pi) using the get_node_table function from the Microeco R package. The network was constructed with co-occurrences between ASVs and visualized using Gephi v0.10.1 (Bastian et al., 2009). Finally, to evaluate the specific influence of each macronutrient on microbial recruitment, we performed a sequential overlap analysis based on presence/absence matrices with the trans_venn function. Unclassified ASVs and Incertae sedis were filtered out before analysis to ensure that all overlapping communities corresponded to identifiable biological entities. Initially, a universal core microbiome was identified by intersecting all experimental conditions (N, P, and K gradients). Nutrient-specific accessory microbiomes were then filtered by excluding ASVs present under the basal condition (N1) or associated with multiple fertilization treatments. We identified nutrient-specific specialists by calculating the intersection of ASVs present exclusively in the modified treatments (e.g., nitrogen gradient) while simultaneously excluding those detected in the basal soil or across parallel treatments (e.g., phosphorus or potassium gradients for nitrogen analysis). The resulting intersection patterns and microbial recruitment dynamics were visualized as UpSet plots using the UpSetR package v1.4.1.

3 Results

3.1 Physicochemical characterization of brown sialitic soils under different fertilization regimes

The physicochemical properties of soils across fertilization treatments are summarized in Table 3. The data showed minimal variation among treatments, suggesting strong homogenization of soil conditions across the experimental site. Soil pH values remained nearly constant, ranging from 6.4 to 6.5 (HCl) and 7.7 to 7.9 (H₂O). This indicates a slightly acidic-to-neutral reaction, typical of sialitic soils. This stability is consistent with the buffering capacity conferred by high levels of exchangeable calcium (58–69 cmol(+) kg⁻¹) and magnesium (0.4–11 cmol(+) kg⁻¹). While statistical differences were detected among plots (one-way ANOVA, p<0.05, Table S1), the sum of exchangeable calcium and magnesium remained consistently high across the entire experimental site, strictly constrained within a highly saturated range (60.7–72.2 cmol(+) kg⁻¹). This persistently high reservoir of basic cations reinforces the notion of a chemically buffered soil matrix that can resist dramatic acidification despite long-term fertilization. Organic matter (OM) content ranged from 2.9 % to 3.9 %, with no clear trends associated with fertilization. Similarly, available phosphorus (P₂O₅) and potassium (K₂O) concentrations exhibited broad ranges (P₂O₅: 63–139 mg per 100 g; K₂O: 25–88 mg per 100 g). Despite statistical differences among plots (one-way ANOVA, p<0.01; Table S1), these variations did not translate into distinct patterns of nutrient enrichment in response to applied fertilizers. This lack of correlation suggests that rapid nutrient turnover or leaching, potentially driven by intense rainfall and the clay-rich matrix of sialitic soils, may diminish the effect of fertilization inputs.

Table 3 Physicochemical characterization of sites.

Soil respiration, i.e., microbial CO₂ release after 9 d of incubation, ranged between 1440 and 2100 mg CO₂, with significant differences among treatments (one-way ANOVA, p<0.05, Table S1). However, soil respiration exhibited erratic spatial variations, with the unplanted bulk soil (blank) showing higher respiration rates than several fertilized rhizosphere plots. This lack of logical correlation strongly suggests that historical field heterogeneity, rapid nutrient turnover, and severe leaching, driven by the intense rainfall characteristic of the tropical climate and the clay-rich matrix of sialitic soils, have overridden the immediate effects of the fertilizer inputs. Although microbial activity was slightly higher in plots with intermediate and high nitrogen and medium phosphorus fertilization (> 2000 mg CO₂ in each case), the inconsistent pattern prevents establishing a direct link between fertilizer dose and microbial respiration. Crop yield also showed minor fluctuations, ranging from 87.7 to 99.8 t ha⁻¹. It appeared to be negatively affected by high nitrogen fertilizer levels, with the high (N150) and very high (N200) treatments yielding the lowest among the fertilized plants (Table S1). Despite the application of nutrient gradients, soil and plant responses remained broadly decoupled from the theoretical fertilizer inputs.

3.2 Characterization of bacterial and fungal communities

Microbial DNA was successfully extracted from a total of 39 samples, comprising 33 rhizosphere samples (11 treatments with three replicates each), three non-fertilized rhizosphere samples (control), and three bulk soil (non-rhizosphere) samples (blank). For the bacterial and fungal analyses, we obtained 3034 and 2707 ASVs identified post-quality filters among all the samples, with a median of 16 883 and 48 883 reads per sample, respectively. Taxonomic assignments were made at the phylum (98.26 % and 91.28 %), class (96.94 % and 85.52 %), order (89.46 % and 84.97 %), family (68.98 % and 73.59 %), and genus (34.46 % and 66.89 %) levels, in bacteria and fungi, respectively. The percentage of classified ASVs diminished substantially for the bacterial sequences. Meanwhile, nearly twice as many fungal ASV sequences were classified at the genus level.

Across all fertilization treatments, bacterial communities associated with the rhizosphere of Saccharum spp. were dominated by a few major phyla. This reflected a relatively stable community structure. The most abundant bacterial phyla were Actinomycetota, Pseudomonadota, and Acidobacteriota, which together accounted for approximately 60.4 % of the total relative abundance. Subdominant groups included Bacillota and Planctomycetota, each representing less than 16.6 % of the total sequences (Fig. 1A). Within these major lineages, members of Bacilli (Bacillota), Alphaproteobacteria (Pseudomonadota), Vicinamibacteria (Acidobacteriota), and Thermoleophilia (Actinomycetota) were present across all treatments at relative abundance exceeding 11 %. This suggests that the fertilization regime did not substantially alter the core microbiota composition (Fig. 1B). Minor variations were observed in the relative abundance of Actinomycetota, which increased in soils receiving higher nitrogen inputs. However, these differences were not statistically significant (Kruskal–Wallis, P>0.05). Although Actinomycetota was the dominant phylum, Bacillota comprised the most represented taxa from class to genus levels, demonstrating that most ASVs in this phylum were successfully classified. At the order level, the number of incertae sedis and unclassified taxa increases, representing the fifth and sixth in relative abundance, respectively, with 7.4 % and 3.1 % of the total relative abundance (Fig. 1C). Relative abundances at the genus level are detailed in Fig. S1 in the Supplement.

Among fungal communities, the phylum Ascomycota was dominant across all treatments, with a relative abundance exceeding 69 %. This was followed by Basidiomycota and Chytridiomycota, which together represented almost 18 % of relative abundance (Fig. 2A). Six minor groups, such as Kickxellomycota, Glomeromycota, Calcarisporiellomycota, Mucoromycota, Neocallimastigomycota, and Zoopagomycota, were present in lower abundances (< 1 %) and did not show significant treatment-related variation (Fig. 2A). At the class level, even when several taxa appeared to be involved in the ecosystem, only three classes represent more than 65.4 % (Sordariomycetes, Eurotiomycetes, and Leotiomycetes), which belong to the Ascomycota phylum (Fig. 2B). This behavior repeats at the order level, where six taxa (Eurotiales, Helotiales, Hypocreales, Sordariales, Chaetothyriales, and Pleosporales) that belong to the phylum Ascomycota represent more than 54 % of the total (Fig. 2C). This result underscores the importance and prevalence of a few representatives of the phylum Ascomycota in sugarcane communities. Relative abundances at the genus level are detailed in Fig. S2.

Figure 1 Bacterial relative abundance at three different taxonomic levels. (A) 10 most abundant bacteria per sample at the phylum level, (B) 10 most abundant bacteria per sample at the class level, and (C) 10 most abundant bacteria per sample at the order level. Sample abbreviations: K1, medium potassium application; K3, high potassium application; K4, very high potassium application; N0, no nitrogen application; N1, low nitrogen application; N2, medium nitrogen application; N3, high nitrogen application; N4, very high N application; N5, only nitrogen application; P1, medium phosphorus application; P2, high phosphorus application.

Figure 2 Fungal relative abundance at three different taxonomic levels. (A) 10 most abundant fungi per sample at the phylum level, (B) 10 most abundant fungi per sample at the class level, and (C) 10 most abundant fungi per sample at the order level. Sample abbreviations: K1, medium potassium application; K3, high potassium application; K4, very high potassium application; N0, no nitrogen application; N1, low nitrogen application; N2, medium nitrogen application; N3, high nitrogen application; N4, very high N application; N5, only nitrogen application; P1, medium phosphorus application; P2, high phosphorus application.

3.3 Alpha and beta diversity analysis

Alpha diversity indices (Observed, Shannon, and Simpson) showed no significant differences among fertilization treatments for either bacterial or fungal communities (Kruskal–Wallis, P>0.05, Figs. S3 and S4, respectively). This pattern indicates that nutrient fertilization did not significantly affect within-sample diversity. Microbial communities remained compositionally stable despite long-term exposure to different NPK regimes (Fig. 3). Even so, bacterial communities exhibited relatively high species richness and evenness, with average Shannon indices ranging from 4.9 to 5.9. Meanwhile, fungal diversity indices were lower, ranging from 2.9 to 4.3, suggesting a more specialized community structure in the rhizosphere.

Figure 3 Alpha diversity for rhizosphere and soil samples based on the Shannon index. (A) Bacterial and (B) Fungal. Dunn's Kruskal–Wallis test, non-significant. Sample abbreviations: K1, medium potassium application; K3, high potassium application; K4, very high potassium application; N0, no nitrogen application; N1, low nitrogen application; N2, medium nitrogen application; N3, high nitrogen application; N4, very high N application; N5, only nitrogen application; P1, medium phosphorus application; P2, high phosphorus application.

Based on six indices (Bray–Curtis, Canberra, Jaccard, Robust Aitchison, Unifrac, and weighted Unifrac matrices), beta diversity analysis showed no clear clustering of samples according to fertilization treatment (PERMANOVA, P>0.05; Fig. 4 shows the Bray–Curtis matrix as a representative). This implies that other environmental factors, such as soil texture, OM, or rainfall regime, may have a stronger influence on microbial assemblages than the applied fertilization treatments. Principal Coordinate Analysis does not reveal clustering among samples across different fertilization treatments. The microbial communities associated with Saccharum spp. in brown sialitic soils are compositionally rich, yet functionally stable, showing limited responsiveness to changes in nutrient fertilization.

Figure 4 Beta diversity for rhizosphere and soil samples. Principal coordinate analysis determined by using the Bray–Curtis method test for (A) Bacterial and (B) Fungal communities; and group distances calculated with Dunn's Kruskal–Wallis test for (C) Bacterial and (D) Fungal communities, with non-significant differences.

To further investigate potential treatment-driven shifts at specific taxonomic levels, differential abundance analyses were conducted using LEfSe and ANCOM-BC2. Notably, neither approach identified any taxa with statistically significant differential abundance among the evaluated fertilization regimes (LDA score < 2.0; adjusted p>0.05). This lack of significant biomarkers further corroborates the overall compositional stability of the sugarcane rhizosphere microbiome across treatments.

3.4 Co-occurrence between bacterial and fungal taxa

The global cross-domain co-occurrence network integrating the bacterial and fungal datasets indicated a complex and heterogeneous microbial structure, encompassing 371 nodes and 747 edges (Fig. 5). Taxonomic annotations demonstrated a clear numerical predominance of fungal taxa within the network topology (218 fungal nodes vs 153 bacterial nodes). The phylum Ascomycota was the most heavily represented in the network, followed closely by Basidiomycota. Among the bacterial domain, the network was primarily supported by members of Pseudomonadota, Bacillota, Actinomycetota, and Planctomycetota. This topological distribution suggests that fungi may play a central structural role within the rhizosphere-associated microbiome in sugarcane cultivated in brown sialitic soils, while diverse bacterial taxa appear to integrate into this framework, potentially providing extensive functional connections. Notably, while the phylum Methylomirabilota exhibited a notable relative abundance (1.8 %), its representation in the network was limited, highlighting the distinction between raw abundance and ecological interaction roles.

Figure 5 Microbial interaction network. Bacterial and fungal ASVs are graphed as nodes. The relative abundance is displayed as the size of the nodes, while the relationship between nodes is displayed as edges (positive relations, lightblue line; negative correlations, red). The networks show: (A) Fungi (orange) and Bacteria (blue); (B) modules; and (C) Ecological role related to general community (grey), nitrogen associated (green), phosphorus associated (red), and potassium associated (purple) ASVS. Significant positive Spearman correlations are graphed as edges (ρ > 0.6, false discovery rate adjusted P < 0.01, two-sided).

The global network was organized into 48 distinct structural modules (Fig. 5B), indicating that independent trans-domain clusters of microorganisms coexist within the rhizosphere. However, the community architecture was heavily dominated by five major modules (M1–M5), which accounted for 48.2 % of the total network nodes. According to the Zi-Pi topological analysis, the vast majority of ASVs (95 %) acted as peripheral nodes, while a small but highly connected subset served as connectors and module hubs, potentially serving as stabilizers of the overall architecture (Fig. S5). Notably, fungi generally exhibited higher connectivity (degree) than bacteria, suggesting that they may serve as key structural components of the network. Among these, one ASV was classified as a potential hub or keystone taxon; both belonged to the phylum Pseudomonadota but with no classified genus (Fig. S5). Among the bacteria, highly connected taxa included members of Pseudomonadota, such as Pseudomonas and Steroidobacter, and Actinomycetota (Gaiella and Streptomyces). Meanwhile, fungi showed greater connectivity than bacteria, with higher degrees within the network. Even though the largest number of nodes belonged to Ascomycota and Basidiomycota, Mortierellomycota also showed higher connectivity than bacteria with higher node numbers. The interactions between bacteria and fungi appear to be even (Fig. S6).

A substantial portion of the network nodes corresponded to unclassified taxa from both fungal and bacterial groups. This indicates that a substantial portion of the interacting “microbial dark matter” in these soils remains taxonomically unresolved but may represent a consistent and relevant structural component of the community. Furthermore, the network topology was overwhelmingly dominated by positive correlations (99.8 %), with only five negative correlations, all between Pseudomonadota and Bacillota ASVs. This result reflects potential mutualistic, commensal, or shared environmental responses among microorganisms rather than competitive exclusion. This strong predominance of positive associations may contribute to the observed ecological resilience of these soils, as community composition remained largely unchanged across fertilization regimes.

3.5 Taxonomic breakdown of shared and exclusive microbial groups across fertilization regimes

To identify shared and exclusive ASVs across different fertilization strategies without the statistical limitations of small-sample networks, trans-domain overlap analysis was generated directly from the standardized abundance matrices (Figs. 6 and S7). To identify nutrient-specific ASVs, we compared the basal condition (N1) with each subsequent nutrient addition and identified ASVs detected exclusively following the addition of that nutrient. From here, we identified microorganisms associated with the presence of NPK (described in Supplementary Table “Shared ASs”). The Universal Resiliency Core: By crossing the conserved taxa across all experimental gradients (NPK), we identified a highly conserved, trans-domain universal core of 76 ASVs (Fig. 6A and Supplement Table “Shared ASVs”). This basal core likely represents a foundational component of the sugarcane rhizosphere in these soils, relying heavily on a dual-domain alliance between robust bacterial phyla, such as Bacillota (17 ASVs) and Pseudomonadota (13 ASVs), while most ASVs belong to the fungal Ascomycota (26 ASVs, including genera like Fusarium and Cladophialophora). The preservation of this core suggests that long-term sugarcane cultivation establishes a stable basal microbiome that resists profound nutritional fluctuations. Nutrient-specific specialists and General Oportunists: Beyond the universal core, we investigated the accessory microbiomes triggered by nutrient addition compared to the basal condition (Fig. 6B–D and Supplement Table “Shared ASVs”). To ensure strict nutrient-specificity, we filtered out overlapping ASVs that responded indiscriminately to multiple fertilizers. When isolating strictly specific specialists, macronutrients acted as profound environmental filters: Nitrogen variation (N2, N3, and N4) specifically recruited a highly specialized group of 18 unique ASVs, predominantly composed of fungal taxa, with Ascomycota predominating (Fig. 6B and E). In contrast, the application of secondary nutrients generated a targeted recruitment response. Specifically, phosphorus fertilization (P1 and P2) recruited 40 unique ASVs that were completely absent in the basal nitrogen plots (Fig. 6C and E). Crucially, this P-driven accessory microbiome was dominated by fungi (27 fungal ASVs vs. 13 bacterial ASVs), primarily by Ascomycota (e.g., genus Talaromyces) and Basidiomycota. On the other hand, potassium gradients (K1, K3, and K4) exhibited the most restricted specific recruitment, with an accessory cohort of only six unique ASVs (Fig. 6D and E). These taxonomic breakdowns align with the topological network findings, suggesting that although bacteria are highly diverse in the soil, nutrient variations – particularly phosphorus – may act as a potent environmental filter that promotes the emergence and recruitment of specialized trans-domain fungal niches rather than fragmenting the overall community structure.

Figure 6 UpSet plots depicting the core and accessory microbiome assemblies across fertilization regimes. Vertical bars indicate the number of shared and unique ASVs for each intersection, derived directly from standardized abundance matrices. Connected dots below the x-axis specify which fertilization treatments are included in each intersection. (A) Global trans-domain overlap analysis across all nine experimental conditions, highlighting the highly conserved Universal Resiliency Core of 76 ASVs shared across all treatments. (B–D) Intersections of the complete microbial communities within individual nutrient variations compared to the basal soil condition (N1): Nitrogen gradient (B), Phosphorus addition (C), and Potassium gradient (D). (E) Accessory microbiome isolating the ASVs strictly recruited by specific nutrient additions. This panel explicitly differentiates strict nutrient-specific specialists (single dots, e.g., 40 ASVs exclusive to Phosphorus) from general opportunistic taxa shared across multiple fertilization strategies (connected dots).

4 Discussion

4.1 Soil buffering capacity and leaching override dose-dependent physicochemical gradients

As Cuba is one of the few countries that uses a nationally developed soil classification system (Hernández-Jiménez et al., 2019), brown sialitic soils are rarely referenced in the literature. Therefore, in this study, we have provided one of the first integrative characterizations of rhizosphere-associated microbial communities of Saccharum spp. cultivated on this soil type. However, the Cuban classification can be linked to international systems. Brown sialitic soils have historically been classified as inceptisols in United States Department of Agriculture (USDA) Soil Taxonomy and Cambisols in the World Reference Base (WRB) and as Food and Agriculture Organization of the United Nations (FAO) classification system (Sanchez, 2019). Inceptisols are young soils with weak horizon development, widely distributed across diverse landscapes (Soil Survey Staff, 2015). They typically have low natural fertility (Syamsiyah et al., 2018), characterized by limited availability of nitrogen, phosphorus, and potassium; low OM content and slightly acidic pH (Sofyan et al., 2024). These conditions contribute to nutrient leaching, with most N losses in inceptisols occurring within the first week after fertilization (Sofyan et al., 2024). When nutrient uptake by plants during that time is insufficient, crops may experience long-term nutrient limitation (Rashid et al., 2024). Slow-release fertilizers are an effective strategy to mitigate nutrient leaching in inceptisols (Hakim et al., 2022). Therefore, this may be relevant for improving fertilizer efficiency in Cuban sugarcane production systems. Our physicochemical analysis revealed significant spatial heterogeneity across the experimental site; however, these variations were broadly decoupled from the theoretical fertilizer doses applied. Rather than forming distinct chemical gradients, the inputs appear to have been rapidly redistributed. This lack of dose-dependency strongly suggests that intense rainfall and the inherent leaching potential of sialitic clay matrices override the immediate local effects of fertilization (Oishy et al., 2025), as well as intrinsic soil properties, agricultural practices, such as tilling of sugarcane, and external climatic influences (McDonald et al., 2024; Peng et al., 2024; Skidmore et al., 2023). Brown sialitic soils have mid-to-high nutrient loss rates (Hernández et al., 2015; Kaur et al., 2024), and under long-term cultivation and repeated fertilizer inputs, nutrient redistribution via rainfall and leaching can lead to chemical homogenization within the soil matrix. Furthermore, the robust stability of soil pH across all plots reflects the massive buffering capacity conferred by the consistently high reservoir of exchangeable basic cations (Ca + Mg). Consequently, microorganisms in this long-term agroecosystem are exposed to a highly buffered, yet spatially erratic, nutritional environment. This phenomenon highlights the importance of soil microhabitats, particularly physically protected aggregates formed by roots, fungal hyphae, and bacterial polymers (Liao et al., 2020; Wilpiszeski et al., 2019), which can maintain biological organization even when bulk soil chemistry is homogenized or heavily leached by tropical environmental processes. This is especially relevant under long-term fertilization, as organic and inorganic fertilizer addition can increase organic and inorganic carbon contents and enhance cation-exchange capacity of soil macroaggregates (Meena et al., 2022; Niu et al., 2022). These aggregates form physically protected microhabitats, and the rhizosphere promotes their formation (Li et al., 2020). Here, plant roots, fungal hyphae, and bacterial extracellular polymers bind particles together and create localized niches with steep gradients of oxygen, carbon availability, moisture, and nutrients (Wilpiszeski et al., 2019; Zhong et al., 2017). The distribution and relative abundance of micro- and macroaggregates also influence the bulk properties of soil, including organic carbon content, water content, and niche availability. The interiors of aggregates can exhibit properties distinct from the surrounding matrix (Wilpiszeski et al., 2019). This is consistent with the results of the present study, in which microbial communities showed nutrient-dependent restructuring despite the absence of substantial changes in community composition.

4.2 Long-term fertilization maintains overall community composition and diversity

Despite variation in synthetic fertilizer application, alpha- and beta-diversity analyses demonstrated that the overall microbial community composition remained remarkably stable (Figs. S3 and S4). This taxonomic resilience was further confirmed by differential abundance analyses (LEfSe and ANCOM-BC2), which detected no significant biological markers associated with any specific treatment. This reflects the adaptation of microbial communities to persistent, long-term management regimes, reaching a resilient equilibrium where taxonomic structure is maintained despite environmental fluctuations, a pattern consistent with other intensively cultivated soils (Kostin et al., 2021; Lori et al., 2023). In such systems, microbial assemblages reach a resilient equilibrium, maintaining taxonomic structure despite environmental or nutritional variations. Similar results have been reported in long-term agricultural trials (Lori et al., 2023; Luo et al., 2023). Furthermore, the results showed no difference relative to non-fertilized soil or bulk soil, suggesting that the microbial consortia associated with sugarcane in brown sialitic soils are resilient to variations in nitrogen, phosphorus, and potassium availability. Such resilience may be linked to the soil physicochemical properties and climatic conditions in eastern Cuba, characterized by high rainfall and seasonal variability (Centella et al., 2007). These conditions could homogenize nutrient availability and microbial niches across treatments. Another possible explanation is that natural conditions, including heavy rain, could mobilize nutrients through leaching and translocation (lateral movement), thereby homogenizing soil physicochemical properties and microbial communities (Stover and Henry, 2018).

When the bacterial and fungal communities were analyzed at the genus level, a considerable fraction of the identified ASVs remained taxonomically unresolved (65 % of bacteria and 33 % of fungi). The unclassified ASVs at the genus level spanned multiple bacterial and fungal phyla, with representation from Actinomycetota, Acidobacteriota, and Bacillota among bacteria, and Ascomycota, Basidiomycota, and Chytridiomycota among fungi. These groups are typically associated with key ecological processes in soil systems, including carbon decomposition (Actinomycetota, Ascomycota, and Basidiomycota) (Singavarapu et al., 2023), nutrient cycling under oligotrophic conditions (Acidobacteriota) (Gonçalves et al., 2024), stress tolerance and sporulation (Bacillota) (Filippidou et al., 2016), lignin decay and soil aggregate formation (Basidiomycota) (Datta et al., 2017), and the degradation of complex biopolymers (Chytridiomycota) (Zeghal et al., 2021). The presence of numerous unclassified taxa within these ecologically important groups suggests that brown sialitic soils may harbor poorly described or potentially novel microbial lineages that contribute to the unique functioning of these tropical, nutrient-leached systems. This taxonomic ambiguity is likely related to two major factors: (i) the limited representation of tropical and Caribbean soil microorganisms in public reference databases, such as SILVA and UNITE, and (ii) the physicochemical characteristics of sialitic soils, which may select for endemic or highly specialized microbial taxa, such as Stretpomyces (Andam et al., 2016; Martiny, 2016). Despite the classification limitations, these unassigned ASVs are an important component of the microbial diversity in this agroecosystem. These findings highlight the need to expand genomic reference databases to better represent tropical soil microbiomes.

Despite the large proportion of poorly classified microorganisms, both bacterial and fungal community structures remained taxonomically stable across the fertilization gradient, suggesting resilience to moderate variations in nutrient inputs. Interestingly, bacteria have a lower percentage of unclassified taxa at the phylum level than fungi at 1.7 % and 8.8 %, respectively. However, at lower taxonomic levels, bacteria could be classified to a lesser extent than fungi (65.5 % and 33.1 %, respectively). This further emphasizes the limited reference coverage for bacterial taxa from tropical soils.

4.3 Trans-domain co-occurrence networks and taxonomic overlap suggest a stable architecture with nutrient-specific recruitment

Although long-term fertilization did not substantially alter overall community diversity, our global topological analysis revealed a highly connected and stable microbial architecture across the evaluated soils. While diversity metrics were relatively stable, the global co-occurrence network analysis indicated a robust web of potential microbial interactions that persists under long-term management. These findings highlight that the soil microbiome maintains a conserved taxonomic composition and a cohesive connectivity. This is consistent with nutrient enrichment, which homogenizes bulk soil conditions over time through leaching, repeated inputs, and microbial convergence. Meanwhile, the rhizosphere remains structured, likely by plant-driven microenvironments. Root exudates and microbial extracellular polymers promote soil aggregation (Liao et al., 2020; Wilpiszeski et al., 2019), generating microsites with distinct nutrient, carbon, and oxygen dynamics. These can harbor highly structured and functionally specialized microbial associations (Wilpiszeski et al., 2019). Such soil aggregates act as ecological refuges that stabilize microbial communities and promote network-level interactions, even when taxonomic profiles remain similar (Akter et al., 2025; Li et al., 2020; Meena et al., 2022; Wilpiszeski et al., 2019). Therefore, the observed bacterial–fungal co-occurrence network may reflect a conserved functional organization within rhizosphere microhabitats rather than profound structural disruptions triggered by fertilization.

The global architecture was highly modular and dominated by positive correlations (99.8 %), reflecting prevalent mutualistic, commensal, or shared environmental responses rather than competitive exclusion (Cheng et al., 2024). Fungal taxa, particularly Ascomycota and Basidiomycota, exhibited higher connectivity than bacterial nodes, suggesting their role as a potential structural backbone of the network. Although co-occurrence network analysis represents the presence of ASVs on the same sample at some point, these findings suggest that bacterial-fungal functional interdependence is critical for nutrient turnover and the maintenance of soil structure within rhizosphere microhabitats (Fabian et al., 2017). Notably, similar patterns have been reported across diverse ecosystems, where functional conservation is often stronger than taxonomic conservation, and ecosystem processes are maintained by different taxa performing equivalent ecological roles (Hou et al., 2021). In this context, microbial taxa with unique or specialized functions may exert an important influence on ecosystem functioning, sustaining the observed microbial interactions and network organization even in the absence of major taxonomic changes (Tang et al., 2023). This functional redundancy and specialization likely underpins the stable microbial network observed across these soils.

Rather than fragmenting the community, our overlapping taxonomic analysis revealed a core-and-satellite assembly strategy in response to fertilization (Yan et al., 2021). We identified a Universal Resiliency Core of 76 trans-domain ASVs shared across all treatments (N, P, and K). This basal core, heavily reliant on alliances between bacterial (Bacillota and Pseudomonadota) and fungal (Ascomycota and Basidiomycota) taxa, represents the stable structural foundation of the sugarcane rhizosphere, capable of withstanding profound nutritional fluctuations.

Interestingly, while the universal core remains structurally stable, each macronutrient actively recruited its own distinct satellite community. By excluding opportunistic taxa that responded generally to any nutrient addition (22 ASVs), we identified the strict specialists for each fertilization strategy. Nitrogen variations stimulated the emergence of 18 exclusive microbial specialists, indicating a highly specific, fungus-driven response to N enrichment. Conversely, secondary macronutrients acted as much stronger environmental filters for accessory recruitment. The application of phosphorus recruited a massive, specialized accessory microbiome of 40 unique ASVs, heavily dominated by fungi (Talaromyces and various Basidiomycota). Potassium variations exhibited a highly restricted recruitment pattern (six strict ASVs). These findings suggest a threshold-dependent ecological response: while the universal core microbiome is highly resilient, variations in specific nutrients dictate the selective recruitment of specialized, trans-domain accessory niches to exploit available resources (Mu et al., 2024).

4.4 Ecological and agronomic implications

The observed resilience in diversity, coupled with the nutrient-specific recruitment of accessory taxa, indicates that microbial ecosystems in brown sialitic soils maintain a stable core architecture while dynamically exploiting available resources under changing fertilization regimes. In natural systems, fungal–bacterial networks emerge from trophic exchanges, metabolite fluxes, and the physical architecture of hyphae, biofilms, and soil aggregates. Their structure and functioning vary with environmental factors, such as pH, soil moisture, aggregate type and size, and management practices (Gao et al., 2024; Jiménez-Venegas et al., 2024, 2026; Matazimov et al., 2025). These networks can coordinate biogeochemical cycling, influence plant nutrition, and contribute to soil structure and system-level resilience. Within biogeochemical cycles, these microbial frameworks regulate nitrogen and phosphorus mobilization, thereby influencing nutrient retention or leaching (Tedersoo and Bahram, 2019). Their enhancement can increase functional enzyme activity, enabling a 25 %–50 % substitution of synthetic nitrogen with organic nitrogen, while also supporting carbon, nitrogen, and phosphorus cycling (Gao et al., 2024). Likewise, biocrusts formed by fungal–bacterial associations improve the connectivity of carbon pools, stimulating or suppressing carbon uptake depending on the dominant guild, and modulating soil organic carbon gains (Naylor et al., 2020; Tedersoo and Bahram, 2019). Common mycorrhizal networks connect plant roots and facilitate the movement of key nutrients, such as carbon and nitrogen (Matazimov et al., 2025). Fungal hyphae bind soil particles, while bacterial extracellular polymeric substances contribute to aggregate stabilization, influencing pore connectivity and water retention. Furthermore, established fungal–bacterial consortia can help suppress pathogens in the rhizosphere through antibiosis, resource competition, and the induction of plant defenses, enhancing system resilience (Matazimov et al., 2025).

From an agronomic perspective, these findings suggest that further intensification of nitrogen and phosphorus fertilization may not significantly alter overall microbial diversity but drives the targeted recruitment of specialized accessory microbes. Understanding how these distinct ecological niches are occupied without disrupting the overarching stable community framework is crucial for functional stability and soil health. Therefore, sustainable fertilization strategies should aim to support these foundational microbial architectures while optimizing nutrient inputs to avoid potential imbalances in the rhizosphere.

5 Conclusions

This study demonstrates that long-term fertilization in brown sialitic soils does not completely restructure the overall rhizosphere community; rather, specific nutrient applications (such as Phosphorus) act as targeted drivers that recruit and organize specialized trans-domain accessory microbiomes within a highly stable core network. Contrary to our initial expectation that fertilization would primarily alter community structure by shifting taxon abundances, our results support that microbial ecosystems adapt through the selective recruitment of specialized accessory taxa, revealing ecological responses not captured by conventional diversity-based approaches.

Using a unified global trans-domain network, we uncovered a robust and modular rhizosphere architecture predominantly driven by cooperative interactions. Within this framework, a stable Nitrogen-Resiliency Core, supported by alliances between bacterial and fungal phyla, persists across varying nitrogen gradients and may represent a structural foundation of the soil microbiome. Conversely, variations in secondary macronutrients triggered distinct ecological responses. Phosphorus, in particular, acted as a potent environmental filter, awakening and recruiting a specialized, fungal-dominated accessory microbiome to exploit available resources without fragmenting the overall basal community.

Importantly, the observed spatial heterogeneity in soil nutrient availability, largely uncoupled from the theoretical fertilization doses, highlights the intense leaching dynamics and rapid nutrient turnover characteristic of tropical sialitic soils. Nevertheless, the massive buffering capacity conferred by exchangeable basic cations (Ca + Mg) provided a stable edaphic baseline. In these highly weathered environments, physically protected rhizosphere microhabitats provide ecological refuges, allowing microbial interaction networks to maintain ecosystem functioning and resilience despite the chemically erratic bulk soil matrix.

Overall, our findings emphasize that trans-domain interaction networks provide crucial insights into soil ecosystem responses to long-term fertilization that are not evident from diversity or abundance data alone. Incorporating network-based perspectives is essential to a deep understanding of the ecological consequences of fertilization practices and for designing more sustainable nutrient management strategies in tropical agroecosystems.